Dimensionality reduction aims at representing high-dimensional data in a lower dimensional space, which can be crucial to reduce the computational complexity of manipulating and processing this data, and is a core task in modern data analysis, machine learning, and related areas. The standard linear dimensionality reduction tool is principal component analysis (PCA), which allows expressing a high-dimensional dataset on the basis formed by the top eigenvectors of its sample covariance matrix. PCA first appeared in the statistics community in the beginning of the 20th century [3] and became one of the workhorses of statistical data analysis, with dimensionality reduction being a notable application. Nowadays, with the ever increasing collection of data by spatially dispersed agents, developing algorithms for distributed PCA constitutes a relevant area of research – see, e.g., [4, 5, 6, 7, 8, 9, 10, 11, 12] (master-slave communication architecture), and [13, 14, 15, 16, 17, 18, 19, 20] (arbitrarily meshed network communication architecture). For a recent and comprehensive review on these works, see, e.g., [21]; for a very recent work see [22].

In the (arbitrarily meshed network) distributed setting, consider a set of $N$ agents linked by an undirected and connected communication graph; the nodes are the agents, and the edges represent the communication channels between the agents. Each agent $n$ holds a finite set of points in $\mathbb{R}^{d}$, $\mathbf{Y}_{n}\subseteq\mathbb{R}^{d}$, and the agents seek to collectively find the top $m$ eigenvectors¹¹1This means the $m$ eigenvectors associated to the largest $m$ eigenvalues. of

(assumed to be positive definite, i.e., $C\succ 0$), where $M=\sum_{n=1}^{N}|\mathbf{Y}_{n}|$, i.e., the sum of the cardinalities of each $\mathbf{Y}_{n}$, and

Consider a collection of spatially distributed sensors monitoring the environment, a common scenario for information processing or decision making tasks see, e.g., [23, 24, 25, 26, 27, 28, 29, 30]. Often, these sensors communicate wirelessly, maybe in a harsh environment, which may result in faulty communications or sensor malfunctions [31]. The setup is modelled as follows: $N$ agents, linked by an undirected and connected communication graph, each holding an independent observation given by

In (1), $\mu^{\star}\in\mathbb{R}^{d}$ is a fixed and unknown parameter, each $h_{n}\in\mathbb{R}^{d}$ is assumed to be known only at agent $n$, $\{W_{j}\}_{n=1}^{N}$ are independent and identically distributed (i.i.d.) zero-mean Gaussian random variables with variance $(\sigma^{\star})^{2}$, and $\{Z_{n}\}_{n=1}^{N}$ are i.i.d. Bernoulli random variables ($Z_{n}\in\{0,1\}$) with $f_{Z_{n}}(z_{n}|p^{\star})=(p^{\star})^{z_{n}}(1-p^{\star})^{1-z_{n}}$. This formulation models a scenario where sensor $n$ measures the parameter $\mu^{\star}$ with probability $p^{\star}$ and, with probability $1-p^{\star}$, it senses only noise, indicating a transducer failure [31]. The agents seek to collectively estimate $\mu^{\star}$, treating $p^{\star}$ and $\sigma^{\star}$ as nuisance (or latent) parameters. Observe that if the binary variables $Z_{n}$ were not random, but fixed and known, the model could be regarded as a (distributed) linear regression problem. However, the randomness introduces an extra layer of difficulty accounting for potential sensor failures.

A decentralized algorithm, rather than one where each sensor sends its data to a central node, is potentially more robust to faulty wireless communications that may render a sensor useless. Moreover, a decentralized algorithm can yield considerable energy savings [23], a very desirable feature.

Our recent work [2] addressed a general distributed setup where $N$ agents, linked by a communication network, collaborate to collectively find an attractor $x^{\star}$ of a map $H$ that can be implicitly represented as an average of local maps, i.e.,

where $H_{n}$ is the map held by agent $n$. As defined in [2], an attractor $x^{\star}$ of $H$ is a fixed point thereof, $H(x^{\star})=x^{\star}$, satisfying

