Non-Bayesian Estimation Framework for Signal Recovery on Graphs

In the rest of this paper, we denote vectors by boldface lowercase letters and matrices by boldface uppercase letters. The operators $(\cdot)^{T}$, $(\cdot)^{-1}$, $(\cdot)^{\dagger}$, and ${\text{Tr}}(\cdot)$ denote the transpose, inverse, Moore-Penrose pseudo-inverse, and trace operators, respectively. For a matrix ${\bf{A}}\in{\mathbb{R}}^{M\times K}$ with a full column rank, ${\bf{P}}_{\bf{A}}^{\bot}={\bf{I}}_{M}-{\bf{A}}({\bf{A}}^{T}{\bf{A}})^{-1}{\bf{A}}^{T}$, where ${\bf{I}}_{M}$ is the identity matrix of order $M$. The matrix ${\text{diag}}({\bf{a}})$ denotes the diagonal matrix with vector ${\bf{a}}$ on the diagonal. The $m$th element of the vector ${\bf{a}}$ and the $(m,q)$th element of the matrix ${\bf{A}}$ are denoted by $a_{m}$ and ${\bf{A}}_{m,q}$, respectively. Similarly, ${\bf{A}}_{\mathcal{S}_{1},\mathcal{S}_{2}}$ is used to denote the submatrix of ${\bf{A}}$ whose rows are indicated by the set $\mathcal{S}_{1}$ and columns are indicated by the set $\mathcal{S}_{2}$. For simplicity, we denote ${\bf{A}}_{\mathcal{S},\mathcal{S}}$ by ${\bf{A}}_{\mathcal{S}}$. For two symmetric matrices of the same size ${\bf{A}}$ and ${\bf{B}}$, ${\bf{A}}\succeq{\bf{B}}$ means that ${\bf{A}}-{\bf{B}}$ is a positive semidefinite matrix. The gradient of a vector function ${\bf{g}}({\mbox{\boldmath$\theta$}})$, $\nabla_{{\mbox{\boldmath{\scriptsize$\theta$}}}}{\bf{g}}({\mbox{\boldmath$\theta$}})$, is a matrix in $\mathbb{R}^{K\times M}$, with the $(k,m)$th element equal to $\frac{\partial g_{k}}{\partial\theta_{m}}$, where ${\bf{g}}=\left[g_{1},\ldots,g_{K}\right]^{T}$ and ${\mbox{\boldmath$\theta$}}=\left[\theta_{1},\ldots,\theta_{M}\right]^{T}$. For any index set, ${\mathcal{S}}\subset\{1,\dots,M\}$, ${\mbox{\boldmath$\theta$}}_{{\mathcal{S}}}$ is a subvector of $\theta$ containing the elements indexed by ${\mathcal{S}}$, where $|{\mathcal{S}}|$ and ${\mathcal{S}}^{c}\stackrel{{\scriptstyle\triangle}}{{=}}\{1,\dots,M\}\backslash{\mathcal{S}}$ denote the set’s cardinality and the complement set, respectively. The vectors ${\bf{1}}$ and ${\bf{0}}$ are column vectors of ones and zeros, respectively, ${\bf{e}}_{m}$ is the $m$th column of the identity matrix, all with appropriate dimensions. The notation $\mathbf{1}_{\mathcal{A}}$ denotes the indicator of an event $\mathcal{A}$. The number of non-zero entries in $\theta$ is denoted by $||{\mbox{\boldmath$\theta$}}||_{0}$. Finally, the notation ${\rm{E}}_{{\mbox{\boldmath{\scriptsize$\theta$}}}}[\cdot]$ represents the expected value parameterized by a deterministic parameter $\theta$.

