1

Distributed Low-Rank Adaptive Algorithms Based on
Alternating Optimization and Applications

arXiv:1710.00617v1 [] 29 Sep 2017

Rodrigo C. de Lamare .

Abstract—This paper presents a novel distributed low-rank
scheme and adaptive algorithms for distributed estimation over
wireless networks. The proposed distributed scheme is based
on a transformation that performs dimensionality reduction
at each agent of the network followed by transmission of a
reduced set of parameters to other agents and reduced-dimension
parameter estimation. Distributed low-rank joint iterative estimation algorithms based on alternating optimization strategies
are developed, which can achieve significantly reduced communication overhead and improved performance when compared
with existing techniques. A computational complexity analysis of
the proposed and existing low-rank algorithms is presented along
with an analysis of the convergence of the proposed techniques.
Simulations illustrate the performance of the proposed strategies
in applications of wireless sensor networks and smart grids.
Index Terms—Dimensionality reduction, distributed estimation
techniques, low-rank algorithms, wireless sensor networks, smart
grids.

I. I NTRODUCTION

D

ISTRIBUTED strategies have become fundamental for
parameter estimation in wireless networks and applications such as sensor networks [1], [2], [3] and smart grids
[4], [5]. Distributed techniques deal with the extraction of
information from data collected at nodes that are distributed
over a geographic area [1]. In this context, a specific sensor
node or agent in the network collects processed data from
its neighbors and combines them with its local information
to generate improved estimates. However, when the unknown
parameter vector to be estimated has a large number of parameters, the network requires a large communication bandwidth
between neighboring nodes to transmit their local estimates.
This problem limits the usefulness of existing algorithms in
applications with large data sets as the convergence speed is
dependent on the number of parameters [2], [6], [7]. Hence,
distributed dimensionality reduction has become an important
tool for distributed inference problems with large data sets.
In order to perform dimensionality reduction or compression, several algorithms have been proposed in the literature
in the context of distributed quantized Kalman filtering [8],
[9], quantized consensus algorithms [10], distributed principal
subspace estimation [11], the single bit strategy [12] and
Krylov subspace optimization techniques [13]. However, these
distributed algorithms [8]-[13] have drawbacks such as high
computational complexity, unsatisfactory performance and
implementation issues. Available distributed approaches for
dimensionality reduction or compression [8]-[13] have tradeoffs between the amount of cooperation, communication and
R. C. de Lamare is with CETUC / PUC-Rio, Brazil and the Department of
Electronics, University of York, U.K. (e-mail: rodrigo.delamare@york.ac.uk)

system performance. This calls for the development of costeffective techniques that can approach the performance of theoretical bounds for parameter estimation, have flexibility and
high-compression capability, and exhibit low computational
complexity. In this context, low-rank techniques are powerful
tools to perform dimensionality reduction, which have been
applied to spread-spectrum systems [14], [15], [16], [17], [18],
[19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29],
[30], [31], multi-input-multi-output (MIMO) systems [32],
[33], [34], [35], [36], [37], [38], [39], [40] and beamforming
[41], [42], [43]. However, limited research has been carried out
on distributed low-rank estimation, in which the distributed
principal subspace estimation [11] and the Krylov subspace
optimization [13] techniques are recent contributions. Related approaches to low–rank techniques include compressive
sensing-based strategies [44],[45], which exploit sparsity to
reduce the number of parameters for estimation, distributed
dictionary learning [46], [47], [48], [49], [50], [51], which employs a bilinear dimensionality-reduction factorization scheme
similar to some low-rank schemes but assumes no regression
vectors, and attribute-distributed learning [52], which employs
agents and a fusion center to meet communication constraints.
Another important tool in recent related work is the principle
of alternating optimization [53], [54], which consists of fixing
a set of parameters, adjusting the remaining parameters and
then proceeding in cycles [18], [35], [55], [56], [57], [58]. Ling
and Ribeiro have studied dynamic decentralized optimization
using the alternating direction method of multipliers [55].
Bai et al. have examined alternating optimization procedures
to design sensing matrices and dictionaries for compressive
sensing. Yan et al. have developed an alternating optimization
for multigraph matching, whereas Magnusson et al. [58] have
studied convergence of nonconvex optimization problems.
In this paper, we propose a scheme for distributed signal
processing and distributed low–rank algorithms for parameter
estimation. In particular, the proposed algorithms are based on
an alternating optimization strategy [53], [54], [18], [35] and
are called the distributed reduced-rank joint iterative optimization normalized least mean squares (DRJIO–NLMS) algorithm
and the distributed reduced-rank joint iterative optimization
recursive least squares (DRJIO–RLS) algorithm. In contrast to
prior work on low-rank techniques [18]-[43] and distributed
methods [8]-[13], [59], [60], [61], [62], [63], [64], [65],
[66], distributed adaptive techniques based on the alternating
optimization strategy are investigated. The proposed low-rank
strategies are distributed and perform dimensionality reduction
without costly decompositions at each agent. The proposed
DRJIO–NLMS and DRJIO–RLS algorithms are flexible with
regards to the amount of information that is exchanged, have

2

low cost and high performance. We also present a computational complexity analysis of the proposed and existing lowrank algorithms along with an analysis of the convergence of
the proposed techniques. Applications to parameter estimation
in wireless sensor networks and smart grids are then studied.
The main contributions of this work can be summarized as:
• Distributed low-rank adaptive algorithms based on alternating optimization .
• An analysis of the convergence and the computational
complexity of the proposed distributed algorithms.
• A study of the proposed and existing distributed algorithms in wireless sensor networks and smart grids.
This paper is organized as follows: In Section II, the system
model and the problem statement are described. In Section
III, the proposed distributed dimensionality reduction and
adaptive processing scheme is presented. Section IV details
the proposed distributed low-rank algorithms. In Section V,
an analysis of the convergence of the proposed algorithms
is carried out along with a study of their computational
complexity. Simulation results are presented and discussed in
Section VI, whereas conclusions are drawn in Section VII.
II. S YSTEM M ODEL AND P ROBLEM S TATEMENT

20

14

12
11

18

16

10
4

14

18

20

J(ω k (i)) =

N
X

2

E dk (i) − ω k H (i)xk (i) ,

(2)

k=1

where E denotes expected value and ωk H (i) is the estimator
at time i. To solve this problem, one suitable technique is the
adapt–then–combine (ATC) diffusion strategy [2] described by

∗

H

 ψ k (i) = ωk (i − 1) + µk xk (i) dk (i) − ωk (i − 1)xk (i) ,
P

ckl ψ l (i),
 ω k (i) =
l∈Nk

(3)
where µk is the step size, Nk indicates the set of neighbors
for node k, ψ k (i) is the local estimator, |Nk | denotes the
cardinality of Nk and ckl > 0 are the combination coefficients,
which are calculated in this work using the Metropolis rule [2]
given by

if k 6= l are linked
c = max(|N1k |,|Nl |) ,

 kl
ckl = 0,
for k and l not linked
P

ckl ,
for k = l
 ckk = 1 −
l∈Nk /k

(4)

12

10

and should satisfy

9

15

7

X

8

19

8
13

4

3

6
17

2

16
0

2

4

6

2

1
8

10

12

5
14

16

18

20

Fig. 1. Network topology with N nodes

A distributed network with N nodes, which have limited
processing capabilities, is considered with a partially connected topology as illustrated in Fig. 1. A diffusion protocol
in which nodes from the same neighborhood communicate
with each other at every iteration is employed [2], although
other strategies such as incremental [1] and consensus-based
[4] could also be used. A partially connected network means
that nodes can exchange information only with their neighbors
determined by the connectivity topology. In contrast, a fully
connected network means that, data broadcast by a node can
be captured by all other nodes in the network [67]. At every
time instant i, each node k takes a scalar measurement dk (i)
according to
dk (i) = ω H
0 xk (i) + nk (i),

ckl = 1, l ∈ Nk ∀k.

(5)

l

6

0

where xk (i) is the M × 1 input signal vector with zero
2
mean and variance σx,k
that is also observed by node k,
nk (i) is the noise sample measured at node k which has
2
zero mean and variance σn,k
. Observing (1), we can see that
the measurements for all nodes are related to an unknown
parameter vector ω0 with size M × 1, that would be estimated
by the network. The aim of such a network is to compute an
estimate of ω 0 in a distributed fashion, which can minimize
the global cost function

i = 1, 2, . . . , I,

(1)

Note that other combination rules can also be employed.
With this adaptation strategy, when the dimension of the
unknown parameter vector ω0 is large, this could lead to a
high communication overhead between each neighbor node
and the learning speed of the network is reduced. In order to
reduce the communication overhead, accelerate the learning
and optimize the distributed processing, we incorporate at the
kth node of the network distributed low-rank strategies based
on alternating optimization techniques.
III. D ISTRIBUTED D IMENSIONALITY R EDUCTION
A DAPTIVE P ROCESSING

AND

The proposed distributed dimensionality reduction scheme,
depicted in Fig.2, employs a transformation matrix S Dk (i)
to process the input signal xk (i) with dimensions M × 1
and projects it onto a lower D × 1 dimensional subspace
x̄k (i), where D ≪ M . Following this procedure, a low-rank
estimator ω̄k (i) is computed, and the ω̄ k (i) is transmitted
by each node. In particular, the transformation matrix S Dk (i)
and low-rank estimator ω̄ k (i) will be jointly optimized in the
proposed scheme according to the mean squared error (MSE)
criterion.

3