where $\rho\big{(}\mathbf{J}_{H}(x^{\star})\big{)}$ is the spectral radius of the Jacobian of $H$ at $x^{\star}$. Moreover, the map $H$ is not assumed to have a symmetric Jacobian and no global structural assumptions (e.g., Lipschitzianity or coercivity) are made.

The main contributions of [2] are a distributed algorithm to find $x^{\star}$ – DBPI – and the proof that it enjoys the local linear convergence of its centralized counterpart, i.e., of the (standard) Banach-Picard iteration: $x^{k+1}=H(x^{k})$.

In this work, we propose addressing the distributed inference problems described in Subsections I-A and I-B using two instantiations of DBPI. More concretely, we propose:

1. 1.

   A distributed algorithm for PCA, which results from considering a map that can be implicitly written as an average of local maps and that has as a fixed point the solution to the PCA problem.

2. 2.

   An algorithm that stems from formulating the problem described in Subsection I-B as a fixed point of a map induced by the stationary equations of the corresponding maximum likelihood estimation (MLE) criterion. This map corresponds to the iterations of a slightly modified EM algorithm for a mixture of linear regressions [32].

The guarantees of local linear convergence for these distributed algorithms involve verifying condition (2) for the maps inducing them, which allows invoking the results from [2]. Consequently, a great portion of this paper is devoted to proving that (2) holds for these maps, which is far from trivial.

The distributed PCA problem (see [23] for a review of distributed PCA) described in Subsection I-A was addressed in [33] , where an algorithm termed accelerated distributed Sanger’s algorithm (ADSA) was proposed. The authors consider a “mini-batch variant” of Sanger’s algorithm (SA, see [34]) and, inspired by [35], arrive at ADSA. Although no proof of convergence was presented in [33], very recent work by the same authors proves convergence of their algorithm [36]. Our contributions in this context are twofold: we show that ADSA is recovered by applying DBPI to SA, and that condition (2) holds for SA, thus, the guarantees of local convergence follow directly as a consequence of the results in [2]. We mention that no computer simulations of ADSA are presented in this work, since these can be found in [33].

The problem presented in Subsection I-B was addressed in [31], where (1) is regarded as a finite mixture model [37]. To estimate $\mu^{\star}$, the authors proposed a distributed version of the expectation-maximization (EM [38]) algorithm, termed diffusion-averaging distributed EM (DA-DEM). However, DA-DEM, very much in the spirit of [39, 40], uses a diminishing step-size to achieve convergence, leading to a sublinear convergence rate. Our contribution is an algorithm for this problem that extends a slightly modified version of the centralized EM algorithm to distributed settings. The key challenge is to show that we can “expect” condition (2) to hold, and we dedicate a considerable amount of effort to this endeavor. We use the term “expect”, since the operator underlying DBPI depends on the observed samples and, therefore, the existence of an attractor is a probabilistic question. Finally, we compare our algorithm with DA-DEM through Monte Carlo simulations, confirming the linear convergence rate of our algorithm and the sublinear convergence of DA-DEM.

There is considerable work on the “probabilistic linear convergence” of EM [41], [42], [43]. However, neither the results in [41], nor those in [42] encompass the mixture model underlying (1). The mixture of regressions presented in [43] bears some similarity with the model underlying (1), but it is not the same: in [43], $p$ is fixed at $1/2$ and $Z_{n}\in\{-1,1\}$ (rather than $\{0,1\}$), thus there are no measurements that are just noise. Furthermore, [43] is primarily concerned with statistical guarantees for the error with respect to the ground truth, while we address the goal of establishing (2).

As mentioned in [31], there are two other relevant works on distributed EM, namely, [44] and [45]. However (see [31]), both these works address a different problem of Gaussian mixture density estimation. Moreover, in the case of [44], the algorithm demands a cyclic network topology, and, in [45] the algorithm requires higher computational load on each node, since it is based on alternating direction method of multipliers.