In this section, we briefly review relevant concepts related to GSP [2, 3] that will be used in this paper. Consider an undirected weighted graph, ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$, where ${\mathcal{M}}=\left\{1,\ldots,M\right\}$ denotes the set of $M$ nodes or vertices and $\xi$ denotes the set of edges with cardinality $|\xi|$. We only consider simple graphs, with no self-loops and multi-edges. The symmetric matrix ${\bf{W}}$ is the weighted adjacency matrix with entry ${\bf{W}}_{m,k}$ denoting the weight of the edge $(m,k)\in\xi$, reflecting the strength of the connection between the nodes $m$ and $k$. This weight may be a physical measure or conceptual, such as a similarity measure. We assume for simplicity that the edge weights in ${\bf{W}}$ are non-negative (${\bf{W}}_{m,k}\geq 0$). When no edge exists between $m$ and $k$, the weight is set to 0, i.e. ${\bf{W}}_{m,k}=0$.

The Laplacian matrix, which contains the information on the graph structure, is defined by ${\bf{L}}\stackrel{{\scriptstyle\triangle}}{{=}}{\bf{D}}-{\bf{W}}$, where ${\bf{D}}$ is a diagonal matrix with ${\bf{D}}_{m,m}=\sum_{k=1}^{M}{\bf{W}}_{m,k}$. The Laplacian matrix, ${\bf{L}}$, is a real, symmetric, and positive semidefinite matrix, which satisfies the null-space property, ${\bf{L}}{\bf{1}}_{M}={\bf{0}}$, and has nonpositive off-diagonal elements. Thus, its associated singular value decomposition (SVD) is given by

where the columns of ${\bf{V}}$ are the eigenvectors of ${\bf{L}}$, ${\bf{V}}^{T}={\bf{V}}^{-1}$, and ${\bf{\Lambda}}\in\mathbb{R}^{M\times M}$ is a diagonal matrix consisting of the distinct eigenvalues of ${\bf{L}}$, $0=\lambda_{1}<\lambda_{2}<\ldots<\lambda_{M}$. Throughout this paper we will focus on the case where the observed graph is connected and, thus, $\lambda_{m}>0$, $m=2,\ldots,M$. If the graph is not connected, the proposed approach can be applied to each connected component separately. The eigenvalues $\lambda_{1},\ldots,\lambda_{M}$ can be interpreted as graph frequencies, and eigenvectors, i.e. the columns of the matrix ${\bf{V}}$, can be interpreted as corresponding graph frequency components. Together they define the graph spectrum for graph ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$.

In this framework, a graph signal is defined as a function ${\mbox{\boldmath$\theta$}}:{\mathcal{M}}\rightarrow{\mathbb{R}}^{M}$, assigning a scalar value to each vertex, where entry $\theta_{m}$ denotes the signal value at node $m\in{\mathcal{M}}$. The graph Fourier transform (GFT) of a graph signal $\theta$ w.r.t. the graph ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$ is defined as the projection onto the orthogonal set of the eigenvectors of ${\bf{L}}$ [2, 3]:

Similarly, the inverse GFT is obtained by left multiplication of a vector by ${\bf{V}}$, i.e. by ${\bf{V}}\tilde{{\mbox{\boldmath$\theta$}}}$. The GFT plays a central role in GSP since it is a natural extension of filtering operations and the notion of the spectrum of the graph signals [3, 2].

We consider the problem of estimating the graph signal, ${\mbox{\boldmath$\theta$}}\in{\mathbb{R}}^{M}$, which is considered in this paper to be a deterministic parameter vector. The estimation is based on a random observation vector, ${\bf{x}}\in\Omega_{\bf{x}}$, where $\Omega_{\bf{x}}$ is a general observation space. We assume that ${\bf{x}}$ is distributed according to a known probability distribution function (pdf), $f({\bf{x}};{\mbox{\boldmath$\theta$}})$, which is parametrized by the graph signal, $\theta$. For example, for the problem of mean estimation in Gaussian graphical models (GGM) [38, 39], $f({\bf{x}};{\mbox{\boldmath$\theta$}})$ is a Gaussian distribution with mean $\theta$ and covariance matrix, $\Sigma$, where ${\mbox{\boldmath$\Sigma$}}^{-1}$ preserves the connectivity pattern of the graph. Our goal is to integrate the side information in the form of graph structure, encoded by the Laplacian matrix, in the estimation approach.

