Network Topology Identification from Spectral Templates
Network topology inference is a cornerstone problem in statistical analyses of complex systems. %Networks and graphs are mathematical objects that encode pairwise relationships. %Under this definition, In this context, the fresh look advocated here permeates benefits from convex optimization and graph signal processing, %A (graph) signal can be understood as a collection of values related by a weighted graph structure $\bbS$ that indicates their level of proximity. Taking this vantage point, this work investigates how to identify the so-termed graph shift operator (encoding the network topology) given only the eigenvectors of the shift. These spectral templates can be obtained, for example, from principal component analysis of a set of graph signals defined on the particular network. The novel idea is to find a graph shift that while being consistent with the provided spectral information, it endows the network structure with certain desired properties such as sparsity. The focus is on developing efficient recovery algorithms along with identifiability conditions for two particular shifts, the adjacency matrix and the normalized graph Laplacian. %In these cases, the problem naturally boils down to network topology identification and we derive conditions under which the shifts are identifiable. Application domains include network topology identification from steady-state signals generated by a diffusion process, and design of a graph filter that facilitates the distributed implementation of a prescribed linear network operator. Numerical tests showcase the effectiveness of the proposed algorithms in recovering synthetic and structural brain networks.
num_experiments.tex

\def \thisplotscale {0.95}
\def \unit {\thisplotscale cm}

{\small \begin{tikzpicture}[thick, x = 1*\unit, y = 1*\unit]

\node (colorbar) at (0.3,-0.05)
{\includegraphics[width=0.025\textwidth]{figures/colorbar.pdf}};

\node (graph_phase_1) at (2.15,1)
{\includegraphics[width=0.15\textwidth]{figures/succ_unique_adj_test.pdf}};

\node (graph_phase_2) at (2.15,-1)
{\includegraphics[width=0.15\textwidth]{figures/succ_unique_lap_test.pdf}};

\node (a) at (2.15,-2.2) {{{\small (a)}}};


\node (graph_phase_3) at (5.2,1)
{\includegraphics[width=0.15\textwidth]{figures/succ_adj_test.pdf}};

\node (graph_phase_4) at (5.2,-1)
{\includegraphics[width=0.15\textwidth]{figures/succ_lap_test.pdf}};

\node (b) at (5.2,-2.2) {{{\small (b)}}};


\node (graph_phase_5) at (9.7, 0)
{\includegraphics[width=0.3\textwidth]{figures/hist_lap_N_10_p_2_100.pdf}};

\node (c) at (9.7,-2.2) {{{\small (c)}}};


\node (comparisons) at (15.2,0)
{\includegraphics[width=0.3\textwidth]{figures/brain_covariance_plot.pdf}};

\node (d) at (15.2,-2.2) {{{\small (d)}}};


\end{tikzpicture}}
tikz_styles.tex
\tikzstyle{phantom vertex} = [ ellipse, 
                               anchor = center, 
                               minimum height = 1*\unit, 
                               minimum width  = 1*\unit,
                               inner sep=0pt,
                               anchor=center]
\tikzstyle{red vertex}   = [black, fill = red!20,   phantom vertex, draw]
\tikzstyle{black vertex} = [black, fill = black!20, phantom vertex, draw]
\tikzstyle{blue vertex}  = [black, fill = blue!20,  phantom vertex, draw]
\tikzstyle{green vertex} = [black, fill = green!20,  phantom vertex, draw]
\tikzstyle{yellow vertex} = [black, fill = yellow!20,  phantom vertex, draw]
\tikzstyle{cyan vertex} = [black, fill = cyan!20,  phantom vertex, draw]
\tikzstyle{vertex}       = [draw, phantom vertex]

\tikzstyle{point} = [ellipse, inner sep=0pt, draw, fill=white, anchor = center,
                     minimum height = 0.05*\unit, minimum width  = 0.05*\unit]

my_sections.tex
\usepackage{needspace}

