Distributed weighted least-squares estimation with fast convergence for large-scale systems

In this paper we study a distributed weighted least-squares estimation problem for a large-scale system consisting of a network of interconnected sub-systems. Each sub-system is concerned with a subset of the unknown parameters and has a measurement linear in the unknown parameters with additive noise. The distributed estimation task is for each sub-system to compute the globally optimal estimate of its own parameters using its own measurement and information shared with the network through neighborhood communication. We first provide a fully distributed iterative algorithm to asymptotically compute the global optimal estimate. The convergence rate of the algorithm will be maximized using a scaling parameter and a preconditioning method. This algorithm works for a general network. For a network without loops, we also provide a different iterative algorithm to compute the global optimal estimate which converges in a finite number of steps. We include numerical experiments to illustrate the performances of the proposed methods.


Introduction
A sensor network is a web of a large number of intelligent sensing and computing nodes connected via a communication network (Dargie & Poellabauer, 2010). The emergence of sensor networks calls for the development of distributed algorithms for a number of tasks to replace the traditional centralized methods. In particular, the development of distributed algorithms for parameter estimation has recently attracted a great deal of attention (Fang & Li, 2008;Kar, Moura, & Ramanan, 2012;Li & AlRegib, 2007Lopes & Sayed, 2008;Ribeiro & Giannakis, 2006a,b;Ribeiro, Schizas, Roumeliotis, & Giannakis, 2010;Xiao, Ribeiro, Luo, & Giannakis, 2006). They find applications in environmental and (2008), Kar et al. (2012), AlRegib (2007, 2009), Lopes and Sayed (2008) and Ribeiro and Giannakis (2006a,b). On the other hand, in dynamic estimation all the nodes track the evolution of a set of parameters for which a dynamic model is available, as in the Refs. Carli, Chiuso, Schenato, and Zampieri (2008), Hlinka, Sluciak, Hlawatsch, Djuric, and Rupp (2012), Khan and Moura (2008) and Ribeiro et al. (2010). Some ''hybrid'' methods exist, which permit tracking a time-varying sequence of parameters, without a dynamic model, by updating a static estimate at each time (Lopes & Sayed, 2008). A final classification is done by whether a distributed estimation method is a small-scale one or a large-scale one. In a small-scale method, all nodes estimate a common set of parameters (Fang & Li, 2008;Kar et al., 2012;Li & AlRegib, 2007Lopes & Sayed, 2008;Ribeiro & Giannakis, 2006a,b). But in a largescale method, each node only reconstructs a subset of the parameters, i.e., the collective knowledge of the parameters is distributed among all the nodes (Conejo et al., 2007;Gómez-Expósito et al., 2011;Huang, Werner, Huang, Kashyap, & Gupta, 2012;Khan & Moura, 2008). Large-scale estimation is in general more challenging than its small-scale counterpart.
In this paper we study distributed static estimation for a large-scale system consisting of a network of interconnected subsystems. Each sub-system is concerned with a subset of the unknown parameters and has measurements, linear in the unknown parameters, corrupted by additive noise. The distributed estimation task is for each sub-system to estimate its local state using its own measurement and information shared with the network through neighborhood communication. We use the weighted least squares (WLS) criterion for optimal estimation. The goal is that the composite estimate of the whole system, consisting of all local estimates, will become globally optimal in the sense that it is the same as the optimal estimate obtained using all the measurements and a centralized estimation method. This problem is motivated by many applications involving a large-scale spatially distributed system. For example, the state estimation problem for a large power network is concerned with estimating the voltages and phases of the power supply at each sub-system, consisting of a number of buses or a substation, using measurements (provided by, for example, a phasor measurement unit (PMU) or a supervisory control and data acquisition (SCADA) system) (Conejo et al., 2007;Jiang, Vittal, & Heydt, 2008). Interactions of sub-systems are reflected by the fact that local measurements available at each sub-system typically involve neighboring sub-systems. For example, a current measurement at a conjunction depends on the voltage difference of two neighboring buses. In a smart grid setting, each sub-system is only interested in the local state, i.e., its own voltages and phases, using local measurements and information acquired from neighboring sub-systems via neighborhood communication (Tai, Marelli, Rohr, & Fu, 2013). For a large power network, it is both impractical and unnecessary for every sub-system to estimate the whole state of the system, thus distributed estimation methods for local state estimation are naturally called for. The second example is the localization problem for a wireless sensor network, involving estimating the locations of all sensors using relative measurements (e.g., relative distances or relative positions) between sensors (Diao, Fu, & Zhang, 2013;Khan, Kar, & Moura, 2009). For a small sensor network with a few sensing nodes, it is possible to aggregate all the measurements at an FC to compute a global estimate of all sensor locations. Such an algorithm is not scalable, and would require massive computing resources when the network becomes large. It is also unnecessary for each sensor to localize other nodes. A distributed method is preferred, where each node aims to estimate its own location using local measurements and neighborhood communication. The third example is a traffic network for a city or a highway system, where each node or sub-system wants to estimate its local state (e.g., traffic flow rates, delays, etc.). Due to the spatial correlations of the traffic flows in different sub-systems, neighboring traffic information is certainly useful in predicting the traffic conditions at each sub-system. Again, distributed estimation methods are preferred over centralized methods. Many other examples in sensor networks can be found in, e.g., Ribeiro et al. (2010), Kar et al. (2012) and Yang et al. (2010).
We first provide a fully distributed iterative algorithm for each node to compute its own local estimate. This algorithm works for a general connected network. Contrary to the method proposed in Conejo et al. (2007), our algorithm guarantees that the composite estimate will converge asymptotically to the global WLS estimate. We then focus on the convergence rate of the algorithm. Since our method is based on Richardson's method for solving systems of linear equations (Bertsekas & Tsitsiklis, 1997), its convergence rate depends on certain scaling parameter and a preconditioning matrix. Choosing the optimum scaling parameter requires the knowledge of the largest and the smallest eigenvalues of certain positive definite matrix (the estimation error covariance). A distributed algorithm for estimating these values can be obtained using the power method (Bertsekas & Tsitsiklis, 1997). However, to prevent numerical instability, this approach requires periodically executing a normalization step, which needs to be carried out in a distributed manner. In Yang et al. (2010), this is done using average consensus. A drawback of this approach is that consensus itself converges asymptotically. This significantly slows down the convergence of the eigenvalue estimation. To avoid this, we propose a different method in which normalization is done locally, at each node, without inter-node communication. In this way, the optimal scaling parameter can be obtained using a distributed method. We then address the problem of designing the preconditioning matrix. Our distributed scenario constrains us to use a block diagonal matrix. Choosing the optimal matrix under this constraint, and using only distributed processing, is a very challenging problem for which we are not able to provide a solution. Nevertheless, we are able to bound the difference between the convergence rate achieved using the optimal matrix, and the convergence rate resulting using a practically feasible matrix design. This bound turns out to have a simple expression which depends on the network connectivity. A shortened version of this method appears in the conference paper .
For an acyclic network (i.e., its communication graph contains no loops), we provide a different iterative algorithm for distributed estimation. As opposed to the previous algorithm, in this one, the composite estimate is guaranteed to converge to the globally optimal estimate in a finite number of steps. Indeed, we show that the convergence time equals the diameter of the aforementioned graph. Numerical experiments show that this method offers a major reduction in convergence time.
The rest of paper is organized as follows. In Section 2 we describe the distributed WLS estimation the problem. In Section 3 we derive the first distributed WLS method, which converges asymptotically. In Section 3.1 we describe distributed methods for finding the value of the scaling parameter which yields the fastest convergence rate. In Section 3.2 we describe a sub-optimal choice for preconditioning matrix, and we bound its sub-optimality. In Section 4 we introduce the second distributed WLS method which converges in finite time. Numerical experiments are presented in Section 5, and concluding remarks are given in Section 6. For the ease of readability, all proofs are contained in the Appendix.
Notation 1. For a vector x, ∥x∥ denotes its 2-norm, and [x] i denotes its ith entry. For a matrix X , ∥X∥ denotes its operator (induced) norm. For square symmetric matrices X and Y , X < Y means that matrix Y − X is positive-definite. For vectors and matrices, the superscript T denotes its transpose, and * denotes its transpose conjugate.