Let $\hat{{\mbox{\boldmath$\theta$}}}:\Omega_{\bf{x}}\rightarrow{\mathbb{R}}^{M}$ be an estimator of $\theta$, based on a random observation vector, ${\bf{x}}\in\Omega_{\bf{x}}$. We consider estimators in the Hilbert space of bounded energy on ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$ [40], i.e. we assume that ${\rm{E}}_{{\mbox{\boldmath{\scriptsize$\theta$}}}}[\hat{{\mbox{\boldmath$\theta$}}}^{T}{\bf{L}}\hat{{\mbox{\boldmath$\theta$}}}]<\infty$. Our goal in this paper is to investigate estimation methods and bounds that are based on defining a structural estimation performance measure that reflects the graph topology and graph signal properties.

In various GSP problems there is side information on the graph signals in the form of linear parametric constraints. These linear constraints describe properties of the graph signals, such as bandlimitedness (see the example in Subsection VI), or properties of the sensing approach, such as the existence of reference nodes (anchors) [20, 21], can serve to obtain a well-posed estimation problem [33]. Formally, in these cases it is known a-priori that $\theta$ satisfies the following linear constraint:

where ${\bf{G}}\in{\mathbb{R}}^{K\times M}$ and ${\bf{a}}\in\mathbb{R}^{K}$ are known, and $0\leq K\leq M$. We assume that the matrix ${\bf{G}}$ has full row rank, i.e. the constraints are not redundant. The constrained set is:

Thus, we are interested in the problem of recovering graph signals that belong to $\Omega_{\mbox{\boldmath{\scriptsize$\theta$}}}$. We define the orthonormal null space matrix, ${\bf{U}}\in{\mathbb{R}}^{M\times(M-K)}$, such that [41]

The case $K=0$ represents an unconstrained estimation problem, in which we use the convention that ${\bf{U}}={\bf{I}}_{M}$ [42, 43, 44]. The matrix ${\bf{U}}$ can be found based on the eigenvector matrix of the orthogonal projection matrix ${\bf{P}}_{\bf{G}}^{\perp}\stackrel{{\scriptstyle\triangle}}{{=}}{\bf{I}}_{M}-{\bf{G}}^{T}({\bf{G}}{\bf{G}}^{T})^{-1}{\bf{G}}$.

This paper deals with a general parameter estimation over networks. In order to demonstrate how our proposal can be applied in practice, we give a physically-motivated example of state estimation in power systems. A power system can be represented as an undirected weighted graph, ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$, where the set of vertices, $\mathcal{M}$, is the set of buses (generators or loads) and the edge set, $\xi$, is the set of transmission lines between these buses. The weight matrix, ${\bf{W}}$, is determined by the branch susceptances [13, 23]. The goal of PSSE, which is at the core of energy management systems for various monitoring and analysis purposes, is to estimate $\theta$ based on system measurements, ${\bf{x}}$, that usually include power measurements [22]. Since $\theta$ is measured over the buses of the electrical network, PSSE can be interpreted as a graph signal recovery problem. The pdf, $f({\bf{x}};{\mbox{\boldmath$\theta$}})$, describes the physical relation between the power and the voltages, as well as the statistical behavior of noises due to, e.g. errors and varying temperatures. In Fig. 1.a, we show the one-line diagram of the IEEE 30-bus system, which is a well-known test grid for power system applications [45]. In Fig. 1.b, we show the graphical representation of this grid. The node color represents the state signal value. It can be seen that neighboring buses have a similar value (similar color) and, thus, the state signal is smooth.