% \nbsubsubsection{} provides a numbered subsection in bold without a line break. The section will contain at least three lines of text before a pagebreak
\newcommand{\nbsubsubsection}[1]{\needspace{1\baselineskip}\color{white}\subsubsection{#1}\vspace{-3\baselineskip}\color{black}\medskip{\noindent \bf \thesubsubsection. #1.}}

% \myparagraph provides a paragraph title in italics. 
\newcommand{\myparagraph}[1]{\needspace{1\baselineskip}\medskip\noindent {\it #1.}}

% \myindenetedparagraph provides an indented paragraph with title in italics.
\newcommand{\myindentedparagraph}[1]{\needspace{1\baselineskip}\medskip \hangindent=11pt \hangafter=0 \noindent{\it #1.}}

% \myparagraphtc provides a paragraph title in italics. It adds an enter to the table of contents
\newcommand{\myparagraphtc}[1]{\needspace{1\baselineskip}\medskip\noindent {\it #1.}\addcontentsline{toc}{subsubsection}{\qquad\qquad\quad#1}}


TopID_SSP.tex
% Template for SSP-2016 paper; to be used with:
%          spconf.sty  - ICASSP/ICIP LaTeX style file, and
%          IEEEbib.bst - IEEE bibliography style file.
% --------------------------------------------------------------------------
\documentclass{article}
\usepackage{spconf}
\usepackage{color}
\usepackage{theorem}
\usepackage{times,amsmath,epsfig}
\usepackage{amssymb}
\usepackage{subfigure}
\usepackage{cite}
%\usepackage{setspace}
%\doublespacing

%\usepackage{jmlr2e}
%\usepackage{amsmath}
%\usepackage [TABBOTCAP]{subfigure}
%\usepackage{psfrag}
\usepackage{algorithmic}
\usepackage{algorithm}
%\usepackage{stackengine}
%\usepackage[font=small,format=plain,up,up]{caption}
%\usepackage{caption}
%\usepackage{subcaption}
\input{my_sections.tex}
\usepackage{tikz}
\usetikzlibrary{shapes,arrows}
\input{mysymbol.sty}
\input{figures/tikz_styles}

%\renewcommand\floatpagefraction{.9}
%\renewcommand\topfraction{.9}
%\renewcommand\bottomfraction{.9}
%\renewcommand\textfraction{.1}

%\newcommand{\dataset}{{\cal D}}
%\newcommand{\fracpartial}[2]{\frac{\partial #1}{\partial  #2}}
%\renewcommand{\Re}{{\mathbb R}}
%\newcommand{\EE}{{\mathbb E}}
%\newcommand{\diag}{{\mathrm{diag}}}
%\newcommand{\bo}{\mathbf{o}}
\newcommand{\trace}{{\text{trace}}}
\newcommand{\bbeps}{{\boldsymbol{\epsilon}}}
\newcommand{\EE}{{\mathbb E}}
\newcommand{\bo}{\mathbf{o}}
\newcommand{\mymod}{\mathrm{mod}_N}
\newcommand{\Ni}{\mathrm{ni}}
\newcommand{\Nv}{\mathrm{nv}}
\newcommand{\QED}{\hfill\ensuremath{\blacksquare}}



\DeclareMathOperator{\Tr}{Trace}
%\newcommand\barbelow[1]{\stackunder[1.2pt]{$#1$}{\rule{.8ex}{.075ex}}}

\newtheorem{mytheorem}{\bf Theorem}
\newtheorem{mydefinition}{\bf Definition}
\newtheorem{mycorollary}{\bf Corollary}
\newtheorem{mylemma}{\bf Lemma}
\newtheorem{myproposition}{\bf Proposition}
\newtheorem{remark}{\bf Remark}
\newtheorem{property}{\hspace{-11pt}\bf Property}


\newenvironment{myproofsketch}[1][$\!\!$]{{\noindent\bf Proof (sketch) #1: }}
{\hfill\QED\medskip}

% Title.
% ------
\title{Network Topology Identification from Spectral Templates}
%
% Single address.
% ---------------
\name{Santiago Segarra$^{\dag}$, Antonio G. Marques$^{\ddagger}$, Gonzalo Mateos$^{*}$ and Alejandro Ribeiro$^{\dag}$\thanks{Work in this paper is supported by the NSF award  CCF-1217963 and Spanish MINECO TEC2013-41604-R. %Emails: ssegarra@seas.upenn.edu, gmateosb@ece.rochester.edu, antonio.garcia.marques@urjc.es, and aribeiro@seas.upenn.edu.
		}}
\address{$^{\dag}$Dept. of ESE, University of Pennsylvania, Philadelphia, PA, USA\\
	$^{\ddag}$Dept. of TSC, King Juan Carlos University, Madrid, Spain\\
	$^{*}$Dept. of ECE, University of Rochester, Rochester, NY, USA}

%
% For example:
% ------------
%\address{School\\
%	Department\\
%	Address}
%
% Two addresses (uncomment and modify for two-address case).
% ----------------------------------------------------------
%\twoauthors
%  {A. Author-one, B. Author-two\sthanks{Thanks to XYZ agency for funding.}}
%	{School A-B\\
%	Department A-B\\
%	Address A-B}
%  {C. Author-three, D. Author-four\sthanks{The fourth author performed the work
%	while at ...}}
%	{School C-D\\
%	Department C-D\\
%	Address C-D}
%
\begin{document}
\ninept
%
\maketitle
%
\begin{abstract}
Network topology inference is a cornerstone problem in statistical analyses of complex systems.
%Networks and graphs are mathematical objects that encode pairwise relationships. 
%Under this definition, 
In this context, the fresh look advocated here permeates benefits from convex optimization and graph signal processing, 
%A (graph) signal can be understood as a collection of values related by a weighted graph structure $\bbS$ that indicates their level of proximity. Taking this vantage point, this work investigates how 
to identify the so-termed graph shift operator (encoding the network topology) given only the eigenvectors of the shift. These \emph{spectral templates} can be obtained, for example, from principal component analysis of a set of graph signals defined on the particular network. 
The novel idea is to find a graph shift that while being consistent with the provided spectral information, it endows the network structure with certain desired properties such as sparsity. The focus is on developing efficient recovery algorithms along with identifiability conditions for two particular shifts, the adjacency matrix and the normalized graph Laplacian. 
%In these cases, the problem naturally boils down to network topology identification and we derive conditions under which the shifts are identifiable. 
Application domains include network topology identification from steady-state signals generated by a diffusion process, and design of a graph filter that facilitates the distributed implementation of a prescribed linear network operator. Numerical tests
showcase the effectiveness of the proposed algorithms in recovering synthetic and structural brain networks.
\end{abstract}
%
\begin{keywords}
Network topology inference, graph signal processing, spectral graph theory, principal component analysis
\end{keywords}
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION: INTRODUCTION - LET'S MOTIVATE THIS THING!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}\label{S:Introduction}
\vspace{-0.05in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Advancing a holistic theory of networks necessitates fundamental breakthorughs in modeling,
identification, and controllability of distributed network processes -- often conceptualized as \emph{signals defined
on the vertices of a graph}~\cite{barrat2012book,kolaczyk2009book}.
Under the assumption that the signal properties are related to the topology of the graph where they are supported, the goal of graph signal processing (GSP) is to develop algorithms that fruitfully leverage this relational structure~\cite{EmergingFieldGSP,SandryMouraSPG_TSP13}.
Instrumental to that end is the so-termed graph-shift operator (GSO)~\cite{SandryMouraSPG_TSP13}, a matrix capturing the graph's local topology and whose eigenbasis is central to defining graph Fourier transforms~\cite{SandryMouraSPG_TSP14Freq}. Most GSP works assume that the GSO (hence the graph) is known, and then analyze how the algebraic and spectral characteristics of the GSO impact the properties of the signals and filters defined on such a graph. Here instead we take the reverse path and investigate how to use information available from graph signals and filters to infer the underlying graph topology; see also~\cite{DongLaplacianLearning,MeiGraphStructure}. In a nutshell, we follow a two-step approach whereby we first leverage results from GSP theory to identify the GSO's eigenbasis using the available information, and then rely on these \emph{spectral templates} to recover the GSO itself.

%Scenarios where $\bbV$ is known but $\bbS$ is not can be more the rule than the exception....

Network topology inference from a set of (graph-signal) observations is a prominent problem in Network Science~\cite{kolaczyk2009book,sporns2012book}. Since networks encode similarities between nodes, most approaches construct graphs whose edge weights correspond to the correlation, or coherence between signal profiles at incident nodes. 
%To find a sparse connectivity structure, small values are typically thresholded. 
This approach is not without merit and widely used in practice, but it exhibits several drawbacks, the main one being that links are formed taking into account only pairwise interactions, ignoring that the observed correlations can be due to latent network effects. Acknowledging these limitations, alternative methods rely on partial correlations~\cite{GLasso2008,kolaczyk2009book}, Gaussian graphical models~\cite{slawski2015estimation,meinshausen06}, or, Granger causality~\cite{Brovelli04Granger,sporns2012book}.  Differently, recent GSP-based network inference frameworks postulate that the network exists as a latent underlying structure, and that observations are generated as a result of a network process defined in such graph. For instance, network structure is estimated in~\cite{MeiGraphStructure}  to unveil unknown relations among nodal time series adhering to an autoregressive model involving graph-filter dynamics. A factor analysis-based approach was put forth in~\cite{DongLaplacianLearning} to infer graph Laplacians, seeking that input graph signals are smooth over the learned topologies. Different from~\cite{DongLaplacianLearning,MeiGraphStructure} that operate on the graph domain, the goal here is to identify graphs that endow the given observations with desired spectral (frequency-domain) characteristics.

After surveying the required GSP background in Section \ref{S:GSP_101}, we formulate the GSO identification problem given spectral templates (Section \ref{S:ShiftInfFullEigen}). The novel idea is to search among all feasible networks for the one that endows the resulting graph-signal transforms with prescribed spectral properties, while the inferred graph also exhibits desirable structural characteristics.  Optimization problems are formulated to recover two particular GSOs, namely the adjacency matrix (Section \ref{Ss:Adjacency}) and the normalized graph Laplacian (Section \ref{Ss:Laplacian}). Conditions under which the feasible set reduces to a singleton are derived, and convex relaxations leading to computationally-efficient algorithms are proposed for the general case. Section \ref{S:ObtainingSpecTemplates} outlines two pragmatic scenarios where acquisition of the GSO eigenbasis is well motivated, making the case that setups where the GSO is unknown but its spectral templates are available can be more the rule than the exception. Computer simulations highlight the effectiveness of the proposed algorithms in identifying both synthetic and real-world networks (Section \ref{S:Simulations}).  Concluding remarks are given in Section \ref{S:Conclusions}. 

%
%% 
%\footnote{\noindent\emph{Notation:} The entries of a matrix $\mathbf{X}$ and a (column) vector $\mathbf{x}$ will be denoted as $X_{ij}$ and $x_i$;
%however, when contributing to avoid confusion, the alternative notation $[\mathbf{X}]_{ij}$ and $[\mathbf{x}]_i$ will be used. The notation $^T$ and $^H$ stands for transpose and transpose conjugate, respectively; $\diag(\mathbf{x})$ is a diagonal matrix satisfying $[\diag(\mathbf{x})]_{ii}=[\mathbf{x}]_i$; $\mathbf{0}$ is the all-zero vector; and $\mymod(x)$ denotes the modulus (remainder) obtained after dividing $x$ by $N$. Finally, $\| \bbX\|_{0,0}:=\| \mathrm{vec}(\bbX)\|_0$ stands for the number of non-zero entries in $\bbX$. \red{Missing: Vector of all ones, identity matrix, Hadamard product, Frobenius norm} }


% The latter can be achieved by imposing that certain norm of the matrix $(\hat{\bbS}-\bbS)$ has to be small. Since the computational burden associated with finding the optimal solution will depend on the specific formulation, different alternatives will be investigated and corresponding relaxations giving rise to efficient solvers will be designed.
 
 
% Setups where only a few noisy eigenbases are known -- and, hence, combine the challenges in a) and b) -- will also be investigated. As already mentioned, the theoretical analysis characterizing the identifiability and recovery performance of the proposed designs will be completed with numerical experiments. Datasets that will be considered include different types of brain signals, gene expression levels and social network data.


%\newline






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION: GRAPH SP 101, INCLUDING DIFFUSION PROCESSES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Graph signals and graph filters}\label{S:GSP_101}
\vspace{-0.05in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here we formally introduce notation and terminology as well as different GSP tools that will be used throughout the paper. %, including the mathematical definition of a graph signal, the graph shift operator, a graph filter, and the Graph Fourier Transform (GFT).  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%x
\noindent\textbf{Graphs.} Let $\ccalG$ denote an undirected graph with a set of nodes $\mathcal{N}$ (with cardinality $N$) and a set of links $\ccalE$, such that if node $i$ is connected to $j$, then both $(i,j)$ and $(j,i)$ belong to $\ccalE$. The set $\mathcal{N}_i:=\{j\:|(j,i)\in\mathcal{E}\}$ stands for the neighborhood of $i$. 
For any given $\mathcal{G}$ we define the adjacency matrix $\bbA\in\mathbb{R}^{N\times N}$ as a sparse matrix with non-zero elements $A_{ij}=A_{ji}$ if and only if $(i,j)\in\ccalE$. The values of $A_{ji}$ can be binary, or real in the weighted case to capture the strength of the connection from $i$ to $j$. Matrix $\bbA$ can be used to define the degree vector $\bbd:=\bbA\bbone$ and the degree matrix $\bbD:=\diag(\bbd)$. Moreover, the normalized Laplacian matrix is $\bbL:=\bbI - \bbD^{-1/2} \bbA\bbD^{-1/2}\in\mathbb{R}^{N\times N}$, which is positive semidefinite and has eigenvector $\sqrt{\bbd}:=\bbD^{1/2}\mathbf{1}$ with associated eigenvalue zero. 

\noindent\textbf{Graph signals and shift operator.} 
Graph signals defined on the nodes of $\ccalG$ are functions $f : \mathcal{N} \to \mathbb{R}$, equivalently represented as vectors $\mathbf{x}=[x_1,...,x_N]^T \in  \mathbb{R}^N$, where $x_i$ denotes the signal value at node $i$. Since $\bbx$ does not account explicitly for the structure of the graph where the signal is defined, $\mathcal{G}$ is endowed with the so-called GSO $\mathbf{S}$ \cite{SandryMouraSPG_TSP13,SandryMouraSPG_TSP14Freq}. The shift $\mathbf{S}\in\mathbb{R}^{N\times N}$ is a matrix whose entry $S_{ij}$ can be non-zero only if $i=j$ or if $(i,j)\in\mathcal{E}$. The sparsity pattern of  $\bbS$ captures the local structure of $\ccalG$, but we make no specific assumptions on the values of its non-zero entries. The shift $\mathbf{S}$ can also be understood as a linear transformation that can be computed locally at the nodes of the graph. More rigorously, if $\mathbf{y}$ is defined as $\mathbf{y}=\mathbf{S}\mathbf{x}$, then node $i$ can compute $y_i$ provided that it has access to the value of $x_j$ at $j\in \mathcal{N}_i$.
%\label{Foot:footnote1}}.
Typical choices for $\mathbf{S}$ are the adjacency matrix $\bbA$ \cite{SandryMouraSPG_TSP13,SandryMouraSPG_TSP14Freq}, the (normalized) Laplacian $\bbL$~\cite{EmergingFieldGSP}, and their respective generalizations \cite{godsil2001algebraic}. We assume henceforth that $\bbS$ is symmetric, so that $\bbS=\bbV\bbLambda\bbV^{T}$ with $\bbLambda\in\mathbb{R}^{N\times N}$ being diagonal, but our results hold for any normal GSO.
%In particular, $\bbS$ is diagonalizable when it is normal, i.e., it satisfies $\bbS\bbS^H=\bbS^H\bbS$. In that case we have that $\bbV$ is unitary, which implies $\bbV^{-1}=\bbV^{H}$, and leads to the decomposition $\bbS=\bbV\bbLambda\bbV^H$.

\noindent\textbf{Graph filters.} The shift $\bbS$ can be used to define linear, shift-invariant graph-signal \emph{operators} of the form
%
\begin{eqnarray}\label{E:Filter_input_output_time}
	&\mathbf{H}:=\sum_{l=0}^{L-1}h_l \mathbf{S}^l&
\end{eqnarray}
%
which are called \emph{graph filters} \cite{SandryMouraSPG_TSP13}. For a given input $\bbx$, the output of the filter is simply $\bby=\bbH\bbx$. The coefficients of the filter are collected into $\mathbf{h}:=[h_0,\ldots,h_{L-1}]^T$, with $L-1$ denoting the filter degree. Graph filters are of particular interest because they represent linear transformations that can be implemented locally~\cite{EUSIPCO_our_interp_2015,ssamar_distfilters_allerton15}.
%; cf. footnote \ref{Foot:footnote1} and

\noindent\textbf{Frequency domain representation.} Leveraging the spectral decomposition of $\bbS$, graph filters and signals can be represented in the frequency domain.
To be precise, let us use the eigenvectors of $\bbS$ to define the $N\times N$ matrix $\bbU:=\bbV^T$, and the eigenvalues of $\bbS$ to define the $N\times L$ Vandermonde matrix $\bbPsi$, where $\Psi_{ij}:=(\Lambda_{ii})^{j-1}$.
Using these conventions, the frequency representations of a \emph{signal} $\bbx$ and of a \emph{filter} $\bbh$ are defined as $\widehat{\bbx}:=\bbU\bbx$ and $\widehat{\bbh}:=\bbPsi\bbh$, respectively~\cite{SandryMouraSPG_TSP14Freq}. 
Exploiting such representations, the filter's output $\bby\!=\!\bbH\bbx$ in the frequency domain is ($\odot$ denotes Hadamard product)
%
\begin{equation}\label{E:Filter_input_output_freq}
	\widehat{\bby}=\diag\big(\bbPsi\bbh\big)\bbU \bbx=\diag\big(\widehat{\bbh}\big)\widehat{\bbx}=\widehat{\bbh}\odot\widehat{\bbx}.
\end{equation}
%
This identity is the counterpart of the convolution theorem for temporal signals%, and follows from $\mathbf{H}=\bbV\big(\sum_{l=0}^{L-1}h_l \boldsymbol{\Lambda}^l\big)\bbU$ [cf. \eqref{E:Filter_input_output_time}] and $\sum_{l=0}^{L-1}h_l \boldsymbol{\Lambda}^l=\diag{(\bbPsi\bbh)}$
; see e.g. \cite{SSGMAMAR_camsap15} for a derivation.
%\blue{\textbf{Optional:} To establish further connections with the time domain, let us consider the directed cycle graph whose adjacency matrix $\mathbf{A}_{dc}$ is zero, except for entries $A_{ij}=1$ whenever $i=\mymod(j)+1$, where
%$\mymod(x)$ denotes the modulus (remainder) obtained after dividing $x$ by $N$. If $\bbS=\bbA_{dc}$, one can verify that: i) $\bby=\bbH\bbx$ can be found as the circular convolution of $\bbh$ and $\bbx$, and ii) both $\bbU$ and $\boldsymbol{\Psi}$ correspond to the Discrete Fourier Transform (DFT) matrix. }
Interestingly, while in the time domain $\bbU=\boldsymbol{\Psi}$ and they correspond to the Discrete Fourier Transform (DFT) matrix, this is not true for general graphs.

\noindent\textbf{Network diffusion processes.} Graph filters can be used to model network diffusion processes. %that depend on the network topology. To be mathematically precise, 
Specifically, the signal at node $i$ during the step $(l+1)$ of a linear diffusion process can be written as
%
\begin{eqnarray}
\label{E:local_diffusion_singlenode}
&x_i^{(l+1)}=\alpha_{i,i} x_{i}^{(l)} + \sum_{j\in \ccalN_i} \alpha_{i,j} x_{j}^{(l)}&
\end{eqnarray}
%
where % we recall that $\ccalN_i$ stands for the neighborhood of $i$ and 
$\alpha_{i,j}$ are the diffusion coefficients; see e.g.,~\cite{SSGMAMAR_camsap15}. Leveraging the GSP framework, \eqref{E:local_diffusion_singlenode} implies that the graph signal $\bbx^{(l+1)}=\bbS\bbx^{(l)}$ at iteration $l+1$ is the shifted version of $\bbx^{(l)}$, for a shift $\bbS$ with entries  $S_{ij}=\alpha_{i,j}$ if either $i=j$ or $(i,j) \in \ccalE$, and $S_{ij}=0$ otherwise.
%where the shift operator $\bbS$ is defined as: $S_{ij}:=\alpha_{ij}$, if $i=j$ or if $(i,j)\in \ccalE$; and $S_{ij}:=0$, otherwise. 
For instance, if we set $\bbS=\bbI - \beta \bbL$ and let the signal of interest be $\bbx := \bbx^{(\infty)}$, then $\bbx$ solves the heat diffusion equation~\cite{segarra2015reconstruction,segarra2015blindCAMSAP}. However, more complex diffusion dynamics such as $\bbx=\Pi_{l=0}^{\infty}(\bbI-\beta_l \bbS)\bbx^{(0)}$ and $\bbx=\sum_{l=0}^{\infty}\gamma_l \bbx^{(l)}=\sum_{l=0}^{\infty}\gamma_l \bbS^l\bbx^{(0)}$, could also be of interest. 

The Cayley��-Hamilton theorem guarantees that the aforementioned infinite-horizon processes can be equivalently described by a filter of degree $N$.
Accordingly, several works have recognized that the steady-state signal $\bbx$ generated by a diffusion process can be modeled as the output of a graph filter $\bbH=\sum_{l=0}^{N-1}h_l \bbS^l$ with input (seed) $\bbx^{(0)}$~\cite{segarra2015reconstruction,segarra2015blindCAMSAP}.
%Note that although the examples given before consider that the diffusion dynamics run infinite steps, 
%Several works have leveraged the fact that steady-state signals generated by diffusion dynamics can be viewed as outputs of graph filters  The fact that steady-state signals can be written as $\bbx=\bbH\bbx^{(0)}$ 
This key insight can be used to relate the statistical and spectral properties of $\bbx$ and $\bbH$, which will be leveraged in the ensuing sections to identify $\bbS$ itself. 

%For the problem of network topology inference, since the steady-state signals can be written as $\bbx=\bbH\bbx^{(0)}$, one 
%would expect the statistical properties of $\bbx$ to be dependent on the seeding signal and the graph shift $\bbS$ that we aim to obtain. 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION: SHIFT INFERENCE BASED ON SPECTRAL TEMPLATES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Shift inference from spectral templates}\label{S:ShiftInfFullEigen}
\vspace{-0.05in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Given a set of eigenvectors $\bbV=[\bbv_1,\ldots,\bbv_N]$, also termed \emph{spectral templates}, our goal is to find a graph shift $\bbS$ that is diagonalized by $\bbV$. As detailed in Section \ref{S:ObtainingSpecTemplates}, the knowledge of $\bbV$ can for instance be obtained from the PCA decomposition of signals like the ones generated by the diffusion processes outlined in Section \ref{S:GSP_101}; and whose covariance depends on the structure of $\bbS$. Since the postulated problem has infinitely many solutions, we further impose conditions on $\bbS$ promoting desirable properties such as sparsity, or a priori information on the graph of interest such as non-negative edge weights.
Note that by definition, $\bbS$ encodes the local structure of the graph it represents, thus, its recovery implies a successful identification of the graph topology of interest.

\vspace{-0.05in}
\subsection{Identifying the adjacency matrix}\label{Ss:Adjacency}
\vspace{-0.05in}

Consider the unknown adjacency matrix $\bbS=\bbA$ of an unweighted and undirected graph, and suppose its spectral templates $\bbV$ are given. With $\bblambda=[\lambda_1,...,\lambda_N]^T$ collecting the unknown eigenvalues of $\bbS$ (hence $\bbLambda=\text{diag}(\bblambda)$), we aim at identifying $\bbS$ by solving
%defining $\ccalS$ as the set of matrices that satisfy $\ccalS:=\{\bbS\in\reals^{n\times n},\;S_{i,j}\in\{0,1\},\; S_{i,i}=0,\; \bbS\bbone \geq \bbone\}$, the problem to solve is 
\begin{align}
	\label{E:SparseAdj_l00_obj} \min_{ \{\bbS, \bblambda\}} \;\;&\| \bbS\|_{0}\\
	\label{E:SparseAdj_l00_const1}\text{s. to } \;\;&\;\bbS=\textstyle\sum_{k =1}^N \lambda_k\bbv_k \bbv_k^T,\\ \label{E:SparseAdj_l00_const2}&\;S_{ij}\in\{0,1\}, \;\; S_{ii}=0, \;\; \bbS\in\ccalM^N,\;\; \bbS \mathbf{1}\geq \mathbf{1}
\end{align}
%\begin{align}
%	\label{E:SparseAdj_l00_obj} \bbS^*\ =\ \argmin_{ \{\bbS, \bblambda\}} \;\;\| \bbS\|_{0,0}\;\;\;\;
%	\text{s. to } \;\;\bbS=\sum_{k =0}^{n-1} \lambda_k \bbv_k \bbv_k^T,\;\;\bbS\in\ccalS,
%\end{align}
where $\ccalM^N$ denotes the set of $N\times N$ real and symmetric matrices. Among the potentially multiple feasible solutions,  the cardinality function $\|\bbS\|_{0}$ in the objective selects as optimum the one that minimizes the number of edges. The role of constraint $\bbS \bbone \geq \bbone$ (entrywise), which reasonably requires each node to have at least one neighbor, is to prevent the trivial solution $\bbS=\bb0$.
As can be seen from \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2}, when all eigenvectors $\{\bbv_k\}_{k=1}^N$ are given the design of $\bbS$ amounts to finding the $N$ eigenvalues in $\bblambda$. The constraint in \eqref{E:SparseAdj_l00_const1} encodes the definition of a general graph shift $\bbS = \bbV \bbLambda \bbV^T$, while those in \eqref{E:SparseAdj_l00_const2} incorporate the supplementary conditions implied by the fact that $\bbS$ is an adjacency matrix.
%For notational convenience we will rewrite the latter as $\bbS \in \ccalS$ where $\ccalS:=\{\bbS\in\ccalM^N\;\;|S_{i,j}\in\{0,1\},\; S_{i,i}=0,\; \bbS\bbone \geq \bbone\}$. Clearly, the set $\ccalS$ could incorporate additional a priori information on $\bbS$, such as some of the entries being known (either zero or one) or information about the nodes' degree. On the other hand, the feasible set could be enlarged to allow for adjacency matrices with negative entries.
Additional information can be incorporated in \eqref{E:SparseAdj_l00_const2}, such as a priori knowledge on particular entries of $\bbS$ or node degrees. Alternatively, the feasible set can be enlarged to accommodate for adjacency matrices with negative entries.

An interesting property of the proposed optimization is that the feasible set described by \eqref{E:SparseAdj_l00_const1}-\eqref{E:SparseAdj_l00_const2} is generally small. To be more precise, we define the matrix $\bbW  := \bbV\odot\bbV \in \reals_+^{N\times N}$ and denote by $Q$ the number of singular values of $\bbW$ that are zero. %Then, the following proposition holds true.

%To have a better understanding of the structure of the problem of recovering the graph shift $\bbS$ based on its eigenvalues it is key to notice that the set of feasible solutions is small. To discuss this formally, we define the nonnegative matrix $\bbW\in \reals^{N\times N}$ as $\bbW := \bbV\odot\bbV$, which collects the values of $\{\bbv_k\bbv_k^T\}_{k=1}^N$ required to compute the diagonal elements of the graph shift in \eqref{E:SparseAdj_l00_const1}. Matrix $\bbW$ can decomposed as $\bbW=\bbQ \bbSigma \bbT^H$ with $\bbSigma$ being diagonal and $\bbQ$ and $\bbT$ unitary. Let $Q$ denote the number of singular values of $\bbW$ that are zero and $\{\bbq_n\}_{n=1}^{Q}$ the corresponding left singular vectors, so that the column nullspace of $\bbW$ has dimension $Q$ and it is spanned by $\{\bbq_n\}_{n=1}^{Q}$. Based on these notations, the following result holds true.
%
\begin{myproposition}\label{P:Feasibility_Adjacency}
	Assume that \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} is feasible, then it holds that:\\
	a) The nullspace of $\bbW$ has dimension at least one, so that $Q\geq 1$.\\
	b) If $Q=1$, the feasible set given by \eqref{E:SparseAdj_l00_const1}-\eqref{E:SparseAdj_l00_const2} is a singleton. 
\end{myproposition}
%
%\begin{myproof}
\vspace{-0.05in}
\noindent {\bf Proof (sketch):} The key of the proof resides in noting that the columns of $\bbW$ correspond to the diagonal entries of $\bbv_k \bbv_k^T$ for $k = 1, \ldots, N$. Then, by combining \eqref{E:SparseAdj_l00_const1} and \eqref{E:SparseAdj_l00_const2} it follows that $\bbW \bblambda = \mathbf{0}$ for all feasible $\bblambda$, thus, feasibility implies $Q\geq 1$. When $Q=1$, $\bblambda$ (and hence $\bbS$) is unique up to a scaling factor. Therefore using the fact that $\bbS$ is binary, statement b) follows. {\hfill\QED\medskip}

%\end{myproof}

\noindent \textbf{Relaxation and algorithmic discussion.}
The two sources of non-convexity that render the solution to \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} challenging are the presence of the $\ell_0$ norm in the objective and the binary constraints $S_{ij}\in\{0,1\}$. We relax the former to an iteratively re-weighted $\ell_1$ norm and the latter by replacing $\{0,1\}$ with its convex hull $[0,1]$. Specifically, with $p$ denoting an iteration index, we aim to solve a sequence $p=1,...,P$ of $\ell_1$-norm penalized problems
%
\begin{align}
	\label{E:SparseAdj_wl11_obj} \min_{ \{\bbS, \bblambda\}} \;\;& { \textstyle\sum_{i=1}^N\sum_{j=1}^N} \omega_{ij}(p)|S_{ij}| \\
	\label{E:SparseAdj_wl11_const1}\text{s. to } \;\;&\;\bbS=\textstyle\sum_{k =1}^N \lambda_k\bbv_k \bbv_k^T,\\ \label{E:SparseAdj_wl11_const2}&\;S_{ij}\in[0,1], \;\; S_{ii}=0, \;\; \bbS\in\ccalM^N,\;\; \bbS \mathbf{1}\geq \mathbf{1}
\end{align}
%
with weights $\omega_{ij}(p):= \left(|S_{ij}(p-1)|+\delta\right)^{-1}$ where $\delta$ is a small and positive constant. Intuitively, the goal of the re-weighted scheme in \eqref{E:SparseAdj_wl11_obj} is that if the value of $|S_{ij}(p-1)|$ is small, in the next iteration the penalization weight $\omega_{ij}(p)$ is large, promoting further shrinkage of $S_{ij}$ towards zero; see, e.g., \cite{candes_l0_surrogate} for technical details. The recovery performance of the proposed scheme is analyzed in Section~\ref{S:Simulations}.

Building on Prop.~\ref{P:Feasibility_Adjacency}, whenever $Q = 1$ the solutions of \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2} and \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} coincide up to a scaling factor. 
%
\begin{myproposition}\label{P:Identifiability_Adjacency_Norm1}
	Let $\bbS_0^*$ and $\bbS_1^*$ denote solutions to \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} and \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2}, respectively. If $Q=1$, then $\bbS_1^*=d_{\min}^{-1} \bbS_0^*$, with $d_{\min}$ being the minimum node degree in the graph to recover.    