Problem description
Consider a network formed by I nodes. For each i = 1, . . . , I, Node i has an associated parameter vector x i ∈ C d i , and measures the vector y i ∈ C p i , which is given by where v i ∼ N (0, R i ) denotes the measurement noise. The noises v i and v j are statistically independent, whenever i ̸ = j.
..,I . Then, we can write the measurement model of the whole network as with v ∼ N (0, R). The WLS estimatorx of x is defined by Kay (1993, Eq. (8.14)) Its solution is given bŷ For the WLS problem to be well defined, we make the following assumption: Assumption 2. Matrix A has full column rank and R is nonsingular.
Computing (3) requires global information, i.e., all the measurements and the information about A and R need to be made available together. Our goal is to derive distributed methods in which Node i computes the componentx i of the estimatex, corresponding to x i , using only the local measurement y i and information received from its neighbors (a formal definition of neighborhood will be given later).
In the rest of the paper we use the following notation: denote the set of nodes whose measurements involve the parameters of Node i, and denote the set of nodes whose parameters are involved in the measurements of Node i.
to be the set of nodes whose neighborhood have a non-empty intersection with that of Node

Asymptotic method for WLS estimation
The distributed WLS method derived in this section uses the following definition of neighbor node.

Definition 4. Node j is a neighbor of Node
Also, the proposed method requires the following connectivity assumption: Assumption 5. For each i = 1, . . . , I, Node i can send/receive information to/from all its neighbors. Also, A i,j , for all j ∈ O i , and R i are available at Node i.
Although our method works regardless of whether the network is sparse or not, it works most efficiently for sparse networks. A network is called sparse if the cardinality of N i is small for all i = 1, 2, . . . , I.