The characteristics of the state estimation problem in power systems are as follows. First, the power flow equations are an up-to-a-constant function w.r.t. the state vector, $\theta$ [23], i.e. relative measurements, as described in Section V. Thus, the MSE cannot be used directly as a performance measure without the use of a reference bus (node). Additionally, it has recently been shown in [12, 11] that the state vector in power systems is a smooth graph signal w.r.t. the associated Laplacian matrix that is determined by the susceptances values. Thus, the GSP framework is well suited for PSSE, as shown, for example, in [11, 12] for addressing the problem of detection of false data injection attacks in power systems based on the GFT of the states.

Lehmann [37] proposed a generalization of the unbiasedness concept based on the chosen cost function for each scenario, which can be used for various cost functions (see, e.g. [48, 17, 17, 43]). The Lehmann unbiasedness definition for a matrix cost functions is defined as follows (p. 13 in [49]):

The following proposition defines a sufficient condition for the graph unbiasedness of estimators in the Lehmann sense w.r.t. the Laplacian-based WMSE and under the linear parametric constraints.

It can be seen that if an estimator has a zero mean bias, i.e. ${\rm{E}}_{{\mbox{\boldmath{\scriptsize$\theta$}}}}[\hat{{\mbox{\boldmath$\theta$}}}]-{\mbox{\boldmath$\theta$}}={\bf{0}}$, $\forall{\mbox{\boldmath$\theta$}}\in\Omega_{\mbox{\boldmath{\scriptsize$\theta$}}}$, then it satisfies (19), but not vice versa. Thus, the uniform graph unbiasedness condition is a weaker condition than requiring the mean-unbiasedness property. For example, since ${\bf{L}}{\bf{1}}={\bf{0}}$, the condition in (19) is oblivious to the estimation error of a constant bias over the graph, $c{\bf{1}}$, for any constant $c\in{\mathbb{R}}$, in contrast with mean unbiasedness. It can be seen that the graph unbiasedness in (19) is a function of the specific graph topology. In addition, the unbiasedness definition can be reformulated by using any matrix that spans the null space of ${\bf{G}}$, ${\mathcal{N}}({\bf{G}})$, instead of ${\bf{U}}$.

Another special case of the graph-unbiasedness for a bandlimited graph signal is discussed in Subsection VI-B.

In non-Bayesian estimation theory, two types of unbiasedness are usually considered: uniform unbiasedness at any point in the parameter space, and local unbiasedness, in which the estimator is assumed to be unbiased only in the vicinity of the parameter ${\mbox{\boldmath$\theta$}}_{0}$. By using the feasible directions of the constraint set, similar to the derivations in [44, 43], it can be shown that the local graph-unbiasedness conditions are as follows:

In the following, a novel graph CRB for the estimation of graph parameters is derived. Various low-complexity, distributive algorithms exist for graph signal recovery. The new bound is a useful tool for assessing their performance. The proposed graph CRB is a lower bound on the WMSE of local graph-unbiased estimators in the vicinity of $\theta$.

Similar to the explanations in Subsection III-1, it can be seen that the first row and column of the proposed matrix bound in (26), ${\bf{B}}({{\mbox{\boldmath$\theta$}}})$, is zero. Thus, in practice, we use the submatrix bound, ${\bf{B}}_{\mathcal{M}\setminus 1}({{\mbox{\boldmath$\theta$}}})$. Discussion of pseudo-inverse matrix bounds in the general case can be found, for example, in [50].

We consider a noisy measurement of the weighted relative state of each edge, as follows:

where $\{\nu_{m,k}\}$, $\forall(m,k)\in\bar{\xi}$, $m>k$, is an i.i.d. Gaussian noise sequence with variance $\sigma^{2}$. The sequence $\{\bar{w}_{m,k}\}$, $\forall(m,k)\in\bar{\xi}$, contains positive weights that are given by the system parameters and is assumed to be known. We assume that the edge weights satisfy $\bar{w}_{k,m}=\bar{w}_{m,k}$. Thus, from symmetry, the measurement and the noise sequences satisfy $h_{k,m}=-h_{m,k}$ and $\nu_{k,m}=-\nu_{m,k}$, respectively, $\forall(m,k)\in\bar{\xi}$. In general, $\bar{\xi}$ and $\{\bar{w}_{m,k}\}$, that are associated with the measurements and determined by the sensing approach, are different from $\xi$ and $\{{w}_{m,k}\}$, that are based on the physical graph. The goal here is to estimate the state vector, ${\mbox{\boldmath$\theta$}}=[\theta_{1},\ldots,\theta_{M}]^{T}$, from the observations in (33). For example, this problem with the model in (33) where $w_{m,k}=1$, $\forall(m,k)\in\bar{\xi}$, is the problem of synchronization of translations from [33, 53] . Another example is PSSE in electrical networks, [10, 11, 13], as described in Section VII.