xk (i)

x̄k (i)

Dimensionality reduction
matrix SDk (i)

dk (i)

Reduced-Rank
estimator ω̄k (i)

−

+
H
H
R̄k (i) = E[S H
Dk (i)xk (i)xk (i)S Dk (i)] = E[x̄k (i)x̄k (i)]
(13)
and the cross-correlation vector is given by

D×1

M ×1

ek (i)

Design Algorithm

∗
p̄k (i) = E[d∗k (i)S H
Dk (i)xk (i)] = E[dk (i)x̄k (i)].

Fig. 2. Proposed dimensionality reduction scheme at each node or agent

Specifically, we start the description of the method with an
M × D matrix S Dk (i), which carries out a dimensionality
reduction on the input signal of each agent as given by
x̄k (i) = S H
Dk (i)xk (i),

(6)

min

S Dk (i),ω̄ k (i)

N
X

H

The full-dimension estimator is then obtained through a rankD approximation:

k=1

+ γE[

D
X

H
eH
d S Dk (i)I M,D ed ]

d=1

(7)
∗

where denotes complex conjugation, ω̄k (i) is the low-rank
estimator, the M × 1 vectors ed contain one in the dth entry
and zeros elsewhere, the M × D matrix I M,D contains a Ddimensional identity matrix on the top and zeros elsewhere,
and the parameters γ and δ are regularization terms that ensure
the solution has rank D. In what follows, we describe the
adaptation step that computes the parameters of ω̄k (i) and
S Dk (i) based on an alternating minimization strategy, which
consists of fixing one set of parameters and then minimizing
the other.
By fixing ω̄ k (i) and minimizing (7) with respect to S Dk (i),
we arrive at the following expression:
−1
R−1
k (i)P Dk (i)R̄ω̄ k (i),

(8)

where the covariance matrix of the input signal vector at node
k xk (i) is assumed to be full-rank and is given by
Rk (i) = E[xk (i)xH
k (i)],

(9)

the cross-correlation matrix is given by
P Dk (i) = E[d∗k (i)(xk (i)ω̄ H
k (i) + γI M,D )]

(10)

and the covariance matrix of the reduced-rank parameter
vector is described by
R̄ω̄k (i) = ω̄k (i)ω̄ H
k (i) + δI D .

(11)

We then fix S Dk (i) and minimize (7) with respect to ω̄k (i),
which results in
−1

ω̄ k (i) = R̄k (i)p̄k (i),

(16)

2
(i)S H
]
Dk (i)xk (i)|which

2
+ δE[||d∗k (i)S H
Dk (i)xk (i)|| ]

S Dk (i) =

The reduced-dimension parameter vector ω̄k (i) computed at
each agent is then transmitted as a local low-rank estimator
ψ̄ k (i) to the other agents according to the network topology.
At the receiver of each agent, there is a combination and reconstruction step in which the received data from neighboring
nodes is combined to obtain a low-rank estimator:
X
(15)
ω̄ k (i) =
ckl ψ̄ l (i),