Consider any block diagonal positive definite matrix
From (6), 0 < γ Υ < 2I. Hence, −I < γ Υ − I < I and therefore ∥I − γ Υ ∥ < 1. In view of this, we can use Richardson's method (Bertsekas & Tsitsiklis, 1997) to computex recursively. This leads tõ Then, by substituting the expressions ofα andx, we obtain straightforwardlŷ We call Π the preconditioning matrix, because, as it will be explained in Section 3.2, it is used to increase the convergence rate of the recursions (8). Let Also, for i, j = 1, . . . , I, let Then, from (8), (11) and (9), we obtain Notice that the matrices Ψ (k) i,j are only available at Node k. That is, Node k acts as an intermediary between Node j, which transmits This means that node j should transmitx j (t) to all nodes k with j ∈ O k , or equivalently, to all nodes in I j . However, Node j does not know which nodes are in I j . Thus, Node j simply transmitsx j (t) to all nodes in N j , and it is up to the receiving Node k to detect whether Node j ∈ O k or not.
Following the discussion above, we obtain the following algorithm to implement (12).

Algorithm 1 -distributed WLS estimation: Initialization:
(1) For each k = 1, · · · , I and i ∈ O k , Node k computes α (k) i and sends it to Node i.
(2) On reception, Node i computes Main loop: At time t ∈ N: (1) For each j = 1, · · · , I and k ∈ N j , Node j sends its current estimatex j (t) to Node k.
(2) On reception, for each k = 1, · · · , I and i ∈ O k , Node k sends to Node ǐ (3) On reception, for each i = 1, · · · , I, Node i computeŝ To implement Algorithm 1, we need to design the scaling factor γ and the preconditioning matrices Π i , for all i = 1, . . . , I. We address these two tasks in Sections 3.1 and 3.2, respectively.

Distributed design of the scaling factor γ
In this section we study two approaches for designing the scaling factor γ . In Section 3.1.1 we describe a distributed algorithm which converges asymptotically to the optimal value of γ , i.e., the resulting γ will achieve the maximum convergence speed. In Section 3.1.2, we give another distributed algorithm which converges in finite time to a sub-optimal value of γ .