Let ${\mathcal{G}}({\mathcal{M}},\bar{\xi},\bar{{\bf{W}}})$ be defined as the measurement graph associated with the model in (33) over the same vertex set as in the “physical” graph, ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$, which is used in (7). We assume that ${\mathcal{G}}({\mathcal{M}},\bar{\xi},\bar{{\bf{W}}})$ is a connected simple graph and define its associated Laplacian matrix, $\bar{{\bf{L}}}{\in{\mathbb{R}}^{M\times M}}$, with the elements

where, similar to Definition 1, $\bar{\mathcal{N}}_{m}=\{k\in{\mathcal{M}}:(m,k)\in\bar{\xi}\}$. In addition, let the matrix $\bar{{\bf{E}}}$ be the ${M\times|\bar{\xi}|}$ oriented incidence matrix of the graph ${\mathcal{G}}({\mathcal{M}},\bar{\xi},\bar{{\bf{W}}})$. Thus, each of the columns of $\bar{{\bf{E}}}$ has two nonzero elements, $1$ in the $m$th row and $-1$ in the $k$th row, representing an edge connecting nodes $m$ and $k$, where the sign is chosen arbitrarily [55]. Then, the model in (33) can be rewritten in a matrix form as follows:

where ${\bf{x}}$, $\bar{{\bf{w}}}$, and $\nu$ include the elements $\{h_{m,k}\}$, $\{\bar{w}_{m,k}\}$, $\{\nu_{m,k}\}$, respectively, $\forall(m,k)\in\bar{\xi}$, $m>k$, in the same order, and we used the symmetry in the model, i.e. the fact that $\bar{w}_{k,m}=\bar{w}_{m,k}$, $h_{k,m}=-h_{m,k}$, and $\nu_{k,m}=-\nu_{m,k}$. By multiplying (35) by $\bar{{\bf{E}}}$ from the left and using the fact that [55]

we obtain

Since $\{\nu_{m,k}\}$, $m>k$, is an i.i.d. Gaussian noise sequence with variance $\sigma^{2}$, then $\bar{{\bf{E}}}{\mbox{\boldmath$\nu$}}$ is a zero-mean Gaussian vector with a covariance matrix $\sigma^{2}\bar{{\bf{E}}}\bar{{\bf{E}}}^{T}$.

The log-likelihood function for the modified model in (37), after removing constant terms and by using $\bar{{\bf{L}}}^{T}=\bar{{\bf{L}}}$, is

Thus, the gradient of (38) satisfies

The matrix $\bar{{\bf{E}}}\bar{{\bf{E}}}^{T}$ is a Laplacian matrix with $M$ nodes and equal unit weights for all edges. Thus, its pseudo-inverse is given by [56]

By substituting (40) in (39) and using the null-space property, $\bar{{\bf{L}}}{\bf{1}}={\bf{0}}$, one obtains

Thus, the FIM from (24) for this case is given by

The FIM in (42) is a function of the graph topology via the Laplacian and the oriented incidence matrices. For the special case of unit weights, i.e. for $\bar{{\bf{L}}}=\bar{{\bf{E}}}\bar{{\bf{E}}}^{T}$, the FIM satisfies ${\bf{J}}({\mbox{\boldmath$\theta$}})=\frac{1}{\sigma^{2}}\bar{{\bf{L}}}$, which coincides with the result in [33].