To summarize, we show that ADSA [33] is an instance of DBPI and propose an algorithm to solve the mixture model underlying (1), also as an instance of DBPI. Consequently, their corresponding guarantees of local linear convergence result from the attractor condition (2) for the map underlying the corresponding centralized counterparts. We compare DA-DEM and our proposed algorithm through numerical Monte Carlo simulations, and the results confirm the linear convergence of our algorithm and the sublinear convergence of DA-DEM.

Section II briefly reviews the DBPI proposed in [2] and the main convergence result therein proved. The characterization of the fixed points of the “mini-batch” variant of Sanger’s algorithm, as well as the attractor condition, are presented in Section III. Section IV describes the centralized variant of EM underlying the proposed distributed algorithm for the problem in Subsection I-B, presents the verification of the attractor condition, and reports the results of simulations comparing DBPI with DA-DEM.

The set of real $n$ dimensional vectors with positive components is denoted by $\mathbb{R}^{n}_{>0}$. Matrices and vectors are denoted by upper and lower case letters, respectively. The spectral radius of a matrix $A$ is denoted by $\rho(A)$ and its Frobenius norm by $\|A\|_{F}$. Given a map $H$, $\mathbf{J}_{H}(x)$, and $dH(x)$ denote, respectively, the Jacobian of $H$ at $x$ and the differential of $H$ at $x$. Given a vector $v$, $v_{s}$ denotes its $s$th component; given a matrix $A$, $A_{st}$ denotes the element on the $s$th row and $t$th column and $A^{T}$ its transpose. The $d$-dimensional identity matrix is denoted by $I_{d}$, $\mathbf{1}_{d}$ is the $d$-dimensional vector of ones, and $\mathbf{0}_{m,n}$ is the $m\times n$ matrix of zeros. Whenever convenient, we will denote a vector with two stacked blocks, $[v^{T},u^{T}]^{T}$, simply as $(v,u)$. Given a square matrix $A$, $\mathcal{U}(A)$ is an upper triangular matrix of the same dimension as $A$ and whose upper triangular part coincides with that of $A$. Given a norm $\|\cdot\|$, $\bar{B}^{\|\cdot\|}_{\delta,\theta}$ denotes the closed ball of center $\theta$ and radius $\delta$ with respect to $\|\cdot\|$. Random variables and vectors are denoted by upper case letters and, for random variable $Y$, the probability density (or mass) function of $Y$ is denoted by $f_{Y}$. The probability density of a Gaussian of mean $\mu$ and variance $\sigma^{2}$ is denoted by $\mathcal{N}(\cdot|\mu,\sigma^{2})$.

We obtain a distributed algorithm for solving the PCA problem described in section I-A as an instantiation of DBPI by introducing a map $H$ with a fixed point at the desired solution. Moreover, the guarantees of local linear convergence follow as a result of verifying (2).

The “mini batch variant” of Sanger’s algorithm (SA) proposed in [33] and inspired by [34] is the Banach-Picard iteration

where $H:\mathbb{R}^{d\times m}\to\mathbb{R}^{d\times m}$ is given by

and $\mathcal{U}$ was defined in Subsection I-F. Observe that $H$ can be written as an average of local maps, i.e.,

where $H_{n}:\mathbb{R}^{d\times m}\to\mathbb{R}^{d\times m}$ is defined by

Let $R$ be the map on $\mathbb{R}^{(d\times m)N}$ defined, for $z=[z_{1}^{T},\ldots,z_{N}^{T}]^{T}$, as in (4), with $H_{n}$ as in (6). The distributed algorithm presented in [33], named ADSA, is exactly the DBPI, i.e., Algorithm 1, with this choice of $R$.

The convergence analysis amounts to verifying the attractor condition (2) for $H$, thus establishing, as a corollary of the results in [2], the local linear convergence of Algorithm 1 with each $H_{n}$ in (4) defined as in (6) (equivalently, ADSA).

We start with the following lemma (proved in Appendix D) showing that the solution sought in the PCA problem is a fixed point of $H$ (as defined in (5)).

The following theorem guarantees that the Banach-Picard iteration of $H$ has local linear convergence to its fixed points.

First, note that $H(X)=I+\eta S(X)$, where