\end{myproposition}
%
\vspace{-0.05in}
\begin{myproofsketch}
%\noindent {\bf Proof (sketch):} 
	Denoting by $\bblambda_0^*$ the unique solution to \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} (cf.~Prop.~\ref{P:Feasibility_Adjacency}), the solution to \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2} is $\bblambda_1^* =  \alpha \bblambda_0^*$  with $\alpha$ being the smallest positive number satisfying $\bbS \mathbf{1}\geq \mathbf{1}$, which is $\alpha=d_{\min}^{-1}$.
\end{myproofsketch}

%If the dimension of the nullspace of $\bbW = \bbV \odot \bbV$ is one, the solution generated by \eqref{E:SparseAdj_wl11_obj} is the same than the solution to \eqref{E:SparseAdj_l00_obj} up to a scaling factor. When the dimension is greater than one, \eqref{E:SparseAdj_wl11_obj} yields a sparse solution, which oftentimes coincides with the solution to \eqref{E:SparseAdj_l00_obj}. This is formalized in the next result. 
%\begin{myproposition}\label{P:Identifiability_Adjacency_Norm1}
%Let $S_0^*$ and $S_1^*$ denote a solution to \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} and \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2}, respectively, and recall that $Q$ denotes the dimension of the nullspace of $\bbW=\bbV\odot\bbV$. Then, it holds that\\
%a) If $Q=1$, then $\bbS_1^*=c_a^{-1} \bbS_0^*$, with $c_a$ being the degree of node $1$.\\    
%b) If $Q\geq 2$, there exist tuning constants in \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2} such that $\bbS_1^*=c_b^{-1}\bbS_0^*$, with $c_b$ being a generic non-zero constant. 
%\end{myproposition}