It is shown in Appendix D that the pseudo-inverse of the FIM in (42) is given by

By substituting (43) and ${\bf{U}}={\bf{I}}_{M}$ (since there are no constraints in the considered model) in (25)-(26), we obtain that the graph CRB for this case is

which is not a function of the specific values of $\theta$. By applying the trace operator on (44), we obtain the bound on the expected Dirichlet energy, which is

The graph CRBs are lower bounds on the WMSE and on the Dirichlet energy of any graph unbiased estimators, as defined for this unconstrained setting in (20) in Special case 1.

For the estimator, by substituting (39), (43), and ${\bf{U}}={\bf{I}}_{M}$ in (27), equality in (44) is obtained iff

By using the properties of the Laplacian and incidence matrices, as well as its pseudo-inverse properties, it can be shown that $\bar{{\bf{E}}}^{T}\bar{{\bf{E}}}\bar{{\bf{L}}}^{\dagger}\bar{{\bf{L}}}(\bar{{\bf{E}}}\bar{{\bf{E}}}^{T})^{\dagger}\bar{{\bf{E}}}=\bar{{\bf{E}}}$ and $\bar{{\bf{E}}}^{T}\bar{{\bf{E}}}\bar{{\bf{L}}}^{\dagger}\bar{{\bf{L}}}(\bar{{\bf{E}}}\bar{{\bf{E}}}^{T})^{\dagger}\bar{{\bf{L}}}=\bar{{\bf{L}}}$. Thus, (46) implies the following condition for achievablity of the bound:

where the second equality is obtained by using the fact that for connected graphs there is only one zero eigenvalue of the Laplacian matrix, which is associated with the eigenvector ${\bf{v}}_{1}=\frac{1}{\sqrt{M}}{\bf{1}}$, and the last equality is obtained by using the Laplacian null-space property, ${\bf{\Lambda}}^{\frac{1}{2}}{{\bf{V}}}^{T}{\bf{1}}={\bf{0}}$. Since the estimator on the r.h.s. of (47) is not a function of $\theta$, then it is also an efficient estimator, in the sense of Definition 4. By using the null-space property ${{\bf{V}}}^{T}{\bf{1}}={\bf{0}}$, it can be verified that any estimator of the form

with an arbitrary scalar $c\in{\mathbb{R}}$, satisfies the equality condition in (47). Thus, the efficient estimator is not unique and the true signals, $\theta$, can be recovered up to a constant vector, $c{\bf{1}}$. This result is reasonable, since, with relative measurements alone, as given in the model in (33), determining the signal is possible only up to an additive constant [20, 21, 23].

In this subsection, we design a sample allocation rule for the sensing model from Subsection V-A based on solving an optimization problem that aims to minimize the graph CRB in (45). In this case, the graph CRB is achievable by the efficient estimator in (48) and is not a function of the unknown graph signal, $\theta$. Thus, minimizing the graph CRB in (45) will result in the minimum WMSE over graph unbiased estimators.

In many cases, the physical topology is known and is the one that is both used in the cost function in (7) and governs the measurement model in (33) for any edge that is measured. Thus, in the following we assume that the edge weights in (33) satisfy $\bar{w}_{m,k}={w}_{m,k}$, $\forall(m,k)\in\xi\bigcap\bar{\xi}$, and the sensors can be located only in existing edges of the physical network, where if $(m,k)\in\bar{\xi}$, then $(k,m)\in\bar{\xi}$. We assume a constrained amount of sensing resources, e.g. due to limited energy and communication budget. We thus state the sensor placement problem as follows:

The requirement of correct graph signal recovery in Problem 1 is defined as recovery up to a constant vector, $c{\bf{1}}$, which is an inherent ambiguity with relative measurements. In order to correctly reconstruct the full graph signal by the estimator from (48) up-to-a-constant, we need to have ${\text{rank}}(\bar{{\bf{L}}})=M-1$. That is, for each node, $m\in M$, we need to have at least one edge in $\bar{\xi}$ which is associated with this node. It can be shown that this condition for unique up-to-a-constant recovery is that $\bar{\xi}$ is a spanning tree of ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$ [57]. Thus, Problem 1 is equivalent to the following optimization problem.