ω k (i) = S Dk (i)ω̄ k (i),
E[|dk (i) − ω̄ k

(14)

l∈Nk

where, in what follows, all D–dimensional quantities are
designated with an overbar. The design of S Dk (i) and ω̄ k (i)
corresponds to the optimization problem given by
 opt opt
=
S Dk , ω̄ k

where the covariance matrix of the reduced-rank input signal
vector that is also assumed to be full-rank is expressed by

(12)

is derived in the Appendix.
The associated low-rank MSE is obtained by substituting
the expressions obtained in (12) and (8) into the cost function
and is described by [18]
−1

MSE = σd2k − p̄H
k (i)R̄k (i)p̄k (i)

(17)

where σd2k = E[|dk (i)|2 ]. Because there is no closed-form
expression for S Dk (i) and ω̄k (i) as they depend on each other,
a strategy to compute the parameters is needed. The proposed
strategy is based on an alternating optimization of S Dk (i)
and ω̄ k (i). The rank D must be set by the designer to ensure
appropriate performance taking into account the bias-variance
tradeoff [15]. Furthermore, for the selection of D the reader is
referred to [68] for rank selection methods. In the next section,
we develop a distributed low-rank algorithm to compute the
parameters of interest.
IV. P ROPOSED D ISTRIBUTED L OW-R ANK A LGORITHMS
In this section, we present the proposed distributed low-rank
adaptive algorithms for distributed estimation, namely DRJIO–
NLMS and DRJIO–RLS. Unlike prior work [11], [12], [13],
the proposed algorithms do not require:
• Additional cost to perform eigen–decompositions [11]
• Extra adaptive processing at the local node [12]
• Multiple matrix-vector multiplications to to build the
Krylov subspace [13]
• Costly convex optimization at the local node, which
introduces extra complexity [13].
The objective of the DRJIO–NLMS and DRJIO–RLS algorithms is to perform compression/decompression and distributed parameter estimation subject to the constraint of
transmitting only D < M parameters. The algorithms are
flexible, have low cost, very fast convergence speed and
assume that the parameter D is given. Alternatively, a modelorder selection algorithm [16], [41], [69] can be employed to
compute D for each node. In particular, the algorithms rely on

4

an alternating optimization strategy which consists of fixing a
set of parameters S Dk (i), updating the other set of parameters
ω̄ k (i), then fixing ω̄ k (i) and updating S Dk (i). This alternating
approach with the recursions for S Dk (i) and ω̄ k (i) is carried
out in cycles until convergence is achieved.
A. Proposed DRJIO–NLMS algorithm
In the DRJIO–NLMS algorithm, the parameters in (7)
are optimized by an alternating procedure that adjusts one
of the parameters while keeping the other parameter fixed
using NLMS recursions. Therefore, the proposed DRJIONLMS algorithm solves the optimization problem in (7) in
an alternating fashion. Using the method of steepest descent,
computing the gradient terms of the cost function in (7)
with respect to S Dk (i), replacing the expected value with
instantaneous estimates and considering the recursions in an
alternating fashion, we arrive at the proposed DRJIO–NLMS
algorithm:
S Dk (i) = S Dk (i − 1) + η(i)e∗k (i)xk (i)ω̄ H
k (i − 1)


+ η(i) γd∗k (i)I M,D − δxk (i)xH
k (i)S Dk (i − 1) ,
(18)
ω̄ k (i) = ω̄ k (i − 1) + µ(i)e∗k (i)x̄k (i),

(19)

varying step sizes. The normalization makes the setting of the
convergence factors easier, improves the convergence speed
and facilitates the comparison with other distributed LMS–
type algorithms. The recursions are computed in an alternating
way with one iteration per time instant at each node.
The proposed DRJIO–NLMS algorithm includes two steps,
namely, adaptation step and combination and reconstruction
step, which are performed using an alternating procedure
which is detailed next.
• Adaptation step
For the adaptation step, at each time instant i=1,2, . . . , I,
each node k=1,2, . . . , N, starts from generating a local lowrank estimator through

H
dk (i)− ω̄ H
k (i−1)S Dk (i)xk (i).

(20)

where ek (i) =
This local lowrank estimator ψ̄ k (i) will be transmitted to all its neighboring
nodes under the network topology structure.
Then, each node k=1,2, . . . , N, will locally update its dimensionality reduction matrix according to (18) and keep it locally.
Note that S Dk (i) only employs the low-rank estimators from
neighboring nodes and the local reference signal dk (i).
• Combination and reconstruction step
At each time instant i=1,2, . . . , I, the combination and reconstruction step starts after the adaptation step. Each node will
combine the local low-rank estimators from its neighboring
nodes and itself through
X
(21)
ω̄ k (i) =
ckl ψ̄ l (i),
l∈Nk

to compute the low-rank estimator ω̄ k (i).

In conclusion, during the distributed processing steps, only the
local low-rank estimator ψ̄ k (i) will be transmitted through the
network. The proposed DRJIO–NLMS algorithm is detailed in
Table I.
TABLE I
T HE DRJIO–NLMS A LGORITHM
Initialize: ω̄ k (0) = 0, S Dk (0) = I M,D
For each time instant i=1,2, . . . , I
For each node k=1,2, . . . , N
ψ̄ k (i) = ω̄ k (i − 1) + µ(i)e∗k (i)x̄k (i)
H
where ek (i) = dk (i) − ω̄ H
k (i − 1)S Dk (i)x k (i)
% ψ̄ k (i) is the local low-rank estimator and will be
% sent to all neighboring nodes of node k under the network
% topology structure.
S Dk (i) = S Dk (i − 1) + η(i)e∗k (i)xk (i)ω̄ k (i − 1) 
+η(i) γd∗k (i)I M,D − δxk (i)xH
k (i)S Dk (i − 1)
% The dimensionality reduction matrix S Dk (i)
% will be updated and kept locally.
end
For each node k=1,2,
P ..., N
ckl ψ̄ l (i)
ω̄ k (i) =
l∈Nk

− 1)S H
where ek (i) = dk (i) −
Dk (i − 1)xk (i), µ(i) =
µ0
η0
are the time–
and
η(i)
=
H
ω̄ H
xH
k (i−1)ω̄ k (i−1)xk (i)xk (i)
k (i)xk (i)
ω̄ H
k (i

ψ̄ k (i) = ω̄ k (i − 1) + µ(i)e∗k (i)x̄k (i),

After the last iteration I, each node will reconstruct a full–
dimensional estimator ω k (I) through the rank-D approximation given by
ωk (I) = S Dk (I)ω̄ k (I).
(22)

% The low-rank estimator ω̄ k (i)
% will be updated and kept locally.

end
end
After the last iteration I
For each node k=1,2, . . . , N
% Reconstruction.
ω k (I) = S Dk (I)ω̄ k (I)
where ω k (I) is the final full–rank estimator.
end

B. Proposed DRJIO–RLS algorithm
In this subsection, we develop the DRJIO–RLS algorithm
for computing S Dk (i) and ω̄ k (i), which is inspired by the
derivation of the standard recursive least squares (RLS) algorithm. The main differences are that we have two sets of
recursions that update the parameters: one for S Dk (i), which
performs compression/decompression, and another for ω̄k (i),
which performs parameter estimation; and the recursions are
distributed and computed in an alternating fashion. Therefore,
we first fix ω̄ k (i) in the derivation and then derive a set of
RLS recursions to compute the parameters for S Dk (i). Subsequently, we fix S Dk (i) and derive a set of RLS recursions to
compute the parameters for ω̄ k (i). The DRJIO–RLS algorithm
consists of an adaptation step, which computes S Dk (i) and
ω̄k (i), and a combination and reconstruction step, which is
identical to that of the DRJIO-NLMS algorithm.
• Adaptation step
In order to derive the proposed algorithm, we first define
P k (i) , R−1
k (i),

(23)

P Dk (i) , λP Dk (i − 1) + d∗k (i)xk (i)ω̄ H
k (i),

(24)

−1

Qω̄ k (i) , R̄ω̄k (i − 1),

(25)

5

and rewrite the expression in (8) as follows
S Dk (i) =

−1
R−1
k (i)P Dk (i)R̄ω̄ k (i

TABLE II
T HE DRJIO-RLS A LGORITHM

− 1)

Initialize: ω̄ k (0)=0
P k (0) = δ−1 I M ×M , Qω̄ k (0) = δ−1 I D×D ,
Φ̄k (0) = δ−1 I D×D and δ = small positive constant
∗
H
λP k (i)P Dk (i − 1)Qω̄k (i) + dk (i)P k (i)xk (i)ω̄ k (i)Qω̄
(i)
k each time instant i=1,2, . . . , I

 For
For each node k=1,2, . . . , N
∗
H
H

= P k (i)P Dk (i)Qω̄ k (i)
=

= S Dk (i − 1) + kk (i) dk (i)tk (i) − xk (i)S Dk (i − 1) ,
(26)
where the D × 1 vector tk (i) = Qω̄k (i)ω̄ k (i) and the M × 1
Kalman gain vector is
λ−1 P k (i − 1)xk (i)
kk (i) =
.
1 + λ−1 xH
k (i)P k (i − 1)xk (i)

(27)

In addition, the update for the M × M matrix P k (i) employs
the matrix inversion lemma [6] as follows:
P k (i) = λ−1 P k (i − 1) − λ−1 kk (i)xH
k (i)P k (i − 1) (28)
and the D × 1 vector tk (i) is updated as
tk (i) =

1+

λ−1 Qω̄k (i − 1)ω̄ k (i − 1)
−1
λ ω̄ H
k (i − 1)Qω̄ k (i − 1)ω̄ k (i

kk (i) =

tk (i) =

λ−1 P k (i−1)xk (i)
1+λ−1 xH
(i)P k (i−1)xk (i)
k
λ−1 Qω̄ (i−1)ω̄ k (i−1)
k

1+λ−1 ω̄ H
(i−1)Qω̄ (i−1)ω̄ k (i−1)
k
k



H (i)S
(i
−
1)
(i)
−
x
S Dk (i) = S Dk (i − 1) + kk (i) d∗k (i)tH
D
k
k
k

P k (i) = λ−1 P k (i − 1) − λ−1 kk (i)xH
k (i)P k (i − 1)
Qω̄ k (i) = λ−1 Qω̄ k (i − 1) − λ−1 tk (i)ω̄ H
k (i − 1)Qω̄ k (i − 1)
k̄k (i) =

λ−1 Φ̄k (i−1)x̄k (i)
1+λ−1 x̄H
(i)Φ̄k (i−1)x̄
k
 k (i)


(i)ω̄
(i
−
1)
ψ̄ k (i) = ω̄ k (i − 1) + k̄k (i) d∗k (i) − x̄H
k
k

Φ̄k (i) = λ−1 Φ̄k (i − 1) − λ−1 k̄k (i)x̄H
k (i)Φ̄k (i − 1)
end
For each node k=1,2,
P ..., N
ckl ψ̄ l (i)
ω̄ k (i) =
l∈Nk

− 1)

.

(29)

The matrix inversion lemma [6] is then used to update the
D × D matrix Qω̄k (i) as described by
Qω̄ k (i) = λ−1 Qω̄k (i − 1) − λ−1 tk (i)ω̄ H
k (i − 1)Qω̄ k (i − 1).
(30)
Equations (23)–(30) constitute the key steps of the proposed
DRJIO-RLS algorithm for computing S Dk (i).
To derive the expression for updating ω̄ k (i), the following
associated quantities are defined
−1

Φ̄k (i) , R̄k (i)

(31)

p̄k (i) = λp̄k (i − 1) + d∗k (i)xk (i).

(32)

end
end
For each node k=1,2, . . . , N
% Reconstruction.
ω k (I) = S Dk (I)ω̄ k (I)
where ω k (I) is the final full–rank estimator.
end

V. A NALYSIS

OF THE PROPOSED ALGORITHMS

In this section, the computational complexity of the proposed algorithms is detailed and an analysis of sufficient conditions for convergence and a convergence proof to the optimal
low-rank estimator are developed. Regarding our convergence
proof, it is worth noting that the proof described in [53] was
performed for a non-adaptive scenario, whereas the work in
[54] was carried out for an adaptive setting. Our proof has
been developed for the distributed adaptive case.

Then, equation (12) will be rewritten as
A. Computational Complexity Analysis

−1

ω̄ k (i) = R̄k (i)p̄k (i)
= Φ̄k (i)p̄k (i)
= λΦ̄k (i)p̄k (i − 1) + d∗k (i)Φ̄k (i)xk (i)


∗
H
= ω̄ k (i − 1) + k̄k (i) dk (i) − x̄k (i)ω̄ k (i − 1) ,
(33)
where the D × 1 Kalman gain vector is given by
k̄k (i) =

λ−1 Φ̄k (i − 1)x̄k (i)
.
1 + λ−1 x̄H
k (i)Φ̄k (i − 1)x̄k (i)

(34)

and the update for the matrix inverse Φ̄k (i) employs the matrix
inversion lemma [6]
Φ̄k (i) = λ−1 Φ̄k (i − 1) − λ−1 k̄k (i)x̄H
k (i)Φ̄k (i − 1). (35)
Equations (31)–(35) are the key steps of the proposed DRJIORLS algorithm for computing ω̄ k (i). Since the combination
and reconstruction step is identical to that of DRJIO-NLMS we
omit it here. The proposed DRJIO–RLS algorithm is detailed
in Table II.

Here, we evaluate the computational complexity of the proposed DRJIO–NLMS and DRJIO–RLS algorithms. The computational complexity of the proposed DRJIO–NLMS algorithm is O(DM ), while the proposed DRJIO–RLS algorithm
has a complexity O(M 2 + D2 ), where O(·) is used to classify
algorithms according to how their requirements in arithmetic
operations grow as the input size grows. The distributed NLMS
algorithm [2] requires O(M ), while the complexity of the
distributed RLS algorithm [70] is O(M 2 ). For the Krylov
Subspace NLMS [13] the complexity reaches O(DM 2 ), while
for the distributed principal subspace estimation algorithms
[11], the complexity is O(M 3 ). Thus, the proposed DRJIO–
NLMS algorithm has a much lower computational complexity,
and because we consider D ≪ M , it has a comparable cost
to the distributed NLMS algorithm [2]. The computational
complexity of the model-order selection algorithm of [69] with
extended filters requires 3(Dmax − Dmin ) + 1 additions and a
sorting algorithm to find the best model order. An additional
and very important aspect of distributed low-rank algorithms
is that the dimensionality reduction results in a decrease in

6

The proposed algorithms can be stated as alternating minimization strategies performed in a distributed fashion and
expressed as


,
ω̄
(i)
for k = 1, 2, . . . , N,
Dk S opt
S Dk (i) ∈ arg optmin
k
Dk

6

10

5

Number of Multiplications

10

DRJIO−NLMS
DRJIO−RLS
Distributed NLMS[2]
Distributed RLS[26]
Krylov Subspace NLMS[13]
Distributed Principal
Subspace Estimation[11]

S D ∈S D (i)

(37)
where Dk (·) is a distance metric and


opt
ω̄k (i) ∈ arg optmin Dk S Dk (i), ω̄ k
for k = 1, 2, . . . , N

4

10

ω̄ k ∈ω̄ k (i)

3

10

2

10

20

30

40

50

60
M

70

80

90

100

Fig. 3. Complexity in terms of multiplications

the number of transmitted parameters from M to D, which
corresponds to a less stringent bandwidth requirement.
The details of the computational complexity of the the
proposed and the existing algorithms, are shown in Table III,
where M is the length of the unknown parameter that needs
to be estimated, D is the reduced dimension and |Nk | is
the cardinality of Nk . To further illustrate the computational
complexity for different algorithms, we present the main trends
in terms of the number of multiplications for the proposed and
existing algorithms in Fig. 3. For the parameters, we consider
a network with N = 20 nodes, take node 14 as an example,
and set D = 5 and |Nk | = 5.

To start the analysis, we assume that the transformation matrix S Dk (i) that performs compression and decompression/reconstruction and the low-rank parameter estimator
ω k (i) aim to estimate the optimum pair S Dk,opt and ω̄ k,opt (i)
containing a common set of parameters of interest in the
network. Then, to develop the analysis and proofs, we need
to define a metric space and the Hausdorff distance that will
extensively be used. A metric space is an ordered pair (M, r),
where M is a nonempty set, and r is a metric on M, i.e., a
function r : M × M → R such that, for any x, y, z, and M,
the following conditions hold.
r(x, y) ≥ 0.
r(x, y) = 0 if x = y.
r(x, y) = r(y, x).
r(x, y) ≤ r(x, y) + r(y, z) (triangle inequality).

i→∞

To present a set of sufficient conditions under which the
proposed algorithms converge, we employ the so–called three–
and four–point properties [53], [54], [71], which are used in
the study of the theory of convex sets. Let us assume that
there is a function f : M × M → R such that the following
conditions are satisfied.
opt
1) Three–point property (S opt
Dk , S̃ Dk , ω̄ k ). For all i ≥
∈ ω̄k (i) and S̃ Dk ∈
1, S opt
∈ S Dk (i), ω̄ opt
Dk
k
opt
opt
, we have
arg minω̄opt ∈ω̄ (i) Dk S Dk , ω̄ k
k