\vspace{-0.05in}
\subsection{Identifying the Laplacian matrix and generalized shifts}\label{Ss:Laplacian}
\vspace{-0.05in}

Consider the unknown normalized Laplacian matrix $\bbS=\bbL$ of an unweighted and undirected graph. With $\ccalM_{+}^N$ denoting the conical set of positive semidefinite $N\times N$ matrices, the problem to solve now is [cf. \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2}]
%
\begin{align}
	\label{E:SparseLap_l00_obj} \min_{ \{\bbS, \bblambda\}} \;\;&\| \bbS\|_{0}\\
	\label{E:SparseLap_l00_const1}\text{s. to } \;\;&\;\bbS=\textstyle\sum_{k =1}^N \lambda_k\bbv_k \bbv_k^T,\;\;\;\; \bbS\! \in \! \ccalM_{+}^N, \;\;\;\; \lambda_1=0,\\ 
	%\;\;\bblambda\succeq \bbzero
	\label{E:SparseLap_l00_const2}&\;\;\,\,\,\, S_{ij} \in [-1, 0]  \,\, \text{for} \,\, i\neq j, \qquad S_{ii}=1 \,\, \forall \,\, i.
\end{align}
%
In the above formulation we assumed that $\bbv_1$ corresponds to the eigenvector $\sqrt{\bbd}$ so that we can leverage the structure of normalized Laplacians and enforce $\lambda_1=0$. This, in turn, prevents the trivial solution $\bbS=\bbI$. Note also that the eigenvector $\sqrt{\bbd}$ can be easily identified from the $N$ eigenvectors in $\bbV$ since it is the only one whose entries have all the same sign \cite{biyikougu2007laplacian}. As was the case for \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2}, if more than one feasible solution exists, the identified Laplacian corresponds to the topology that minimizes the number of edges. 
%Notice that for connected shifts one of the eigenvectors, say the first one, is $\bbv_1 = \frac{1}{\sqrt{N}}\bbone$, so that the condition $\bbS\bbone=\bbzero$ can be equivalently written as $\lambda_1=0$. 