Asymptotic algorithm for γ
In view of (7), the value of γ that maximizes the convergence rate is because this is the value that minimizes ∥I − γ Υ ∥. In order for each node to find the value of γ given in (13), we need distributed methods for finding ∥Υ ∥ and   Υ −1   −1 . We give these methods below. These methods yield, at Node i and time step t, estimates Υ i (t) and Υ i (t), of ∥Υ ∥ and   Υ −1   −1 , respectively. Then, at the same node and time step, the estimate γ i (t) of γ is obtained by Distributed method for finding ∥Υ ∥: Choose any vector b(0) ̸ = 0 and let b(t) = Υ t b(0). Then, using (11) with Π 1/2 b(t) in place ofx(t), we obtain at Node i, where b i (t) denotes the ith block component of b(t). Then, using the power method (Bertsekas & Tsitsiklis, 1997), Node i can asymptotically compute ∥Υ ∥ as follows An inconvenience with the approach above is that b(t) either increases or decreases indefinitely. To avoid this, the vector b(t) can be periodically normalized. In Yang et al. (2010), this was done using average consensus (in the continuous-time case). As we mentioned in Section 1, we avoid the drawbacks of that method by providing an alternative approach in which b(t) is normalized at each node, and each iteration, without inter-node communication. This algorithm is given below: (2) On reception, for each k = 1, · · · , I and i ∈ O k , Node k sends b (3) On reception, for each i = 1, · · · , I, Node i computes The convergence of Algorithm 2 to ∥Υ ∥ is guaranteed by the next theorem.
Theorem 6. Consider the network (1) together with Assumptions 2 and 5. Then, for each i ∈ {1, . . . , I}, Distributed method for finding where, for a symmetric matrix X , eig (X) and eig (X) denote the smallest and largest eigenvalues of X , respectively. Hence, we can find   Υ −1   −1 by applying Algorithm 2 on Φ, to find ∥Φ∥. To this end, at Node i and time t, we choose c = Υ i (t). This leads to the following algorithm: and (17) replaced by

Finite-time algorithm for γ
A sub-optimal design of γ can be achieved using the following result.
Theorem 7. Condition (6) is satisfied by choosing γ so that The design of γ using Theorem 7 requires the global information max i φ i . For each, i = 1, . . . , I, Node i can obtain φ i from an initialization stage in which it receives υ i,k , from each Node k, with k ∈ I i . Then, max i φ i can be obtained by running the maxconsensus algorithm (Olfati-Saber & Murray, 2004), in parallel with the estimation Algorithm 1. Notice that the max-consensus algorithm converges in finite time.

Design of the preconditioning matrix Π
As mentioned above, for a given choice of Υ , the fastest convergence rate of Algorithm 1 is achieved when γ is chosen as in (13). Under this choice of γ , we have that where we recall thatx denotes the global estimate of x, given by (3). Then, we define the time constant τ (Υ ) of the distributed WLS algorithm by Hence, a natural question is whether the preconditioning matrices Π i , i = 1, . . . , I, can be chosen so that τ (Υ ) is minimized. While we are not able to answer this question, we have the following result, which follows using an argument similar to the one in Demmel (1983, Theorem 2).
with P denoting the set of positive definite block diagonal matrices of the form (4).
Theorem 8 states that, if the preconditioning matrices Π i , i = 1, . . . , I, are chosen as then κ(Υ ) is at most β times bigger than the smallest possible value κ ⋆ achievable using block diagonal preconditioning matrices.

Remark 9.
In view of (20) and (10), Hence, its computation requires the matrices Ψ (k) i,i , k ∈ I i , to be transmitted from Node k to Node i during an initialization stage.

Finite-time method for WLS estimation
In this method we replace the definition of neighborhood by the following one: Consequently, we replace the connectivity Assumption 5 by the following one: Assumption 11. For each i = 1, . . . , I, Node i can send/receive information to/from all its neighbors. Also, Ψ j,i , for all j ∈ B i , and α i are available at Node i.
To illustrate the idea behind the proposed algorithm, we consider a network with two nodes. The next lemma states how to obtain the global optimal solution, at each node, in this simple case.
Lemma 12. Consider the network (1) together with Assumption 2. If there are only two nodes, labeled by a and b, then Ψ a,a − Ψ a,bΣb Ψ b,a is an invertible matrix and the global estimatex a of the components x a associated to Node a is given bŷ Our next result is an immediate generalization of the one above, to a network with a star topology, i.e., in which all nodes are only possibly connected to a single one.

Main loop:
For each i = 1, · · · , I, and time t ∈ N (1) Node i computes, for each j ∈ B i \ {i}, to Node j.
(2) Node i computeš Our next step is to show that Algorithm 4 converges in finite time to the global WLS solution.
Definition 14. Each pair (i, j), i, j ∈ {1, . . . , I}, is called an edge if Ψ i,j ̸ = 0. A path is a concatenation of contiguous edges, and its length is the number of edges forming it. For each i, j ∈ {1, . . . , I}, the distance d i,j between Nodes i and j is defined as the minimum length of a path joining these two nodes. The radius ρ i of Node i is defined as the maximum distance between Node i and any other node in the network. The diameter of the network is the maximum radius between all its nodes. A network is called acyclic if it does not contain a path forming a cycle.
The next theorem states that, if the network is acyclic, then the algorithm above yields the global estimate at each node in finite time.

