\section{Introduction}
\label{sec:intro}

\cmt{Motivation}
\begin{myitemize}%
\cmt{monitoring spatial fields vs (possibly W)SNs}A frequent
inference problem in signal processing involves the estimation of a
spatial field using measurements collected by a (possibly wireless)
sensor
network~\cite{raghavendra2006,chen2015ubiquitous,nevat2015heterogeneous,asensio2015energy,nowak2004inhomogeneous}.
The field of interest may quantify magnitudes such as temperature,
electromagnetic radiation, concentration of airborne or liquid
pollutants, flows of gas or liquids in porous soils and rocks such as
oil reservoirs, acoustic pressure, and radioactivity to name a
few. Instead of spatial fields, one may  be alternatively interested in fields
defined on the nodes or edges of a network; see
e.g.~\cite{forero2014dictionary}.
\cmt{Inference$\rightarrow$~ projection}In either case, a number of
common inference tasks such as least squares estimation, denoising,
(weighted) consensus, and decentralized detection can be cast as projecting
the observations onto a given signal subspace; see
e.g.~\cite{kay1,harsanyi1994hyperspectral,behrens1994oblique}.
\cmt{decentralized}
\begin{myitemize}%
\cmt{centralized}Such a fundamental task can be implemented in
a centralized fashion,
\begin{myitemize}%
\cmt{description}where a fusion center gathers and processes
the measurements collected by all sensors.
\cmt{limitations}Unfortunately, this approach
\begin{myitemize}%
\cmt{comp}gives rise to i) communication bottlenecks, since those nodes
near the fusion center are required to forward data packets from many
sensors, 
\cmt{comm}ii) computational challenges, since the load is
concentrated in the fusion center,
\cmt{robustness}and iii)  vulnerability to attacks or failure of the fusion center.
\end{myitemize}%
\end{myitemize}%
\cmt{description}For these reasons, the decentralized
paradigm, where there is no central processor and all nodes share the
computational load, is oftentimes preferred~\cite{nedic2018decentralized}.
\cmt{strengths}These implementations are therefore scalable,
robust, and balance the communication and processing requirements across nodes.
\end{myitemize}%
\cmt{fast}The present paper capitalizes on the notion of graph
filter~\cite{shuman2013emerging,sandryhaila2014bigdata} to develop
algorithms for computing projections in a decentralized fashion with
an approximately minimal number of iterations.
\end{myitemize}%




\cmt{literature: approaches for subspace proj}
\begin{myitemize}%
\cmt{general inference}To tackle problems involving
decentralized processing, it is common to define a \emph{communication
graph} where each node represents a sensor and there exists an edge
between two nodes if the corresponding sensors can communicate,
e.g. via a radio link.  One could therefore feel inclined to address
the problem at hand using standard inference tools for data defined on
graphs; see
e.g.~\cite{belkin2004regularization,romero2017spacetime,romero2017multikernel,lauritzen1996graphical,isufi2019forecasting,scarselli2008graph}. However,
these methods are based on exploiting a certain relation between the
data and the graph;
e.g. smoothness~\cite{romero2017multikernel}. Therefore, this
framework fundamentally differs from the one at hand since here the
graph only provides information about how the sensors communicate,
i.e., it does not generally provide information about the spatial
field of interest.\footnote{\arev{If the field is smooth over space
and the graph is also smooth over space, meaning that nodes that are
spatially close have a low geodesic distance in the graph, then both
notions may approximately coincide. However, even in this case,
exploiting spatial information would still be more
accurate. Fortunately, as described in this paper, graph signal
processing tools can still be applied to exploit spatial smoothness
rather than graph smoothness. }}
\cmt{decentralized optimization}To obtain a decentralized
    subspace projection algorithm, one could instead adopt a
    decentralized optimization standpoint, e.g. via
\begin{myitemize}%
\cmt{works}
\begin{myitemize}%
    \cmt{DLMS}the \emph{distributed least mean squares} (DLMS)
    method in~\cite{mateos2012drls}, based on the \emph{alternating
    direction method of multipliers} (ADMM)~\cite{giannakis2016decentralized,boyd2011distributed},
    \cmt{DGD}or the \emph{decentralized gradient descent} (DGD)
    method in~\cite{shi2017distributed}, which builds upon gradient descent.
\end{myitemize}%
\cmt{limitations}Although these algorithms can accommodate
    general objective functions, their convergence is only
    asymptotic and can be significantly improved by exploiting the
    structure of the subspace projection problem; see
    Sec.~\ref{sec:sim}. 
\end{myitemize}%
\cmt{general projections}%
\begin{myitemize}%
\cmt{works}%
\begin{myitemize}%
\cmt{barbarossa and insausti}For this reason, a method tailored to computing projections in a
decentralized fashion was proposed in~\cite{barbarossa2009projection},
later extended in~\cite{insausti2012in-network} and~\cite{dilorenzo2020innetwork}, where every node
obtains each iterate by linearly combining its previous iterate with
the previous iterate of its neighbors.  The combination weights are
adjusted to achieve a fast asymptotic convergence.
\end{myitemize}%
\cmt{limitations $\rightarrow$~ asymptotic}The main strength of this 
approach is its simplicity, since each node simply repeats the same
operation over and over. The price to be paid is that convergence is
asymptotic, which means that a large number of data packages need to
be exchanged to attain a prescribed projection accuracy. In addition,
these algorithms can only accommodate a limited set of
topologies~\cite{camaro2013reducing}.
\end{myitemize}%
\cmt{finte-time consensus}A special case of the subspace
projection problem is average consensus, where the signal subspace
comprises the vectors whose entries are all equal.
\begin{myitemize}%
\cmt{works no GSP}For this special case, the
algorithms in~\cite{kibangou2012graph,safavi2015nulling} produce
projections with a finite number of communication rounds. 
\cmt{works GSP}A graph signal
processing~\cite{shuman2013emerging,sandryhaila2014bigdata}
perspective to tackle this special case was adopted
in~\cite{sandryhaila2014finitetime}.
\cmt{Limitations}Unfortunately, these schemes cannot be applied
to the general subspace projection problem.
\end{myitemize}%
\end{myitemize}%
\begin{myitemize}%
\cmt{GSP}%
\begin{myitemize}%
\cmt{works}%
\begin{myitemize}%
\cmt{segarra}A more general framework is proposed
in~\cite{segarra2017graphfilterdesign}, which would allow
implementation of a projection with a graph filter if one were given
a \emph{shift matrix} such that a subset of its eigenvectors spans the
subspace of interest. \cmt{limitations}Unfortunately, this framework
does not include any method to find such a shift matrix in a general case unless the
target subspace is of dimension 1. 
\cmt{coutinho}An even more general setup is presented
in~\cite{coutino2019dgf}, which approximates an arbitrary linear
 transformation at the expense of more complex node computations
 through the notion of ``edge-variant'' graph filters.
\cmt{limitations}%
\begin{myitemize}%
\cmt{Non-convex}Unfortunately, the non-convex~\cite{boyd}
 nature of the optimization problem involved therein
\begin{myitemize}%
\cmt{guarantees}yields no guarantees that a projection filter can be
found even if it exists
\cmt{unpredictable behavior}and may lead to unpredictable
behavior if the network topology is modified and the filter weights
need to be updated.
\cmt{complexity}
\end{myitemize}%
\cmt{higher complexity $\rightarrow$~ Each node needs to store more coefficients}
\cmt{speed}Besides, even in the unlikely event that  the optimization algorithm finds
a global optimum, the resulting filter is not necessarily
implementable in a small number of iterations. 
\end{myitemize}%
\cmt{freq response}Further graph-filter design  schemes abound,
but they typically seek implementing a given frequency
    response \cite{isufi2017autoregressive,tay2014near,narang2012perfect}
    relative to a given shift matrix.
\end{myitemize}%
\end{myitemize}%
\end{myitemize}%


\cmt{contribution}
\begin{myitemize}%
\cmt{gap to fill}To sum up, there is no decentralized algorithm for computing
    general subspace projections in a finite number of iterations. 
\cmt{exact}The present paper fills this gap by suitably
designing a shift matrix and the graph filter coefficients. The sought
filter is of approximately minimal order, which implies that the
number of data exchanges among nodes is approximately minimized. 
\cmt{convex surrogate for the multiplicity}The minimal order is
seen to depend on the  multiplicity of the eigenvalues of the shift
matrix. Since maximizing this multiplicity would lead to a non-convex
problem, a novel convex relaxation technique is developed relying on
a nuclear norm functional of the shift matrix. 
\cmt{convexity}To solve the resulting optimization problem, a
solver based on ADMM is also developed. 
\cmt{approximate}For those scenarios where there exists no
graph filter  that can implement the desired
projection on the given topology, a method is proposed to approximate such a projection
while trading approximation error for filter order.
\cmt{solvers}Another ADMM solver is developed for this case.

\end{myitemize}%

\cmt{conference}\arev{
The conference precursor~\cite{weerasinghe2018fast} of this work
  contains~\thref{prop:disjointevals} and the key steps
  leading to \eqref{prob:exactrelaxed}. Most of the analysis,
  simulations, the approximate projection method, and the ADMM solvers
  are presented here anew. We also published a related subgradient
  method in~\cite{mollaebrahim2018large} but it is not contained in the
  present manuscript. }

\cmt{paper structure}The paper is structured as
follows. Sec.~\ref{prob} formulates the problem, reviews common
applications, and lies some background on graph
filters. Secs.~\ref{sec:exactproj} and~\ref{sec:approximate}
respectively propose methods for exact and approximate projection
implementation. Finally, Sec.~\ref{sec:sim} \arev{presents the simulations
and Sec.~\ref{sec:conclusions} summarizes the main conclusions and provides a discussion. One proof
and the derivations of the ADMM methods are provided in the supplementary
material.}