This implies $\mathbf{J}_{H}(X^{\star})=I+\eta\mathbf{J}_{S}(X^{\star})$, and, as a consequence, each eigenvalue of $\mathbf{J}_{H}(X^{\star})$ is of the form $1+\eta\beta$, with $\beta$ being an eigenvalue of $\mathbf{J}_{S}(X^{\star})$. The idea is to show that these eigenvalues $\beta$ of $\mathbf{J}_{S}(X^{\star})$ enjoy a key property: they are real-valued and negative, $\beta<0$. Such property means that, for sufficiently small $\eta>0$, we have $|1+\eta\beta|<1$. To establish this key property, we divide the proof of Theorem 1 in two lemmas: Lemma 2 and Lemma 3.

Lemma 2 will show that the eigenvalues of $\mathbf{J}_{S}(X^{\star})$ coincide with those of the linear map from $\mathbb{R}^{d\times m}$ to $\mathbb{R}^{d\times m}$ given by

where

and

Lemma 3 will show that the eigenvalues of (8) are real and negative.

In the proof of the following lemma it is crucial that the eigenvalues are in decreasing order.

This is a rather long section, hence the need for a road map. The analysis of (1) is simplified if the measurements are identically distributed besides being just independent and, therefore, we start by introducing a probability distribution on the vectors $h_{n}$ and a joint model on $(Y,H)$.

To estimate $\mu^{\star}$ ($p^{\star}$ and $\sigma^{\star}$ are treated as nuisance parameters), we consider the stationary equations imposed by equating to zero the gradient of the log-likelihood function, a necessary condition satisfied by the maximum likelihood estimator (MLE). Once the particular form of the stationary equations is realized, we reformulate them as a fixed point equation of the form $g_{1}\circ g_{2}(\theta^{\star})=\theta^{\star}$ that naturally suggests the Banach-Picard iteration

Observing that the map $g_{1}\circ g_{2}$ cannot be written as an average of local maps, we switch to the map $H=g_{2}\circ g_{1}$, which can be implicitly written as an average of local maps. With this map, we arrive at a distributed algorithm by considering the map $R$ (see section II) arising from $H$ and appealing to Algorithm 1.

Finally, we observe that the existence of a fixed point of $H=g_{2}\circ g_{1}$ satisfying (2) follows from the existence of a fixed point of $g_{1}\circ g_{2}$ satisfying (2). The final part of the section is thus devoted to verifying (2) for the map $g_{1}\circ g_{2}$, and to a numerical simulation comparing Algorithm 1 and DA-DEM from [31].

Let $\theta^{\star}=(\mu^{\star},p^{\star},(\sigma^{\star})^{2})\in\Theta=\mathbb{R}^{d}\times(0,1)\times(0,+\infty)$ be an unknown and fixed vector which we term the ground truth.

The agents’ measurements are assumed to be independent (see (1)); however, they are not identically distributed, given the presence of the vectors $h_{n}$ in (1). To address this issue, let $Z\in\{0,1\}$, $H\in\mathbb{R}^{d}$, and $Y\in\mathbb{R}$ be, respectively, a binary random variable, a random vector, and a real random variable. Suppose the joint density on $(Y,H,Z)$ factors as

where

and

Instead of assuming that the $h_{n}$ are fixed as in [31], we assume that each sensor $n$ has a measurement $(y_{n},h_{n})$, where $(y_{n},h_{n},z_{n})$ was drawn from (20), but agent $n$ has no knowledge of $z_{n}$. After marginalization, the joint density of $Y,H$ is given by

which is a mixture model [49].

To estimate $\mu^{\star}$, the agents seek $\theta\in\Theta$ such that

where $\phi$ is the log-likelihood of $(Y,H)$, i.e.,

Since $f_{H}(h)$ does not depend on $\theta$,

in other words, (23) is a necessary condition satisfied by the MLE corresponding to the log-likelihood function $\log f_{Y|H}(y|h,\theta)$, thus independent of $f_{H}$.

Before explicitly writing the stationary equations corresponding to (23), we introduce the responsibility functions [49],