Cyclic network topology
In the first simulation we cluster the buses into eight nodes, as shown in Table 1. From the definition of N i , it follows that j ∈ N i Table 1. if there is a bus in either, Node i or j, with a PMU installed, having an attached line coming from a bus inside the other node. Fig. 2 shows the topology of the communication network induced by the clustering given in Table 1. Fig. 3 shows the convergence of the asymptotic method without preconditioning. To this end, we show the modulus of the estimated voltage of each bus at each step. We see that the convergence is very slow, with a relative difference of r  10 6  = −37 dB, between the global estimate and the one obtained by the distributed algorithm at t = 10 6 . The reason for the slow convergence is that the condition number of Ψ is 478 972. The preconditioning matrix in (20) gives a condition number of 700, which leads to a much faster convergence. This is shown in Fig. 4, where r  2 × 10 3  = −52.47 dB. Fig. 5 shows that the convergence of the estimation of ∥Υ ∥ and   Υ −1   −1 , at each node, is much faster than that of the WLS estimation algorithm. Finally, Table 2 shows the complexity at each node. To this end, we measure the number of multiplications in a whole cycle of Algorithms 1-3.

Fig. 7.
Convergence of the asymptotic method (with preconditioning) in an acyclic network.

Acyclic network topology
In the second simulation we do the clustering such that the induced topology is acyclic. From the definition of B i , it follows that j ∈ B i if there is a bus (possibly neither in Node i nor in j), with a PMU installed, having one neighbor bus (i.e., a bus connected to it via an attached line), including possibly itself, in each node, i and j. The clustering and its induced topology are shown in Table 3 and Fig. 6, respectively. The convergence of the asymptotic method (with preconditioning) is shown in Fig. 7, with a final relative difference with the global estimate of r  20 × 10 3  = −69.67 dB. In this case, we wait for 100 steps before starting with the distributed estimation. The gap caused by this delay can be seen at the beginning of the graph. We introduced this late start so as to give time for Algorithm 2 and 3 to obtain reasonable approximations of ∥Υ ∥ and   Υ −1   −1 , respectively. We see that the asymptotic method presents an oscillating behavior between time steps 1500 and 3500. This is because the transients in the estimation of the scaling factor γ (t) cause the recursions (8) to become temporarily unstable. We also see that the asymptotic method requires about 20 × 10 3 steps to converge. This is because in this case, preconditioning leads to a  condition number of 5264. On the other hand, the convergence of the finite-time method does not depend on the condition number, but on the network diameter, which in this case is four. Fig. 8 shows the convergence of this method in four steps, with a final error of r (4) = −223.7 dB, caused by numerical inaccuracy. Table 4 shows the complexity at each node. To this end, we consider that solving the positive-definite linear system for computingx i,j (t), using the Cholesky decomposition, requires n 3 /3 + 2n 2 multiplications (n 3 /3 for the decomposition and 2n 2 for solving two triangular linear systems). Also, computing the inverse of the matrixΣ i (t), using also the Cholesky Decomposition, requires n 3 /2 multiplications (Krishnamoorthy & Menon, 2011).

Sensor localization
Sensor localization refers to the problem of obtaining the locations of each node in a network, based on the knowledge of the locations of a few anchor nodes, a well as the mutual distances between neighbor nodes. A distributed method for carrying out this task is proposed in Khan et al. (2009). This method requires that, for each i = 1, . . . , I, Node i lies inside of at least one triangle defined by three of its neighbors N i = {j, k, l}. Then, the coordinates x i of Node i can be written as where the barycentric coordinates c i,j are given by with S(i, j, k) denoting the area of the triangle formed by Nodes i, j and k. The latter can be computed using the Cayley-Menger determinant as follows Fig. 9. Node positions and estimates. be written as or equivalently, with y = (D ⊗ I 2 ) a and A = I − C ⊗ I 2 . Due to inaccuracy in distance measurements, (26) can be approximately expressed as in (2). In that case, we can use our proposed distributed method to obtain, at each node, a WLS estimation of its coordinates. The experiment setup is shown in Fig. 9. It includes three anchor nodes, defining a triangle containing I = 20 randomly placed nodes. We use a noise covariance matrix R = σ 2 I d , where I d denotes the identity matrix, and σ = √ 10 −3 ≃ 31.62 centimeters.
With this setup, the global estimatex yields a relative localization error of e = −33.39 dB, defined as in (23). The convergence of the coordinate estimates at each node, using the proposed method, with preconditioning, is shown in Fig. 10. As before, we wait for 10 steps before starting the iterations, to give time for Algorithms 2 and 3 to obtain reasonable approximations of ∥Υ ∥ and   Υ −1   −1 , respectively. The convergences of these estimates are shown in Fig. 11. Finally, the complexity at each node is shown in Table 5. For comparison, we also consider the distributed iterative localization algorithm (DILOC) proposed in Khan et al. (2009). This Table 5 Complexity at each sensor node, in a number of multiplications per iteration .   Node  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20   Complexity  170  282  282  426  170  282  602  602  282  426  282  426  282  170  170  426  426  426  426