\emph{Notation}:
\begin{myitemize}%
\cmt{}Symbol $:=$ denotes equality by definition.
\cmt{Set notation}%
\begin{myitemize
\cmt{}For sets $\mathcal{A}$ and $\mathcal{B}$, the cardinality
of $\mathcal{A}$ is denoted as $|\mathcal{A}|$ whereas
\cmt{}$\mathcal{A}\subsetneq\mathcal{B}$ indicates that
$\mathcal{A}$ is a proper subset of $\mathcal{B}$. 
\end{myitemize
\cmt{matrix}
\begin{myitemize
\cmt{}Boldface lowercase (uppercase) letters represent column
vectors (matrices).
\cmt{}The $\ell_n$ norm of vector $\bm
v$ is denoted as $\|\bm v\|_n$.
\cmt{}With $\bm A$ and $\bm B$ matrices of appropriate dimensions, 
$[\bm A;\bm B]$ and $[\bm A,\bm B]$ respectively denote
their vertical and horizontal concatenation, 
 \cmt{}$\bm A^\top$ the transpose of $\bm A$,
 \cmt{}$\mathsf{cols}(\bm A)$ the set
 of the  columns of $\bm A$,
 \cmt{}$\mathop{\rm diag}(\bm A)$ a vector comprising the diagonal
 entries of $\bm A$,
 \cmt{}${\hc{\mathcal{R}}}(\bm A)$ the span of the columns of $\bm A$, 
 \cmt{}$\bm A\otimes\bm B$  the Kronecker product of $\bm A$
 and $\bm B$,
  \cmt{}$\mathsf{evals}(\bm A)$ the set of  eigenvalues of $\bm A$,
 \cmt{}$\lambda_i(\bm A)$ the $i$-th largest eigenvalue of $\bm
 A$,
 \cmt{}$\sigma_{i}(\bm A)$ the
$i$-th largest singular value of $\bm A$,
\cmt{norms}%
\begin{myitemize
 \cmt{}$
 ||\bm A||_2 := \sigma_1(\bm A) $
 the $2$-norm of $\bm A$, and
\cmt{}$||\bm A||_\star:=\sum_i\sigma_i(\bm A) $
 the nuclear norm of $\bm A$.
\end{myitemize
\cmt{}For a subspace $\mathcal{A}$, notation
$\mathcal{A}^\perp$ represents the orthogonal complement.
\end{myitemize
\cmt{Probability}Finally,
\begin{myitemize
\cmt{}$\mathop{\textrm{E}}\nolimits$ denotes expectation and
\cmt{}$\mathcal{N}$  the normal distribution.
\end{myitemize
\end{myitemize}%





\section{Preliminaries} \label{prob}
\subsection{The Subspace Projection Problem}
\label{sec:projformulation}
\cmt{graph}
\begin{myitemize}%
\cmt{def}Let ${\hc{\mathcal{G}}}({\hc{\mathcal{V}}},{\hc{\mathcal{E}}})$ denote a
 graph with vertex set ${\hc{\mathcal{V}}}=\{1,\ldots,\nodenum\}$, where  each vertex 
 corresponds to a sensor or \emph{node}, and edge set
 ${\hc{\mathcal{E}}}\subset {\hc{\mathcal{V}}}^2$.  Let there be an edge
 $(n,n')$ in $ {\hc{\mathcal{E}}}$ if and only if (iff) the nodes $n$ and
 $n'$ can communicate directly, e.g. through their radio
 interface. Thus, it is natural to assume (i) that  ${\hc{\mathcal{E}}}$ contains all
 self loops, i.e., $(n,n)\in{\hc{\mathcal{E}}}~\forall n\in{\hc{\mathcal{V}}}$, and (ii) that ${\hc{\mathcal{G}}}$ is undirected, which means
 that $(n,n')\in {\hc{\mathcal{E}}}$ implies that  $(n',n)\in {\hc{\mathcal{E}}}$.
\cmt{neighborhood}The neighborhood of the $n$-th node is defined as $\neighborhood{\nodeind}=\{n'\mid(n,n')\in{\hc{\mathcal{E}}}\}$.
\end{myitemize}%

\cmt{signal}
\begin{myitemize}%
\cmt{observations}Given  $\obsvec=[\obsentry_1,\ldots,
\obsentry_\nodenum]^\top$, where $ \obsentry_n\in \mathbb{R}$ denotes the observation or measurement
acquired by the $n$-th node,
\cmt{goal}the goal is to estimate the
signal vector ${\hc{\bm \xi}} \in \mathbb{R}^\nodenum$, which quantifies
the phenomenon of interest (e.g. temperature field).
\cmt{model}
\begin{myitemize}%
\cmt{observation model}The latter is
related to $\obsvec$ via
\begin{equation}
\label{eq:obsvecdef}
{\obsvec}={\hc{\bm \xi}}+{{\hc{\bm v}}},
\end{equation}
where  ${\hc{\bm v}} \in \mathbb{R}^\nodenum$ stands for additive noise.
\cmt{subspace}Vector ${\hc{\bm \xi}}$ is known to
lie in a given subspace $\range{{\hc{\bm U}}_{\parallel}}$ of dimension
$\subspacedim < \nodenum$, where the columns of   ${\hc{\bm U}}_{\parallel}
\in \mathbb{R}^{\nodenum \times \subspacedim} $ are assumed orthonormal without loss of
generality (w.l.o.g.). Hence, ${\hc{\bm \xi}}$ can be expressed as $
{\hc{\bm \xi}}={\hc{\bm U}}_{\parallel} {\bm{\coord}} $ for some
${\bm{\coord}}\in \mathbb{R}^\subspacedim$. 
\end{myitemize}%
\end{myitemize}%

\cmt{subspace projection problem}
\begin{myitemize}%
\cmt{projection}The
orthogonal projection of $\obsvec$ onto
$\range{{\hc{\bm U}}_{\parallel}}$, also known as the least-squares
estimate\footnote{Also the \emph{best linear unbiased
estimator}, \emph{minimum variance unbiased estimator},
and \emph{maximum likelihood estimator} \cite{kay1} under appropriate
assumptions. } of ${\hc{\bm \xi}}$, is given
by: \begin{equation}\label{projection}
{\hat{{\hc{\bm \xi}}}}:=
{{\hc{\bm U}}_{\parallel}\uparmat^{\top}}{{\obsvec}}\mathop{=}^{\Delta} {{\hc{\bm P}}\obsvec}\hspace{0.5mm},
\end{equation}
\noindent where ${\hc{\bm P}} \in \mathbb{R}^{\nodenum \times \nodenum} $ is the
projection matrix onto $\range{{\hc{\bm U}}_{\parallel}}$.
\cmt{problem}The subspace projection problem is to find $ \hat{\hc{\bm \xi}}$ given
$\obsvec$ and ${\hc{\bm U}}_{\parallel}$. Vector ${\hat{\hc{\bm \xi}}}$ is expected to be a better
estimate of $ {\hc{\bm \xi}}$ than $\obsvec$ since the noise 
 is annihilated along $\nodenum-\subspacedim$ dimensions.

\end{myitemize}%


\subsection{The Choice of the Basis}
\label{sec:basis} 

\cmt{Overview}This section discusses specific choices of the basis
${\hc{\mathcal{U}}}:=\{{\hc{\bm u}}_1,\ldots, {\hc{\bm u}}_\subspacedim\}$ formed by the
columns of ${\hc{\bm U}}_{\parallel}$ in different application scenarios where a
spatial field needs to be monitored. 
\cmt{Notation}
\begin{myitemize}%
\cmt{space}To this end, let
${\hc{\bm x}}_\nodeind\in \mathbb{R}^{\hc{d}}$ denote the spatial
location of the $\nodeind$-th sensor, where ${\hc{d}}=2$ or
3. Similarly, let ${\hc{\bm x}}\in \mathbb{R}^{\hc{d}}$ denote an arbitrary
location in the area of interest.
\cmt{field}Suppose that the goal is to estimate a spatial field
${\hc{\xi}}:\mathbb{R}^{\hc{d}}\rightarrow\mathbb{R}$ given the
measurements in \eqref{eq:obsvecdef}, where ${\hc{\bm \xi}}:=
[{\hc{\xi}}({\hc{\bm x}}_1),\ldots,{\hc{\xi}}({\hc{\bm x}}_\nodenum)]^\top$.
\end{myitemize}%

\cmt{Parametric model}Oftentimes, the physics of the problem directly
provides a linear parametric expansion for ${\hc{\xi}}$. 
\begin{myitemize}%
\cmt{examples}
\begin{myitemize
\cmt{diffusion}For example,
in the case of a diffusion field, such as a temperature field, one has
\begin{align}
\label{eq:diffusion}
{\hc{\xi}}({\hc{\bm x}})
=&\sum_{i=1}^{\subspacedim} \frac{ \exp\left\{-  \|{\hc{\bm x}} -
\sourcelocvec{i}\|_2^2 /(2 \sigma_i^2)\right\}}{2 \pi \sigma_i^2}
\arev{\tilde{\hc{\alpha}}_i}
\end{align}
\arev{for some coefficients $\tilde{\hc{\alpha}}_i$}, where $\sourcelocvec{i}$ is the location of the $i$-th source and
the parameters $\{\sigma_i^2\}_{i=1}^r$ are related to the diffusivity
of the medium.
\cmt{RSS}In some cases governed by a wave
equation, as occurs in wireless communications (see e.g.~\cite{romero2017spectrummaps}),
 ${\hc{\xi}}$ may admit an expansion in terms of  Cauchy bells~\cite{barbarossa2009projection}: 
\begin{align}
\label{eq:cauchy}
{\hc{\xi}}({\hc{\bm x}}) = \sum_{i=1}^{\subspacedim} \frac{1}{1 + { \|{\hc{\bm x}} -
\sourcelocvec{i}\|_2^2}/{\sigma_i^2}}
\arev{\tilde{\hc{\alpha}}_i}.
\end{align}
\end{myitemize}%
\cmt{Find basis}\arev{To obtain ${\hc{\mathcal{U}}}$, %
evaluate \eqref{eq:diffusion} or \eqref{eq:cauchy} at the
sensor locations and collect the coefficients that multiply each
$\tilde{\hc{\alpha}}_i$ to form the vector
$\tilde {\hc{\bm u}}_i\in \mathbb{R}^\nodenum$, $i=1,\ldots,r$. 
This yields the expansion
${\hc{\bm \xi}}=\sum_{i=1}^\subspacedim \tilde{\hc{\bm u}}_i\tilde{\hc{\alpha}}_i$. Finally, 
orthonormalize  $\{\tilde {\hc{\bm u}}_i\}_{i=1}^{r}$.}

\end{myitemize}%


\cmt{Approximations model}\arev{This approach applies when   ${\hc{\xi}}$
satisfies a parametric expansion as in \eqref{eq:diffusion}
or \eqref{eq:cauchy} and this expansion is known. However, it is
often the case that the form of the expansion is known but it contains
unknown parameters, the form of the expansion is not even known,
or the field does not even admit a linear expansion but it
approximately does. In these situations, one may still pursue a linear
inference approach by capitalizing on some form of smoothness that the
target field exhibits across space.}
\begin{myitemize}%
\cmt{approximations}%
\arev{For instance, ${\hc{\xi}}$} can be approximately bandlimited, which means
that ${\hc{\xi}}$ can be reasonably approximated by a reduced number
$\subspacedim$ of Fourier or \emph{discrete cosine transform} (DCT) basis functions.
\begin{myitemize}%
\cmt{bandlimited}In the latter case, upon letting
${\hc{\bm x}}:=[{\hc{x}}_1,{\hc{x}}_2]^\top$, one can write
\begin{align}
 {\hc{\xi}}({\hc{\bm x}}) \approx& \sum_{i_1=0}^{\subspacedim_1 -1} \sum_{i_2=0}^{\subspacedim_2 -1} {\hc{\alpha}}_{i_1,i_2}  & \\
& \times \cos \left ( \frac{\pi}{{\hc{X}}_1} i_1 \left ( {\hc{x}}_1 + \frac{1}{2} \right )   \right ) \; 
\cos \left ( \frac{\pi}{{\hc{X}}_2} i_2 \left ( {\hc{x}}_2 + \frac{1}{2} \right )   \right ) \nonumber.
\end{align}
Here, ${\hc{X}}_1$ and
 ${\hc{X}}_2$ denote the length  along the 1st and 2nd dimensions of the region where ${\hc{\xi}}$ is
 defined. The vector ${\bm{\coord}}$ defined in
 Sec.~\ref{sec:projformulation} can be recovered by stacking the
 $\subspacedim_1\subspacedim_2$ coefficients $\{{\hc{\alpha}}_{i_1,i_2}\}$,
 whereas ${\hc{\mathcal{U}}}$ can be found as described earlier in this section.

\cmt{Others}Besides Fourier or DCT bases, one may pursue
 approximations based on any other collection of 
 basis functions such as conventional polynomials, discrete prolate
 spheroidal functions, and wavelets.
\end{myitemize}%
\cmt{tradeoff}Note that in any approximation of this kind there
 is a fundamental variance-bias trade-off; see
 e.g.~\cite[Ch. 3.4]{cherkassky2007}. To see this, note that the \emph{signal-to-noise ratio} after
 projection in the model
 \eqref{eq:obsvecdef} is given by
 $\|{\hc{\bm P}}{\hc{\bm \xi}}\|^2_2/\mathop{\textrm{E}}\nolimits[\|{\hc{\bm P}}{\hc{\bm v}}\|_2^2]$. If
 ${\hc{\bm v}}$ has zero mean and  covariance matrix
 $\mathop{\textrm{E}}\nolimits[{\hc{\bm v}} {\hc{\bm v}}^\top]={\hc{\sigma^2}}\bm I_\nodenum$, then 
$\mathop{\textrm{E}}\nolimits[\|{\hc{\bm P}}{\hc{\bm v}}\|_2^2]=\mathop{\textrm{E}}\nolimits[\|{\hc{\bm U}}_{\parallel}\uparmat^\top{\hc{\bm v}}\|_2^2]={\hc{\sigma^2}}\mathop{\rm Tr}[{\hc{\bm P}}]={\hc{\sigma^2}}\subspacedim$. Thus, although a basis with a larger $\subspacedim$ may  capture more signal
 energy $\|{\hc{\bm P}}{\hc{\bm \xi}}\|^2_2$, the power of the noise
 component in $\hat{\hc{\bm \xi}}$ is also
 increased. 





\end{myitemize}%

\subsection{Graph Filters}
\label{sec:graphfilters}

\cmt{overview}This section briefly reviews the notion of graph
 filters~\cite{shuman2013emerging,sandryhaila2014bigdata}, which
 constitute a central part of the proposed algorithms.
\cmt{Graph signal}In this context, vector $\bm
z:=[z_1,\ldots,z_N]^\top$ is referred to as a \emph{graph signal},
which emphasizes the fact that the entry $z_n$ is stored at the 
$n$-th node.

\cmt{graph filter}A graph filter involves two steps, as described
next.
\begin{myitemize}%
\cmt{first: local aggregations}In the first step,
\begin{myitemize}%
\cmt{overview}a finite
sequence of graph signals
$\{\obsvec\itnot{\itind}\}_{\itind=0}^{\itnum}$, where $\obsvec\itnot{\itind}:=[\obsentry_1\itnot{\itind},\ldots,
\obsentry_\nodenum\itnot{\itind}]^\top$, is collaboratively
obtained by the network through a sequence of $\itnum$ \emph{local
data exchange rounds}, or just \emph{local exchanges} for short, where
$\obsvec\itnot{0}:=\obsvec$ is the graph signal to filter.
\cmt{l-th round}At the $\itind$-th round,  each node sends
its $\obsentry\itnodenot{\itind-1}{\nodeind}$ to its neighbors and
computes a linear combination of the entries
$\{\obsentry\itnodenot{\itind-1}{\nodeind'}\}_{\nodeind'\in \neighborhood{\nodeind}}$
that it receives from them. Specifically, the next graph signal is
obtained as $\obsentry\itnodenot{\itind}{\nodeind}
= \sum_{\nodeind'\in \neighborhood{\nodeind}}{\hc{s}}\nodenodenot{\nodeind}{\nodeind'} \obsentry\itnodenot{\itind-1}{\nodeind'},~\nodeind=1,\ldots,\nodenum,$
where  ${\hc{s}}\nodenodenot{\nodeind}{\nodeind'}$ 
is the
 coefficient corresponding to the linear aggregation that takes place
 between nodes ${\nodeind}$ and ${\nodeind'}$.
 \cmt{shift operator}
\begin{myitemize}%
\cmt{matrix form}By letting ${\hc{s}}\nodenodenot{\nodeind}{\nodeind'}=0$
 whenever $\nodeind'\notin\neighborhood{\nodeind}$, one can
 equivalently write  $\obsentry\itnodenot{\itind}{\nodeind}
 = \sum_{\nodeind'=1}^\nodenum{\hc{s}}\nodenodenot{\nodeind}{\nodeind'} \obsentry\itnodenot{\itind-1}{\nodeind'}$
 or, in matrix form, $\obsvec\itnot{\itind}
 = {\hc{\bm S}}\obsvec\itnot{\itind-1}$, where
 ${\hc{\bm S}}\in \mathbb{R}^{\nodenum\times \nodenum}$ is given by
 $({\hc{\bm S}})_{\nodeind,\nodeind'}={\hc{s}}_{\nodeind,\nodeind'},~\nodeind,\nodeind'=1,\ldots,\nodenum$.
\cmt{term}In the graph signal processing literature, 
the matrix ${\hc{\bm S}}$ is usually referred to as  \emph{shift
 matrix}~\cite{segarra2017graphfilterdesign}.
\cmt{def}More generally, an $\nodenum\times\nodenum$ matrix
 ${\hc{\bm S}}$
 is said to be a shift matrix over the graph
 ${\hc{\mathcal{G}}}:=({\hc{\mathcal{V}}},{\hc{\mathcal{E}}})$ if $({\hc{\bm S}})_{\nodeind,\nodeind'}=0$
 for all $(\nodeind,\nodeind')\notin {\hc{\mathcal{E}}}$. The set of all possible shift
 matrices over   ${\hc{\mathcal{G}}}$
 will be denoted as $\shiftmatset$.
\cmt{examples}Examples of  matrices in  $\shiftmatset$ include the adjacency
 and Laplacian matrices of ${\hc{\mathcal{G}}}$~\cite{segarra2017graphfilterdesign}.

  \cmt{operator}Associated with the shift matrix is the \emph{shift
  operator}, defined as the function
  $\obsvec \mapsto {\hc{\bm S}} \obsvec$.
\end{myitemize}%
\cmt{explicit form}Notice that, upon recursively applying  $\obsvec\itnot{\itind}
 = {\hc{\bm S}}\obsvec\itnot{\itind-1}$,  one can write\footnote{Throughout the paper,  $\bm A^0$ for a
 square matrix $\bm A$ denotes the identity matrix of the same size as
 $\bm A$, regardless of whether $\bm A$ is invertible.}
 $\obsvec\itnot{\itind}={\hc{\bm S}}^\itind \obsvec$,
 $\itind=0,\ldots, \itnum$.
\end{myitemize}%


\cmt{second: linear combination}In the second step of the graph filter,
all nodes linearly combine the iterates in the first step.
 Specifically, the following graph signal is computed:
\begin{align}
\label{eq:outputgfilter}
\projobsvec
= \sum_{\itind=0}^{\itnum}\filtercoef_\itind \obsvec\itnot{\itind}
= \sum_{\itind=0}^{\itnum}\filtercoef_\itind
 {\hc{\bm S}}^{\itind}\obsvec
\end{align}
where $\filtercoef_\itind\in\mathbb{R},~{\itind=0},\ldots,{\itnum},$ are
the so-called filter coefficients.

\cmt{Operator}%
\begin{myitemize}%
\cmt{graph filter}The operation in \eqref{eq:outputgfilter} can be
generically expressed as $\obsvec \mapsto {\hc{\bm H}}\obsvec$, where
\begin{equation} \label{filterH}
{{\hc{\bm H}}}:=\sum_{\itind=0}^{\itnum}\filtercoef_\itind{{\hc{\bm S}}}^{\itind}, 
\end{equation}
and is commonly referred to as an order-$\itnum$ \emph{graph filter}.
\cmt{cayley-hamilton}An important implication of  the
Cayley-Hamilton Theorem \cite{horn} is that for any order-$\itnum$
graph filter ${\hc{\bm H}}$ with $\itnum\geq \nodenum$, there exists an
order-($\nodenum-1$) graph filter ${\hc{\bm H}}'$ with shift matrix ${\hc{\bm S}}$
and coefficients $\filtercoef_\itind'$ such that
${\hc{\bm H}}={\hc{\bm H}}'$. This establishes an upper bound on the order
and, therefore, the number of local
exchanges required to apply a graph filter. Thus,   one can
assume w.l.o.g. that
$\itnum\leq\nodenum-1$.
\end{myitemize}%
\end{myitemize}%

\subsection{Asymptotic Decentralized Projections}
\label{sec:barbarossa}

\begin{myitemize}%
\cmt{description}A decentralized scheme for subspace projection
was proposed in~\cite{barbarossa2009projection}. There, a matrix
${\hc{\bm S}}$ is found such that (i) ${\hc{\bm S}}\in \shiftmatset$ and
(ii) $\lim_{\itind\rightarrow\infty} {\hc{\bm S}}^\itind
= {\hc{\bm P}}$. Then, the nodes compute the sequence
$\{\obsvec\itnot{\itind}, \itind=0, 1,\ldots\}$, where
$\obsvec\itnot{\itind}={\hc{\bm S}}^\itind \obsvec$. This
constitutes the infinite counterpart of the first step in a graph
filter; cf.  Sec.~\ref{sec:graphfilters}. Due to (ii), it follows that
$\lim_{\itind\rightarrow\infty}  \obsvec\itnot{\itind}
= {\hc{\bm P}}\obsvec$, as desired.
\cmt{strengths}The main strength of this method is its
simplicity, since each node just needs to store one coefficient for
each neighbor and the same operation is repeated over and
over. 
\cmt{limitations}A limitation is that 
the number of local exchanges required to attain a target error
$||\obsvec\itnot{\itind}- {\hc{\bm P}} \obsvec||$ is generally high since this
approach only provides asymptotic convergence. Furthermore, the set of
graphs for which (i) and (ii) can be simultaneously satisfied is
considerably limited; see Sec.~\ref{sec:sim}
and \cite{camaro2013reducing}.



\end{myitemize}%




\section{Exact Projection Filters} \label{sec:exactproj}

\cmt{overview}This section proposes an algorithm to find
a graph filter that yields a subspace projection in an approximately
minimal number of iterations. To this end, Secs.~\ref{sec:minorder}
and~\ref{sec:polfeas}
formalize the problem and 
characterize the set of feasible shift matrices for a given
$\mathcal{E}$ and ${\hc{\bm U}}_{\parallel}$. Subsequent sections introduce an
optimization methodology  to approximately minimize the order of the
filter, i.e. the number of communication steps needed to obtain the
projection via graph filtering.

\subsection{Minimum-order Projection Filters}
\label{sec:minorder}

\cmt{feasibility problem}
\begin{myitemize}%
\cmt{goal}To solve the subspace projection problem formulated
in Sec.~\ref{sec:projformulation} with a graph filter,  one could think of
finding a shift matrix ${\hc{\bm S}}\in \shiftmatset$ and a set of
coefficients $\{\filtercoef_\itind\}_{\itind=0}^{\itnum}$ such that
${\hc{\bm P}} \obsvec
= \sum_{\itind=0}^{\itnum}\filtercoef_\itind {\hc{\bm S}}^{\itind}\obsvec$
for all $\obsvec\in\mathbb{R}^\nodenum$ or, equivalently, such that
${\hc{\bm P}}
= \sum_{\itind=0}^{\itnum}\filtercoef_\itind {\hc{\bm S}}^{\itind}$. Since
${\hc{\bm P}}$ is symmetric, it will be assumed that ${\hc{\bm S}}$ is also
symmetric.
\cmt{poly feasible}To assist in this quest, consider the following
definition:
\begin{mydef}
\thlabel{def:polfeas}
Let ${\hc{\bm Q}}\in\mathbb{R}^{\nodenum\times\nodenum}$ be an arbitrary
(not necessarily a projection) matrix.  A symmetric matrix (not
necessarily a shift matrix) ${\hc{\bm S}}\in{\mathbb{R}}^{\nodenum\times{\nodenum}}$
is polynomially feasible to implement the operator
$\obsvec \mapsto{\hc{\bm Q}}\obsvec$ if there exist $\itnum$ and
$\filtercoefvec:=[\filtercoef_0,\ldots,\filtercoef_\itnum]^\top$
such that $\sum_{\itind=0}^{\itnum}\filtercoef_\itind{\hc{\bm S}}^{\itind}={\hc{\bm Q}}$.
\end{mydef}
For a given ${\hc{\bm Q}}$, the set of all polynomially feasible matrices
${\hc{\bm S}}\in \mathbb{R}^{\nodenum\times \nodenum}$  will be denoted as
$\polfeasmatset{{\hc{\bm Q}}}$.
Except for the tip, this set is a cone since $ {\hc{\bm S}}\in\polfeasmatset{{\hc{\bm Q}}}
$ implies
${\hc{\kappa}} {\hc{\bm S}}\in\polfeasmatset{{\hc{\bm Q}}} $ for all ${\hc{\kappa}}\neq 0$.

\cmt{coefficients straightfwd if shift known}Note that given a matrix in
$\polfeasmatset{{\hc{\bm Q}}}$, it is straightforward to obtain 
$\filtercoefvec$
such that $\sum_{\itind=0}^{\itnum}\filtercoef_{\itind}{\hc{\bm S}}^{l}={\hc{\bm Q}}$; see e.g.~\cite{segarra2017graphfilterdesign} and 
Sec.~\ref{sec:polfeas}.
\cmt{formulation}One may then consider the following \emph{feasibility problem}:
\begin{myitemize}%
\cmt{given}Given a graph ${\hc{\mathcal{G}}}:=({\hc{\mathcal{V}}},{\hc{\mathcal{E}}})$ and a
projection matrix ${\hc{\bm P}}$, 
\cmt{find}find a polynomially feasible shift matrix, i.e.,
find ${\hc{\bm S}} \in \shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}$.
\end{myitemize}%
\cmt{feasibility}When this problem admits a solution, we say
that \emph{there exists an exact projection filter} for implementing
${\hc{\bm P}}$ on ${\hc{\mathcal{G}}}$. Whether this is the case depends on ${\hc{\bm P}}$
and ${\hc{\mathcal{G}}}$. For example, if ${\hc{\mathcal{G}}}$ is too sparse, then ${\hc{\bm P}}$
will not be computable as a graph filter. In the extreme case where
${\hc{\mathcal{G}}}$ is fully disconnected, then the only computable projection is
${\hc{\bm P}} = {\hc{\bm I}}_\nodenum$. Conversely, when ${\hc{\mathcal{G}}}$ is fully
connected, then all projections can be computed as a graph
filter. Furthermore, given  that $ \shiftmatset $ is a
subspace and that $ {\hc{\bm S}}\in\polfeasmatset{{\hc{\bm Q}}} $ implies
${\hc{\kappa}} {\hc{\bm S}}\in\polfeasmatset{{\hc{\bm Q}}} $ for all
${\hc{\kappa}}\neq 0$, it is easy to see that such a feasibility problem
either has no solution or has infinitely many.

\end{myitemize}%
\cmt{fast filter}
\begin{myitemize}%
\cmt{informally}When feasible shifts exist, it is reasonable to
seek the ${\hc{\bm S}}$ that minimizes the number of local exchanges
$\itnum$.
\cmt{formally}For arbitrary  matrices ${\hc{\bm S}}$ and ${\hc{\bm Q}}$, define $\orderfun{{\hc{\bm Q}}}({\hc{\bm S}})$ 
as the minimum $\itnum$ such that ${\hc{\bm Q}}
= \sum_{\itind=0}^{\itnum}\filtercoef_\itind {\hc{\bm S}}^{\itind}$ for
some $\{\filtercoef_\itind\}_{\itind=0}^\itnum$. In view of the bound
dictated by the Cayley-Hamilton Theorem (Sec.~\ref{sec:graphfilters}),   $\orderfun{{\hc{\bm Q}}}({\hc{\bm S}})$ can
be viewed as a function
$\orderfun{{\hc{\bm Q}}}:\polfeasmatset{{\hc{\bm Q}}}\rightarrow\{0,\ldots,\nodenum-1\}$.
\cmt{formulation}Given a projection matrix ${\hc{\bm P}}$, the
problem of finding the shift matrix associated with the minimum-order
filter can therefore be  formulated as: 

\begin{align}
\tag{P1}\label{prob:exact}~~\mathop{\text{minimize}}_{{\hc{\bm S}}}~~&\orderfun{{\hc{\bm P}}}({\hc{\bm S}})\\
 \mathop{\text{s.t.}}~~& {\hc{\bm S}} \in \shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}.
 \nonumber
\end{align}
\tagsleft@false 
\end{myitemize}%
\cmt{approximate subspace proj}Conversely, when
$\shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}=\emptyset$, there exists
no graph filter capable of implementing ${\hc{\bm P}}$. For those cases,
Sec.~\ref{sec:approximate} describes how to find a graph filter that
approximates ${\hc{\bm P}}$.






\subsection{Polynomially Feasible Matrices}
\label{sec:polfeas}

\cmt{overview} To assist in solving \eqref{prob:exact}, this section presents an
algebraic characterization of the set $\polfeasmatset{{\hc{\bm P}}}$ of
polynomially feasible matrices. Recall that the matrices in this set
need not be shift matrices, that is, they need not satisfy the
topology constraints determined by the edge set $\mathcal{E}$.

\cmt{Prefeasible shifts}
\begin{myitemize}%
\cmt{condition}

\begin{mylemma}\thlabel{prop:prefeasible}
Let ${\hc{\bm U}}_{\parallel}\in\mathbb{R}^{\nodenum\times\subspacedim}$ with orthonormal
columns be given and let ${\hc{\bm P}} ={\hc{\bm U}}_{\parallel}\uparmat^{\top}$.  Let
also ${\hc{\bm U}}_{\perp}\in\mathbb{R}^{\nodenum\times{\nodenum-\subspacedim}}$ with
orthonormal columns satisfy
${{\hc{\mathcal{R}}}({\hc{\bm U}}_{\perp})={\hc{\mathcal{R}}}^{\perp}({\hc{\bm U}}_{\parallel})}$. If
${\hc{\bm S}}\in\polfeasmatset{{\hc{\bm P}}}$, then there exist symmetric
matrices $ {\hc{\bm F}}_{\parallel}\in\mathbb{R}^{{\subspacedim}\times{\subspacedim}}$
and
$ {\hc{\bm F}}_{\perp}\in\mathbb{R}^{\nodenum-\subspacedim\times{\nodenum-\subspacedim}}$
such that:
\begin{align}\label{eq:def:cpref}
{\hc{\bm S}}=\begin {bmatrix}{\hc{\bm U}}_{\parallel}& {\hc{\bm U}}_{\perp}\end {bmatrix}\left[	\begin{matrix} {\hc{\bm F}}_{\parallel}& \bm{0}\\ \bm{0}&{\hc{\bm F}}_{\perp}\end{matrix} \right]\begin {bmatrix}{\hc{\bm U}}_{\parallel}^{\top}\\ {\hc{\bm U}}_{\perp}^{\top}\end {bmatrix}.
\end{align} 
\end{mylemma}
\begin{proof}
See Appendix~\ref{proof:prefeasible}.
\end{proof}

\cmt{indep. choice basis}Note that matrices  ${\hc{\bm F}}_{\parallel}$
and ${\hc{\bm F}}_{\perp}$ satisfying  \eqref{eq:def:cpref}
exist regardless\footnote{Note that the algorithms in this paper
produce the same filters for all matrices ${\hc{\bm U}}_{\parallel}$ that span a given
signal subspace $\range{{\hc{\bm U}}_{\parallel}}$; likewise for ${\hc{\bm U}}_{\perp}$. Thus, the obtained
filters only depend on the signal subspace and not on the specific
choice of the basis. This property is what bypasses the difficulty
encountered in~\cite[eq. (20)]{segarra2017graphfilterdesign}, which
will typically be infeasible for graphs with more than $\nodenum$
missing edges.  } of the
choice of ${\hc{\bm U}}_{\parallel}$ and ${\hc{\bm U}}_{\perp}$ as long as the columns of ${\hc{\bm U}}_{\parallel}$ and ${\hc{\bm U}}_{\perp}$ 
respectively constitute an orthonormal basis for the signal subspace  $\range{{\hc{\bm U}}_{\parallel}}$  and its
orthogonal complement ${\hc{\mathcal{R}}}^{\perp}({\hc{\bm U}}_{\parallel})$. 
\cmt{converse}As seen later, the converse of \thref{prop:prefeasible} does not hold. 

\cmt{implications: eigenstructure}To understand the implications of \thref{prop:prefeasible},
\begin{myitemize
\cmt{orth decomposition}rewrite  \eqref{eq:def:cpref} as
\begin{align}\label{eq:def:prefparperp}
{\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^{\top}
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^{\top} = {{\hc{\bm S}}}_{\parallel} + {{\hc{\bm S}}}_{\perp}, 
\end{align}
where ${{\hc{\bm S}}}_{\parallel}:= {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^\top$ and
${{\hc{\bm S}}}_{\perp} := {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^\top$ are 
symmetric matrices whose column spans are respectively contained in the signal
subspace and its orthogonal complement. Thus, they clearly satisfy
${{\hc{\bm S}}}_{\parallel}^\top{{\hc{\bm S}}}_{\perp}={{\hc{\bm S}}}_{\perp}^\top{{\hc{\bm S}}}_{\parallel}=\bm 0$.
\cmt{not nec. shifts}Note that even if ${\hc{\bm S}}$ is a shift
matrix, i.e. ${\hc{\bm S}}\in\shiftmatset$, matrices
${{\hc{\bm S}}}_{\parallel}$ and ${{\hc{\bm S}}}_{\perp}$ may not be shift matrices. 
\cmt{eigendecomposition}
Consider now the eigendecompositions
${\hc{\bm F}}_{\parallel}={\hc{\bm Q}}_{\parallel}{\hc{\bm \Lambda}}_{\parallel}{\hc{\bm Q}}_{\parallel}^\top,{\hc{\bm F}}_{\perp}={\hc{\bm Q}}_{\perp}{\hc{\bm \Lambda}}_{\perp}{\hc{\bm Q}}_{\perp}^\top$ 
for orthogonal $ {\hc{\bm Q}}_{\parallel} \in \mathbb{R}^{\subspacedim \times
\subspacedim}, {\hc{\bm Q}}_{\perp}\in \mathbb{R}^{\nodenum-\subspacedim \times \nodenum-\subspacedim}$ and diagonal ${\hc{\bm \Lambda}}_{\parallel} \in \mathbb{R}^{\subspacedim \times \subspacedim},{\hc{\bm \Lambda}}_{\perp} \in \mathbb{R}^{\nodenum-\subspacedim \times\nodenum-\subspacedim}$.
Then, \eqref{eq:def:cpref} can be rewritten as:
\begin{align}\label{eq:def:prefevd}
{\hc{\bm S}}=\begin
{bmatrix}{\hc{\bm U}}_{\parallel}{\hc{\bm Q}}_{\parallel}& {\hc{\bm U}}_{\perp}{\hc{\bm Q}}_{\perp}\end
{bmatrix}
\left[	\begin{matrix} {\hc{\bm \Lambda}}_{\parallel}& \bm{0}\\ \bm{0}&{\hc{\bm \Lambda}}_{\perp}\end{matrix} \right]
\begin {bmatrix}{\hc{\bm Q}}_{\parallel}^\top{\hc{\bm U}}_{\parallel}^{\top}\\ {\hc{\bm Q}}_{\perp}^\top{\hc{\bm U}}_{\perp}^{\top}\end {bmatrix}.
\end{align} 
This is clearly an eigenvalue decomposition of ${\hc{\bm S}}$. It further
shows that
$\mathrm{evals}({\hc{\bm S}})=\mathrm{evals}({\hc{\bm F}}_{\parallel})\cup\mathrm{evals}({\hc{\bm F}}_{\perp})$. In
view of \eqref{eq:def:prefevd},
\thref{prop:prefeasible} establishes that  any 
${\hc{\bm S}}\in \polfeasmatset{{\hc{\bm P}}}$
has exactly $\subspacedim$ orthogonal eigenvectors (the
columns of ${\hc{\bm U}}_{\parallel}{\hc{\bm Q}}_{\parallel}$) in the
signal subspace and $\nodenum-\subspacedim$ (the
columns of ${\hc{\bm U}}_{\perp}{\hc{\bm Q}}_{\perp}$) in its orthogonal complement.

\end{myitemize}%

\cmt{definition}The following definition builds
upon  \thref{prop:prefeasible} to  introduce 
a necessary condition  for feasibility of a shift matrix 
that will prove instrumental in subsequent sections. 
\begin{mydef}
\thlabel{def:prefeas}
Let ${\hc{\bm U}}_{\parallel}$ and ${\hc{\bm U}}_{\perp}$ be given. 
If  ${\hc{\bm S}}\in\mathbb{R}^{\nodenum\times{\nodenum}}$ is such that it
satisfies \eqref{eq:def:cpref} for some symmetric
${\hc{\bm F}}_{\parallel}\in{\mathbb{R}}^{{\subspacedim}\times{\subspacedim}}$ and
${\hc{\bm F}}_{\perp}\in{\mathbb{R}}^{\nodenum-\subspacedim\times{\nodenum-\subspacedim}}$, then ${\hc{\bm S}}$
is said to be pre-feasible.
\end{mydef}
\cmt{indep choice bases}Note again that this definition is independent
of the choice of  ${\hc{\bm U}}_{\parallel}$ and ${\hc{\bm U}}_{\perp}$ so long as their columns
respectively form a basis for the signal subspace and its
orthogonal complement.
\cmt{notation}Given
${\hc{\bm P}}\in \mathbb{R}^{\nodenum\times\nodenum}$, the set of all
pre-feasible matrices in $\mathbb{R}^{\nodenum\times\nodenum}$ will be
denoted as $\prefeasmatset{{\hc{\bm P}}}$.
\end{myitemize}%

\cmt{Feasible shifts}%
\begin{myitemize
\cmt{feas $\subset$ prefeas}Observe that all
matrices in  $\polfeasmatset{{\hc{\bm P}}}$ are also 
in $\prefeasmatset{{\hc{\bm P}}}$.
\cmt{prefeas $\subsetneq$ feas}However, 
not all matrices in $\prefeasmatset{{\hc{\bm P}}}$ are in
$\polfeasmatset{{\hc{\bm P}}}$. Trivial examples include ${\hc{\bm S}}
= {\hc{\bm I}}_\nodenum$ (recall that $\subspacedim<\nodenum$) and
${\hc{\bm S}}=\bm{0}$. The rest of this section will characterize the
matrices in $\prefeasmatset{{\hc{\bm P}}}$ that are in
$\polfeasmatset{{\hc{\bm P}}}$.
\cmt{which prefeasible mats are feasible}In particular, it will be seen that any pre-feasible
matrix ${\hc{\bm S}}$ where ${\hc{\bm \Lambda}}_{\parallel}$ and ${\hc{\bm \Lambda}}_{\perp}$ share at
least an
eigenvalue is not polynomially feasible.  To this end, note
from \thref{def:polfeas} and \eqref{eq:def:prefevd} that any
pre-feasible matrix must satisfy
\begin{align}
\label{eq:polfeasprefeas}
{\hc{\bm P}}= {\hc{\bm U}}_{\parallel}{\hc{\bm Q}}_{\parallel} \Big[\sum_{\itind=0}^{\itnum}\filtercoef_\itind {\hc{\bm \Lambda}}_{\parallel}^{l}\Big]{\hc{\bm Q}}_{\parallel}^{\top} {\hc{\bm U}}_{\parallel}^{\top} +{\hc{\bm U}}_{\perp} {\hc{\bm Q}}_{\perp} \Big[\sum_{\itind=0}^{\itnum}\filtercoef_\itind {\hc{\bm \Lambda}}_{\perp}^{l}\Big]{\hc{\bm Q}}_{\perp}^{\top} {\hc{\bm U}}_{\perp}^{\top}
\end{align}
for some  $\{\filtercoef_\itind\}_{\itind=0}^{\itnum}$  to be
polynomially feasible.  Multiplying both sides of \eqref{eq:polfeasprefeas} on the left by 
${\hc{\bm U}}_{\parallel}^\top$ and on the right by
${\hc{\bm U}}_{\parallel}$, it follows that ${\hc{\bm Q}}_{\parallel}\big[{\sum_{\itind=0}^{\itnum}\filtercoef_\itind\ {\hc{\bm \Lambda}}_{\parallel}^{\itind}}\big]{\hc{\bm Q}}_{\parallel}^\top = {\hc{\bm I}}_\subspacedim$ or, equivalently,
$\sum_{\itind=0}^{\itnum}\filtercoef_\itind\ {\hc{\bm \Lambda}}_{\parallel}^{\itind}={\hc{\bm I}}_\subspacedim$. Likewise, multiplying both sides of  \eqref{eq:polfeasprefeas} on the left by ${\hc{\bm U}}_{\perp}^\top$ and on the right by ${\hc{\bm U}}_{\perp}$, it follows that
${\sum_{\itind=0}^{\itnum}\filtercoef_\itind\ {\hc{\bm \Lambda}}_{\perp}^{\itind}=\bm
0}$. Arranging these two conditions in matrix form yields:
 \begin{equation} \label{eq:vander}
  \begin{bmatrix}{\hc{\bm 1}}_\subspacedim \\  \bm 0_{\nodenum-\subspacedim}\end{bmatrix} = 
\begin{bmatrix}
    1       & \lambda_{1}  &\dots & \lambda_{1}^{\itnum} \\
    \vdots  & \vdots      &\ddots & \vdots \\
    1       & \lambda_{\subspacedim} &\dots
    & \lambda_{\subspacedim}^{\itnum} \\
    1       & \lambda_{\subspacedim+1} &\dots & \lambda_{\subspacedim+1}^{\itnum} \\ 
    \vdots & \vdots  & \ddots & \vdots \\
    1       & \lambda_{\nodenum}  &\dots & \lambda_{\nodenum}^{\itnum}
\end{bmatrix}
 \begin{bmatrix}
  \filtercoef_0 \\ \filtercoef_1\\ \vdots\\ c_{\itnum}
 \end{bmatrix},
 \end{equation}
\noindent where $\lambda_1,\ldots,\lambda_\nodenum $ are such that
  ${\hc{\bm \Lambda}}_{\parallel}\triangleq \text{diag}\{\lambda_1,\ldots,\lambda_\subspacedim \}$
  and
  ${\hc{\bm \Lambda}}_{\perp}\hspace{-2mm}\triangleq \text{diag}\{\lambda_{\subspacedim+1},\ldots,\lambda_\nodenum \}.$
  Vandermonde systems such as  \eqref{eq:vander}  frequently arise when designing graph
  filters; see
  e.g.~\cite{segarra2017graphfilterdesign,sandryhaila2014finitetime}. With
  the appropriate definitions, it can be expressed in matrix
  form as:
\begin{equation} \label{eq:vander2}
\boldsymbol \lambda_{\hc{\bm P}} = {\hc{\bm \Psi}} \filtercoefvec ,
\end{equation}
\noindent which provides a means to obtain the coefficients
  $\{\filtercoef_\itind\}_{\itind=0}^{\itnum}$ when \eqref{eq:vander2} admits a solution.
  
\begin{myitemize}%
\cmt{}To understand when the latter is the case, assume w.l.o.g. that $\itnum=\nodenum-1$ since the existence
of a solution to \eqref{eq:vander} for some $\itnum$ implies its
existence for $\itnum=\nodenum-1$.
\cmt{}Since ${\hc{\bm \Psi}}$ is a square Vandermonde matrix, any two rows
corresponding to distinct eigenvalues are linearly independent. 
\cmt{}Looking at the left-hand side of \eqref{eq:vander}, it is easy to
see that the system \eqref{eq:vander} admits a solution iff
${\hc{\bm \Lambda}}_{\parallel}$ and ${\hc{\bm \Lambda}}_{\perp}$ do not share eigenvalues.
\end{myitemize}%



\cmt{statement}This conclusion can be combined
with \thref{prop:prefeasible} as follows:
\begin{mytheorem}
\thlabel{prop:disjointevals}
Let ${\hc{\bm S}}\in\mathbb{R}^{\nodenum\times{\nodenum}}$ be
symmetric. Then, ${\hc{\bm S}}\in\polfeasmatset{{\hc{\bm P}}}$ iff both
the following two conditions hold:

\begin{align}
\tag{C1}&{\hc{\bm S}}\in\prefeasmatset{{\hc{\bm P}}}\text{, i.e., it
satisfies \eqref{eq:def:cpref}} \\&\text{ for some symmetric  }{\hc{\bm F}}_{\perp}\text{ and
  }{\hc{\bm F}}_{\parallel},\nonumber\\ \tag{C2} & \mathrm{evals}({\hc{\bm F}}_{\parallel})\cap\mathrm{evals}({\hc{\bm F}}_{\perp})=\emptyset.
\end{align}
\tagsleft@false

\end{mytheorem}

\end{myitemize}%


\subsection{Filter Order Minimization}
\label{sec:filterorderminimization}


\cmt{overview}\thref{prop:disjointevals} implies
  that \eqref{prob:exact} can be reformulated as the minimization of
$\orderfun{{\hc{\bm P}}}({\hc{\bm S}})$ subject to
${\hc{\bm S}} \in \shiftmatset$, (C1), and (C2). This section and the
next develop a reformulation more amenable to application of a
numerical solver.

\cmt{objective}The first step is to express the objective function
$\orderfun{{\hc{\bm P}}}({\hc{\bm S}})$  more explicitly. To that end,
consider the following result:
\begin{myitemize}%
\cmt{recall def $\orderfun{{\hc{\bm P}}}$}
\cmt{Order = \# distinct evals}
\begin{mylemma}
\thlabel{prop:minimalorder}If ${\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^\top
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^\top\in \polfeasmatset{{\hc{\bm P}}}$, then
\begin{align}
\label{eq:minimalorder}
\orderfun{{\hc{\bm P}}}({\hc{\bm S}}) \leq \noevalsfun({\hc{\bm F}}_{\parallel})
+\noevalsfun({\hc{\bm F}}_{\perp}) -1,
\end{align}
where $\noevalsfun(\cdot)$ is the number of distinct eigenvalues of
its argument.
\end{mylemma}

\begin{IEEEproof}
Computing $\orderfun{{\hc{\bm P}}}({\hc{\bm S}})$ for a given ${\hc{\bm S}}$
amounts to determining the minimum $\itnum$ for
which \eqref{eq:vander}, or equivalently its compact
version \eqref{eq:vander2}, admits a solution.  Since
${\hc{\bm S}} \in \polfeasmatset{{\hc{\bm P}}}$ by hypothesis, one has
that \eqref{eq:vander2} is satisfied for at least one value of
$\itnum$. Let $\itnum_0$ denote the smallest value for
which \eqref{eq:vander2} holds. Since ${\hc{\bm \Psi}}$ is
Vandermonde, it has at most $\noevalsfun({\hc{\bm S}})$ linearly
independent rows, which implies that 
$\itnum_0+1\leq \noevalsfun({\hc{\bm S}})$. The proof is completed by noting
from \thref{prop:disjointevals} that $\noevalsfun({\hc{\bm S}})
= \noevalsfun({\hc{\bm F}}_{\parallel})
+\noevalsfun({\hc{\bm F}}_{\perp})$.
\end{IEEEproof}
\cmt{tightness}In practice, the bound in \eqref{eq:minimalorder}
will hold with equality unless in degenerate cases. For example, when
$\nodenum=4$, $\subspacedim=2$, $\lambda_1=-\lambda_2$ and
$\lambda_3=-\lambda_4\neq |\lambda_1|$, it can be seen that
$\orderfun{{\hc{\bm P}}}({\hc{\bm S}})=2$ whereas $ \noevalsfun({\hc{\bm F}}_{\parallel})
+\noevalsfun({\hc{\bm F}}_{\perp})-1=3$. However, in general, such an ${\hc{\bm S}}$
will only be in $\shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}$ if
${\hc{\bm P}}$ and ${\hc{\mathcal{G}}}$ are jointly selected to achieve this end,
which will not occur in a practical application.
\cmt{interpretation}In words, 
\thref{prop:minimalorder} implies that one may seek the shift matrix
of an approximately minimal-order projection filter as the matrix with
the smallest number of distinct eigenvalues among all matrices in the
feasible set of \eqref{prob:exact}.
\end{myitemize}%

\cmt{Feasible set}Next, this feasible set is
rewritten more explicitly.
\begin{myitemize}%
\cmt{split into two constraints}To this end, split the constraint
${\hc{\bm S}} \in \shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}$ 
into the two constraints
${\hc{\bm S}} \in \shiftmatset$ and
${\hc{\bm S}} \in   \polfeasmatset{{\hc{\bm P}}}$.
\begin{myitemize}%
\cmt{topology constraint}Regarding ${\hc{\bm S}} \in \shiftmatset$,
\begin{myitemize}%
\cmt{def shift}recall from the definition of $\shiftmatset$
in Sec.~\ref{sec:graphfilters} that
${\hc{\bm S}} \in \shiftmatset$ iff  $({\hc{\bm S}})_{\nodeind,\nodeind'}=0$
for all $(\nodeind,\nodeind') $ such that
 $(\nodeind,\nodeind')\notin {\hc{\mathcal{E}}}$.
 \cmt{symmetry}Because any 
 ${\hc{\bm S}}$ that is feasible for  \eqref{prob:exact} has to be 
 in $\polfeasmatset{{\hc{\bm P}}}$ and all the 
 matrices in this set are symmetric, any feasible ${\hc{\bm S}}$ is necessarily
 symmetric. Thus, one can just require that  $({\hc{\bm S}})_{\nodeind,\nodeind'}=0$
 only  for those $(\nodeind,\nodeind') $ such that
 $(\nodeind,\nodeind')\notin {\hc{\mathcal{E}}}$  and 
 $\nodeind<\nodeind'$.
 \cmt{canonicalbasis}With $\canbasisvec_\nodeind$ the
 $\nodeind$-th column of the identity matrix ${\hc{\bm I}}_\nodenum$, it
 follows that   $({\hc{\bm S}})_{\nodeind,\nodeind'}=
\canbasisvec_\nodeind^\top{\hc{\bm S}}\canbasisvec_{\nodeind'}=
(\canbasisvec_{\nodeind'}\otimes\canbasisvec_\nodeind)^\top\mathop{\rm vec}({\hc{\bm S}})
$. Thus, the constraint ${\hc{\bm S}} \in \shiftmatset$ can be expressed
as ${\hc{\bm{W}}}\mathop{\rm vec}({\hc{\bm S}}) = \bm 0$, where
${\hc{\bm{W}}}$ is a matrix whose rows are given by the
vectors
$\{(\canbasisvec_{\nodeind'}\otimes\canbasisvec_\nodeind)^\top, 
\forall(\nodeind,\nodeind')$ such that
$(\nodeind,\nodeind')\notin {\hc{\mathcal{E}}}$ and
 $\nodeind<\nodeind'\}$. As expected, the fewer edges in the graph,
 the more rows ${\hc{\bm{W}}}$ has and, consequently, the
 smaller the feasible set. In the extreme case of a fully disconnected
 graph, only the diagonal matrices satisfy
 ${\hc{\bm S}} \in \shiftmatset$ (recall that  ${\hc{\mathcal{E}}}$ contains all self-loops;
 cf. Sec.~\ref{sec:projformulation}). 
 
\end{myitemize}%

\cmt{polfeas constraint}On the other hand, the constraint
${\hc{\bm S}} \in   \polfeasmatset{{\hc{\bm P}}}$ can be easily expressed invoking
\thref{prop:disjointevals} and introducing two auxiliary optimization
variables ${\hc{\bm F}}_{\parallel}$ and ${\hc{\bm F}}_{\perp}$. 

\end{myitemize}%


\cmt{resulting problem}In view of these observations and \thref{prop:minimalorder},
problem \eqref{prob:exact} becomes 

\begin{align}
\tag{P1'}\label{prob:exactwithconstraints}
\begin{aligned}
~~\mathop{\text{minimize}}_{{\hc{\bm S}},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp}}~~&
 \noevalsfun({\hc{\bm F}}_{\parallel})
+\noevalsfun({\hc{\bm F}}_{\perp})\\
 \mathop{\text{s.t.}}~~&  {\hc{\bm{W}}}\mathop{\rm vec}({\hc{\bm S}}) = \bm 0\\
& {\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^\top
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^\top\\
&{\hc{\bm F}}_{\parallel}={\hc{\bm F}}_{\parallel}^\top,~{\hc{\bm F}}_{\perp}={\hc{\bm F}}_{\perp}^\top\\
&\lambda_\nodeind({\hc{\bm F}}_{\parallel})\neq\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})~\forall\nodeind,\nodeind'
\end{aligned}
\end{align}
\tagsleft@false 
for an arbitrary choice of
${\hc{\bm U}}_{\perp}\in \mathbb{R}^{\nodenum\times \nodenum-\subspacedim}$ with
orthonormal columns spanning ${\hc{\mathcal{R}}}^\perp({\hc{\bm U}}_{\parallel})$.
\end{myitemize}%

\cmt{two modifications}Two further modifications are in order.
\begin{myitemize}%
\cmt{invariance to scalings}First, note that
\eqref{prob:exactwithconstraints} is invariant to scalings in the
sense that if $({\hc{\bm S}},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp})$ is feasible, then
 $({\hc{\kappa}}{\hc{\bm S}},{\hc{\kappa}}{\hc{\bm F}}_{\parallel},{\hc{\kappa}}{\hc{\bm F}}_{\perp})$ is also
 feasible and attains the same objective value $\forall {\hc{\kappa}}\neq
 0$. For this reason, the constraint $\mathop{\rm Tr}({\hc{\bm F}}_{\parallel}) = \subspacedim$ will be
 introduced w.l.o.g.  to eliminate this ambiguity.

\cmt{feasible set is not closed}Second, the feasible set
 of \eqref{prob:exactwithconstraints} is not a closed set due to
 the constraint
 $\lambda_\nodeind({\hc{\bm F}}_{\parallel})\neq\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})~\forall\nodeind,\nodeind'$,
 implying that the optimum may not be attained through an iterative algorithm. In practice, this
 constraint holds so long as the eigenvalues of ${\hc{\bm F}}_{\parallel}$ differ from
 those of ${\hc{\bm F}}_{\perp}$, even if there is an eigenvalue of ${\hc{\bm F}}_{\parallel}$
 arbitrarily close to an eigenvalue of ${\hc{\bm F}}_{\perp}$. But in the latter case,
 the numerical conditioning of \eqref{eq:vander} would be poor, implying
 that the target projection cannot be implemented as a graph filter
 using finite-precision arithmetic. Thus, the aforementioned
 constraint should be replaced with another one that ensures that (i)
 the feasible set is closed, and (ii) the eigenvalues of ${\hc{\bm F}}_{\parallel}$
 are sufficiently different from those of ${\hc{\bm F}}_{\perp}$. One natural possibility
 is
 $|\lambda_\nodeind({\hc{\bm F}}_{\parallel})-\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})|\geq {\hc{\epsilon}}~\forall\nodeind,\nodeind'$,
 where ${\hc{\epsilon}}>0$ is a user-selected parameter. Unfortunately, this
 constraint is not convex, but an effective  relaxation will be
 presented in the next section.

\end{myitemize}%

\cmt{resulting problem} To sum up, the optimization problem 
 to be solved is:

\begin{align}
\tag{P1''}\label{prob:exactclosedfs}
\begin{aligned}
~~\mathop{\text{minimize}}_{{\hc{\bm S}},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp}}~~&
 \noevalsfun({\hc{\bm F}}_{\parallel})
+\noevalsfun({\hc{\bm F}}_{\perp})\\
 \mathop{\text{s.t.}}~~&  {\hc{\bm{W}}}\mathop{\rm vec}({\hc{\bm S}}) = \bm 0,~\mathop{\rm Tr}({\hc{\bm F}}_{\parallel}) = \subspacedim\\
& {\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^\top
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^\top\\
&{\hc{\bm F}}_{\parallel}={\hc{\bm F}}_{\parallel}^\top,~{\hc{\bm F}}_{\perp}={\hc{\bm F}}_{\perp}^\top\\
&|\lambda_\nodeind({\hc{\bm F}}_{\parallel})-\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})|\geq {\hc{\epsilon}}~\forall\nodeind,\nodeind'.
\end{aligned}
\end{align}
\tagsleft@false 



\cmt{relation to barbarossa's feasible set}Remarkably, the set of
 topologies for which \eqref{prob:exactclosedfs} is feasible is
strictly larger than the set of topologies for which the
method proposed in~\cite{barbarossa2009projection} (and reviewed in
Sec.~\ref{sec:barbarossa}) can be applied. It can be easily seen that
any feasible matrix for the problem in~\cite{barbarossa2009projection}
must be in $ \shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}$ and,
additionally, $\nodenum-\subspacedim$ of its eigenvalues must be equal
to 1 whereas the rest must be less than 1. It follows that the
feasible set therein is strictly contained in the feasible set of \eqref{prob:exactwithconstraints}.




\subsection{Exact Projection Filters via Convex Relaxation}
\label{sec:convex}

\cmt{overview}Problem \eqref{prob:exactclosedfs}  is non-convex because
\begin{myitemize}%
\cmt{}(i) the objective function is not convex and
\cmt{}(ii) the last constraint, which does not define a convex set.
\end{myitemize}%
This section proposes a convex problem to approximate the solution
to \eqref{prob:exactclosedfs}.

\cmt{objective}To address (i), recall  that $ \noevalsfun({\hc{\bm F}}_{\parallel})$
equals the number of distinct eigenvalues of ${\hc{\bm F}}_{\parallel}$ and 
\begin{myitemize}%
\cmt{transforming objective}note that the larger  $ \noevalsfun({\hc{\bm F}}_{\parallel})$, the
larger the number of non-zero elements of the vector
$[\lambda_1-\lambda_2,\lambda_1-\lambda_3,\ldots,\lambda_1-\lambda_r,\lambda_2-\lambda_3,\ldots,\lambda_{r-1}-\lambda_r]^\top$. A
similar observation applies to $ \noevalsfun({\hc{\bm F}}_{\perp})$.  This
suggests replacing the objective in \eqref{prob:exactclosedfs} with
\begin{align}
\label{eq:surrobj}
 \sum_{\nodeind=1}^\subspacedim
\sum_{\nodeind'=1}^\subspacedim||\lambda_\nodeind({\hc{\bm F}}_{\parallel})-\lambda_{\nodeind'}({\hc{\bm F}}_{\parallel})||_0
+\\
\sum_{\nodeind=1}^{\nodenum-\subspacedim}
\sum_{\nodeind'=1}^{\nodenum-\subspacedim}||\lambda_\nodeind({\hc{\bm F}}_{\perp})-\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})||_0\nonumber,
\end{align}
where $\|{\hc{\bm x}}\|_0$ is the so-called \emph{zero norm} or number of
non-zero elements of vector ${\hc{\bm x}}$.
\cmt{obj. surrogate}A
typical convex surrogate of the zero norm is the
$\itind_1$-norm \cite{candes2008introduction}. However, just replacing
$||\cdot ||_0$ in \eqref{eq:surrobj} with an $l_1$-norm would still
give rise to a non-convex objective since it involves the functions
$\lambda_\nodeind(\cdot)$. 
\cmt{Kronecker}One of the key ideas in this paper is to apply
the following result:
\begin{mylemma}
\thlabel{prop:evalsnucleardif}
Let ${\hc{\bm A}}$ be an $\nodenum\times \nodenum$ matrix with eigenvalues $\lambda_1, \lambda_2, \ldots, \lambda_\nodenum$. Then,
 \begin{align}
 {{\left\| {\hc{\bm A}}\otimes {{\hc{\bm I}}_{\nodenum}}-{{\hc{\bm I}}_{\nodenum}}\otimes {\hc{\bm A}}\right\|}_{*}}=\sum\limits_{\nodeind=1}^{\nodenum}{\sum\limits_{\nodeind'=1}^{\nodenum}{\left| {{\lambda }_{\nodeind}}-{{\lambda }_{\nodeind'}} \right|}}.
 \end{align}
\end{mylemma}
\begin{IEEEproof}
See Appendix~\ref{proof:evalsnucleardif}.
\end{IEEEproof}
Applying \thref{prop:evalsnucleardif} to \eqref{eq:surrobj} suggests
replacing the objective in \eqref{prob:exactclosedfs} with $ ||{\hc{\bm F}}_{\parallel}\otimes {\hc{\bm I}}_{\subspacedim}
-{\hc{\bm I}}_{\subspacedim} \otimes {\hc{\bm F}}_{\parallel} ||_\star + 
||{\hc{\bm F}}_{\perp}\otimes {\hc{\bm I}}_{\nodenum-\subspacedim}-{\hc{\bm I}}_{\nodenum-\subspacedim} \otimes
{\hc{\bm F}}_{\perp}||_\star$. \arev{The implications of this relaxation are
further discussed in Appendix~\ref{sec:tightness}.}


\end{myitemize}%


\cmt{feasible set}Regarding (ii), note that 
\begin{myitemize}%
\cmt{relax feasible set}constraint
$|\lambda_\nodeind({\hc{\bm F}}_{\parallel})-\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})|\geq {\hc{\epsilon}}~\forall\nodeind,\nodeind'$
renders the feasible set non-convex due to (ii-1) the functions
$\lambda_\nodeind(\cdot)$ and (ii-2) the absolute value.
\begin{myitemize}%
\cmt{relax lambdas}To deal with (ii-1), a sensible approach is to relax the constraints $|\lambda_\nodeind({\hc{\bm F}}_{\parallel})-\lambda_{\nodeind'}({\hc{\bm F}}_{\perp})|\geq {\hc{\epsilon}}~\forall\nodeind,\nodeind'$
into a single constraint which requires only that the average of the
eigenvalues of ${\hc{\bm F}}_{\parallel}$ must differ sufficiently from the average of
the eigenvalues of ${\hc{\bm F}}_{\perp}$, that is
\begin{align}
\left|\frac{1}{\subspacedim}\sum_{\nodeind=1}^\subspacedim\lambda_\nodeind({\hc{\bm F}}_{\parallel})-
\frac{1}{\nodenum-\subspacedim}\sum_{\nodeind=1}^{\nodenum-\subspacedim}
\lambda_{\nodeind}({\hc{\bm F}}_{\perp})\right|\geq {\hc{\epsilon}}\nonumber
\end{align}
for some user-selected ${\hc{\epsilon}}>0$. Clearly, this constraint is
equivalent to    $\left|\mathop{\rm Tr}({\hc{\bm F}}_{\parallel})/{\subspacedim}-
\mathop{\rm Tr}({\hc{\bm F}}_{\perp})/({\nodenum-\subspacedim})\right|\geq {\hc{\epsilon}}$
and, using the constraint $\mathop{\rm Tr}({\hc{\bm F}}_{\perp})=\subspacedim$ introduced
in \eqref{prob:exactclosedfs},  becomes equivalent to 
  $\left|1-
\mathop{\rm Tr}({\hc{\bm F}}_{\perp})/({\nodenum-\subspacedim})\right|\geq {\hc{\epsilon}}$.
\cmt{absolute value}To deal with (ii-2), one can now use the definition of absolute
value to conclude that this constraint is satisfied if either
$\mathop{\rm Tr}({\hc{\bm F}}_{\perp}) \leq (1-{\hc{\epsilon}})({\nodenum-\subspacedim})$ or
$
\mathop{\rm Tr}({\hc{\bm F}}_{\perp})\geq (1+{\hc{\epsilon}})({\nodenum-\subspacedim})$. Both of
these inequalities are affine and therefore convex, however the OR
condition renders the feasible set 
non-convex. However, a solution can be readily found by 
first solving the problem with the constraint $\mathop{\rm Tr}({\hc{\bm F}}_{\perp}) \leq
(1-{\hc{\epsilon}})({\nodenum-\subspacedim})$, then solving it 
with the constraint
$\mathop{\rm Tr}({\hc{\bm F}}_{\perp})\geq (1+{\hc{\epsilon}})({\nodenum-\subspacedim})$ instead, and
finally comparing the objective values achieved at the optimum of both these
sub-problems.
\cmt{ineq $\rightarrow$~ eq}To simplify the task of solving the resulting
problem, note that the solution in both subproblems will
satisfy the corresponding constraint with equality. To see this, note
that if the inequality constraint is removed, the optimum of the
resulting convex problem becomes
$({\hc{\bm S}},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp})=({\hc{\bm I}}_\nodenum,{\hc{\bm I}}_\subspacedim,
{\hc{\bm I}}_{\nodenum-\subspacedim})$. Due to convexity and since this optimum does not satisfy the removed constraint, such a
constraint will become necessarily active when introduced.

\cmt{plus minus}In view of these observations, the last
constraint in \eqref{prob:exactclosedfs} will be replaced with  $ 
\mathop{\rm Tr}({\hc{\bm F}}_{\perp})= (1\pm{\hc{\epsilon}})({\nodenum-\subspacedim})$, where the
$\pm$ sign indicates that both subproblems must be solved separately.
\cmt{Regularizer}Except for degenerate cases, one expects that
these relaxed constraints will suffice to ensure that ${\hc{\bm F}}_{\perp}$ and
${\hc{\bm F}}_{\parallel}$ do not share eigenvalues. However, although this was indeed
observed in our experiments, sometimes an eigenvalue of ${\hc{\bm F}}_{\perp}$ may
become sufficiently close to an eigenvalue of ${\hc{\bm F}}_{\parallel}$ in such a way
that the Vandermonde system \eqref{eq:vander} becomes poorly
conditioned. To alleviate this situation, a penalty proportional to
$||{\hc{\bm F}}_{\perp} ||_\text{F}^{2}$ will be added to the objective to
``push'' the eigenvalues of ${\hc{\bm F}}_{\perp}$ towards zero.

\end{myitemize}%


\cmt{Convex problem}With all these considerations, the resulting relaxed problem becomes:

\begin{align}
\nonumber~~\mathop{\text{minimize}}_{{\hc{\bm S}},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp}}~~&
 {\hc{ \eta}}_{\parallel}||{\hc{\bm F}}_{\parallel}\otimes {\hc{\bm I}}_\subspacedim
-{\hc{\bm I}}_\subspacedim \otimes {\hc{\bm F}}_{\parallel} ||_\star \\&+ {\hc{ \eta}}_{\perp}
||{\hc{\bm F}}_{\perp}\otimes {\hc{\bm I}}_{\nodenum-\subspacedim}-{\hc{\bm I}}_{\nodenum-\subspacedim}\otimes
{\hc{\bm F}}_{\perp}||_\star+||{\hc{\bm F}}_{\perp} ||_\text{F}^{2}\nonumber\\
 \mathop{\text{s.t.}}~~&  {\hc{\bm{W}}}\mathop{\rm vec}({\hc{\bm S}}) = \bm 0,\nonumber\\
\tag{P1-R}\label{prob:exactrelaxed}& {\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^\top
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^\top,\\
&{\hc{\bm F}}_{\parallel}={\hc{\bm F}}_{\parallel}^\top,~{\hc{\bm F}}_{\perp}={\hc{\bm F}}_{\perp}^\top\nonumber,\\
&\mathop{\rm Tr}({\hc{\bm F}}_{\perp})= (1\pm{\hc{\epsilon}})({\nodenum-\subspacedim}),
\mathop{\rm Tr}({\hc{\bm F}}_{\parallel}) = \subspacedim\nonumber,
\end{align}
\tagsleft@false 
where ${\hc{ \eta}}_{\parallel}>0$ and ${\hc{ \eta}}_{\perp}>0$ are user-specified parameters
to control the relative weight of each term in the objective. Having
two parameters rather than just one that scales the term $||{\hc{\bm F}}_{\perp}
||_\text{F}^{2}$ empowers  the user with more flexibility to alleviate
possible numerical issues, as described earlier.
\cmt{relaxation}\arev{A study of the ability
of \eqref{prob:exactrelaxed} to approximate the solution
of \eqref{prob:exact} is presented in Appendix~\ref{sec:tightness}.
}

\cmt{solver}An iterative solver for \eqref{prob:exactrelaxed} is
proposed in Appendix~\ref{sec:exactsolver} in the supplementary
 material based on ADMM, therefore inheriting its solid convergence
 guarantees.


\end{myitemize}%


\section{Approximate Projection Filters}
\label{sec:approximate}

\cmt{Overview}When the graph is too sparse for \eqref{prob:exact} to
admit a solution, one may  instead seek  a low-order graph filter that
approximates ${\hc{\bm P}}$ reasonably well. Before presenting an
optimization criterion to obtain the corresponding shift matrix, the next
section analyzes in which situations there exists an exact projection
filter.



\subsection{Feasibility of Exact Projection Filters}
\label{sec:feasibility}

\cmt{Relate topology an ${\hc{\mathcal{U}}}$ to existence}The following result provides a necessary condition for the existence
of an exact projection shift matrix. 
\begin{myitemize}%
\cmt{reduced edge set}
Let $\reducededgeset:=\{
(\nodeind,\nodeind')\in{\hc{\mathcal{E}}}:~\nodeind\leq\nodeind'
\}=\{(\nodeind_1,\nodeind'_1),\ldots,(\nodeind_\reducededgenum,\nodeind'_\reducededgenum)\}$
denote the reduced edge set, where each of the $\reducededgenum$
undirected  edges shows up only once.
\cmt{single edge shift}Consider also the symmetric shift matrices
$\singleedgeshiftmat_i:= \canbasisvec_{n_{i}}\canbasisvec_{n_{i}'}^{\top}+\canbasisvec_{n_{i}'}\canbasisvec_{n_{i}},~i=1,\ldots,\reducededgenum$,  where only a single
edge is used.
\cmt{basis for $\shiftmatset$}Clearly, all feasible shift matrices are
linear combinations of these matrices.
\end{myitemize}%
\begin{mytheorem}
\thlabel{prop:necessaryfeasibility}
Let
$\singleedgeshiftmat :=[\mathop{\rm vec}(\singleedgeshiftmat_1),\ldots,\mathop{\rm vec}(\singleedgeshiftmat_\reducededgenum)]$. It
holds that
\begin{align}
\label{eq:necessaryfeasibility}
\dim(\shiftmatset\cap \prefeasmatset{{\hc{\bm P}}})
=\reducededgenum -
 \mathrm{rank}((
{\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat ).
\end{align}

\end{mytheorem}
\begin{IEEEproof}
See Appendix~\ref{proof:necessaryfeasibility}. 
\end{IEEEproof}





\cmt{Implications}
\begin{myitemize}%
\cmt{rank condition}Due to the presence of the scaled identity matrices
 in $\shiftmatset\cap \prefeasmatset{{\hc{\bm P}}}$ (see Sec.~\ref{sec:polfeas}), 
a necessary condition for a feasible projection filter to exist,
i.e. 
$\shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}\neq \emptyset$, is that
$\dim(\shiftmatset\cap \prefeasmatset{{\hc{\bm P}}})>1$. From \eqref{eq:necessaryfeasibility},
this condition becomes $ \mathrm{rank}((
{\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat )\leq
\reducededgenum-2$.
For future reference, this is summarized as follows:
\begin{mycorollary}
\thlabel{prop:necessaryfullfeas}
If there exists an exact projection filter, i.e. $\shiftmatset \cap \polfeasmatset{{\hc{\bm P}}}\neq \emptyset$, 
then
\begin{align}
\label{eq:necessaryfeasibility2}
 \mathrm{rank}((
{\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat )\leq
\reducededgenum-2.
\end{align}
\end{mycorollary}
Note that this provides a condition that can be easily checked before
attempting to solve \eqref{prob:exactrelaxed}.
\cmt{guiding relation}It also provides a guideline
on the minimum number of edges required to exactly implement a
projection with a given $\subspacedim$.  Since
$({\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat \in{\mathbb{R}}^{\subspacedim(\nodenum-\subspacedim)\times\reducededgenum}$,
it follows that $\mathrm{rank}((
{\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat)\leq
\min(\subspacedim(\nodenum-\subspacedim),\reducededgenum)
$. To obtain a reference for the  number of
edges required for the existence of exact projection filters, suppose that $(
{\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat $ is full rank,
as often is the case. Then, \eqref{eq:necessaryfeasibility2}
becomes
$\min(\subspacedim(\nodenum-\subspacedim),\reducededgenum)
\leq
\reducededgenum-2$, which is equivalent to 
$\subspacedim(\nodenum-\subspacedim)\leq
\reducededgenum-2$. Thus, the number of edges required for an exact
projection filter to exist is in the order of
$\subspacedim(\nodenum-\subspacedim)$.
This agrees with intuition because the difficulty to implement a
projection must depend equally on $\subspacedim$ and
$\nodenum-\subspacedim$. This follows by noting that (i)
due to the term
$\filtercoef_0\bm I_\nodenum$ in \eqref{filterH},
 implementing ${\hc{\bm P}}$
with a graph filter is
equally difficult as implementing $\bm I_{\nodenum}-{\hc{\bm P}}$; and
(ii)  ${\hc{\bm P}}$ has rank  $\subspacedim$ whereas $\bm
I_{\nodenum}-{\hc{\bm P}}$ has rank
$\nodenum-\subspacedim$. \arev{In any case, note that this is just a
guideline on the number of edges. The existence of an exact projection
filter will not only depend on the number of edges but also on their
 locations.}


\cmt{when POLF is empty}Note that \thref{prop:necessaryfullfeas}
provides a necessary condition for the existence of an exact
projection filter. Obtaining a sufficient condition, in turn, is much
more challenging. However, one expects that 
$\shiftmatset\cap\polfeasmatset{{\hc{\bm P}}}\neq\emptyset$ whenever
$\dim(\shiftmatset\cap \prefeasmatset{{\hc{\bm P}}})>1$. In fact we
conjecture that, given a topology and generating ${\hc{\bm U}}_{\parallel}$ by
orthonormalizing $\subspacedim$ vectors in $\mathbb{R}^\nodenum$ drawn
i.i.d. from a continuous distribution, the event where
$\dim(\shiftmatset\cap \prefeasmatset{{\hc{\bm P}}})>1$ and
$\shiftmatset\cap\polfeasmatset{{\hc{\bm P}}}=\emptyset$ has zero
probability. In practice, however, whether
$\shiftmatset\cap\polfeasmatset{{\hc{\bm P}}}$ is empty or not is not as
relevant as it may seem.
\cmt{not as important}To see this, remember
from \thref{prop:disjointevals} that
$\shiftmatset\cap\polfeasmatset{{\hc{\bm P}}}\neq \emptyset$ when there
exists a pre-feasible shift matrix for which the eigenvalues of ${\hc{\bm F}}_{\parallel}$ differ
from those of ${\hc{\bm F}}_{\perp}$. However, even when this is the case, if one
of the eigenvalues of ${\hc{\bm F}}_{\parallel}$ lies too close to an eigenvalue of
${\hc{\bm F}}_{\perp}$, the corresponding filter cannot be  implemented  due to poor conditioning
of \eqref{eq:vander}. In short, even when an exact projection filter
may exist from a theoretical perspective, such a filter may not be
implementable in practice. Instead, having a sufficiently large
$\dim(\shiftmatset\cap \polfeasmatset{{\hc{\bm P}}})$ seems more important
since that would allow the user to choose a shift that
leads to a good conditioning of \eqref{eq:vander}. In this sense,
\thref{prop:necessaryfeasibility} suggests that the \emph{margin}
$\reducededgenum-
 \mathrm{rank}((
{\hc{\bm U}}_{\parallel}^{\top}\otimes{\hc{\bm U}}_{\perp}^{\top})\singleedgeshiftmat )$ would be
 a reasonable indicator of how easy it  is to obtain an exact
 projection filter.



\end{myitemize}%

\subsection{Approximate Projection Criterion}
\label{sec:approximatecriterion}

\cmt{overview}The method proposed in Sec.~\ref{sec:convex} to obtain an exact
 projection filter relies on solving \eqref{prob:exactrelaxed}. Using
 this problem as a starting point, the present section develops an
 optimization criterion that yields low-order graph
 filters that approximate a given projection operator.


\cmt{rewrite \eqref{prob:exactrelaxed} in terms of ${\hc{\bm S}}$}
\begin{myitemize}%
\cmt{change vars}To this end,  note from the second
 constraint in \eqref{prob:exactrelaxed} and the orthogonality of
 ${\hc{\bm U}}_{\parallel}$ and ${\hc{\bm U}}_{\perp}$ that 
${\hc{\bm F}}_{\parallel}= {\hc{\bm U}}_{\parallel}^{\top}{\hc{\bm S}}{\hc{\bm U}}_{\parallel}$ and
${\hc{\bm F}}_{\perp}= {\hc{\bm U}}_{\perp}^{\top}{\hc{\bm S}}{\hc{\bm U}}_{\perp}$.
\cmt{symmetry constraints}It is also easy to see that 
  ${\hc{\bm S}}$ is symmetric iff ${\hc{\bm F}}_{\parallel}$ and ${\hc{\bm F}}_{\perp}$ are
  symmetric\arev{; cf. Appendix~\ref{proof:prefeasible}}.
\cmt{separate eigenspaces}Additionally, one can also easily
  prove that any  ${\hc{\bm S}}$ can be expressed as ${\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}{\hc{\bm U}}_{\parallel}^\top
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}{\hc{\bm U}}_{\perp}^\top$ for some ${\hc{\bm F}}_{\parallel}$ and
${\hc{\bm F}}_{\perp}$ iff ${\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel} = \bm 0$. In
words, the latter condition states that each eigenvector of
${\hc{\bm S}}$ must be either in the signal subspace or in its orthogonal
complement. 


\cmt{problem}Then, problem \eqref{prob:exactrelaxed} can be equivalently
expressed in terms of ${\hc{\bm S}}$ as

\begin{align}
~~\mathop{\text{minimize}}_{{\hc{\bm S}}}~~&
 {\hc{ \eta}}_{\parallel}||{\hc{\bm U}}_{\parallel}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel}\otimes {\hc{\bm I}}_\subspacedim
-{\hc{\bm I}}_\subspacedim \otimes {\hc{\bm U}}_{\parallel}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel} ||_\star +\nonumber\\& {\hc{ \eta}}_{\perp}
||{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp}\otimes {\hc{\bm I}}_{\nodenum-\subspacedim}-{\hc{\bm I}}_{\nodenum-\subspacedim} \otimes
{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp}||_\star+\nonumber\\&||{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp} ||_\text{F}^{2}
\tag{P1-R'}\label{prob:exactrelaxedS}
\\\nonumber
 \mathop{\text{s.t.}}~~
&{\hc{\bm{W}}}\mathop{\rm vec}({\hc{\bm S}}) = \bm 0,~{\hc{\bm S}} = {\hc{\bm S}}^\top,~{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel} = \bm 0,\\\nonumber
&\mathop{\rm Tr}({\hc{\bm U}}_{\parallel}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel}) = \subspacedim, \\&\mathop{\rm Tr}({\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp})= (1\pm{\hc{\epsilon}})({\nodenum-\subspacedim}).\nonumber
\end{align}
\tagsleft@false
\end{myitemize}%

\cmt{now relax}The constraint that renders \eqref{prob:exactrelaxedS}
infeasible when an exact projection filter does not exist is precisely
${\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel} = \bm 0$. This suggests finding
a shift for an approximate projection filter by solving

\begin{align}
~~\mathop{\text{minimize}}_{{\hc{\bm S}}}~~&
 {\hc{ \eta}}_{\parallel}||{\hc{\bm U}}_{\parallel}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel}\otimes {\hc{\bm I}}_\subspacedim
-{\hc{\bm I}}_\subspacedim \otimes {\hc{\bm U}}_{\parallel}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel} ||_\star +\nonumber\\& {\hc{ \eta}}_{\perp}
||{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp}\otimes {\hc{\bm I}}_{\nodenum-\subspacedim}-{\hc{\bm I}}_{\nodenum-\subspacedim} \otimes
{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp}||_\star+\nonumber\\&||{\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp} ||_\text{F}^{2}+\regpar|| {\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel}||_\text{F}^2
\tag{P2-R}\label{prob:approxrelaxed}
\\\nonumber
 \mathop{\text{s.t.}}~~
&{\hc{\bm{W}}}\mathop{\rm vec}({\hc{\bm S}}) = \bm 0,~{\hc{\bm S}} = {\hc{\bm S}}^\top,\\\nonumber
&\mathop{\rm Tr}({\hc{\bm U}}_{\parallel}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel}) = \subspacedim, \\&\mathop{\rm Tr}({\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\perp})= (1\pm{\hc{\epsilon}})({\nodenum-\subspacedim}),
\nonumber
\end{align}
\tagsleft@false
where the separation of eigenspaces \emph{enforced} by the constraint
${\hc{\bm U}}_{\perp}^\top{\hc{\bm S}}{\hc{\bm U}}_{\parallel} = \bm 0$
in \eqref{prob:exactrelaxedS} is now simply \emph{promoted} through
the last term in the objective. The parameter $\regpar>0$ is selected by
the user to balance the trade off between approximation error and
filter order. 
\cmt{solver}\arev{An iterative solver for \eqref{prob:approxrelaxed} is
proposed in Appendix~\ref{sec:approxsolver} of the supplementary
 material based on ADMM, therefore inheriting its solid convergence
 guarantees.}







\section{Numerical Experiments}
\label{sec:sim}



     \begin{figure}[t]
     \centering
         \includegraphics[width=0.4\textwidth]{./figures/rank-one-eps-converted-to.pdf}
   \caption{\small{NMPE as a function of the number of
  communications performed per node  for the Erdős–Rényi
         networks ($ \nodenum=30, \subspacedim=1$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$).}}
         \label{rank-one}
     \end{figure}
     

     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/NMPE2-eps-converted-to.pdf}
         \caption{\small{NMPE as a function of the number of
  communications performed per node  for the Erdős–Rényi
         networks ($ \nodenum=20$, $p_{\text{miss}}=0.6$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$).}}
\label{NMPE-2}
     \end{figure}

     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/fig-1d-eps-converted-to.pdf}
         \caption{\small{NMPE as a function of the number of
  communications performed per node for WSN ($ d_{max}=.6, \rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$).}}
\label{NMPE-3}
     \end{figure}


     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/NMSE1-new-eps-converted-to.pdf}
         \caption{\small{NMSE as a function of the number of
  communications performed per node ($ \nodenum=20$, $\subspacedim=3, \beta=5,$ Erdős–Rényi graph with $p_{\text{miss}}=0.6$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$).}}
         \label{NMSE-4}
     \end{figure}

     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/NMSE2-new-eps-converted-to.pdf}
         \caption{\small{NMSE as a function of the number of
  communications performed per node ($ \nodenum=30, \subspacedim=5$, $\beta=5,$ Erdős–Rényi graph with $p_{\text{miss}}=0.7$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$).}}
         \label{NMSE-5}
     \end{figure}

     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/NMSE3-eps-converted-to.pdf}
     \caption{\small{NMSE as a function of the subspace dimension (Erdős–Rényi network
     with $p_{\text{miss}}=.7$, $\beta=5,  L=20$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$).}}
        \label{NMSE-6}
     \end{figure}




\begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/MNDE-eps-converted-to.pdf}
    \caption{\small{MNDE for the exact projection method as a function of $p_{\text{miss}}$ ($\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$, $\tau=0.005$, Erdős–Rényi network).}}
         \label{MNDE}
     \end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.4\textwidth]{./figures/need-eps-converted-to.pdf}
       \caption{\small{Average number of  local exchanges required to satisfy
       $\text{NMPE}<\gamma_\text{target}=.1$ vs. the WSN graph  parameter
       $d_{\text{max}}$  ($\rho=0.1$, $I_{\text{max}}=1000$, ${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1$.)}}
\label{WSN-1}
\end{figure}


\begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/NMPEapp1-eps-converted-to.pdf}
         \caption{\small{NMPE as a function of $p_{\text{miss}}$ for the
         Erdős–Rényi networks ($\itnum=\nodenum-1$, $\nodenum=20, \subspacedim=3$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.2, \regpar=10$).}}
         \label{approx-new1}
\end{figure}


     \begin{figure}
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/approximation1-eps-converted-to.pdf}
     \caption{\small{NMPE as a function of the number of
 communications performed per node  for the Erdős–Rényi networks ($\nodenum=25, \subspacedim=3$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.2, \regpar=10$).}}
        \label{NMPE-approx-er}
     \end{figure}

     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/approxim2-eps-converted-to.pdf}
    \caption{\small{NMPE as a function of the number of
  communications performed per node for the WSN networks ($\nodenum=35, \subspacedim=5$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.2, \regpar=10$).}}
         \label{NMPE-approx-wsn}
     \end{figure}




     \begin{figure}[t]
         \centering
         \includegraphics[width=0.4\textwidth]{./figures/extra-eps-converted-to.pdf}
    \caption{\small{NMPE as a function of $d_{\text{max}}$ for the WSN networks ($\nodenum=20$, $\rho=0.1$, $I_{\text{max}}=1000$,
${\hc{ \eta}}_{\perp}=0.9, {\hc{ \eta}}_{\parallel}=0.1, {\hc{\epsilon}}=0.1, \regpar=10$).}}
         \label{extra}
     \end{figure}






\cmt{overview}This section validates the performance of the proposed algorithms by
means of numerical experiments.\footnote{The code necessary to reproduce the
experiments is available at \texttt{github.com/uiano/fast\_decentralized\_projections}.}

\cmt{Data Generation}The data generation process is as follows.
\begin{myitemize}%
\cmt{subspace}Matrix ${\hc{\bm U}}_{\parallel}$ is obtained by orthonormalizing an
$\nodenum \times {\subspacedim}$ matrix with i.i.d. standard Gaussian
entries. \arev{This basis is assumed exactly known, meaning that the
error due to the selection of the basis is disregarded; see
Sec.~\ref{sec:basis}.}
\cmt{observations}To generate the observations
$\obsvec={\hc{\bm \xi}}+{\hc{\bm v}}$, the noise is drawn as
${\hc{\bm v}}\sim\mathcal{N}(\bm 0_\nodenum,\bm I_\nodenum)$ and the signal
as
${\hc{\bm \xi}}={\hc{\bm U}}_{\parallel}{\bm{\coord}}$, where  ${\bm{\coord}}$ is obtained as 
${\bm{\coord}}=\sqrt{\beta\nodenum/\subspacedim}~{\bm{\coord}}_0$
with ${\bm{\coord}}_0\sim\mathcal{N}(\bm 0_\subspacedim,\bm
I_\subspacedim)$ and with  $\beta:= \mathop{\textrm{E}}\nolimits{\|{\hc{\bm \xi}}\|_2^2}/\mathop{\textrm{E}}\nolimits{\|{\hc{\bm v}}\|_2^2}$ the \emph{signal-to-noise ratio} (SNR).

\end{myitemize}%



\cmt{Graph Models}The following kinds of networks are generated:
\begin{myitemize
\cmt{Erdos-Renyi}(i) Erdős–Rényi  graphs, where the presence of each undirected edge is an
i.i.d. Bernoulli($1-p_{\text{miss}}$) random variable with
$p_{\text{miss}}$ the missing edge probability. 
\cmt{WSN}(ii) wireless sensor network (WSN)  graphs,  generated  by deploying the nodes uniformly at random
over a square area of unit side length and connecting them with an
edge if the internode distance is smaller than $d_{\text{max}}$.
\end{myitemize}%
\cmt{ensure connectivity}If the generated graph in (i) and (ii) is not connected, then additional edges are 
introduced to ensure connectivity. More specifically, if the graph has $c$
components, $c-1$ eges are added between nodes chosen at random from
each component.

\cmt{The compared methods}%
\begin{myitemize
\cmt{least squares}Among the compared methods, those for
decentralized optimization 
obtain the projection by solving the least squares problem
$\argmin_{{\bm{\coord}}}\|\obsvec-{\hc{\bm U}}_{\parallel}{\bm{\coord}}\|^2$. This includes
    \begin{myitemize
    \cmt{DLMS}(ii) the \emph{distributed least mean squares} (DLMS)
    method in~\cite{mateos2012drls} with augmented Lagrangian
    parameter $\dlmsauglagpar$ and step size
    $ \dlmsstepsize$, which builds upon ADMM,
    \cmt{DGD}and (ii) the \emph{decentralized gradient descent} (DGD)
    method in~\cite{shi2017distributed} with step size
    $\mu_{\text{DGD}}$, which is based on gradient descent.
    \end{myitemize
    \cmt{Shifts}Other methods iteratively apply a shift matrix:
    \begin{myitemize
    \cmt{barbarossa}(iii) The \emph{gossiping} scheme
    in~\cite{barbarossa2009projection} obtains the shift matrix
    ${\hc{\bm S}}$ that provides fastest convergence of
    ${\hc{\bm S}}^\itind$ to ${\hc{\bm P}}$ as $l\rightarrow \infty$
    according to a certain criterion. Then, the nodes collaboratively
    obtain $\obsvec\itnot{\itind}={\hc{\bm H}}_\itind \obsvec$ with
    ${\hc{\bm H}}_\itind={\hc{\bm S}}^\itind$ at the $\itind$-th exchange
    round; see also Sec.~\ref{sec:graphfilters}. Signal
    $\obsvec\itnot{\itind}$ asymptotically converges to the desired
    projection.
    \cmt{Segarra}(iv) The \emph{rank-1} method
    in~\cite{segarra2017graphfilterdesign} obtains a shift matrix
    ${\hc{\bm S}}$ when ${\hc{\bm P}}$ is of rank 1 and then applies a graph
    filter. The resulting graph filter is generally of
    maximum order 
    $\nodenum-1$.
    \end{myitemize}%
    \end{myitemize}%


    \cmt{Implementation details}Some implementation details follow.
    \begin{myitemize
    \cmt{filter coefficients}To alleviate numerical issues, the method
    in~\cite{segarra2017graphfilterdesign} and the proposed algorithms
    use node-dependent
    coefficients~\cite[Sec.~II-B]{segarra2017graphfilterdesign}. For
    comparison purposes, an order $\itind$ filter ${\hc{\bm H}}_\itind$
    is obtained for each number $\itind$ of local exchanges by fitting
    the node-dependent coefficients to minimize
    $\|{\hc{\bm P}}-{\hc{\bm H}}_\itind\|_\text{F}^2$. 
    \cmt{feasibility}Regarding feasibility (see Sec.~\ref{sec:approximate}), the solver
     for \eqref{prob:exactrelaxed} in Appendix~\ref{sec:exactsolver}
    declares the problem as infeasible if the convergence criterion is
    not met after $I_{\text{max}}$ updates. 
    \cmt{Parameters}Both  proposed  methods use the
    same parameters in most experiments. This illustrates that a
     single set of parameters works reasonably well in a wide range of
     scenarios, which facilitates parameter tuning. For fairness,
     the competing methods use the same parameters in all experiments:
     $\dlmsauglagpar=0.001, \dlmsstepsize=1$,
     $\mu_{\text{DGD}}=0.1$. These values were adjusted to
     approximately yield the best performance in these scenarios.
    \end{myitemize}%


    \cmt{Metrics}Two main performance metrics will be adopted.
    \begin{myitemize
    \cmt{NMSE}To quantify estimation error, some experiments
    obtain the \emph{normalized mean square error} (NMSE), defined as
    $\text{NMSE}({\hc{\bm H}}_\itind):= \mathbb{E}\left[||\boldsymbol
    {\xi}-{\hc{\bm H}}_\itind \obsvec ||_2^2 \right]
    / \mathbb{E}\left[||{\hc{\bm \xi}} ||_2^2\right] $, where the
    expectation is taken over ${\hc{\mathcal{G}}}$,
    ${\hc{\bm U}}_{\parallel}, {\bm{\coord}},$ and ${\hc{\bm v}}$.
    \cmt{NMPE}For the methods based on graph shift operators, the
     error between the obtained ${\hc{\bm H}}_\itind$ and the target
     ${\hc{\bm P}}$ is measured through the \emph{normalized mean
     projection error}
     $\text{NMPE}({\hc{\bm H}}_{\itind})\triangleq
    \mathbb{E}\|{\hc{\bm P}} -{\hc{\bm H}}_\itind\|_\text{F}^2/
       \subspacedim$, where $\mathbb{E}$ averages over ${\hc{\mathcal{G}}}$,
     and ${\hc{\bm U}}_{\parallel}$. The normalization factor $\subspacedim$ was chosen so
     that $\text{NMPE}({\hc{\bm H}}_{\itind})$ equals the NMSE when
     $\bm \xi={\hc{\bm P}} \bm z'$ with $\bm z'\sim\mathcal{N}(\bm
     0_\nodenum,\bm I_\nodenum)$.
     \cmt{MC}\arev{Both the NMSE and NMPE will be estimated by
    averaging across 500 Monte Carlo realizations.}
    \cmt{infeasibility}The optimization problems solved by some of these methods, namely the exact
     projection and the gossiping method, may be infeasible for
     certain realizations of ${\hc{\mathcal{G}}}$ and ${\hc{\bm U}}_{\parallel}$; see
     Sec.~\ref{sec:approximate}. In that case, ${\hc{\bm H}}_\itind$ is
 set to $\bm 0$, which would penalize the tested algorithm by setting its  NMSE or NMPE closer to 1. 


\end{myitemize}%

\subsection{Exact Projection Filters}

\cmt{Shift-based approaches}Figs.~\ref{rank-one}-\ref{NMPE-3}
 depict the NMPE for those algorithms that rely on shift
 matrices. Whereas Figs.~\ref{rank-one}-\ref{NMPE-2} adopt an
 Erdős–Rényi random graph, a WSN is used in Fig.~\ref{NMPE-3}. The
 method in~\cite{segarra2017graphfilterdesign} is not displayed in
 Figs.~\ref{NMPE-2}-\ref{NMPE-3} because it cannot be applied when
 $\subspacedim>1$.
\begin{myitemize}%
\cmt{overall performance}Relative to the competing
 alternatives, the proposed exact projection method is seen to
 generally require a significantly smaller number of local exchange
 rounds to obtain an NMPE close to 0.
\cmt{finite iterations}%
\begin{myitemize
\cmt{proposed}The aforementioned figures also reveal that the proposed
 method provides an exact projection in a finite number of
 iterations. As predicted by the Cayley-Hamilton Theorem, the order is
 never greater than $\nodenum-1$. However, most of the times, the
 actual order is much smaller than this upper bound. This phenomenon
 may not be clear at first glance from Figs.~\ref{rank-one}-\ref{NMPE-3} because they reflect \emph{average} behavior across
 a large number of Monte Carlo iterations. For this reason,
     $\text{NMPE}({\hc{\bm H}}_{\itind})$ is generally positive for all
     $\itind<\nodenum-1$ when one or more realizations yield a
     filter with maximum order $\itind=\nodenum-1$.
\cmt{rank-1}On the other hand, the rank-1 method generally yields exact
 projection filters with order $\nodenum-1$ since it is not designed
 to minimize the order.  \cmt{gossiping}The gossiping method,
 however, does not necessarily produce an exact projection in a finite
 number of iterations. In exchange
 its implementation is simpler; cf. Sec.~\ref{sec:barbarossa}.
\end{myitemize}%

\end{myitemize}%



\cmt{Non-shift-based approaches}To compare with the competing
algorithms that are not based on graph filtering,
\begin{myitemize
\cmt{vs. num exchanges}Figs.~\ref{NMSE-4}-\ref{NMSE-5} depict     $\text{NMSE}({\hc{\bm H}}_\itind)$
vs.  the number $\itind$  of local exchanges for different
scenarios. 
\begin{myitemize}%
\cmt{NMSE}Observe that the NMSE of the exact projection  method does not converge to 0 due to
the observation noise ${\hc{\bm v}}$.
\cmt{scale}The miniatures demonstrate the different convergence time scales of
 DGD, DLMS, and the exact projection method. The cause for the slower
 convergence of DGD and DLMS is twofold: first, DGD and DLMS are
 general-purpose methods whereas the proposed method is tailored to
 the subspace projection problem. Second, DGD and DLMS require the
 selection of step size parameters, which in the simulations here were
 set to ensure convergence in virtually all Monte Carlo
 realizations. For certain specific realizations, though, one could
 find step sizes that yield a faster convergence.
\end{myitemize}%
\cmt{vs. rank}An alternative perspective to the comparison in
 Figs.~\ref{NMSE-4}-\ref{NMSE-5} is offered by Fig.~\ref{NMSE-6},
which depicts the NMSE vs. the subspace dimension $\subspacedim$ when
the number of local exchange rounds is fixed to $\itnum$. This
analysis would be necessary e.g. in real-time applications. It is
interesting to observe that the sensitivity to $\subspacedim$ is higher for the
gossiping algorithm as compared to the other methods.

\end{myitemize


\cmt{graph sparsity}
\begin{myitemize
\cmt{MNDE}To analyze the impact of the graph sparsity in the
required number of local exchanges to implement a projection exactly
with the proposed method, Fig.~\ref{MNDE} shows the \emph{mean number
of distinct eigenvalues} (MNDE) yielded by the exact projection method
vs.  $p_{\text{miss}}$. The NMDE is defined as the mean of the number
of distinct eigenvalues of ${\hc{\bm F}}_{\parallel}$ plus the number of distinct
eigenvalues of ${\hc{\bm F}}_{\perp}$, which equals the order of the filter;
cf. Sec.~\ref{sec:filterorderminimization}. A threshold $\tau$ is used
to determine whether two eigenvalues are different. As expected, the
filter order increases as the graph becomes sparser; see
Sec.~\ref{sec:feasibility}. Remarkably, the increase is more
pronounced for larger networks and higher $\subspacedim$.

\cmt{required local exchanges}In the case of WSN graphs, the
impact of sparsity is studied in Fig.~\ref{WSN-1}, which displays the
number of local exchange rounds required to reach an NMPE below a
target value $\gamma_\text{target}$.  The gossiping method is also
shown for comparison purposes, yet it is not capable of attaining this
target for certain values of $d_\text{max}$ due to the infeasibility
of the optimization problem that it solves in a significant fraction
of the Monte Carlo realizations; see explanation after the definition
of NMSE earlier in Sec.~\ref{sec:sim}. As $d_{\text{max}}$ increases,
the network becomes more densely connected and therefore an exact
projection can be obtained in a smaller number of local
exchanges. This agrees with intuition and with Fig.~\ref{MNDE}.


\end{myitemize



\subsection{Approximate Projection Filters}
\cmt{fixed num. local exchanges}As explained in
Sec.~\ref{sec:approximate}, a graph filter capable of exactly
 implementing a projection may not exist when the graph is highly
 sparsely connected. The approximate projection method in that section provides a
 graph filter capable of approximating such a projection. To
 illustrate how this algorithm complements the exact projection
 method from Sec.~\ref{sec:exactproj}, Fig.~\ref{approx-new1} shows
 NMPE after $\itnum=\nodenum-1$ local exchanges for both methods along
 with the gossiping algorithm.
 \begin{myitemize
\cmt{when exact}The exact projection method is seen to provide a filter that
exactly implements the target projection when $p_\text{miss}$ is below
approximately 0.7. For larger $p_\text{miss}$, the NMPE becomes
strictly positive. It is important to emphasize that the reason
is \emph{not} that this method provides filters that do not exactly
implement the target projection. When the problem is feasible, the
filters yield an exact projection. However, as $p_\text{miss}$
increases, a smaller fraction of Monte Carlo realizations give rise to
a feasible problem and this is penalized in the NMSE computation; see
explanation after the definition of NMSE earlier in
Sec.~\ref{sec:sim}.
\cmt{relative to approx}As expected, 
the exact projection method performs better than its approximate
counterpart when exact projection filters exist. However, for
sufficiently sparse graphs, the NMPE of the exact projection method
explodes, whereas the approximate projection method remains low. 
\end{myitemize


\cmt{vs. local exchanges}Figs.~\ref{NMPE-approx-er}
and~\ref{NMPE-approx-wsn} present the evolution of the NMPE vs. the
number of local exchange rounds for Erdős–Rényi and WSN graphs
respectively with several degrees of sparsity. These figures showcase
that the proposed approximate projection method can reasonably
approximate a projection with a small number of local exchanges even
if the graph is highly sparse. Note that for some of the sparsity
levels used in these figures, the problems solved by the gossiping
algorithm and the exact projection method become infeasible.


\cmt{Robustness}\arev{To analyze the robustness of the approximate projection
method, Fig.~\ref{extra} depicts its NMPE for different numbers of
node failures. Note that the exact projection method cannot be applied
since the removal of one or more sensors renders the exact projection
problem infeasible except for trivial cases where their measurements
bear no information about the estimated signal. Whenever the network
detects that a node has failed, the graph filter is recomputed with
the updated topology.  As expected, the degradation is more pronounced
when the network is sparser (smaller $d_{\max}$). Due to the
distribution used to generate ${\hc{\bm U}}_{\parallel}$, each measurement contains,
on average, a fraction $1/\nodenum$ of the signal energy. This informally
indicates that $n/\nodenum$ lower bounds the error attainable in the
presence of $n$ node failures. Observe that the NMPE is indeed close
to such fundamental limits, thereby establishing that the approximate
projection method is reasonably robust to node failures.  }


\section{\arev{Conclusions and Discussion}}
\label{sec:conclusions}

This paper develops methods to obtain graph filters that can be used
 to compute subspace projections in a decentralized fashion. The
 approach relies on transforming the filter design task into a
 shift-matrix design problem. The first method addresses the latter by
 exploiting the eigenstructure of feasible shift matrices. The second
 method builds upon the first to approximate projections when no
 feasible projection filter exists. An exhaustive simulation analysis
 demonstrates the ability of the designed filters to effectively
 produce subspace projections in a small number of iterations. This
 contrasts with existing methods, whose convergence is generally
 asymptotic.

 \arev{The main strength of graph filters lies in their simplicity. As
 described before, the inference operations considered in this paper can be
 implemented with $\itnum\nodenum$ transmissions. Lower and
 higher values can be attained by alternative communication
 strategies, including centralized processing. The approach adopted in
 practice will depend on factors such as energy and hardware
 constraints.  }

\arev{The filters obtained through the proposed methods also inherit
 the general numerical limitations of graph filters. In turn, these
 limitations stem from the well-known  conditioning issues of
 Vandermonde systems, e.g. \eqref{eq:vander}. This limits the
 application of graph filters to networks whose number of nodes is
 comparable to the ones in Sec.~\ref{sec:sim}. Further work is
 required by the graph signal processing community to enable the
 implementation of graph filters in significantly larger networks.}

\arev{Another direction along which the proposed schemes can be
 extended is by lifting the symmetry constraint on ${\hc{\bm S}}$. Since
 the target matrix ${\hc{\bm P}}$ is symmetric, the gain may not be
 significant. Future work will investigate this possibility.}


\balance
\if\editmode1 
\onecolumn
\printbibliography
\twocolumn
\else

\section{ADMM for Exact Projection Filters}

\cmt{Objective}To rewrite the objective, let $\canbasisvec_i$ and $\canbasisperp_j$ be respectively the $i$-th and $j$-th columns of $\bm I_\subspacedim$ and $\bm I_{\nodenum-\subspacedim}$. 
The first and  second norms in the objective of \eqref{prob:exactrelaxed} can be expressed as:
\begin{align}
{{\left\| {\hc{\bm F}}_{\parallel}\otimes\bm{I}_{\subspacedim}-\bm{I}_\subspacedim\otimes{\hc{\bm F}}_{\parallel} \right\|}_{*}}\nonumber &= {{\left\| \mathrm{vec}^{-1}({\hc{\bm A}}\mathrm{vec} ({\hc{\bm F}}_{\parallel})) \right\|}_{*}}\\
{{\left\| {\hc{\bm F}}_{\perp}\otimes\bm{I}_{\nodenum-\subspacedim}-\bm{I}_{\nodenum-\subspacedim}\otimes{\hc{\bm F}}_{\perp} \right\|}_{*}}\nonumber & =  {{\left\| \mathrm{vec}^{-1}({\hc{\bm B}}\mathrm{vec} ({\hc{\bm F}}_{\perp}) )\right\|}_{*}},
\end{align}
  where ${\hc{\bm A}}:= [{\hc{\bm a}}_{11}, {\hc{\bm a}}_{21},...,{\hc{\bm a}}_{\subspacedim\subspacedim}]$, ${\hc{\bm B}}:=[{\hc{\bm b}}_{11}, {\hc{\bm b}}_{21},...,{\hc{\bm b}}_{\nodenum-\subspacedim,\nodenum-\subspacedim}]$, ${\hc{\bm a}}_{ij}:=\mathrm{vec}({\canbasisvec_{i}\canbasisvec_{j}^{\top}\otimes\bm{I}_{\subspacedim}-\bm{I}_{\subspacedim}\otimes\canbasisvec_{i}\canbasisvec_{j}^{\top}})$, ${\hc{\bm b}}_{ij}:=\mathrm{vec}({\canbasisperp_{i}\canbasisperp_{j}^{\top}\otimes\bm{I}_{\nodenum-\subspacedim}-\bm{I}_{\nodenum-\subspacedim}\otimes\canbasisperp_{i}\canbasisperp_{j}^{\top}})$.


\cmt{constraints}To rewrite the  constraints of \eqref{prob:exactrelaxed},
\begin{myitemize}%
\cmt{1-2}invoke  the properties of the Kronecker product to
combine  the  first and  second constraints as ${\hc{\bm{W}}}(({\hc{\bm U}}_{\parallel}\otimes{\hc{\bm U}}_{\parallel})\mathrm{vec}({\hc{\bm F}}_{\parallel})+({\hc{\bm U}}_{\perp}\otimes{\hc{\bm U}}_{\perp})\mathrm{vec}({\hc{\bm F}}_{\perp}))=\bm{0}$.
\cmt{3-4}The third and fourth constraints can be rewritten as
${\hc{\bm G}}_{\parallel}\mathrm{vec}({\hc{\bm F}}_{\parallel})=\bm{0}$,
${\hc{\bm G}}_{\perp}\mathrm{vec}({\hc{\bm F}}_{\perp})=\bm{0}$ where the rows of
${\hc{\bm G}}_{\parallel}\in\mathbb{R}^{{(\subspacedim^{2}-\subspacedim)/2\times{\subspacedim}^{2}}}$
and ${\hc{\bm G}}_{\perp}\in\mathbb{R}^{{((\nodenum-\subspacedim)^{2}-(\nodenum-\subspacedim))/2\times{(\nodenum-\subspacedim)}^{2}}}$  are respectively given by $(\canbasisvec_{j}^{\top}\otimes\canbasisvec_{i}^{\top}-\canbasisvec_{i}^{\top}\otimes\canbasisvec_{j}^{\top})$ and $(\canbasisperp_{j}^{\top}\otimes\canbasisperp_{i}^{\top}-\canbasisperp_{i}^{\top}\otimes\canbasisperp_{j}^{\top})$ for $i<j$.
\cmt{trace}Regarding the trace constraints, rewrite $\mathrm{tr}({\hc{\bm F}}_{\parallel})=\subspacedim$ as $\mathrm{tr}(\bm{I}_\subspacedim{\hc{\bm F}}_{\parallel})=\subspacedim$,  which in turn can be expressed as $\mathrm{vec}^{\top}(\bm{I}_\subspacedim)\mathrm{vec}({\hc{\bm F}}_{\parallel})=\subspacedim$. Similarly, one can rewrite $\mathrm{tr}({\hc{\bm F}}_{\perp})=(1\pm\epsilon)(\nodenum-\subspacedim)$ as $\mathrm{vec}^{\top}(\bm{I}_{\nodenum-\subspacedim})\mathrm{vec}({\hc{\bm F}}_{\perp})=(1\pm\epsilon)(\nodenum-\subspacedim)$.
\end{myitemize}%

\cmt{rewritten problem}Therefore, \eqref{prob:exactrelaxed} can be written as:
 \begin{subequations} \label{eq:app} 
\begin {align}
\nonumber \underset{{\hc{\bm Y}}_{\perp},{\hc{\bm Y}}_{\parallel},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp}}{\text{min}} ~~&
          {\hc{ \eta}}_{\parallel} {{\left\| {\hc{\bm Y}}_{\parallel} \right\|}_{*}}
           +{\hc{ \eta}}_{\perp}{{\left\|{\hc{\bm Y}}_{\perp}\right\|}_{*}} +{{\left\|{\hc{\bm F}}_{\perp} \right\|}_\text{F}^{2}}\label{optt1}\\
\text{s. t.}   \;\;\;\;    	&  {\hc{\bm T}}_{\parallel}\mathrm{vec}({\hc{\bm F}}_{\parallel})+{\hc{\bm T}}_{\perp}\mathrm{vec}({\hc{\bm F}}_{\perp})={\hc{\bm b}}\\
&\mathrm{vec}({\hc{\bm Y}}_{\parallel})={\hc{\bm A}}\mathrm{vec} ({\hc{\bm F}}_{\parallel})\\
&\mathrm{vec}({\hc{\bm Y}}_{\perp})={\hc{\bm B}}\mathrm{vec} ({\hc{\bm F}}_{\perp}),
\end{align}
\end{subequations}
where
\begin{subequations}
\begin{align}
{\hc{\bm T}}_{\parallel}&\overset{\Delta}{=}[ {\hc{\bm{W}}}({\hc{\bm U}}_{\parallel}\otimes{\hc{\bm U}}_{\parallel});
									{\hc{\bm G}}_{\parallel};
                                                           \bm{0};
	                                                    \mathrm{vec^{\top}}(\bm{I}_\subspacedim);
                                                         
                                                           \bm{0}]
                                                           \\
{\hc{\bm T}}_{\perp}&\overset{\Delta}{=}[
                                                        {\hc{\bm{W}}}({\hc{\bm U}}_{\perp}\otimes{\hc{\bm U}}_{\perp});
                                                         \bm{0};
									{\hc{\bm G}}_{\perp};
	                                                    \bm{0};
                                                           \mathrm{vec^{\top}}(\bm{I}_{\nodenum-\subspacedim})]\\
{\hc{\bm b}}&\overset{\Delta}{=}[
                                                          \bm{0};
                                                         \bm{0};
									\bm{0};
	                                                    \subspacedim;
                                                             (1\pm\epsilon)(\nodenum-\subspacedim)].
\end{align}
\end{subequations}


 
    
       

\cmt{ADMM}The ADMM method
in  scaled form \cite[Sec. 3.1.1]{boyd2011distributed} applied to \eqref{eq:app} obtains the $k$-th iterate as follows:
 \begin{subequations} \label{ADMM1}
\begin{align}
&({\hc{\bm Y}}_{\parallel}^{k+1},{\hc{\bm Y}}_{\perp}^{k+1})\label{ADMM1a}\\& \quad\quad :=\underset{{\hc{\bm Y}}_{\parallel},{\hc{\bm Y}}_{\perp}}{\text{argmin}}~~L_{\rho}({\hc{\bm Y}}_{\parallel},{\hc{\bm Y}}_{\perp},{\hc{\bm F}}_{\parallel}^{k},{\hc{\bm F}}_{\perp}^{k},{\hc{\bm q}}_{1}^{k},{\hc{\bm q}}_{2}^{k},{\hc{\bm q}}_{3}^{k})\nonumber\\%
\label{ADMM1ab}&({\hc{\bm F}}_{\parallel}^{k+1},{\hc{\bm F}}_{\perp}^{k+1})\\& \quad\quad :=  \underset{{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp}}{\text{argmin}}~~ L_{\rho}({\hc{\bm Y}}_{\parallel}^{k+1},{\hc{\bm Y}}_{\perp}^{k+1},{\hc{\bm F}}_{\parallel}, {\hc{\bm F}}_{\perp},{\hc{\bm q}}_{1}^{k},{\hc{\bm q}}_{2}^k, {\hc{\bm q}}_{3}^{k}) \nonumber\\%\\
&{\hc{\bm q}}_{1}^{k+1}:={\hc{\bm q}}_{1}^{k}+{\hc{\bm T}}_{\parallel}\mathrm{vec}({\hc{\bm F}}_{\parallel}^{k+1})+{\hc{\bm T}}_{\perp}\mathrm{vec}({\hc{\bm F}}_{\perp}^{k+1})-{\hc{\bm b}}\\
&{\hc{\bm q}}_{2}^{k+1}:={\hc{\bm q}}_{2}^{k}+\mathrm{vec}({\hc{\bm Y}}_{\parallel}^{k+1})-{\hc{\bm A}}\mathrm{vec} ({\hc{\bm F}}_{\parallel}^{k+1}) \\%\label{ADMM_con1}\\
&{\hc{\bm q}}_{3}^{k+1}:={\hc{\bm q}}_{3}^{k}+\mathrm{vec}({\hc{\bm Y}}_{\perp}^{k+1})-{\hc{\bm B}}\mathrm{vec} ({\hc{\bm F}}_{\perp}^{k+1}),
\end{align}
\end{subequations}
where
\begin{align}
 L_{\rho}&({\hc{\bm Y}}_{\parallel}, {\hc{\bm Y}}_{\perp},{\hc{\bm F}}_{\parallel},{\hc{\bm F}}_{\perp},{\hc{\bm q}}_{1},{\hc{\bm q}}_{2},{\hc{\bm q}}_{3}) \nonumber\\:=&
{\hc{ \eta}}_{\parallel}{{\left\| {\hc{\bm Y}}_{\parallel} \right\|}_{*}}+{\hc{ \eta}}_{\perp}{{\left\| {\hc{\bm Y}}_{\perp} \right\|}_{*}}+{{\left\| {\hc{\bm F}}_{\perp} \right\|}_\text{F}^{2}}\nonumber\\
 &+(\rho/2){{\left\| {\hc{\bm T}}_{\parallel}\mathrm{vec}({\hc{\bm F}}_{\parallel})+{\hc{\bm T}}_{\perp}\mathrm{vec}({\hc{\bm F}}_{\perp})-{\hc{\bm b}}+{\hc{\bm q}}_{1} \right\|}_{2}^{2}}\nonumber\\
&+(\rho/2){{\left\| \mathrm{vec}({\hc{\bm Y}}_{\parallel})-{\hc{\bm A}}\mathrm{vec}
 ({\hc{\bm F}}_{\parallel})+{\hc{\bm q}}_{2}\right\|}_{2}^{2}}
 \nonumber\\
 & +(\rho/2){{\left\| \mathrm{vec}({\hc{\bm Y}}_{\perp})-{\hc{\bm B}}\mathrm{vec}
 ({\hc{\bm F}}_{\perp})+{\hc{\bm q}}_{3}\right\|}_{2}^{2}}\nonumber
\end{align}
is the so-called augmented Lagrangian with user-defined \emph{penalty parameter} $\rho>0$. Variables ${\hc{\bm q}}_{1}$, ${\hc{\bm q}}_{2}$ and ${\hc{\bm q}}_{3}$ correspond to the Lagrange multipliers of \eqref{eq:app}. 

To evaluate \eqref{ADMM1a}, one can leverage the proximal
operator of the nuclear norm~\cite[Th. 2.1]{cai2010singular}
\begin{align*}
\mathrm{prox}_\tau(\bm Z) := \argmin_{\bm Y}~||\bm Y||_* + \frac{1}{2\tau}||\bm Y-\bm Z||_\text{F}^2
= D_\tau(\bm Z),
\end{align*}
where the \emph{singular value shrinkage} operator $D_\tau$ is defined for $\bm Z$ with SVD $\bm Z = \bm V_{\bm Z} \bm \Sigma_{\bm Z}\bm V_{\bm Z}^\top$ as $D_\tau(\bm Z) :=
\bm V_{\bm Z} D_\tau(\bm \Sigma_{\bm Z})\bm V_{\bm Z}^\top$ with
$ D_\tau ({\bm \Sigma}_{\bm Z})$ a diagonal matrix whose $(i,i)$-th entry
equals $\max(0, (\bm \Sigma_{\bm Z})_{i,i}-\tau )$.  
 With this operator,  \eqref{ADMM1a} becomes
\begin{subequations}
\begin{align}
{\hc{\bm Y}}_{\parallel}^{k+1}=\mathrm{prox}_{{\hc{ \eta}}_{\parallel}/\rho}(\mathrm{vec}^{-1}({\hc{\bm A}}\mathrm{vec} ({\hc{\bm F}}_{\parallel}^{k})-{\hc{\bm q}}_{2}^{k})) \label{eq:Y_par}\\
{\hc{\bm Y}}_{\perp}^{k+1}=\mathrm{prox}_{{\hc{ \eta}}_{\perp}/\rho}(\mathrm{vec}^{-1}({\hc{\bm B}}\mathrm{vec} ({\hc{\bm F}}_{\perp}^{k})-{\hc{\bm q}}_{3}^{k})) \label{eq:Y_per}.
\end{align}
\end{subequations}
On the other hand,  \eqref{ADMM1ab} can be obtained in closed-form as
\begin{subequations}
\begin{align}
\label{eq:F_par}
&\mathrm{vec}({\hc{\bm F}}_{\parallel}^{k+1})=[{\hc{\bm T}}_{\parallel}^{\top}{\hc{\bm T}}_{\parallel}+{\hc{\bm A}}^{\top}{\hc{\bm A}}]^{-1}\\&[{\hc{\bm T}}_{\parallel}^{\top}({\hc{\bm b}}-{\hc{\bm q}}_{1}^{k}-{\hc{\bm T}}_{\perp}\mathrm{vec}({\hc{\bm F}}_{\perp}^{k}))
+{\hc{\bm A}}^{\top}({\hc{\bm q}}_{2}^{k}+\mathrm{vec}({\hc{\bm Y}}_{\parallel}^{k+1}))],\nonumber\\
\label{eq:F_per}
&\mathrm{vec}({\hc{\bm F}}_{\perp}^{k+1})=[2\bm{I}+\rho{\hc{\bm T}}_{\perp}^{\top}{\hc{\bm T}}_{\perp}+\rho{\hc{\bm B}}^{\top}{\hc{\bm B}}]^{-1}\\&[\rho{\hc{\bm T}}_{\perp}^{\top}({\hc{\bm b}}-{\hc{\bm q}}_{1}^{k}-{\hc{\bm T}}_{\parallel}\mathrm{vec}({\hc{\bm F}}_{\parallel}^{k}))
+\rho{\hc{\bm B}}^{\top}({\hc{\bm q}}_{3}^{k}+\mathrm{vec}({\hc{\bm Y}}_{\perp}^{k+1}))]\nonumber.
\end{align}
\end{subequations}
The overall ADMM algorithm is summarized as
Algorithm~\ref{exact-algorithm} and drastically improves the converegence rate of
our previous algorithm, reported in~\cite{mollaebrahim2018large}, since the latter
is based on subgradient descent.

\cmt{complexity}\arev{The computational complexity is dominated either  by \eqref{eq:Y_par} or by \eqref{eq:Y_per}, whose SVDs respectively require $O(\subspacedim^6)$  and  $O((\nodenum-\subspacedim)^6)$  arithmetic operations, depending on which quantity is larger. This complexity is much lower than the one of  general-purpose convex solvers and does not limit the values of $\nodenum$ and $\subspacedim$ to be used in practice given the intrinsic limitations of graph filters; see Sec.~\ref{sec:conclusions}.
}

\begin{algorithm}[t]
\caption{Iterative solver for \eqref{prob:exactrelaxed}.}
\label{exact-algorithm}
\begin{algorithmic}[1]
\REQUIRE ${\hc{\mathcal{E}}}$, ${\hc{\bm U}}_{\parallel}$, 
 $\ADMMreg$, ${\hc{ \eta}}_{\parallel}$, ${\hc{ \eta}}_{\perp}$, ${\hc{\epsilon}}$.
\STATE{Obtain ${\hc{\bm U}}_{\perp}$ s.t. ${\hc{\bm U}}_{\perp}^\top{\hc{\bm U}}_{\perp}=\bm
I_{\nodenum-\subspacedim}$ and ${\hc{\bm U}}_{\parallel}^\top{\hc{\bm U}}_{\perp}=\bm 0_{\subspacedim\times\nodenum}$.}
\STATE{Initialize
${\hc{\bm F}}_{\parallel}^0,~{\hc{\bm F}}_{\perp}^0,~ {\hc{\bm q}}_{1}^{0},~{\hc{\bm q}}_{2}^{0},~{\hc{\bm q}}_{3}^{0}, k=0
$.}

\WHILE{stopping\_criterion\_not\_met()}
	

	
\STATE{obtain ${\hc{\bm Y}}_{\parallel}^{{k+1}}$ via \eqref{eq:Y_par}.}	
\STATE{obtain ${\hc{\bm Y}}_{\perp}^{{k+1}}$ via \eqref{eq:Y_per}.}	
\STATE{obtain ${\hc{\bm F}}_{\parallel}^{{k+1}}$ via \eqref{eq:F_par}.}	
\STATE{obtain ${\hc{\bm F}}_{\perp}^{{k+1}}$ via \eqref{eq:F_per}.}	
\STATE{ obtain ${\hc{\bm q}}_{1}^{k+1},{\hc{\bm q}}_{2}^{k+1},{\hc{\bm q}}_{3}^{k+1}$ via \eqref{ADMM1}.}	
\STATE{ $k\leftarrow k+1$.}	
	\ENDWHILE
\RETURN	$ {\hc{\bm S}} = {\hc{\bm U}}_{\parallel}{\hc{\bm F}}_{\parallel}^{k}{\hc{\bm U}}_{\parallel}^\top
+ {\hc{\bm U}}_{\perp}{\hc{\bm F}}_{\perp}^{k}{\hc{\bm U}}_{\perp}^\top$.
\end{algorithmic}
\end{algorithm}