Interestingly, problem \eqref{E:SparseLap_l00_obj}-\eqref{E:SparseLap_l00_const2} also has a small feasible set. To state this formally, consider the matrix $\tilde{\bbV}=[\bbone, \bbv_2, \bbv_3, \ldots, \bbv_N ]$, define $\tilde{\bbW}:=\tilde{\bbV}\odot\tilde{\bbV}$ and denote by $\tilde{Q}$ the number of singular values of $\tilde{\bbW}$ that are zero. Then the following result, which is the counterpart of Prop.~\ref{P:Feasibility_Adjacency} for Laplacian identification, holds.
%
\begin{myproposition}\label{P:Feasibility_normLaplacian}
	Assume that \eqref{E:SparseLap_l00_obj}-\eqref{E:SparseLap_l00_const2} is feasible, then it holds that:\\
	a) The nullspace of $\tilde{\bbW}$ has at least dimension one, so that $\tilde{Q}\geq 1$.\\
	b) If $\tilde{Q}=1$, the feasible set given by \eqref{E:SparseLap_l00_const1}-\eqref{E:SparseLap_l00_const2} is a singleton.
	%\red{and the entries of vector $\tilde{\bbq}_1$ satisfy that $\tilde{q}_{1,1}\neq 0$ and $\sign{(\tilde{q}_{1,n})}\neq\sign{(\tilde{q}_{1,1})}$ for all $n\neq 1$}