If $\bar{{\bf{L}}}$ is a spanning tree of a connected (loopless) graph, then it can be seen that its oriented incidence matrix, $\bar{{\bf{E}}}$, has $M-1$ linearly independent columns. Thus, $\bar{{\bf{E}}}^{T}(\bar{{\bf{E}}}\bar{{\bf{E}}}^{T})^{\dagger}\bar{{\bf{E}}}=\bar{{\bf{E}}}^{\dagger}\bar{{\bf{E}}}={\bf{I}}_{M-1}$ and by substituting (36), we obtain that for spanning trees

That is, $\bar{\bar{{\bf{L}}}}$ is also a Laplacian matrix of a spanning tree of ${\bf{L}}$ with the weights $\{{w}_{m,k}^{2}\}$, $\forall(m,k)\in\bar{\xi}$, and the optimization from (49) is reduced to

It is shown in [58] that

where $st_{{\bar{{\bf{L}}}}}({\bf{L}})$ is the total stretch of the graph ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$ w.r.t. the spanning tree represented by ${\bar{{\bf{L}}}}$. Total stretch is a parameter used to measure the quality of a spanning tree in terms of distance preservation and can represent the average effective resistance [58, 59]. The objective function in Problem 3 is the $\ell_{p=\frac{1}{2}}$ total stretch of the graph ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}})$, i.e. the total stretch from (52) after taking the $p=\frac{1}{2}$ power of the weights of the edges [60]. Thus, the problem of finding the minimum graph CRB is equivalent to the problem of finding a spanning tree with the minimum average $\ell_{p=\frac{1}{2}}$ stretch, which is a classical problem in network design, graph theory, and discrete mathematics. Although this problem in general is an NP-hard problem, it was shown that a standard maximum weight spanning tree algorithm [61] yields good results in practice [62]. Thus, in this paper we solve the sensor allocation problem in Problem 3 by finding the maximum weight spanning tree of ${\mathcal{G}}({\mathcal{M}},\xi,{\bf{W}}^{2})$, which is one of the fundamental problems of algorithmic graph theory and can be solved by existing algorithms in near-linear time [61, 63].

We assume the following measurement model:

where $\theta$ is an unknown signal, ${\bf{w}}$ is a noise vector with zero mean and a known distribution, and ${\bf{M}}\in[0,1]^{D\times M}$ is a binary mask matrix with $D\leq M$. The mask matrix, ${\bf{M}}$, may indicate missing data, i.e. the indicator matrix for the missing values [66], removing data which is irrelevant for a specific task [67, 68], as well as representing any interpolation, filtering, and normalization approaches [64]. The task is to recover the true signal, $\theta$, based on the accessible measurement vector, ${\bf{x}}$. Signal recovery from inaccessible and corrupted measurements requires additional knowledge of signal properties. In GSP, a widely-used assumption is that the signal of interest, $\theta$, is a graph-bandlimited signal [3, 2, 4], i.e. its GFT, $\tilde{{\mbox{\boldmath$\theta$}}}$, defined in (2), is a $R$-sparse vector. Here we assume a graph low-frequency signal defined as follows:

The condition in (54) can be represented as the linear constraint from (3) with

where ${\bf{Q}}$ is an $(M-R)\times M$ matrix of the form $\ {\bf{Q}}=\left[{\bf{0}}_{(M-R)\times R},{\bf{I}}_{M-R}\right]$. Thus, a null-space matrix ${\bf{U}}$, defined in (5), can be chosen to be ${\bf{U}}={\bf{V}}\bar{{\bf{Q}}}$, where $\bar{{\bf{Q}}}=\left[\begin{array}[]{c}{\bf{I}}_{R}\\ {\bf{0}}_{(M-R)\times R}\end{array}\right]$.