Notice that $r(y,h,\theta)=\mathbb{P}(z=1|y,h,\theta)$, the posterior probability that the observation $y$ was not a result of measuring only noise.

For reasons that will be clear later, the following set of equalities, which can be easily verified, will be convenient:

Using (26), (23) can be explicitly written as

where

If the matrix $\frac{1}{N}\sum_{n=1}^{N}\Gamma(y_{n},h_{n},\theta)$ is invertible, then (27)-(29) can be written as a fixed point equation.²²2The invertibility of this matrix is assumed throughout the rest of the paper. In fact, if $N$ is sufficiently large - greater than $d$ - this happens with probability one. This constitutes the motivation for the centralized algorithm that we suggest next (see Algorithm 2) and from which we will derive the distributed version; observe that it is the Banach-Picard iteration motivated by (27)-(29).

Another way to write (30)-(32) (see Algorithm 2 below) is $\theta^{k+1}=g_{1}\circ g_{2}(\theta^{k})$, where

and

Although, $g_{2}$ is an average of local maps, the map $g_{1}\circ g_{2}$ is not, due to the matrix inversion in (30). As a consequence, (30)-(32) cannot be directly extended to a distributed algorithm. However, switching the order of $g_{1}$ and $g_{2}$ results in a map that can be implicitly written as an average of local maps. In fact, instead of the iteration $\theta^{k+1}=g_{1}\circ g_{2}(\theta^{k})$, consider the iteration

where $H=g_{2}\circ g_{1}$ and $z=(\Gamma,\psi,p,\sigma^{2})$.

Let

and it follows that $H=\frac{1}{N}\sum_{n=1}^{N}H_{n}$, where

To conclude, the distributed algorithm we suggest amounts to Algorithm 1, with $R:\mathbb{R}^{(d^{2}+d+2)N}\to\mathbb{R}^{(d^{2}+d+2)N}$ defined, for $z=(z_{1}^{T},\ldots,z_{N}^{T})^{T}$, as in (4), and $H_{n}$ as in (33). Additionally, following [31], we suggest the initialization

Some remarks are due:

a)

The existence of a fixed point of $g_{1}\circ g_{2}$ satisfying (2) is addressed in section IV-E;

b)

The existence of a fixed point of $g_{2}\circ g_{1}$ satisfying (2) follows from the existence of a fixed point of $g_{1}\circ g_{2}$ satisfying (2), by the chain rule;

c)

Expanding (32) yields

	$\displaystyle(\sigma^{k+1})^{2}$	$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}r(y_{n},h_{n},\theta^{k})(y_{n}-h_{n}^{T}\mu^{k})^{2}$
		$\displaystyle+\big{(}1-r(y_{n},h_{n},\theta^{k})\big{)}y_{n}^{2},$

and, if the update rule is modified according to

	$\displaystyle(\sigma^{k+1})^{2}$	$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}r(y_{n},h_{n},\theta^{k})(y_{n}-h_{n}^{T}\mu^{k+1})^{2}$
		$\displaystyle+\big{(}1-r(y_{n},h_{n},\theta^{k})\big{)}y_{n}^{2},$

then, a straightforward manipulation recovers the EM algorithm presented in [31]. Moreover, the EM algorithm derived in [31] is still amenable to a distributed implementation using Algorithm 1. However, we found it easier to prove (2) for this variant of EM, than for the standard EM.

The proof of local linear convergence of the centralized variant of EM, i.e., Algorithm 2, is not trivial. In fact, this question is probabilistic in nature, because updates (30)-(32) depend on observations that are, in turn, samples from a probability distribution. Before presenting the main convergence result (Theorem 3), we need to introduce some definitions and only one mild assumption that is instrumental in the proof of Lemma 4 below: the Fisher information at $\theta^{\star}$, given by,

is non-singular.

Let $T_{N}=g_{1}\circ g_{2}$ denote the map underlying the Banach-Picard iteration (30)-(32).³³3The subscript $N$ emphasizes that $T_{N}$ depends on $N$ observations. A straightforward manipulation, using (26), shows that

where