\end{myproposition}

The result follows using arguments similar to those in the proof of Prop.~\ref{P:Feasibility_Adjacency}. Regarding the design of efficient algorithms to solve \eqref{E:SparseLap_l00_obj}-\eqref{E:SparseLap_l00_const2}, since \eqref{E:SparseLap_l00_obj} is the only source of non-convexity, the iteratively re-weighted scheme in \eqref{E:SparseAdj_wl11_obj} can also be applied here. Details are omitted due to space limitations.

Modifications to the proposed formulations can be made to accommodate for shifts that exhibit properties other than sparsity. For example, one could go after Laplacians with good mixing conditions, related to having a large second smallest eigenvalue \cite{chung1997spectral}. This can be accomplished by introducing the optimization variable $\lambda_{\min}$, adding the constraints $\lambda_{\min}\leq\lambda_k$ for all $k\geq 2$ and augmenting the cost with the regularization $-\eta \lambda_{\min}$, with $\eta$ being a tuning constant.


The identification of graph shifts different from the adjacency and the normalized Laplacian can be of interest, including the combinatorial Laplacian $\bbL=\bbD - \bbA$ and the random walk Laplacian $\bbL=\bbD^{-1}\bbA$. These would require minor modifications in the constraints \eqref{E:SparseLap_l00_const1}-\eqref{E:SparseLap_l00_const2}. For the particular case of the combinatorial Laplacian, results in Prop.~\ref{P:Feasibility_normLaplacian} can be extended directly provided that the degree vector $\bbd$, which contains the same information than the first eigenvector of the normalized Laplacian, is known. Specifically, the constraints $S_{ii}=1$ in \eqref{E:SparseLap_l00_const2} need to be replaced with $S_{ii}=d_i$ and matrix $\tilde{\bbV}$, which is used in the definition of $\tilde{\bbW}=\tilde{\bbV}\odot\tilde{\bbV}$, has to be redefined as $\tilde{\bbV}=[\bbd, \bbv_2, \bbv_3, \ldots, \bbv_N ]$. 
%\subsection{Identifiability conditions}
%Besides exploring different alternative formulations for \eqref{E:SparseAdj_l00_obj}, our research will also encompass the design of convex approximations that can be solved efficiently as well as theoretical and experimental analysis to characterize the identifiability and recovery performance.    

%With regard to the algorithm to solve the non-convex optimization in \eqref{E:SparseAdj_l00_obj}, a reasonable approach is to replace the zero norm in the objective with a (possible reweighed) one norm and the binary set $\{0,1\}$ with its convex hull $[0,1]$. Identifying conditions under which the convex relaxation yields the same solution than the original non-convex optimization is of interest. Theoretical and experimental analysis of the recovery properties of the proposed approach will also be pursued.

%%%%%%%%%%%%%%%   F   I   G   U   R    E   :   R    A    N    D        N   E   T   W   O   R   K   %%%%%%%%%%%%%%%%%%%%%
\begin{figure*}[t]
	\centering
	\input{figures/num_experiments.tex}	
	\vspace{-0.175in}
	\caption{(a) Proportion of topology identification problems with unique solution for Erd\H{o}s-R\'enyi graphs as a function of their size $N$ and the probability of edge appearance $p$ for adjacency (top) and normalized Laplacian (bottom) matrices. (b) Recovery rate for the same set of graphs used in (a) when implementing the iteratively re-weighted approach. (c) Histogram of the rank of matrix $\tilde{\bbW}$ for $N=10$ and $p=0.2$. (d) Brain graph recovery error for three patients as a function of the number of signals observed in the estimation of the spectral templates. }
		\vspace{-0.15in}
	\label{F:num_exp}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION: OBTAINING THE GRAPH SPECTRAL TEMPLATES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Obtaining the graph spectral templates}\label{S:ObtainingSpecTemplates}