method solves (26) using Richardson's recursions to invert matrix
A. This requires that N = 1, i.e., only one equation of the form (25) is considered for each node. In this case, the recursions are guaranteed to converge because, as the authors show, ∥I − A∥ < 1 holds in this problem. Fig. 12 shows the evolution of the relative difference r(t) (defined as in (24)) between the estimates of each method, and the global estimatex. We see that the DILOC method has a faster convergence. This is because the condition number of A is smaller than that of A * R −1 A, which is the matrix inverted by our proposed method. However, at t = 500, the DILOC method yields r(500) = −29.44 dB, while the proposed one gives r (500) = −72.71 dB. This difference results from the fact that the DILOC method does not produce the WLS solution on the limit. 1

Conclusion
We proposed two methods for weighted least squares estimation in large-scale systems. Both methods converge to the global solution and aim to maximize the convergence speed. The first method converges asymptotically and involves a distributed estimation of the scaling parameter upon which the convergence 1 Notice that, in the scenario considered in this work, noisy inter-node distances are only measured once, and they remain unchanged during the whole iteration process. This is in contrast to the scenario considered in Khan et al. (2009), where these distances are re-measured at each iteration. speed depends. To further speed up the convergence, we also use a practically feasible preconditioning method, for which we bounded the speed difference with respect to the fastest theoretically achievable. The second proposed method has an even faster convergence, as it achieves the global optimal in finite time. However, it is only suitable for applications where the graph produced by the communication network contains no loops.

Appendix A. Proofs of Section 3
A.1. Proof of Theorem 6 Let k i (0) = 1 and We need to design k i (t + 1), or equivalently ς i (t + 1), to avoid the indefinite increase or decrease of b(t). In principle, this could be achieved by choosing However, the question then arises as to whether some of the scalars υ i,j (t) would grow to infinity. Notice that Hence, this could only happen if some vector in the eigenspace associated to the largest eigenvalue of Υ has zero components in the entries corresponding to b j (t). We call a matrix satisfying this property, ill-posed. Although the set of ill-posed matrices is nowhere dense, (i.e., it is unlikely to have an ill-posed matrix Υ ), we can avoid the indefinite growth of υ i,j (t) by choosing ς i (t + 1) so that This leads to From (15) and (A.1), the estimateῩ i (t) of ∥Υ ∥ at t is However, if Υ is ill-posed,  b i (t)   will tend to zero. In such case, Υ i (t) can be computed bȳ for some neighbor node j for which  b j (t)   does not tend to zero.
Notice that such a neighbor always exists, for otherwise Node i would be isolated from all other nodes.

A.2. Proof of Theorem 7
We need use the following result.
and the result follows.
Proof of Theorem 7. Using Lemma 16 we have Hence, and the result follows.

A.3. Proof of Theorem 8
We need the following lemma.
for any y 0 ̸ = 0. Let x 0 and y 0 have unit norm and be such that D −1/2 x 0 = σ  D −1/2  x 0 and Υ 1/2 y 0 = σ  Υ 1/2  y 0 . Then, Now, since D −1/2 is block diagonal, x 0 can be chosen so that its nonzero components correspond to only one block of D −1/2 . Let x 0,b denote the entries of x 0 in that block, and Υ 1/2 b denote the columns of Υ 1/2 corresponding to the same block. Then, The result follows since the inequality above holds for anyΠ .