k

(40)
opt
˜
,
S̃
,
ω̄
).
For
all
,
ω̄
2) Four–point property (S opt
D
k
k
k
Dk
opt
i ≥ 1, S opt
,
S̃
∈
S
(i),
ω̄
∈
ω̄
(i)
and
Dk
Dk
k
Dk
k 

opt
˜ k ∈ arg minω̄opt ∈ω̄ (i) Dk S̃ Dk , ω̄ k
ω̄
, we have
k

k







opt
opt
opt
˜
S
,
ω̄
+
f
S
,
S̃
≤
D
,
ω̄
Dk S opt
Dk .
k
k
Dk
Dk
Dk
k

(41)

opt
Theorem: Let {S Dk (i)}i≥0 , {ω̄k (i)}i≥0 , S opt
be
Dk , ω̄ k
compact subsets of the compact metric space (M, r) such
that
r

h
S opt
S Dk (i) →
Dk

The Hausdorff distance measures how far two subsets of a
metric space are from each other and is defined by


rH (X, Y ) = max sup inf r(x, y), sup inf r(x, y) .
x∈X y∈Y

(38)
opt
and
ω̄
correspond
to
the
optimal
values
of
where S opt
Dk
k
S Dk (i) and ω̄ k (i), respectively, and the sequences of compact sets {S Dk (i)}i≥0 and {ω̄ k (i)}i≥0 converge to the sets
S Dk ,opt and ω̄ k,opt , respectively.
The sets S Dk ,opt and ω̄ k,opt are not directly given, but
we observe the sequence of compact sets {S Dk (i)}i≥0 and
{ω̄k (i)}i≥0 . The goal of the proposed algorithms is to find a
sequence of S Dk (i) and ω̄ k (i) in a distributed way such that




opt
,
ω̄
.
(39)
lim Dk S Dk (i), ω̄ k (i) = Dk S opt
Dk
k







opt
opt
opt
S
,
ω̄
.
S̃
,
ω̄
≤
D
+
D
f S opt
,
S̃
D
k
k
D
k
k
Dk
Dk
k
k

B. Sufficient Conditions for Convergence

1)
2)
3)
4)

k

k

r

h
ω̄ opt
ω̄ k (i) →
k

and let Dk : M × M → R be a continuous function.
Now, let conditions 1) and 2) hold. Then, for the proposed
algorithms, we have




opt
,
ω̄
.
(43)
lim Dk S Dk (i), ω̄ k (i) = Dk S opt
Dk
k
i→∞

y∈Y x∈X

(36)

(42)

A general proof of this theorem is detailed in [53], [54].

7

TABLE III
C OMPUTATIONAL C OMPLEXITY OF D IFFERENT A LGORITHMS
Algorithm
DRJIO–NLMS
DRJIO–RLS
Distributed NLMS [2]
Distributed RLS [70]
Krylov Subspace
NLMS [13]
Distributed principal
subspace estimation [11]

Multiplications
2(D + 1)M + (3 + |Nk |)D + 5
2M 2 + (3 + 2D)M + 4D 2
+(9 + |Nk |)D
(4 + |Nk |)M + 1
4M 2 + (12 + |Nk |)M − 1
6DM 2 + 4M + (5 + |Nk |)D

Additions
(2D + 1)M + (2 + |Nk |)D − 2
2M 2 + 2DM + 4D 2
+(2 + |Nk |)D
(5 + |Nk |)M − 1
4M 2 + (16 + |Nk |)M + 1
6DM 2 + 2M + (2 + |Nk |)D

M 3 + 2(D + 2)M
+(3 + |Nk |)D + 4

M 3 + (D + 1)M
+(2 + |Nk |)D − 1