\vspace{-0.05in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The methods discussed in Section~\ref{S:ShiftInfFullEigen} are relevant in practice when we have access to the eigenvectors (i.e., the spectral templates) of $\bbS$, but not to the GSO itself. A couple of scenarios where that is indeed the case are presented next.

\vspace{0.05in}
\noindent\textbf{Signals generated through diffusion processes.}
Recall that graph filters can be used to model diffusion processes that depend on the topology of the graph (cf.~Section~\ref{S:GSP_101}). More specifically, a diffused or steady-state signal $\bbx$ can written as $\bbx=\bbH \bbx^{(0)}$ where $\bbx^{(0)}$ is a seed signal and $\bbH=\sum_{l=0}^{N-1}h_l \bbS^l$ is a graph filter.
Under the assumption that $\bbx^{(0)}$ is white (with identity covariance matrix) and zero mean, the covariance matrix of the output $\bbC_\bbx:=\E{\bbx\bbx^T}$ is given by
%
\vspace{-0.03in}
\begin{align}\label{E:cov_output_filter}
\hspace{-0.28cm}\bbC_\bbx=\bbH \E{\bbx^{(0)}({\bbx^{(0)}})^T}\bbH^T=\bbH\bbH^T\stackrel{\star}{=}\bbV\diag(|\widehat{\bbh}|^2)\bbV^T,
\end{align}
\vspace{-0.03in}
%
where for equality $\stackrel{\star}{=}$ we leveraged the frequency interpretation of a graph filter. Expression \eqref{E:cov_output_filter} reveals that the eigenvectors of the covariance matrix $\bbC_\bbx$ and those of $\bbS$ are the same. Thus, if $\bbC_\bbx$ is known, the spectral templates $\bbV$ can be readily obtained. More commonly, if $\bbC_\bbx$ is unknown but we have access to a set of $M$ diffused signals $\{\bbx_m\}_{m=1}^M$, we may approximate $\bbC_\bbx$ with the sample covariance $\hat{\bbC}_\bbx=1/M\sum_{m=1}^M \bbx_m \bbx_m^T$ and estimate the eigenvectors of $\bbS$. Applying the methods in Section~\ref{S:ShiftInfFullEigen}, we can then use the estimated eigenvectors $\hat{\bbV}$ to recover $\bbS$ as illustrated in Section~\ref{S:Simulations}. Detailed analysis on the sensibility of the topology identification performance to the errors in estimating $\hat{\bbV}$ is left as future work.

\vspace{0.05in}
\noindent\textbf{Implementation of linear network operators.}
Matrices $\bbB \in \reals^{N \times N}$ can be viewed as linear network operators since they are linear maps between graph signals. Conceivably, one might be interested in representing (or approximating) a pre-specified operator $\bbB$ as a graph filter $\bbH$, to facilitate its distributed implementation in the network. This problem is thoroughly investigated in \cite{segarra2015graphfilteringTSP15}, and one of the conditions derived for perfect implementation states that all the eigenvectors of $\bbB$ and $\bbS$ must coincide \cite[Prop. 1]{segarra2015graphfilteringTSP15}. Consequently, we can use the eigenvectors of a desired $\bbB$ as spectral templates to generate a graph shift $\bbS$. This way, $\bbB$ can be expressed as a graph filter and implemented via local agent interactions.
%
%\vspace{0.05in}
%\noindent\textbf{Bandlimited graph signals.}
%Considerable attention has been given to the sampling and interpolation of bandlimited graph signals \red{REF}, i.e. signals with sparse frequency representation. From a reasoning similar to that in \eqref{E:cov_output_filter}, it follows that from the covariance matrix of a $K$-bandlimited signal, i.e. with only $K$ non-zero frequencies, the first $K$ eigenvectors of the shift $\bbS$ can be extracted. More specifically, expressing the bandlimited signal as $\bbx = \bbV_K \widehat{\bbx}_K$ (cf.~Section~\ref{S:GSP_101}) it follows that $\bbC_\bbx=\bbV_K \E{\widehat{\bbx}_K\widehat{\bbx}_K^H} \bbV_K^H$. Furthermore, assuming that the active frequency coefficients are independent, $\E{\widehat{\bbx}_K\widehat{\bbx}_K^H}$ is a diagonal matrix implying that $\bbV_K$ can be recovered as the set of eigenvectors associated with the non-zero eigenvalues of $\bbC_\bbx$. Methods similar to those in Section~\ref{S:ShiftInfFullEigen} but where only $K$ eigenvectors are given can be designed but are left as future work. Notice that, in such a case, the additional degrees of freedom coming from the unspecified spectral templates can be leveraged to endow the recovered graph with desired properties.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION: NUMERICAL EXPERIMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{-0.05in}
\section{Numerical experiments}\label{S:Simulations}
\vspace{-0.05in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We illustrate the performance of the topology identification algorithm by solving \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2}, as well as its counterpart for Laplacian identification, for different synthetic and real-world graphs. 

\vspace{0.05in}
\noindent\textbf{Random graphs.}
Consider Erd\H{o}s-R\'enyi graphs \cite{bollobas1998random} of varying size $N \in \{10, 20, \dots, 50\}$ and different edge-formation probabilities $p \in \{0.1, 0.2, \ldots, 0.9\}$. For each combination of $N$ and $p$ we generate 100 graphs and try to recover their adjacency $\bbA$ and normalized Laplacian $\bbL$ matrices from the corresponding spectral templates $\bbV$. In Fig.~\ref{F:num_exp}(a) we plot the proportion of instances where the corresponding optimization problems -- \eqref{E:SparseAdj_l00_obj}-\eqref{E:SparseAdj_l00_const2} for $\bbA$ and \eqref{E:SparseLap_l00_obj}-\eqref{E:SparseLap_l00_const2} for $\bbL$ -- have unique solutions. Notice that multiple solutions are more frequent when the expected number of neighbors of a given node is close to either $1$ or $N$. For intermediate values of $p$, the rank of both $\bbW$ and $\tilde{\bbW}$ is typically $N-1$, guaranteeing that the solution to our optimization is unique (cf. Props.~\ref{P:Feasibility_Adjacency} and~\ref{P:Feasibility_normLaplacian}). Using the same set of graphs that those in Fig.~\ref{F:num_exp}(a), Fig.~\ref{F:num_exp}(b) shows the recovery rate when solving the iteratively re-weighted problem in \eqref{E:SparseAdj_wl11_obj}-\eqref{E:SparseAdj_wl11_const2} and its counterpart for the Laplacian. As expected, the rates in Fig.~\ref{F:num_exp}(b) dominate those in Fig.~\ref{F:num_exp}(a) since every instance with a unique solution is recovered successfully. Moreover, the improvement in the rates observed in Fig.~\ref{F:num_exp}(b) reflect the effect of the weighted $\ell_1$ norm objective [cf.~\eqref{E:SparseAdj_wl11_obj}] in recovering the true underlying graph.

As indicated by Props.~\ref{P:Feasibility_Adjacency} and~\ref{P:Feasibility_normLaplacian}, the rate of recovery is intimately related to the ranks of $\bbW$ and $\tilde{\bbW}$ for the adjacency and normalized Laplacian case, respectively. Fig.~\ref{F:num_exp}(c) further illustrates this relation via a histogram of the rank of $\tilde{\bbW}$ for the 100 graphs with $N=10$ and $p=0.2$. For more than half of the instances, the rank of $\tilde{\bbW}$ was equal to 9 (blue bar) and, as stated in Prop.~\ref{P:Feasibility_normLaplacian}, for all these graphs the solution was unique (yellow bar) and successfully recovered (cyan bar). We see that, as the rank of $\tilde{\bbW}$ degrades, uniqueness is no longer guaranteed but for most cases the true graph can still be recovered following the iteratively re-weighted scheme proposed. Only in 8 of the cases where $\mathrm{rank}(\tilde{\bbW})<9$ the recovery was not successful, entailing a recovery rate of 0.92, as reported in the corresponding entry ($N=10$, $p=0.2$) of the lower plot in Fig.~\ref{F:num_exp}(b).

\vspace{0.05in}
\noindent\textbf{Brain graphs.}
We consider the identification of unweighted and undirected graphs corresponding to human brains~\cite{hagmann2008mapping}, consisting of $N=66$ nodes or regions of interest (ROIs) and whose edges link ROIs with density of anatomical connections greater than a threshold. The threshold is chosen as the largest one that entails a connected graph. We test the recovery performance for noisy spectral templates $\hat{\bbV}$ obtained from sample covariances of signals generated through diffusion processes (cf.~Section~\ref{S:ObtainingSpecTemplates}). 
%More specifically, we generate random white signals $\bbx^{(0)}$, pass them through a filter $\bbH$ whose coefficients are chosen at random, and observe the output $\bbx = \bbH \bbx^{(0)}$.
%More specifically, with $\bbA_i = \bbV_i \bbLambda_i \bbV_i^{-1}$ denoting the adjacency matrix (to be recovered) for patient $i \in \{1, 2, 3\}$, we generate perturbed versions of $\bbV_i$ by normalizing the columns of $\hat{\bbV}_i= \bbV_i + \sigma \, \bbV_i \odot \bbZ$, where the elements of $\bbZ$ are realizations of independent normal random variables and the scalar $\sigma$ controls the magnitude of the perturbation. 
Our framework can handle the recovery of $\bbS$ from noisy versions of $\bbV$ by rewriting the equality constraint \eqref{E:SparseAdj_wl11_const1} as an element-wise inequality bounding the absolute difference between $S_{ij}$ and $\textstyle\sum_{k =1}^N \lambda_k V_{ki} V_{kj}$ for all $i,j$. Denoting by $\hat{\bbV}_i$ the noisy spectral templates of patient $i \! \in \! \{1, 2, 3\}$ and by $\hat{\bbA}_i$ the adjacency matrices recovered from them, Fig.~\ref{F:num_exp}(d) plots the recovery error as a function of the number of signals observed in the computation of the sample covariance. The error is quantified as the proportion of edges misidentified, i.e., $\|\bbA_i - \hat{\bbA}_i\|_0 / \| \bbA_i \|_0$, and each point in Fig.~\ref{F:num_exp}(d) is the average across 50 random realizations. 
%First notice that in the noiseless case ($\sigma = 0$), perfect recovery is achieved for the three patients. 
First notice that for an increasing number of observed signals we see a monotonous decrease in recovery error. For example, when going from $10^4$ to $10^5$ observations the error is (approximately) divided by seven, when averaged across patients. This is reasonable since a larger number of observations gives rise to a more reliable estimate of the covariance matrix entailing a less noisy version of the spectral templates. Notice that, even though the performance increases for all patients with the number of observations, the brain of patient 1 is consistently the hardest to identify. For instance, consider the errors for $10^5$ observations, where for patient 1 we successfully recover $96.4\%$ of the \textit{edges}, whereas for patients 2 and 3 the average recovery rate is $99.5\%$ and $99.4\%$, respectively. This points towards the fact that some graphs are inherently more robust for identification when given noisy spectral templates. A formal analysis of this phenomenon is left as future work.

Traditional methods like graphical lasso \cite{GLasso2008} fail to recover $\bbS$ from the sample covariance of filtered white signals. This occurs because, in general, the filter $\bbH$ introduces conditional dependence between signal values more than one hop apart. Even when focusing on filters of the form $\bbH = h_0 \bbI + h_1 \bbS$, graphical lasso underperforms compared to the method presented. More precisely, based on $10^5$ observations, the recovery error of graphical lasso averaged over 50 realizations and with optimal tuning parameters is 0.303, 0.350, and 0.270 for patients 1, 2, and 3, respectively [cf. Fig.~\ref{F:num_exp}(d)].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION: CONCLUDING REMARKS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{-0.05in}
\section{Conclusions}\label{S:Conclusions}
\vspace{-0.05in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We studied the problem of identifying a graph shift operator $\bbS$, i.e., the topology of a graph $\ccalG$ of interest, given its eigenbasis $\bbV$. The focus was on setups where $\bbS$ represents the adjacency or the normalized Laplacian of $\ccalG$. We formulated optimization problems to recover $\bbS$, presented their convex relaxations, and derived conditions under which the relaxed problem is guaranteed to recover the desired shift. To highlight the practical relevance of the proposed schemes we outlined scenarios where $\bbV$ is known but $\bbS$ is not, and carried out numerical tests showing the effectiveness in recovering synthetic and brain graphs, even from imperfect spectral templates.


%\red{Estimated length of material to be included:}

%\red{1.25 pages of title, abstract, intro and preliminaries}

%\red{1.25 pages of ``Shift inference using the full eigenbasis"}

%\red{0.5 pages of ``Obtaining the graph spectral templates"}

%\red{1 page of numerical experiments and conclusion}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliographystyle{IEEEbib}
\bibliography{citations}

\end{document}

TopID_SSP.bbl
\begin{thebibliography}{10}

\bibitem{barrat2012book}
A.~Barrat, M.~Barth\'{e}lemy, and A.~Vespignani,
\newblock {\em Dynamic Processes on Complex Networks},
\newblock Cambridge University Press, Cambridge, UK, 2012.

\bibitem{kolaczyk2009book}
E.~D. Kolaczyk,
\newblock {\em Statistical Analysis of Network Data: Methods and Models},
\newblock Springer, New York, NY, 2009.

\bibitem{EmergingFieldGSP}
D.I. Shuman, S.K. Narang, P.~Frossard, A.~Ortega, and P.~Vandergheynst,
\newblock ``The emerging field of signal processing on graphs: Extending
  high-dimensional data analysis to networks and other irregular domains,''
\newblock {\em IEEE Signal Process. Mag.}, vol. 30, no. 3, pp. 83--98, Mar.
  2013.

\bibitem{SandryMouraSPG_TSP13}
A.~Sandryhaila and J.M.F. Moura,
\newblock ``Discrete signal processing on graphs,''
\newblock {\em IEEE Trans. Signal Process.}, vol. 61, no. 7, pp. 1644--1656,
  Apr. 2013.

\bibitem{SandryMouraSPG_TSP14Freq}
A.~Sandryhaila and J.M.F. Moura,
\newblock ``Discrete signal processing on graphs: Frequency analysis,''
\newblock {\em IEEE Trans. Signal Process.}, vol. 62, no. 12, pp. 3042--3054,
  June 2014.

\bibitem{DongLaplacianLearning}
X.~Dong, D.~Thanou, P.~Frossard, and P.~Vandergheynst,
\newblock ``Learning laplacian matrix in smooth graph signal representations,''
\newblock {\em arXiv preprint arXiv:1406.7842v2 [cs.LG]}, 2015.

\bibitem{MeiGraphStructure}
J.~Mei and J.~M.~F. Moura,
\newblock ``Signal processing on graphs: Modeling (causal) relations in big
  data,''
\newblock {\em arXiv preprint arXiv:1503.00173v1 [cs.IT]}, 2015.

\bibitem{sporns2012book}
O.~Sporns,
\newblock {\em Discovering the Human Connectome},
\newblock MIT Press, Boston, MA, 2012.

\bibitem{GLasso2008}
J.~Friedman, T.~Hastie, and R.~Tibshirani,
\newblock ``Sparse inverse covariance estimation with the graphical lasso,''
\newblock {\em Biostatistics}, vol. 9, no. 3, pp. 432--441, 2008.

\bibitem{slawski2015estimation}
M.~Slawski and M.~Hein,
\newblock ``Estimation of positive definite m-matrices and structure learning
  for attractive gaussian markov random fields,''
\newblock {\em Linear Algebra and its Applications}, vol. 473, pp. 145--179,
  2015.

\bibitem{meinshausen06}
N.~Meinshausen and P.~Buhlmann,
\newblock ``High-dimensional graphs and variable selection with the lasso,''
\newblock {\em Ann. Stat.}, vol. 34, pp. 1436--1462, 2006.

\bibitem{Brovelli04Granger}
A.~Brovelli, M.~Ding, A.~Ledberg, Y.~Chen, R.~Nakamura, and S.~L. Bressler,
\newblock ``Beta oscillations in a large-scale sensorimotor cortical network:
  directional influences revealed by granger causality,''
\newblock {\em PNAS}, vol. 101, pp. 9849–9854, 2004.

\bibitem{godsil2001algebraic}
C.~Godsil and G.~Royle,
\newblock {\em Algebraic Graph Theory}, vol. 207,
\newblock Springer-Verlag, Graduate Texts in Mathematics, 2001.

\bibitem{EUSIPCO_our_interp_2015}
S.~Segarra, A.~G. Marques, G.~Leus, and A.~Ribeiro,
\newblock ``Interpolation of graph signals using shift-invariant graph
  filters,''
\newblock in {\em European Signal Process. Conf. (EUSIPCO)}, Nice, France, Aug.
  31 - Sept. 4 2015.

\bibitem{ssamar_distfilters_allerton15}
S.~Segarra, A.~G. Marques, and A.~Ribeiro,
\newblock ``Distributed implementation of linear network operators using graph
  filters,''
\newblock in {\em 53rd Allerton Conf. on Commun. Control and Computing}, Univ.
  of Illinois at U-C, Monticello, IL, Sept. 30- Oct. 2, 2015.

\bibitem{SSGMAMAR_camsap15}
S.~Segarra, G.~Mateos, A.~G. Marques, and A.~Ribeiro,
\newblock ``Blind identification of graph filters with sparse inputs,''
\newblock in {\em IEEE Intl. Wrksp. Computat. Advances Multi-Sensor Adaptive
  Process. (CAMSAP)}, Cancun, Mexico, Dec. 13-16, 2015.

\bibitem{segarra2015reconstruction}
S.~Segarra, A.~G. Marques, G.~Leus, and A.~Ribeiro,
\newblock ``Reconstruction of graph signals through percolation from seeding
  nodes,''
\newblock {\em arXiv preprint arXiv:1507.08364}, 2015.

\bibitem{segarra2015blindCAMSAP}
S.~Segarra, G.~Mateos, A.~G. Marques, and A.~Ribeiro,
\newblock ``Blind identification of graph filters with sparse inputs,''
\newblock in {\em Computational Advances in Multi-Sensor Adaptive Processing
  (CAMSAP), 2015 IEEE 6th International Workshop on}. IEEE, 2015.

\bibitem{candes_l0_surrogate}
E.~J. Candes, M.~B. Wakin, and S.~Boyd,
\newblock ``Enhancing sparsity by reweighted $\ell_1$ minimzation,''
\newblock {\em Journal of Fourier Analysis and Applications}, vol. 14, pp.
  877--905, Dec. 2008.

\bibitem{biyikougu2007laplacian}
T.~Biyiko{\u{g}}u, J.~Leydold, and P.~Stadler,
\newblock {\em Laplacian Eigenvectors of Graphs: Perron-Frobenius and
  Faber-Krahn Type Theorems},
\newblock Springer-Verlag Berlin Heidelberg, 2007.

\bibitem{chung1997spectral}
F.~Chung,
\newblock {\em Spectral Graph Theory}, vol.~92,
\newblock American Mathematical Soc., 1997.

\bibitem{segarra2015graphfilteringTSP15}
S.~Segarra, A.~G. Marques, and A.~Ribeiro,
\newblock ``Distributed linear network operators using graph filters,''
\newblock {\em arXiv preprint arXiv:1510.03947}, 2015.

\bibitem{bollobas1998random}
B.~Bollob{\'a}s,
\newblock {\em Random Graphs},
\newblock Springer, 1998.

\bibitem{hagmann2008mapping}
P.~Hagmann, L.~Cammoun, X.~Gigandet, R.~Meuli, C.~J. Honey, V.~J Wedeen, and
  O.~Sporns,
\newblock ``Mapping the structural core of human cerebral cortex,''
\newblock {\em PLoS Biol}, vol. 6, no. 7, pp. e159, 2008.

\end{thebibliography}