C. Convergence to the Optimal Low-Rank Estimator
In this section, we show that the proposed low–rank algorithm globally and exponentially converges to the optimal
low-rank estimator [15], [72]. This result is applicable to leastsquares type algorithms with forgetting factor λ = 1. We
remark that for stochastic gradient (or LMS) algorithms and
least-squares algorithms with forgetting factor λ 6= 1, there
will be a misadjustment or loss in MSE due to the adaptation
with the step size when an LMS algorithms is adopted, or due
to the forgetting factor when an RLS algorithm is chosen. To
proceed with our proof, let us rewrite the expressions in (8)
and (10 for time instant zero as follows:
Rk (0)S Dk (0)R̄ω̄ k (0) = P Dk (0) = pk (0)ω̄H
k (0) + δΥ,
(44)
h
i
ID
where Υ =
, Υ is an M × D matrix containing
0M−D×D
an identity matrix I D with size D and an M − D × D matrix
with zeros 0M−D×D , δ is a small positive scalar used to
regularize the recursion at initialization and ensure that a rankD matrix P Dk (0) is obtained, and the D-dimensional set of
normal equations that must be solved to compute ω̄ k (1) is
given by
R̄k (0)ω̄ k (1) = S H
Dk (0)Rk (0)S Dk (0)ω̄ k (1)
= S Dk (0)pk (0) = p̄k (0),

(45)

where pk (i) = E[dk (i)∗ xk (i) is the cross-correlation vector.
Using (44), we can obtain the following relation

−1
2
(0)R
(0)S
(0)
R̄ω̄k (0) = S H
Dk
Dk
k
× S Dk (0)Rk (0)P Dk (0).

(46)

Substituting the aforementioned result for R̄ω̄ k (0) into the
expression in (44), we get a recursive expression for S Dk (0)
as

−1
H
S Dk (0) = R−1
(0)P
(0)
S
(0)R
(0)P
(0)
Dk
k
Dk
Dk
k
−1

2
.
(47)
× SH
Dk (0)Rk (0)S Dk (0)
Using (44), we can express ω̄k (1) as

−1
(0)R
(0)S
ω̄ k (1) = S H
(0)
SH
k
D
Dk
Dk (0)pk (0).
k

(48)

For the proposed DRJIO–NLMS and DRJIO–RLS, the relation
is given by
X
ω k (1) = S Dk (1)
ckl ω̄ l (1).
(49)
l∈Nk

Substituting S Dk (1) and ω̄ l (1) into (49), we obtain

−1
−1
H
ωk (1) = Rk (1)P Dk (1) S Dk (1)Rk (1)P Dk (1)

−1
H
2
× S Dk (1)Rk (1)S Dk (1)

−1
X
ckl S H
×
(0)R
(0)S
(0)
SH
l
Dl
Dl
Dl (0)pl (0).
l∈Nk

(50)
More generally, we can express the proposed distributed
algorithms by the following recursion:
X
ω k (i) = S Dk (i)
ckl ω̄ l (i)
l∈Nk


−1
H
=
S Dk (i)Rk (i)P Dk (i)

−1
H
2
× S Dk (i)Rk (i)S Dk (i)

−1
X
(i
−
1)R
(i
−
1)S
(i
−
1)
×
ckl S H
l
Dl
Dl
R−1
k (i)P Dk (i)

×

l∈Nk
SH
Dl (i

− 1)pl (i − 1).

(51)

At this point, we resort to the assumption that the matrices
S Dk (i) for each node k must converge to the same values,
which correspond to the optimal transformation matrix. Because the optimal low-rank filter can be described by the
−1/2
eigenvalue decomposition of Rk (i)pk (i) [20], [21], where
−1/2
Rk (i) is the square root of the matrix Rk (i), and pk (i) is
the cross-correlation vector, we have
−1/2

Rk

(i)pk (i) = Φk Λk ΦH
k pk (i),

(52)

where Λk is an M ×M diagonal matrix with the eigenvalues of
Rk , and Φk is a M × M unitary matrix with the eigenvectors
of Rk .
Let us also assume that there exists some ω k (0) such that
the randomly selected S Dk (0) can be written as [21]
−1/2

S Dk (0) = Rk

(i)Φk ω k (0),

(53)

8

Using (52) and (53) in (51) together with the assumption
and some manipulation of the algebraic expressions, we can
express (51) in a more compact way that is suitable for
analysis, as given by
X
−1
2
ω k (i) =
ckl Λ2l ω l (i − 1) ωH
l (i − 1)Λl ω l (i − 1)
l∈Nk
× ωH
l (i

− 1)ω l (i − 1).

(54)

20

14

12
11

18

16

10
4

14

18

20
12

10

9

15

7

8

The aforementioned expression can be decomposed as follows:
X
ω k (i) =
ckl Ql (i)Ql (i − 1) . . . Ql (1)ωl (0),
(55)
l∈Nk

19

8

6

13

4

17
2

16

where

0

Ql (i) =

Λ2i
l ω l (0)

−1 H
4i−2
ωH
ωl (0)
ω l (0)Λl2i−2 .
l (0)Λl

(56)
At this point, we need to establish that the norm of S Dk (i), for
all i, is both lower and upper bounded, i.e., 0P
<k S Dk (i) k<
∞, for all i, and that ω k (i) = S Dk (i) l∈Nk ckl ω̄ l (i)
exponentially approaches ωk,opt (i) as i increases. Due to the
linear mapping, the boundedness of S Dk (i) is equivalent to the
boundedness of ωk (i). Therefore, we have upon convergence
when i → ∞ that
H
ωH
k (i)ω k (i − 1) = ω k (i − 1)ω k (i − 1).

k ω k (∞) k≥k ω k (i) k≥k ω k (0) k .

Using the matrix identity (A + B)−1 = A−1 − A−1 B(I +
A−1 B)−1 A−1 in the decomposed Ql (i) in (59) and making
i large, we get

Ql (i) = diag(1
. . 1} 0| .{z
(60)
. . 0}) + Ol ǫl (i) ,
| .{z
M−D

where ǫl (i) = (σr+1 /σr )2i , in which σr+1 and σr are the
−1/2
(r + 1)th and the rth largest singular values of Rl
(i)pl (i),
respectively, and O(·) denotes the order of the argument.
Based on (60), it follows that, for some positive constant g,
we have k ω l (i) k≤ 1 + gǫl (i). Based on (55), we obtain
X
k ωk (∞) k ≤
ckl (k Ql (∞) k . . . k Ql (1) kk Ql (0) k)
l∈Nk

2

4

2

1

6

8

10

12

5
14

16

18

20

Fig. 4. Network topology with N = 20 nodes

≤

X

ckl



ckl



ckl



l∈Nk

=

X

l∈Nk

≤

X

l∈Nk

=k

X

l∈Nk

ckl

k ω l (0) k

∞
Y

i=0



1 + gǫl (i)

k ω l (0) k exp

X
∞


log 1 + gǫl (i)

k ω l (0) k exp

X
∞


gǫl (i)

i=1



i=1

k ωl (0) k exp






.
1 − (σr+1 /σr )2
g

(61)

(58)

To show that the upper bound k ω k (∞) k is finite, let
us express the M × M matrix Qk (i) as a function of the
ω l,1 (i)
M × 1 vector ω l (i) =
and the M × M matrix
ω l,2 (i)


Λl,1
Λ =
. Substituting the previous expressions
Λl,2
of ω l (i) and Λk l into Ql (i) as given in (56), we obtain

 2i
Λl,1 ωl,1 (0)
4i−2
Ql (i) =
ωH
l,1 (0)Λl,1 ω l,1 (0)
Λ2i
l,2 ω l,2 (0)

2i−2 
−1 ω H
4i−2
l,1 (0)Λl,1
+ ωH
(0)Λ
ω
(0)
(59)
l,2
2i−2
l,2
l,2
ωH
l,2 (0)Λl,2

D

0

(57)

H
Because k ωH
k (i)ω k (i − 1) k≤k ω k (i − 1) kk ω k (i) k
H
and k ωk (i − 1)ωk (i − 1) k=k ω k (i − 1) k2 , the relaH
tion ω H
k (i)ω k (i − 1) = ω k (i − 1)ω k (i − 1) implies that
k ωk (i) k≥k ω k (i − 1) k, and hence we have

3

6

With the previous development, the norm of ω k (i) is proven
to be both lower and upper bounded. Once this case has been
established, the expression in (51) converges for a sufficiently
large i to the low-rank Wiener filter. This condition is verified
by equating the terms of (54), which yields


X

−1/2
H
(i)Φl,1 Λl,1 Φl,1 pl (i) + Ol ǫl (i)
ckl Rl
ω k (i) =
l∈Nk

(62)
where Φl,1 is a M × D matrix with the D largest eigenvectors
of Rl (i), and Λl,1 is a D × D matrix with the largest
eigenvalues of Rl (i).
VI. S IMULATION

RESULTS

In this section, we investigate the performance of the
proposed DRJIO–NLMS and DRJIO–RLS algorithms for distributed estimation in two scenarios: wireless sensor networks
and smart grids.
A. Wireless Sensor Networks
In this subsection, we compare the proposed DRJIO–NLMS
and DRJIO–RLS algorithms with the distributed NLMS algorithm (normalized version of [2]), distributed RLS algorithm
[70], Krylov subspace NLMS [13] and distributed principal
subspace estimation [11], based on their MSE performance.

9

With the network topology structure outlined in Fig. 4 with
N = 20 nodes, we consider numerical simulations under three
scenarios for the parameter vector ωo :

•

The input signal is generated as xk (i) = [xk (i) xk (i −
1) ... xk (i − M + 1)] and xk (i) = uk (i) + αk xk (i − 1),
where αk is a correlation coefficient and uk (i) is a white noise
2
process with variance σu,k
= 1 − |αk |2 , to ensure the variance
2
of xk (i) is σx,k = 1. In particular, this application requires
the estimation of a set of parameters that could be modeled
as a finite-impulse response (FIR) filter (related to a moving
average (MA) model). Furthermore, the algorithms would also
work with an input generated by an auto-regressive (AR)
model but their performance would depend on the condition
number and the rank of the correlation matrix of the input
data. The noise samples are modeled as complex Gaussian
2
= 0.001. We have adopted the
noise with variance σn,k
regularization parameters γ = 0.02 and δ = 0.01 in all
examples. We have also evaluated the impact of different
values of regularization parameters and the results indicate
that the performance of the algorithms degrades when the
parameters are not well chosen. Moreover, the optimized
values work very well for a wide range of scenarios and values
of noise variance. We assume that the network has error–free
transmission between linked nodes.
The step size µ0 for the distributed NLMS algorithm,
Krylov subspace NLMS, distributed principal subspace estimation and DRJIO–NLMS is set to 0.15 and η0 is set to 0.5.
For the distributed RLS algorithm and DRJIO–RLS algorithm,
the forgetting factor λ is equal to 0.99 and δ is set to 0.11. In
Fig. 5, we compare the proposed DRJIO–NLMS and DRJIO–
RLS algorithms with the existing strategies using the full–rank
system with M =20 and D=5. The dimensionality reduction
matrix S Dk (0) is initialized as [I D 0D,M−D ]T .
We observe that the proposed DRJIO–RLS algorithm has
the best performance when compared with other algorithms,
while the proposed DRJIO–NLMS algorithm also has a better
performance, which is very close to the distributed RLS
algorithm. The superior performance of DRJIO-RLS can be
explained by the fact that the convergence rate or learning
speed of adaptive algorithms depends on the number of
parameters that need to be estimated. This is well known in
adaptive signal processing [6]. For instance, when we compare
the DRJIO-RLS and the full-rank RLS, the difference is that
the proposed DRJIO-RLS estimates the unknown parameters
using a reduced dimension and retaining the most relevant features of the data. As a result, the DRJIO-RLS converges faster
than the standard RLS algorithm. However, its complexity is
an order of magnitude lower than those of the distributed RLS
algorithm and the DRJIO–RLS algorithm.
When the full–rank system M increases to 60, Fig. 6
illustrates that, the proposed DRJIO–RLS algorithm still has
the best performance, while DRJIO–NLMS algorithm also
shows a fast convergence rate, which is comparable to the

0

−5
MSE(dB)

•

DRJIO−RLS(D=5)
Distributed Principal
Subspace Estimation(D=5)[11]
Krylov Subspace NLMS(D=5)[13]
Distributed RLS(M=20)[26]
Distributed NLMS(M=20)[2]
DRJIO−NLMS(D=5)

5

Full–rank system with M =20
Sparse system with M =20 (D non-zero coefficients and
M − D zeros coefficients)
Full–rank system with M =60

−10

−15

−20

−25

−30
0

100

200

300

400
500
Time instant, i

600

700

800

900

Fig. 5. Full–rank system with M =20

15

DRJIO−RLS(D=5)
Distributed Principal
Subspace Estimation(D=5)[11]
Krylov Subspace NLMS(D=5)[13]
Distributed RLS(M=20)[26]
Distributed NLMS(M=20)[2]
DRJIO−NLMS(D=5)

10

5

0

MSE(dB)

•

10

−5

−10

−15

−20

−25

−30
0

100

200

300

400
500
Time instant, i

600

700

800

900

Fig. 6. Full–rank system with M =60

distributed RLS algorithm. For the distributed NLMS, Krylov
subspace NLMS and distributed principal subspace estimation
algorithms, their convergence speed is much lower.
In a sparse system scenario with M = 100, we first
evaluate the MSE performance versus the rank D and then we
assess the MSE performance versus the number of iterations,
as shown in Figs. 7 and 8, respectively. In particular, the
curves illustrating the MSE performance versus the rank D
are obtained after 500 iterations for a range of D between
1 and 10. The results depicted in Fig. 7 indicate that the
best rank for both DJRIO–NLMS and DJRIO–RLS algorithms
corresponds to D = 5 and that the MSE performance gradually
degrades for other values. The rank D should be carefully
selected as it affects the performance of these algorithms
and determines the number of parameters that should be

10

−15

MSE (dB)

−20

−25

DJIO−NLMS
DJIO−RLS
−30

B. Smart Grids
1

2

3

4

5
6
Rank (D)

7

8

9

10

Fig. 7. MSE performance versus rank D for a sparse system with M =100

10
DRJIO−RLS(D=5)
Distributed Principal
Subspace Estimation(D=5)[11]
Krylov Subspace NLMS(D=5)[13]
Distributed RLS(M=100)[26]
Distributed NLMS(M=100)[2]
DRJIO−NLMS(D=5)

5

0

−5

MSE(dB)

the proposed DCE scheme outperforms the DRJIO–NLMS
algorithm. With the decrease of the system sparsity level,
the proposed DRJIO–NLMS algorithm outperforms the DCE
scheme. The results of 9 and 10 indicate that the proposed
DRJIO–NLMS algorithm is superior to the DCE scheme when
the level of sparsity is not very high. The computational
complexity of DRJIO-NLMS is much lower than the DCE
scheme because the latter requires a basis pursuit algorithm
to reconstruct the full-dimension estimator and DRJIO-NLMS
employs a rank-D approximation based on a simple matrixvector multiplication.

−10

−15

−20

−25

−30
0

100

200

300

400

500

600

700

800

900

Time instant, i

Fig. 8. Sparse system with M =100

exchanged between nodes. Moreover, we have considered
D = 5 for assessing the MSE performance versus the number
of iterations for the proposed and other existing algorithms, as
shown in Fig. 8. The results indicate that the proposed DRJIO–
RLS and DRJIO–NLMS algorithms have a more pronounced
performance advantage over the distributed NLMS, the Krylov
subspace NLMS and the distributed principal subspace estimation algorithms. Specifically, the proposed DRJIO–NLMS
algorithm performs very close to the distributed RLS algorithm
and outperforms the other analyzed algorithms.
In the last example on wireless sensor networks, we compare the performance between the proposed DRJIO–NLMS
and the DCE scheme in [45], under different sparsity level
scenarios. The step size for both algorithms is set to 0.3 and the
η0 for DRJIO–NLMS is set to 0.5. The length of the unknown
parameter ω 0 is 20 and D = 10. For the first scenario, the
number of non–zero coefficients in the unknown parameter
is 3 and for the second scenario, the number of non–zero
coefficients is set to 10. The comparison results are shown
in Fig. 9 and 10. It is clear that in a very sparse system,

In order to test the proposed algorithms in a possible
smart grid scenario, we consider the Hierarchical IEEE 14–
bus system which has been proposed in [73], where 14 is
the number of substations. At every time instant i, each
bus k, k = 1, 2, . . . , 14, takes a scalar measurement dk (i)
according to

(63)
dk (i) = Xk ω0 (i) + nk (i), k = 1, 2, . . . , 14,

where ω0 (i) is the state vector of the entire interconnected
system, Xk (ω 0 (i)) is a nonlinear measurement function of
bus k. The quantity nk (i) is the measurement error with mean
equal to zero and which corresponds to bus k.
We focus on the linearized DC state estimation problem.
We assume that each bus connects and measures the state of
three users. As a result, for the IEEE–14 bus system, there will
be 42 users in the system. The system is built with 1.0 per
unit (p.u) voltage magnitudes at all users and j1.0 p.u. branch
impedance. Then, the state vector ω 0 (i) is taken as the voltage
phase angle vector ω 0 for all users. Initially, each bus only
knows the voltage phase angle of the three users connected to
it. With the help of distributed estimation algorithms, each bus
is supposed to estimate the state of the voltage phase angles for
all users in the system. Therefore, the nonlinear measurement
model for state estimation (63) is approximated by
dk (i) = xH
k (i)ω 0 + nk (i),

k = 1, 2, . . . , 14,

(64)

where xk (i) is the measurement Jacobian vector for bus k.
Then, the aim of the distributed estimation algorithm is to
compute an estimate of ω 0 , which can minimize the cost
function given by
2
Jω k (i) (ω k (i)) = E|dk (i) − xH
k (i)ω k (i)| .

(65)

and the global network cost function is described by
Jω (ω) =

N
X

2
E|dk (i) − xH
k (i)ω| .

(66)

k=1

We compare the proposed algorithms with the M–CSE algorithm [4], the distributed RLS algorithm [70], the distributed
NLMS algorithm (normalized version of [2]) and distributed
principal subspace estimation [11] in terms of MSE performance. The MSE comparison is used to determine the
accuracy of the algorithms and the rate of convergence. We
define the Hierarchical IEEE–14 bus system as in Fig. 11.

11

Sparsity Level S=3

10

10

5

5

0

0

−5
MSE(dB)

−5
MSE(dB)

DRJIO−RLS(M=10)
M−CSE[4]
Distributed RLS[26]
Distributed NLMS[2]
DRJIO−NLMS(M=10)
Distributed Principal
Subspace Estimation(D=10)[11]

DRJIO−NLMS
DCE Scheme

−10

−10

−15

−15

−20

−20

−25

−25

−30
0

50

100

150

200
250
Time instant, i

300

350

400

−30
0

100

200

300

400
500
Time instant, i

600

700

800

900

Fig. 12. MSE performance for smart grids

Fig. 9. DRJIO–NLMS vs DCE scheme with sparsity level S=3

Sparsity Level S=10
10

DRJIO−NLMS
DCE Scheme

5

0

MSE(dB)

−5

−10

−15

−20

−25

−30
0

50

100

150

200
250
Time instant, i

300

350

Fig. 10. DRJIO–NLMS vs DCE scheme with sparsity level S=10

14

13

400

All buses are corrupted by additive white Gaussian noise
2
with variance σn,k
= 0.001. The step size for the distributed
NLMS [2] and the proposed DRJIO–NLMS algorithms is µ =
0.15 and η0 is set to 0.5. The parameter vector ω0 is set to
an all–one vector with size 42 × 1. For the distributed RLS,
DRJIO–RLS algorithms the forgetting factor λ is set to 0.99
and δ is equal to 0.11. The reduced dimension D is set to
10 for both DRJIO–RLS and DRJIO–NLMS algorithm. The
results are averaged over 100 independent runs. We simulate
the proposed algorithms for smart grids under a static scenario.
From Fig. 12, it can be seen that the proposed DRJIO–
RLS algorithm has the best performance, and significantly
outperforms the distributed NLMS [2] and the M–CSE [4] algorithms. The DRJIO–NLMS is slightly worse than distributed
RLS algorithm [70], but better than the distributed NLMS and
M–CSE algorithms. In addition, the proposed DRJIO–NLMS
and DRJIO–RLS algorithms can compress the data to be
transmitted from each node from M to D, resulting in reduced
bandwidth requirements. These algorithms are also important
tools for dealing with large sets of data which exhibit some
form of redundancy, sparsity and are compressible.

12

VII. C ONCLUSIONS

10
11
6
9
8
7
1
5
4
3
2

Fig. 11. Hierarchical IEEE 14–bus system

In this paper, we have proposed a novel distributed lowrank scheme along with efficient algorithms for distributed
estimation in wireless sensor networks and smart grids. Simulation results have shown that the proposed DRJIO–RLS has
the best performance, while DRJIO–NLMS algorithm has a
better performance and lower cost than existing algorithms
in all the three scenarios considered. We have also compared
the proposed algorithms with the DCE scheme, which was
presented in [45], for systems with different levels of sparsity.
Furthermore, the proposed scheme requires the transmission of
only D parameters instead of M , resulting in higher bandwidth
efficiency than standard schemes.

12

A PPENDIX
R ECONSTRUCTION

USING A RANK -D APPROXIMATION

In this appendix, we show how the reconstruction of the fulldimension estimator ω k (i) can be carried out using a rank-D
approximation with the low-rank estimator ω¯k (i), i. e.,
ωk (D) (i) = S Dk (i)ω̄ k (i),

(67)

In order to show the above relation, we consider the expression
of the low-rank estimator given by
−1

ω¯k (i) = R̄k (i)p̄k (i)

−1
SH
= SH
(i)R
(i)S
(i)
k
D
Dk (i)pk (i).
Dk
k

(68)

and the Wiener filter given by
ω k (i) = R−1
k (i)pk (i).

(69)

Using the fact that the low-rank estimator converges to the
low-rank Wiener filter, S Dk (i) converges to a M × D matrix
with the eigenvectors ΦDP
and an eigenvalue decomposition of
N
H
Rk (i) = ΦN ΛN ΦH
=
N
n=1 λn (i)φn (i)φn (i), where ΛN
and ΦN are the N × N diagonal matrix with the eigenvalues
and the N × N unitary matrix with the eigenvectors of Rk (i),
respectively, λn (i) is the nth eigenvalue and φn (i) is the nth
eigenvector of Rk (i), we have
N

−1
X
ω¯k (i) = S H
(i)
λn (i)φn (i)φH
SH
Dk
n (i)S Dk (i)
D (i)pk (i)
n=1

H
−1 H
= (ΦH
ΦD pk (i)
D ΦN Λ N ΦN ΦD )
H
= Λ−1
D ΦD pk (i),

(70)
then, multiplying S Dk (i) = ΦD on both sides, we obtain
H
S Dk (i)ω¯k (i) = S Dk (i)Λ−1
D ΦD pk (i)
−1 H
ΦD pk (i)
= ΦD Λ D

= ωk

(D)

(71)

(i),

(D)

where Rk (i) = ΦD ΛD ΦH
D is a rank-D approximation
H
of Rk (i) and ω k (D) (i) = ΦD Λ−1
D ΦD pk (i) is the rank-D
approximation of ωk (i), which gives us the relation in (67).
Note that when D = M , the D-rank approximation yields the
full-rank Wiener filter.
R EFERENCES
[1] C. G. Lopes and A. H. Sayed, “Incremental adaptive strategies over
distributed networks,” IEEE Trans. Signal Process., vol. 48, no. 8, pp.
223–229, August 2007.
[2] ——, “Diffusion least–mean squares over adaptive networks: Formulation and performance analysis,” IEEE Trans. Signal Process., vol. 56,
no. 7, pp. 3122–3136, July 2008.
[3] S. Xu, R. C. de Lamare, and H. V. Poor, “Adaptive link selection
algorithms for distributed estimation,” EURASIP Journal on Advances
in Signal Processing, vol. 2015, no. 1, p. 86, 2015.
[4] L. Xie, D. H. Choi, S. Kar, and H. V. Poor, “Fully distributed state
estimation for wide-area monitoring systems,” IEEE Trans. Smart Grid,
vol. 3, no. 3, pp. 1154–1169, September 2012.
[5] S. Xu, R. C. de Lamare, and H. V. Poor, “Dynamic topology adaptation
for distributed estimation in smart grids,” in Proc. IEEE 5th International Workshop on Computational Advances in Multi-Sensor Adaptive
Processing, December 2013, pp. 420–423.

[6] S. Haykin, Adaptive Filter Theory, 4th ed. Upper Saddle River, NJ,
USA: Prentice Hall, 2002.
[7] A. H. Sayed, Fundamentals of Adaptive Filtering. Hoboken, NJ, USA:
John Wiley&Sons, 2003.
[8] E. Msechu, S. Roumeliotis, A. Ribeiro, and G. Giannakis, “Decentralized quantized kalman filtering with scalable communication cost,” IEEE
Trans. Signal Process., vol. 56, no. 8, pp. 3727–3741, October 2008.
[9] J. J. Xiao, A. Ribeiro, Z. Q. Luo, and G. Giannakis, “Decentralized
quantized Kalman filtering with scalable communication cost,” IEEE
Signal Process. Mag., vol. 23, no. 4, pp. 27–41, July 2006.
[10] S. Pereira and A. Pages-Zamora, “Distributed consensus in wireless sensor networks with quantized information exchange,” in Proc. IEEE 9th
Workshop on Signal Processing Advances in Wireless Communications,
July 2008, pp. 241–245.
[11] L. Li, A. Scaglione, and J. H. Manton, “Distributed principal subspace
estimation in wireless sensor networks,” IEEE J. Sel. Topics Signal
Process., vol. 5, no. 4, pp. 725–738, August 2011.
[12] M. O. Sayin and S. Kozat, “Single bit and reduced dimension diffusion
strategies over distributed networks,” IEEE Signal Process. Lett., vol. 20,
no. 10, pp. 976–979, October 2013.
[13] S. Chouvardas, K. Slavakis, and S. Theodoridis, “Trading off complexity
with communication costs in distributed adaptive learning via Krylov
subspaces for dimensionality reduction,” IEEE J. Sel. Topics Signal
Process., vol. 7, no. 2, pp. 257–273, April 2013.
[14] L. L. Scharf and D. W. Tufts, “Rank reduction for modeling stationary
signals,” IEEE Transactions on Acoustics Speech and Signal Processing,
vol. 35, pp. 350–355, 1987.
[15] L. L. Scharf, “The SVD and reduced-rank signal processing,” Signal
Processing, vol. 25, no. 2, pp. 113–133, November 1991.
[16] M. L. Honig and J. S. Goldstein, “Adaptive reduced–rank interference
suppression based on the multistage Wiener filter,” IEEE Trans. Commun., vol. 6, no. 50, pp. 986–994, June 2002.
[17] R. C. de Lamare and R. Sampaio-Neto, “Adaptive reduced-rank mmse
filtering with interpolated fir filters and adaptive interpolators,” IEEE
Signal Processing Letters, vol. 12, no. 3, pp. 177–180, March 2005.
[18] ——, “Reduced–rank adaptive filtering based on joint iterative optimization of adaptive filters,” IEEE Signal Process. Lett., vol. 14, no. 12, pp.
980–983, December 2007.
[19] M. Yukawa, R. C. de Lamare, and R. Sampaio-Neto, “Efficient acoustic
echo cancellation with reduced-rank adaptive filtering based on selective
decimation and adaptive interpolation,” IEEE Transactions on Audio,
Speech, and Language Processing, vol. 16, no. 4, pp. 696–710, May
2008.
[20] R. C. de Lamare and R. Sampaio-Neto, “Adaptive reduced-rank processing based on joint and iterative interpolation, decimation, and filtering,”
IEEE Transactions on Signal Processing, vol. 57, no. 7, pp. 2503–2514,
2009.
[21] ——, “Reduced-rank space-time adaptive interference suppression with
joint iterative least squares algorithms for spread-spectrum systems,”
IEEE Transactions on Vehicular Technology, vol. 59, no. 3, pp. 1217–
1228, March 2010.
[22] R. Fa, R. C. de Lamare, and L. Wang, “Reduced-rank stap schemes
for airborne radar based on switched joint interpolation, decimation and
filtering algorithm,” IEEE Transactions on Signal Processing, vol. 58,
no. 8, pp. 4182–4194, Aug 2010.
[23] L. Wang, R. C. de Lamare, and M. Haardt, “Direction finding algorithms
based on joint iterative subspace optimization,” IEEE Transactions on
Aerospace and Electronic Systems, vol. 50, no. 4, pp. 2541–2553,
October 2014.
[24] N. Song, W. U. Alokozai, R. C. de Lamare, and M. Haardt, “Adaptive
widely linear reduced-rank beamforming based on joint iterative optimization,” IEEE Signal Processing Letters, vol. 21, no. 3, pp. 265–269,
March 2014.
[25] N. Song, R. C. de Lamare, M. Haardt, and M. Wolf, “Adaptive widely
linear reduced-rank interference suppression based on the multistage
wiener filter,” IEEE Transactions on Signal Processing, vol. 60, no. 8,
pp. 4003–4016, Aug 2012.
[26] P. Clarke and R. C. de Lamare, “Low–complexity reduced-rank linear
interference suppression based on set-membership joint iterative optimization for ds–cdma systems,” IEEE Trans. Veh. Technol., vol. 60,
no. 9, pp. 4324–4337, November 2011.
[27] R. C. de Lamare, M. Haardt, and R. Sampaio-Neto, “Blind adaptive
constrained reduced-rank parameter estimation based on constant modulus design for cdma interference suppression,” IEEE Transactions on
Signal Processing, vol. 56, no. 6, pp. 2470–2482, June 2008.

13

[28] H. Ruan and R. C. de Lamare, “Robust adaptive beamforming using a
low-complexity shrinkage-based mismatch estimation algorithm,” IEEE
Signal Processing Letters, vol. 21, no. 1, pp. 60–64, Jan 2014.
[29] L. Qiu, Y. Cai, R. C. de Lamare, and M. Zhao, “Reduced-rank doa
estimation algorithms based on alternating low-rank decomposition,”
IEEE Signal Processing Letters, vol. 23, no. 5, pp. 565–569, May 2016.
[30] H. Ruan and R. C. de Lamare, “Robust adaptive beamforming based
on low-rank and cross-correlation techniques,” IEEE Transactions on
Signal Processing, vol. 64, no. 15, pp. 3919–3932, Aug 2016.
[31] X. Wu, Y. Cai, M. Zhao, R. C. de Lamare, and B. Champagne, “Adaptive
widely linear constrained constant modulus reduced-rank beamforming,”
IEEE Transactions on Aerospace and Electronic Systems, vol. 53, no. 1,
pp. 477–492, Feb 2017.
[32] R. C. de Lamare, “Massive mimo systems: Signal processing challenges
and future trends,” URSI Radio Science Bulletin, vol. 86, no. 4, pp. 8–20,
Dec 2013.
[33] Y. Sun, V. Tripathi, and M. L. Honig, “Adaptive, iterative, reduced–rank
(turbo) equalization,” IEEE Trans. Wireless Commun., vol. 4, no. 6, pp.
2789–2800, March 2005.
[34] W. Zhang, H. Ren, C. Pan, M. Chen, R. C. de Lamare, B. Du, and J. Dai,
“Large-scale antenna systems with ul/dl hardware mismatch: Achievable
rates analysis and calibration,” IEEE Transactions on Communications,
vol. 63, no. 4, pp. 1216–1229, April 2015.
[35] R. C. de Lamare and R. Sampaio-Neto, “Adaptive reduced–rank equalization algorithms based on alternating optimization design techniques
for MIMO systems,” IEEE Trans. Veh. Technol., vol. 60, no. 6, pp.
2482–2494, July 2011.
[36] K. Zu, R. C. de Lamare, and M. Haardt, “Generalized design of lowcomplexity block diagonalization type precoding algorithms for multiuser mimo systems,” IEEE Transactions on Communications, vol. 61,
no. 10, pp. 4232–4242, October 2013.
[37] ——, “Multi-branch tomlinson-harashima precoding design for mumimo systems: Theory and algorithms,” IEEE Transactions on Communications, vol. 62, no. 3, pp. 939–951, March 2014.
[38] L. Zhang, Y. Cai, R. C. de Lamare, and M. Zhao, “Robust multibranch
tomlinson-harashima precoding design in amplify-and-forward mimo
relay systems,” IEEE Transactions on Communications, vol. 62, no. 10,
pp. 3476–3490, Oct 2014.
[39] Y. Cai, R. C. de Lamare, B. Champagne, B. Qin, and M. Zhao, “Adaptive
reduced-rank receive processing based on minimum symbol-error-rate
criterion for large-scale multiple-antenna systems,” IEEE Transactions
on Communications, vol. 63, no. 11, pp. 4185–4201, Nov 2015.
[40] W. Zhang, R. C. de Lamare, C. Pan, M. Chen, J. Dai, B. Wu, and X. Bao,
“Widely linear precoding for large-scale mimo with iqi: Algorithms and
performance analysis,” IEEE Transactions on Wireless Communications,
vol. 16, no. 5, pp. 3298–3312, May 2017.
[41] H. Qian and S. N. Batalama, “Data–record–based criteria for the
selection of an auxiliary vector estimator of the MMSE/MVDR filter,”
IEEE Trans. Commun., vol. 51, no. 10, pp. 1700–1708, October 2003.
[42] S. Nuan, W. U. Alokozai, R. C. de Lamare, and M. Haardt, “Adaptive widely linear reduced–rank beamforming based on joint iterative
optimization,” IEEE Signal Process. Lett., vol. 21, no. 3, pp. 265–269,
March 2014.
[43] L. Wang and R. C. de Lamare, “Low–complexity constrained adaptive
reduced-rank beamforming algorithms,” IEEE Trans. Aerosp. Electron.
Syst., vol. 49, no. 4, pp. 2114–2128, October 2013.
[44] Y. Yao, A. P. Petropulu, and H. V. Poor, “MIMO radar using compressive
sampling,” IEEE J. Sel. Top. Sign. Proces., vol. 4, no. 1, pp. 146–163,
February 2010.
[45] S. Xu, R. C. de Lamare, and H. V. Poor, “Distributed compressed
estimation based on compressive sensing,” IEEE Signal Process. Lett.,
vol. 22, no. 9, pp. 1311–1315, September 2015.
[46] J. Chen, Z. J. Towfic, and A. H. Sayed, “Dictionary learning over
distributed models,” IEEE Transactions on Signal Processing, vol. 63,
no. 4, pp. 1001–1016, February 2015.
[47] P. Chainais and C. Richard, “Learning a common dictionary over a
sensor network,” in Proceedings of IEEE International Workshop on
Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), December 2013, pp. 1–4.
[48] J. Liang, M. Zhang, X. Zeng, and G. Yu, “Distributed dictionary learning
for sparse representation in sensor networks,” IEEE Transactions on
Image Processing, vol. 23, no. 4, pp. 2528–2541, June 2014.
[49] S. Chouvardas, Y. Kopsinis, and S. Theodoridis, “An online algorithm
for distributed dictionary learning,” in Proceedings of the IEEE Internactional Conference on Acoustics, Speech and Signal Processing, April
2015, pp. 1–5.

[50] Z. Yang, R. C. de Lamare, and X. Li, “l1 -regularized stap algorithms
with a generalized sidelobe canceler architecture for airborne radar,”
IEEE Transactions on Signal Processing, vol. 60, no. 2, pp. 674–686,
Feb 2012.
[51] T. G. Miller, S. Xu, R. C. de Lamare, and H. V. Poor, “Distributed
spectrum estimation based on alternating mixed discrete-continuous
adaptation,” IEEE Signal Processing Letters, vol. 23, no. 4, pp. 551–
555, April 2016.
[52] H. Zheng, S. R. Kulkarni, and H. V. Poor, “Attribute-distributed learning:
Models, limits, and algorithms,” IEEE Trans. Signal Process., vol. 59,
no. 1, pp. 386–398, January 2011.
[53] I. Csiszar and G. Tusnady, “Information geometry and alternating
minimization procedures,” Statist. Decision Supplement Issue., no. 1,
pp. 205–237, 1984.
[54] U. Niesen, D. Shah, and G. W. Wornell, “Adaptive alternating minimization algorithms,” IEEE Trans. Inf. Theory., vol. 55, no. 3, pp. 1423–1429,
March 2009.
[55] Q. Ling and A. Ribeiro, “Decentralized dynamic optimization through
the alternating direction method of multipliers,” IEEE Transactions on
Signal Processing, vol. 62, no. 5, pp. 1185–1197, March 2014.
[56] H. Bai, G. Li, S. Li, Q. Li, Q. Jiang, and L. Chang, “Alternating
optimization of sensing matrix and sparsifying dictionary for compressed
sensing,” IEEE Transactions on Signal Processing,, vol. 63, no. 6, pp.
1581–1594, March 2015.
[57] J. Yan, J. Wang, H. Zha, X. Yang, and S. Chu, “Consistency-driven
alternating optimization for multigraph matching: A unified approach,”
IEEE Transactions on Image Processing, vol. 24, no. 3, pp. 994–1009,
March 2015.
[58] S. Magnusson, P. Chathuranga, M. Rabbat, and C. Fischione, “On the
convergence of alternating direction lagrangian methods for nonconvex
structured optimization problems,” IEEE Transactions on Control of
Network Systems, vol. PP, no. 99, pp. 1–1, 2015.
[59] R. C. D. Lamare and R. Sampaio-Neto, “Minimum mean-squared error
iterative successive parallel arbitrated decision feedback detectors for dscdma systems,” IEEE Transactions on Communications, vol. 56, no. 5,
pp. 778–789, May 2008.
[60] P. Li, R. C. de Lamare, and R. Fa, “Multiple feedback successive
interference cancellation detection for multiuser mimo systems,” IEEE
Transactions on Wireless Communications, vol. 10, no. 8, pp. 2434–
2439, August 2011.
[61] R. C. de Lamare, “Adaptive and iterative multi-branch mmse decision
feedback detection algorithms for multi-antenna systems,” IEEE Transactions on Wireless Communications, vol. 12, no. 10, pp. 5294–5308,
October 2013.
[62] P. Clarke and R. C. de Lamare, “Transmit diversity and relay selection
algorithms for multirelay cooperative mimo systems,” IEEE Transactions on Vehicular Technology, vol. 61, no. 3, pp. 1084–1098, March
2012.
[63] T. Peng, R. C. de Lamare, and A. Schmeink, “Adaptive distributed spacetime coding based on adjustable code matrices for cooperative mimo
relaying systems,” IEEE Transactions on Communications, vol. 61,
no. 7, pp. 2692–2703, July 2013.
[64] T. Peng and R. C. de Lamare, “Adaptive buffer-aided distributed spacetime coding for cooperative wireless networks,” IEEE Transactions on
Communications, vol. 64, no. 5, pp. 1888–1900, May 2016.
[65] P. Li and R. C. de Lamare, “Distributed iterative detection with reduced
message passing for networked mimo cellular systems,” IEEE Transactions on Vehicular Technology, vol. 63, no. 6, pp. 2947–2954, July
2014.
[66] J. Gu, R. C. de Lamare, and M. Huemer, “Buffer-aided physicallayer network coding with optimal linear code designs for cooperative
networks,” IEEE Transactions on Communications, vol. PP, no. 99, pp.
1–1, 2017.
[67] A. Bertrand and M. Moonen, “Distributed adaptive node–specific signal
estimation in fully connected sensor networks–part II: Simultaneous and
asynchronous node updating,” IEEE Trans. Signal Process., vol. 58,
no. 10, pp. 5292–5306, October 2010.
[68] H. Qian and S. N. Batalama, “Data record–based criteria for the selection
of an auxiliary vector estimator of the MMSE/MVDR filter,” IEEE
Trans. Commun., vol. 51, no. 10, pp. 1700–1708, October 2003.
[69] R. C. de Lamare and R. Sampaio-Neto, “Reduced–rank space–time
adaptive interference suppression with joint iterative least squares algorithms for spread-spectrum systems,” IEEE Trans. Veh. Technol., vol. 59,
no. 3, pp. 1217–1228, March 2010.
[70] F. S. Cattivelli, C. G. Lopes, and A. H. Sayed, “Diffusion recursive
least-squares for distributed estimation over adaptive networks,” IEEE
Trans. Signal Process., vol. 56, no. 5, pp. 1865–1877, May 2008.

14

[71] M. L. J. van de Vel, Theory of Convex Structures. Elsevier, 1993.
[72] Y. Hua and M. Nikpour, “Computing the reduced-rank Wiener filter by
iqmd,” IEEE Signal Process. Lett., vol. 6, no. 9, pp. 240–242, September
1999.
[73] H. Ma, Y. Yang, Y. Chen, and K. J. R. Liu, “Distributed state estimation
in smart grid with communication constraints,” in Asia–Pacific Signal
Information Processing Association Annual Summit and Conference,
December 2012, pp. 1–4.

5

DRJIO−RLS(D=5)
Distributed Principal
Subspace Estimation(D=5)[11]
Krylov Subspace NLMS(D=5)[13]
Distributed RLS(M=20)[26]
Distributed NLMS(M=20)[2]
DRJIO−NLMS(D=5)

0

MSE(dB)

−5

−10

−15

−20

−25

−30
0

100

200

300

400
500
Time instant, i

600

700

800

900

