Adaptation and learning over networks under subspace constraints -- Part II: Performance Analysis

TL;DR

This paper analyzes the steady-state performance of a distributed learning algorithm over networks with subspace constraints, showing it approaches centralized performance under small step-sizes and highlighting the effects of noise and data characteristics.

Contribution

It provides explicit analysis of how gradient noise, data, and subspace constraints affect network performance, extending previous stability results to steady-state behavior.

Findings

Distributed algorithm achieves centralized steady-state performance for small step-sizes.

Gradient noise and data characteristics significantly influence network estimation errors.

Explicit relationships between subspace constraints and network performance are established.

Abstract

Part I of this paper considered optimization problems over networks where agents have individual objectives to meet, or individual parameter vectors to estimate, subject to subspace constraints that require the objectives across the network to lie in low-dimensional subspaces. Starting from the centralized projected gradient descent, an iterative and distributed solution was proposed that responds to streaming data and employs stochastic approximations in place of actual gradient vectors, which are generally unavailable. We examined the second-order stability of the learning algorithm and we showed that, for small step-sizes $\mu$, the proposed strategy leads to small estimation errors on the order of $\mu$. This Part II examines steady-state performance. The results reveal explicitly the influence of the gradient noise, data characteristics, and subspace constraints, on the network performance. The results also show that in the small step-size regime, the iterates generated by the distributed algorithm achieve the centralized steady-state performance.

Tables

Table 1. TABLE I: Definition of some variables used throughout the analysis. ℐ ℐ \mathcal{I} is a permutation matrix defined by ( 16 ).

Variable	Real data case	Complex data case
Data-type variable $h$	1	2
Gradient vector	${\nabla }_{w_k^{\top }}{J}_k\left(w_k\right)$	${\nabla }_{w_k^{\ast }}{J}_k\left(w_k\right)$
Error vector ${\widetilde{𝒘}}_{k,i}^e$	${\widetilde{𝒘}}_{k,i}$ from (39)	$\left[\begin{array}{c}{\widetilde{𝒘}}_{k,i}\\ {\left({\widetilde{𝒘}}_{k,i}^{\ast }\right)}^{\top }\end{array}\right]$
Gradient noise ${𝒔}_{k,i}^e\left(w\right)$	${𝒔}_{k,i}\left(w\right)$ from (4)	$\left[\begin{array}{c}{𝒔}_{k,i}\left(w\right)\\ {\left({𝒔}_{k,i}^{\ast }\left(w\right)\right)}^{\top }\end{array}\right]$
Bias vector $b_k^e$	$b_k$ from (40)	$\left[\begin{array}{c}{b}_k\\ {\left(b_k^{\ast }\right)}^{\top }\end{array}\right]$
$(k,\ell )$-th block of ${𝒜}^e$	$A_{k\ell }$
Matrix ${𝒰}^e$	$𝒰$
Matrix ${𝒥}_{ϵ}^e$	${𝒥}_{ϵ}$ from (28)
Matrix ${𝒱}_{R,ϵ}^e$	${𝒱}_{R,ϵ}$ from (28)
Matrix ${\left({𝒱}_{L,ϵ}^e\right)}^{\ast }$	${𝒱}_{L,ϵ}^{\ast }$ from (28)
Noise covariance $R_k^o$	$R_{q,k}$ from (51),(52)

Table 2. TABLE II: Performance of strategy ( 2 ) w.r.t. 𝒲 o superscript 𝒲 𝑜 {\scriptstyle\mathcal{W}}^{o} in ( 1 ) for 2 2 2 different choices of μ 𝜇 \mu .

		MSD
Step-size $\mu$	Solution	Exp. (67)	Exp. (64)	Simulation
	Centralized	$-29.66$dB	$-29.74$dB	$-29.603$dB
${10}^{-2}$	Distributed	$-29.66$dB	$-27.6$dB	$-27.298$dB
	Centralized	$-39.66$dB	$-39.67$dB	$-39.691$dB
${10}^{-3}$	Distributed	$-39.66$dB	$-39.26$dB	$-39.196$dB
	Centralized	$-49.66$dB	$-50.19$dB	$-49.772$dB
${10}^{-4}$	Distributed	$-49.66$dB	$-50.14$dB	$-49.335$dB

Equations

$$\begin{array}{rl}{\mathcal{W}}^o=& \arg \underset{\mathcal{W}}{\min }{J}^{\text{glob}}(\mathcal{W})\triangleq \sum _{k=1}^N{J}_k(w_k),\\ & \mathrm{subject}\mathrm{to}\mathcal{W}\in \mathcal{R}(\mathcal{U}),\end{array}$$

$$\left\{\begin{array}[]{rl}\boldsymbol{\psi}_{k,i}=&\boldsymbol{w}_{k,i-1}-\mu\widehat{\nabla_{w_{k}^{*}}J_{k}}(\boldsymbol{w}_{k,i-1}),\\ \boldsymbol{w}_{k,i}=&\sum\limits_{\ell\in\mathcal{N}_{k}}A_{k\ell}\boldsymbol{\psi}_{\ell,i},\end{array}\right.$$

$$\widehat{{\nabla }_{w_k^{\ast }}{J}_k}(w_k)={\nabla }_{w_k^{\ast }}{Q}_k(w_k;{\mathbf{x}}_{k,i}),$$

$${\mathbf{s}}_{k,i}(w)\triangleq {\nabla }_{w_k^{\ast }}{J}_k(w)-\widehat{{\nabla }_{w_k^{\ast }}{J}_k}(w).$$

$$\underset{i\to \infty }{\lim }{\mathcal{A}}^i={\mathcal{P}}_{\mathcal{U}},$$

$$A_{k\ell }=[\mathcal{A}{]}_{k\ell }=0,\text{if }\ell \notin {\mathcal{N}}_k\text{ and }k\ne \ell ,$$

$${\mathcal{P}}_{\mathcal{U}}=\mathcal{U}({\mathcal{U}}^{\ast }\mathcal{U}{)}^{-1}{\mathcal{U}}^{\ast }.$$

$$\mathcal{A}\mathcal{U}$$

$${\mathcal{U}}^{\ast }\mathcal{A}$$

$$\rho (\mathcal{A}-{\mathcal{P}}_{\mathcal{U}})$$

$$\begin{array}[]{cl}\text{find}&\mathcal{A}\\ \text{such that}&\mathcal{A}\,\mathcal{U}=\mathcal{U},~{}~{}\mathcal{U}^{*}\mathcal{A}=\mathcal{U}^{*},\\ &\rho(\mathcal{A}-\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}})<1,\\ &[\mathcal{A}]_{k\ell}=0,~{}\text{if }\ell\notin\mathcal{N}_{k}\text{ and }\ell\neq k.\end{array}$$

$$\text{MSD}\triangleq \mu \underset{\mu \to 0}{\lim }\left(\underset{i\to \infty }{\mathrm{limsup}}\frac{1}{\mu }\mathbb{E}\left(\frac{1}{N}\Vert {\mathcal{W}}^o-{\mathcal{W}}_i{\Vert }^2\right)\right),$$

$${\mathcal{W}}_i^c={\mathcal{P}}_{\mathcal{U}}\left({\mathcal{W}}_{i-1}^c-\mu \text{col}{\left\{\widehat{{\nabla }_{w_k^{\ast }}{J}_k}({\mathbf{w}}_{k,i-1}^c)\right\}}_{k=1}^N\right),i\ge 0,$$

$$\underset{i\to \infty }{\mathrm{limsup}}\mathbb{E}\Vert w_k^o-{\mathbf{w}}_{k,i}{\Vert }^2=O(\mu ),k=1,\dots ,N,$$

$$\text{MSD}=\frac{\mu }{2hN}\text{Tr}\left({\left(({\mathcal{U}}^e{)}^{\ast }{\mathcal{H}}^o{\mathcal{U}}^e\right)}^{-1}\left(({\mathcal{U}}^e{)}^{\ast }\mathcal{S}{\mathcal{U}}^e\right)\right).$$

$$[\mathcal{I}]_{mn}\triangleq\left\{\begin{array}[]{ll}I_{M_{k}},&\text{if }m=k,n=2(k-1)+1\\ I_{M_{k}},&\text{if }m=k+N,n=2k\\ 0,&\text{otherwise}\\ \end{array}\right.$$

$$H_k(w_k)$$

$$\mathcal{H}(\mathcal{W})$$

$$\frac{{\nu }_k}{h}{I}_{h{M}_k}\le H_k(w_k)\le \frac{{\delta }_k}{h}{I}_{h{M}_k},$$

$$0<\frac{\nu }{h}{I}_{hP}\le ({\mathcal{U}}^e{)}^{\ast }\mathcal{H}(\mathcal{W}){\mathcal{U}}^e\le \frac{\delta }{h}{I}_{hP},$$

$$\mathbb{E}[{\mathbf{s}}_{k,i}(\mathbf{w})|{\mathcal{F}}_{i-1}]$$

$$\mathbb{E}[{\mathbf{s}}_{k,i}(\mathbf{w}){\mathbf{s}}_{\ell ,i}^{\ast }(\mathbf{w})|{\mathcal{F}}_{i-1}]$$

$$\mathbb{E}[{\mathbf{s}}_{k,i}(\mathbf{w}){\mathbf{s}}_{\ell ,i}^{\top }(\mathbf{w})|{\mathcal{F}}_{i-1}]$$

$$\mathbb{E}[\Vert {\mathbf{s}}_{k,i}(\mathbf{w}){\Vert }^2|{\mathcal{F}}_{i-1}]$$

$${\mathcal{A}}^e\triangleq {\mathcal{V}}_{ϵ}^e{\Lambda }_{ϵ}^e({\mathcal{V}}_{ϵ}^e{)}^{-1},$$

$$\Lambda_{\epsilon}^{e}=\left[\begin{array}[]{cc}I_{hP}&0\\ 0&\mathcal{J}^{e}_{\epsilon}\end{array}\right],~{}\mathcal{V}^{e}_{\epsilon}=[\mathcal{U}^{e}~{}\mathcal{V}^{e}_{R,\epsilon}],~{}(\mathcal{V}^{e}_{\epsilon})^{-1}=\left[\begin{array}[]{c}(\mathcal{U}^{e})^{*}\\ (\mathcal{V}_{L,\epsilon}^{e})^{*}\end{array}\right]$$

$$\mathcal{A}\triangleq {\mathcal{V}}_{ϵ}{\Lambda }_{ϵ}{\mathcal{V}}_{ϵ},$$

$$\Lambda_{\epsilon}=\left[\begin{array}[]{cc}I_{P}&0\\ 0&\mathcal{J}_{\epsilon}\end{array}\right],~{}\mathcal{V}_{\epsilon}=\left[\begin{array}[]{cc}\mathcal{U}&\mathcal{V}_{R,\epsilon}\end{array}\right],~{}\mathcal{V}_{\epsilon}^{-1}=\left[\begin{array}[]{c}\mathcal{U}^{*}\\ \mathcal{V}_{L,\epsilon}^{*}\end{array}\right],$$

$$\mathbb{E}\left[\Vert {\mathbf{s}}_{k,i}(\mathbf{w}){\Vert }^4|{\mathcal{F}}_{i-1}\right]\le ({\beta }_{4,k}/h{)}^4\Vert \mathbf{w}{\Vert }^4+{\sigma }_{s4,k}^4,$$

$$\underset{i\to \infty }{\mathrm{limsup}}\mathbb{E}\Vert w_k^o-{\mathbf{w}}_{k,i}{\Vert }^4=O({\mu }^2),k=1,\dots ,N.$$

$$\Vert {\nabla }_{w_k}{J}_k(w_k^o+\Delta w_k)-{\nabla }_{w_k}{J}_k(w_k^o)\Vert \le {\kappa }_d\Vert \Delta w_k\Vert ,$$

$$\framebox{$\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}_{i}^{e}=\boldsymbol{\cal{B}}_{i-1}\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e}_{i-1}-\mu\mathcal{A}^{e}\boldsymbol{s}_{i}^{e}+\mu\mathcal{A}^{e}b^{e}$}$$

$${\tilde{\mathcal{W}}}_i^e$$

$${\mathcal{H}}_{i-1}$$

Full text

I Introduction

As pointed out in Part I [2] of this work, most prior literature on distributed inference over networks focuses on consensus problems, where agents with separate objective functions need to agree on a common parameter vector corresponding to the minimizer of the aggregate sum of individual costs [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. In this paper, and its accompanying Part I [2], we focus instead on multitask networks where the agents need to estimate and track multiple objectives simultaneously [13, 14, 15, 16, 17, 18, 19, 20]. Based on the type of prior information that may be available about how the tasks are related to each other, multitask algorithms can be derived by translating the prior information into constraints on the parameter vectors to be inferred.

In this paper, and the accompanying Part I [2], we consider multitask inference problems where each agent seeks to minimize an individual cost, and where the collection of parameter vectors to be estimated across the network is required to lie in a low-dimensional subspace. That is, we let $w_{k}\in{\mathbb{C}^{M_{k}}}$ denote some parameter vector at node $k$ and let ${\scriptstyle\mathcal{W}}=\text{col}\{w_{1},\ldots,w_{N}\}$ denote the collection of parameter vectors from across the network ($N$ is the number of agents in the network). We associate with each agent $k$ a differentiable convex cost $J_{k}(w_{k}):{\mathbb{C}^{M_{k}}}\rightarrow\mathbb{R}$, which is expressed as the expectation of some loss function $Q_{k}(\cdot)$ and written as $J_{k}(w_{k})=\mathbb{E}Q_{k}(w_{k};\boldsymbol{x}_{k})$, where $\boldsymbol{x}_{k}$ denotes the random data. The expectation is computed over the distribution of the data. Let $M=\sum_{k=1}^{N}M_{k}$. We consider constrained problems of the form:

[TABLE]

where $\mathcal{R}(\cdot)$ denotes the range space operator, and $\mathcal{U}$ is an $M\times P$ full-column rank matrix with $P\ll M$. Each agent $k$ is interested in estimating the $k$-th $M_{k}\times 1$ subvector $w^{o}_{k}$ of ${\scriptstyle\mathcal{W}}^{o}=\text{col}\{w^{o}_{1},\ldots,w^{o}_{N}\}$.

In order to solve problem (1), we proposed in Part I [2] the following adaptive and distributed strategy:

[TABLE]

where $\mu>0$ is a small step-size parameter, $\boldsymbol{\psi}_{k,i}$ is an intermediate estimate, $\boldsymbol{w}_{k,i}$ is the estimate of $w^{o}_{k}$ at agent $k$ and iteration $i$, $\mathcal{N}_{k}$ denotes the neighborhood of agent $k$, and ${\nabla_{w_{k}^{*}}J_{k}}(\cdot)$ is the (Wirtinger) complex gradient [4, Appendix A] of $J_{k}(\cdot)$ relative to $w^{*}_{k}$ (complex conjugate of $w_{k}$). Notice that approximate gradient vectors $\widehat{\nabla_{w_{k}^{*}}J_{k}}(\cdot)$ are employed in (2) instead of true gradient vectors ${\nabla_{w_{k}^{*}}J_{k}}(\cdot)$ since we are interested in solving (1) in the stochastic setting when the distribution of the data $\boldsymbol{x}_{k}$ is unknown. A common construction in stochastic approximation theory is to employ the following approximation at iteration $i$:

[TABLE]

where $\boldsymbol{x}_{k,i}$ represents the data observed at iteration $i$. The difference between the true gradient and its approximation is called the gradient noise $\boldsymbol{s}_{k,i}(\cdot)$:

[TABLE]

This noise will seep into the operation of the algorithm and one main challenge is to show that despite its presence, agent $k$ is still able to approach $w^{o}_{k}$ asymptotically. The matrix $A_{k\ell}$ appearing in (2) is of size $M_{k}\times M_{\ell}$. It multiplies the intermediate estimate $\boldsymbol{\psi}_{\ell,i}$ arriving from neighboring agent $\ell$ to agent $k$. Let $\mathcal{A}\triangleq[A_{k\ell}]\in\mathbb{C}^{M\times M}$ denote the matrix that collects all these blocks. This $N\times N$ block matrix is chosen by the designer to satisfy the following two conditions:

[TABLE]

where $[\mathcal{A}]_{k\ell}$ denotes the $(k,\ell)$-th block of $\mathcal{A}$ and $\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}$ is the projector onto the $P$-dimensional subspace of $\mathbb{C}^{M}$ spanned by the columns of $\mathcal{U}$:

[TABLE]

The sparsity condition (6) characterizes the network topology and ensures local exchange of information at each time instant $i$. It is shown in Part I [2] that the matrix equation (5) holds, if and only if, the following conditions on the projector $\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}$ and the matrix $\mathcal{A}$ are satisfied:

[TABLE]

where $\rho(\cdot)$ denotes the spectral radius of its matrix argument. Conditions (8) and (9) state that the $P$ columns of $\mathcal{U}$ are right and left eigenvectors of $\mathcal{A}$ associated with the eigenvalue $1$. Together with these two conditions, condition (10) means that $\mathcal{A}$ has $P$ eigenvalues at one, and that all other eigenvalues are strictly less than one in magnitude. Combining conditions (8)–(10) with the sparsity condition (6), the design of a matrix $\mathcal{A}$ to run (2) can be written as the following feasibility problem:

[TABLE]

Not all network topologies satisfying (6) guarantee the existence of an $\mathcal{A}$ satisfying condition (5). The higher the dimension of the signal subspace is, the greater the graph connectivity has to be. In the works [1, 20], it is assumed that the sparsity constraints (6) and the signal subspace lead to a feasible problem. That is, it is assumed that problem (11) admits at least one solution. As a remedy for the violation of such assumption, one may increase the network connectivity by increasing the transmit power of each node, i.e., adding more links [20]. In Section III of this part, we shall relax the feasibility assumption by considering the problem of finding an $\mathcal{A}$ that minimizes the number of edges to be added to the original topology while satisfying the constraints (8), (9), and (10). In this case, if the original topology leads to a feasible solution, then no links will be added. Otherwise, we assume that the designer is able to add some links to make the problem feasible.

When studying the performance of algorithm (2) relative to ${\scriptstyle\mathcal{W}}^{o}$, we assume that a feasible $\mathcal{A}$ (topology) is computed by the designer and that its blocks $\{A_{k\ell}\}_{\ell\in\mathcal{N}_{k}}$ are provided to agent $k$ in order to run (2). We carried out in Part I [2] a detailed stability analysis of the proposed strategy (2). We showed that, despite the gradient noise, the distributed strategy (2) is able to converge in the mean-square-error sense within $O(\mu)$ from the solution of the constrained problem (1), for sufficiently small step-sizes $\mu$. We particularly established that, for each agent $k$, the error variance relative to $w^{o}_{k}$ enters a bounded region whose size is in the order of $\mu$, namely, $\limsup_{i\rightarrow\infty}\mathbb{E}\|w^{o}_{k}-\boldsymbol{w}_{k,i}\|^{2}=O(\mu)$. In Section II of this Part II, we will assess the size of this mean-square-error by deriving closed-form expression for the network mean-square-deviation (MSD) defined by [4]:

[TABLE]

where $\boldsymbol{{\scriptstyle\mathcal{W}}}_{i}\triangleq\text{col}\{\boldsymbol{w}_{k,i}\}_{k=1}^{N}$. In other words, we will assess the size of the constant multiplying $\mu$ in the $O(\mu)-$term. This closed form expression will reveal explicitly the influence of the data characteristics (captured by the second-order properties of the costs and second-order moments of the gradient noises) and subspace constraints (captured by $\mathcal{U}$), on the network performance. In this way, we will be able to conclude that distributed strategies of the form (2) with small step-sizes are able to lead to reliable performance even in the presence of gradient noise. We will be able also to conclude that the iterates generated by the distributed implementation achieve the centralized steady-state performance. Particularly, we compare the performance of strategy (2) to the following centralized stochastic gradient projection algorithm [21]:

[TABLE]

where $\boldsymbol{{\scriptstyle\mathcal{W}}}^{{c}}_{i}=\text{col}\{\boldsymbol{w}^{c}_{1,i},\ldots,\boldsymbol{w}^{c}_{N,i}\}$ is the estimate of ${\scriptstyle\mathcal{W}}^{o}$ at iteration $i$. Observe that each agent at each iteration needs to send its data to a fusion center, which performs the projection in (13), and then sends the resulting estimates $\boldsymbol{w}^{c}_{k,i}$ back to the agents. Finally, simulations will be provided in Section IV to verify the theoretical findings.

II Stochastic performance analysis

In Part I [2], we carried out a detailed stability analysis of the proposed strategy (2). We showed, under some Assumptions on the risks $\{J_{k}(\cdot)\}$ and on the gradient noise processes $\{\boldsymbol{s}_{k,i}(\cdot)\}$ defined by (4), that a network running strategy (2) with a matrix $\mathcal{A}$ satisfying conditions (6), (8), (9), and (10) is mean-square-error stable for sufficiently small step-sizes, namely, it holds that:

[TABLE]

for small enough $\mu$–see [2, Theorem 1]. Expression (14) indicates that the mean-square error $\mathbb{E}\|{\scriptstyle\mathcal{W}}^{o}-\boldsymbol{{\scriptstyle\mathcal{W}}}_{i}\|^{2}$ is on the order of $\mu$. In this section, we are interested in characterizing how close the $\boldsymbol{{\scriptstyle\mathcal{W}}}_{i}$ gets to the network limit point ${\scriptstyle\mathcal{W}}^{o}$. In particular, we will be able to characterize the network mean-square deviation (MSD) (defined by (12)) value in terms of the step-size $\mu$, the data-type variable $h$ defined in Table I, the second-order properties of the costs (captured by $\mathcal{H}^{o}$ defined below in (45)), the second-order moments of the gradient noises (captured by $\mathcal{S}$ defined below in (66)), and the subspace constraints (captured by $\mathcal{U}^{e}$ defined in Table I) as follows:

[TABLE]

As explained in Part I [2], in the general complex data case, extended vectors and matrices need to be introduced in order to analyze the network evolution. The arguments and results presented in this section are applicable to both cases of real and complex data through the use of data-type variable $h$. Table I lists a couple of variables and symbols that will be used in the sequel for both real and complex data cases. The matrix $\mathcal{I}$ in Table I is a permutation matrix of $2N\times 2N$ blocks with $(m,n)$-th block given by:

[TABLE]

for $m,n=1,\ldots,2N$ and $k=1,\ldots,N$.

II-A Modeling assumptions from Part I [2]

In this section, we recall the assumptions used in Part I [2] to establish the network mean-square error stability (14). We first introduce the Hermitian Hessian matrix functions (see [2, Sec. II-A]):

[TABLE]

Assumption 1**.**

(Conditions on aggregate and individual costs).*

The individual costs $J_{k}(w_{k})\in\mathbb{R}$ are assumed to be twice differentiable and convex such that:*

[TABLE]

where $\nu_{k}\geq 0$ for $k=1,\ldots,N$. It is further assumed that, for any ${\scriptstyle\mathcal{W}}$, $\mathcal{H}({\scriptstyle\mathcal{W}})$ satisfies:

[TABLE]

for some positive parameters $\nu\leq\delta$. The data-type variable $h$ and the matrix $\mathcal{U}^{e}$ are defined in Table I.

As explained in [2], condition (20) ensures that problem (1) has a unique minimizer ${\scriptstyle\mathcal{W}}^{o}$.

Assumption 2**.**

(Conditions on gradient noise).*

The gradient noise process defined in (4) satisfies for any $\boldsymbol{w}\in\boldsymbol{\cal{F}}_{i-1}$ and for all $k,\ell=1,\ldots,N$:*

[TABLE]

for some $\beta_{k}^{2}\geq 0$, $\sigma^{2}_{s,k}\geq 0$, and where $\boldsymbol{\cal{F}}_{i-1}$ denotes the filtration generated by the random processes $\{\boldsymbol{w}_{\ell,j}\}$ for all $\ell=1,\ldots,N$ and $j\leq i-1$.

Assumption 3**.**

(Condition on $\mathcal{U}$).*

The full-column rank matrix $\mathcal{U}$ is assumed to be semi-unitary, i.e., its column vectors are orthonormal and $\mathcal{U}^{*}\mathcal{U}=I_{P}$.*

Consider the $N\times N$ block matrix $\mathcal{A}^{e}$ whose $(k,\ell)$-th block is defined in Table I. This matrix will appear in our subsequent study. In [2, Lemma 2], we showed that this $hM\times hM$ matrix $\mathcal{A}^{e}$ admits a Jordan decomposition of the form:

[TABLE]

with

[TABLE]

where $\mathcal{U}^{e},\mathcal{J}^{e}_{\epsilon},\mathcal{V}^{e}_{R,\epsilon}$, and $(\mathcal{V}_{L,\epsilon}^{e})^{*}$ are defined in Table I with the matrices $\mathcal{J}_{\epsilon}$, $\mathcal{V}_{R,\epsilon}$, and $\mathcal{V}_{L,\epsilon}^{*}$ originating from the eigen-structure of $\mathcal{A}$. Under Assumption 3, the $M\times M$ combination matrix $\mathcal{A}$ satisfying conditions (8), (9), and (10) admits a Jordan canonical decomposition of the form:

[TABLE]

with:

[TABLE]

where $\mathcal{J}_{\epsilon}$ is a Jordan matrix with the eigenvalues (which may be complex but have magnitude less than one) on the diagonal and $\epsilon>0$ on the super-diagonal. The eigen-decomposition (25) will be useful for establishing the mean-square performance.

The results in Part I [2] established that the iterates $\boldsymbol{w}_{k,i}$ converge in the mean-square error sense to a small $O(\mu)-$neighborhood around the solution $w^{o}$. In this part, we will be more precise and determine the size of this neighborhood, i.e., assess the size of the constant multiplying $\mu$ in the $O(\mu)-$term. To do so, we shall derive an accurate first-order expression for the mean-square error (14); the expression will be accurate to first-order in $\mu$.

To arrive at the desired expression, we start by motivating a long-term model for the evolution of the network error vector after sufficient iterations have passed, i.e., for $i\gg 1$. It turns out that the performance expressions obtained from analyzing the long-term model provide accurate expressions for the performance of the original network model to first order in $\mu$. To derive the long-term model, we follow the approach developed in [4]. The first step is to establish the asymptotic stability of the fourth-order moment of the error vector, $\mathbb{E}\|w^{o}_{k}-\boldsymbol{w}_{k,i}\|^{4}$. Under the same settings of Theorem 1 in [2] with the second-order moment condition (24) replaced by the fourth-order moment condition:

[TABLE]

with $\beta_{4,k}^{4}\geq 0$, $\sigma^{4}_{s4,k}\geq 0$, and using similar arguments as in [4, Theorem 9.2], we can show that the fourth-order moments of the network error vectors are stable for sufficiently small $\mu$, namely, it holds that (see Appendix F)

[TABLE]

As explained in [4], condition (29) implies (24). We analyze the long-term model under the same settings of Theorem 1 in [2] and the following smoothness assumption on the individual costs.

Assumption 4**.**

(Smoothness condition on individual costs).*

It is assumed that each $J_{k}(w_{k})$ satisfies the following smoothness condition close to the limit point $w^{o}_{k}$:*

[TABLE]

for small perturbations $\|\Delta w_{k}\|$ and $\kappa_{d}\geq 0$.

II-B Long-term-error model

To introduce the long-term model, we reconsider the network error recursion from Part I [2], namely,

[TABLE]

where:

[TABLE]

where:

[TABLE]

and $\widetilde{\boldsymbol{w}}_{k,i}^{e}$, $\boldsymbol{s}_{k,i}^{e}(\boldsymbol{w}_{k,i-1})$, and $b_{k}^{e}$ are defined in Table I with:

[TABLE]

We rewrite (32) as:

[TABLE]

in terms of the constant matrix $\mathcal{B}$ and the random perturbation sequence $\boldsymbol{c}_{i-1}$:

[TABLE]

where $\mathcal{H}^{o}$ and $\widetilde{\boldsymbol{\cal{H}}}_{i-1}$ are given by:

[TABLE]

with each $H_{k}^{o}$ given by the value of the Hessian matrix at the limit point, namely,

[TABLE]

By exploiting the smoothness condition (31), and following an argument similar to [4, pp. 554], we can show from Theorem 1 in [2] that, for $i\gg 1$, $\|\boldsymbol{c}_{i-1}\|=O(\mu)$ with high probability. Motivated by this observation, we introduce the following approximate model, where the last term $\mu\mathcal{A}^{e}\boldsymbol{c}_{i-1}$ that appears in (41), which is $O(\mu^{2})$, is removed:

[TABLE]

for $i\gg 1$. Obviously, the iterates $\{\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}\}$ generated by (47) are generally different from the iterates generated by the original recursion (32). To highlight this fact, we are using the prime notation for the state of the long-term model. Note that the driving process $\boldsymbol{s}_{i}^{e}(\boldsymbol{{\scriptstyle\mathcal{W}}}_{i-1})$ in (47) is the same gradient noise process from the original recursion (32).

We start by showing that the mean-square difference between $\{\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}_{i}^{e^{\prime}},\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}_{i}^{e}\}$ is asymptotically bounded by $O(\mu^{2})$ and that the mean-square-error of the long term model (47) is within $O(\mu^{\frac{3}{2}})$ from the one of the original recursion (32). Working with (47) is much more tractable for performance analysis because its dynamics is driven by the constant matrix $\mathcal{B}$ as opposed to the random matrix $\boldsymbol{\cal{B}}_{i-1}$ in (32). Therefore, we shall work with model (47) and evaluate its performance, which will provide an accurate representation for the performance of (2) to first order in $\mu$.

Theorem 1**.**

(Size of approximation error).*

Consider a network of $N$ agents running the distributed strategy (2) with a matrix $\mathcal{A}$ satisfying conditions (8), (9), and (10) and $\mathcal{U}$ satisfying Assumption 3. Assume the individual costs, $J_{k}(w_{k})$, satisfy the conditions in Assumptions 1 and 4. Assume further that the gradient noise processes satisfy the conditions in Assumption 2 with the second-order moment condition (24) replaced by the fourth-order moment condition (29). Then, it holds that, for sufficiently small step-sizes:*

[TABLE]

Proof.

See Appendix A. ∎

Using similar eigenvalue perturbation arguments as in [4, Theorem 9.3], we can show that, under the same settings of Theorem 1, the constant matrix $\mathcal{B}$ defined by (42) is stable for sufficiently small step-sizes (see Appendix G).

II-C Mean-square-error performance

We showed in Theorem 1 in Part I [2] that a network running the distributed strategy (2) is mean-square stable for sufficiently small $\mu$. Particularly, we showed that $\limsup_{i\rightarrow\infty}\mathbb{E}\|w^{o}_{k}-\boldsymbol{w}_{k,i}\|^{2}=O(\mu)$. In this section, we assess the size of the mean-square error by measuring the network MSD defined by (12).

We refer to the individual gradient noise process in (4) and denote its conditional covariance matrix by:

[TABLE]

We assume that, in the limit, the following moment matrices tend to constant values when evaluated at $w^{o}_{k}$:

[TABLE]

Assumption 5**.**

(Smoothness condition on noise covariance).*

It is assumed that the conditional second-order moments of the individual noise processes satisfy the following smoothness condition,*

[TABLE]

for small perturbations $\|\Delta w_{k}\|$, and for some constant $\kappa_{d}\geq 0$ and exponent $0<\gamma\leq 4$.

One useful conclusion that follows from (53) is that, for $i\gg 1$ and for sufficiently small step-size, we can express the covariance matrix of $\boldsymbol{s}_{k,i}^{e}(\boldsymbol{w})$ in terms of the limiting matrix $R^{o}_{k}$ defined in Table I as follows (see [4, Lemma 11.1]):

[TABLE]

Before proceeding, we introduce the $(hM)^{2}\times(hM)^{2}$ matrix $\mathcal{F}$ that will play a critical role in characterizing the performance:

[TABLE]

This matrix is defined in in terms of the block Kronecker operation. In the derivation that follows, we shall use the block Kronecker product $\otimes_{b}$ operator [22] and the block vectorization operator $\text{bvec}(\cdot)$111In our derivations, the block Kronecker product and the block vectorization operations are applied to $2\times 2$ block matrices $\mathcal{C}=[C_{k\ell}]$ and $\mathcal{D}=[D_{k\ell}]$ with blocks $\{C_{11},D_{11}\}$ of size $hP\times hP$, blocks $\{C_{12},D_{12}\}$ of size $hP\times h(M-P)$, blocks $\{C_{21},D_{21}\}$ of size $h(M-P)\times hP$, and blocks $\{C_{22},D_{22}\}$ of size $h(M-P)\times h(M-P)$.. As explained in [4], these operations preserve the locality of the blocks in the original matrix arguments. The matrix $\mathcal{F}$ will sometimes appear transformed under the similarity transformation:

[TABLE]

Lemma 1**.**

(Low-rank approximation).* Assume the matrix $\mathcal{A}$ satisfies conditions (8), (9), and (10) with $\mathcal{U}$ satisfying Assumption 3. For sufficiently small step-sizes, it holds that $\mathcal{F}$ is stable and that:*

[TABLE]

where the leading $(hP)^{2}\times(hP)^{2}$ block in $(I-\overline{\mathcal{F}})^{-1}$ is $O(1/\mu)$. Moreover, we can also write:

[TABLE]

in terms of the block Kronecker operation, where the matrix $Z$ has dimension $(hP)^{2}\times(hP)^{2}$:

[TABLE]

with $\mathcal{D}_{11}=\mu\,(\mathcal{U}^{e})^{*}\mathcal{H}^{o}\mathcal{U}^{e}$ which is positive definite under Assumption 1.

Proof.

See Appendix B. ∎

Theorem 2**.**

(Mean-square-error performance).*

Consider the same settings of Theorem 1. Assume further that Assumption 5 holds. Let $\gamma_{m}\triangleq\frac{1}{2}\min\{1,\gamma\}>0$ with $\gamma\in(0,4]$ from (53). Then, it holds that:*

[TABLE]

where:

[TABLE]

Furthermore, it holds that:

[TABLE]

Proof.

See Appendix C. ∎

Since $(I-\mathcal{F})$ is of size $(hM)^{2}\times(hM)^{2}$, the first term on the R.H.S. of expression (63) may be hard to evaluate due to numerical reasons. In comparison, the first term in expression (64) only requires manipulations of matrices of size $hM\times hM$. In practice, a reasonable number of terms can be used instead of $n\rightarrow\infty$ to obtain accurate evaluation.

Note that the MSD of the centralized solution is equal to (67) since the centralized implementation can be obtained from (2) by replacing $\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}$ by $\mathcal{A}$ and by assuming fully-connected network. We therefore conclude, for sufficiently small step-sizes (i.e., in the slow adaptation regime), that the distributed strategy (2) is able to attain the same MSD performance as the centralized solution.

III Finding a combination matrix $\mathcal{A}$

In the following, we consider the problem of finding an $\mathcal{A}$ that minimizes the number of edges to be added to the original topology while satisfying the constraints (10), (8), and (9). That is, we consider the following optimization problem:

[TABLE]

where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A_{k\ell}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{1}\triangleq\sum_{m=1}^{M_{k}}\sum_{n=1}^{M_{\ell}}|[A_{k\ell}]_{mn}|\in\mathbb{R}$, $\|\mathcal{A}\|_{\text{F}}=\sqrt{\text{Tr}(\mathcal{A}^{*}\mathcal{A})}\in\mathbb{R}$ is the Frobenius norm of $\mathcal{A}$, $\gamma\geq 0$ is a regularization parameter, and $\epsilon\in(0,1]$ is a small positive number. In general, the spectral radius of a matrix is not convex over the matrix space. We therefore restrict our search to the class of Hermitian matrices, since their spectral radius coincides with their spectral norm (maximum singular value), which is a convex function. Problem (68) is convex since the objective is convex, the equality constraints are linear, and the inequality constraint function is convex [23]. The parameter $\epsilon$ controls the convergence rate of $\mathcal{A}^{i}$ towards the projector $\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}$. That is, small $\epsilon$ leads to slow convergence and large $\epsilon$ gives fast convergence. The convex $\ell_{1}$-norm based function $\sum_{k=1}^{N}\sum_{\ell\notin\mathcal{N}_{k}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A_{k\ell}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{1}$ is used as a relaxation of the pseudo $\ell_{0}$-norm $h(\mathcal{A})=\sum_{k=1}^{N}\text{card}\{\ell|A_{k\ell}\neq 0,\ell\notin\mathcal{N}_{k}\}$, which is a non-convex function that leads to computational challenges. Among the potentially multiple feasible solutions, the cardinality function $h(\mathcal{A})$ in the objective in (68) selects as optimum the one that minimizes the number of edges to be added to the network topology in order to satisfy constraint (5). The quadratic term $\|\mathcal{A}\|^{2}_{\text{F}}=\sum_{m=1}^{M}\sum_{n=1}^{M}|a_{mn}|^{2}$ in (68) makes the objective strictly convex, and therefore problem (68) has a unique minimum. Problem (68) can be solved using general convex optimization solvers such as CVX [24]. These solvers generally implement second-order methods that require calculation of Hessian matrices. Therefore, problems with more than few thousand entries are probably beyond the capabilities of these solvers. The Douglas-Rachford algorithm can also be employed to solve problem (68). As we shall see in the following, the required proximal operators for implementing this algorithm can be computed efficiently using closed form expressions.

In the following, we shall assume that $\mathcal{U}\in\mathbb{R}^{M\times P}$ and, therefore, we shall solve (68) over real-valued matrices $\mathcal{A}\in\mathbb{R}^{M\times M}$. In order to solve the constrained problem (68), we shall apply the Douglas-Rachford splitting algorithm [25], which is used to solve problems of the form:

[TABLE]

where $g_{1}(\cdot)$ and $g_{2}(\cdot)$ are functions in $\Gamma_{0}(\mathbb{R}^{N})$ such that $(\text{ri dom}g_{1})\cap(\text{ri dom}g_{2})\neq 0$ and $g_{1}(x)+g_{2}(x)\rightarrow+\infty$ as $\|x\|\rightarrow+\infty$. By selecting $g_{1}(\cdot)$ as $f(\cdot)$ in (68) and $g_{2}(\cdot)$ as the indicator function $\mathcal{I}_{\Omega}(\cdot)$ of the closed nonempty convex set:

[TABLE]

defined as:

[TABLE]

the Douglas-Rachford algorithm to solve (68) has the following form:

[TABLE]

where $\eta>0$ and $\text{prox}_{\eta g}:\mathbb{R}^{M\times M}\rightarrow\mathbb{R}^{M\times M}$ is the proximal operator of $\eta g(\cdot)$ ($g:\mathbb{R}^{M\times M}\rightarrow\mathbb{R}\cup\{+\infty\}$) defined as [25, 26]:

[TABLE]

Every sequence $(\mathcal{A}_{i})_{i\in\mathbb{N}}$ generated by algorithm (72) converges to a solution of problem (68) [25, Proposition 4.3]. The Douglas-Rachford algorithm operates by splitting since it employs the functions $f(\cdot)$ and $\mathcal{I}_{\Omega}(\cdot)$ separately. It requires the implementation of two proximal steps at each iteration, which can be computed efficiently as explained in the following.

The function $f(\mathcal{A})$ is an entrywise matrix function that treats the matrix $\mathcal{A}\in\mathbb{R}^{M\times M}$ as a vector in $\mathbb{R}^{M^{2}}$ and then uses a corresponding vector function; the proximal operator is then the same as that of the vector function. Let $C_{k\ell}$ denote the $(k,\ell)$-th block of an $N\times N$ block matrix $\mathcal{C}$ and let $[C_{k\ell}]_{mn}$ denote the $(m,n)$-th entry of $C_{k\ell}$. The $(k,\ell)$-th block of the proximal operator of $\eta f(\cdot)$ is given by:

[TABLE]

where the matrix $C^{s}_{k\ell}$ is of size $M_{k}\times M_{\ell}$ with $(m,n)$-th entry given by:

[TABLE]

Since $\mathcal{I}_{\Omega}$ is the indicator function of the closed convex set $\Omega$, its proximal operator reduces to the projection onto $\Omega$ defined as:

[TABLE]

where the parameter $\eta$ does not appear since the proximal operator is a projection. The set $\Omega$ in (70) can be written alternatively as $\Omega=\Omega_{1}\cap\Omega_{2}$ where $\Omega_{1}$ and $\Omega_{2}$ are two closed convex sets defined as:

[TABLE]

As we shall explain in the following, the projection onto the intersection $\Omega$ can be obtained by properly projecting onto the individual sets $\Omega_{1}$ and $\Omega_{2}$ according to:

[TABLE]

The projection onto $\Omega_{1}$ is given by (see Appendix D):

[TABLE]

and the projection of the symmetric matrix $\Pi_{\Omega_{1}}(\mathcal{D})$ onto $\Omega_{2}$ is given by (see Appendix D):

[TABLE]

where:

[TABLE]

where $\{\lambda_{m},v_{m}\}_{m=1}^{M}$ are the eigenvalues and eigenvectors of the symmetric matrix $(I-\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}})\left(\frac{\mathcal{D}+\mathcal{D}^{\top}}{2}\right)(I-\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}})$. In order to establish (79), we introduce the following lemma.

Lemma 2**.**

(Characterization of the projection).*

If $\Omega_{1}$ is an affine set, $\Omega_{2}$ is a closed convex set, and $\Pi_{\Omega_{2}}(\Pi_{\Omega_{1}}(\mathcal{C}))\in\Omega_{1}$, then $\Pi_{{\Omega_{1}}\cap{\Omega_{2}}}(\mathcal{C})=\Pi_{\Omega_{2}}(\Pi_{\Omega_{1}}(\mathcal{C}))$.*

Proof.

See Appendix E. ∎

Since the projection onto $\Omega_{2}$ (given by (81)) changes only the eigenvalues of a matrix without affecting the eigenvectors, we have $\Pi_{\Omega_{2}}(\Pi_{\Omega_{1}}(\mathcal{C}))\in\Omega_{1}$. We then conclude from Lemma 2 that (79) holds.

IV Simulation results

We apply strategy (2) to solve distributed inference under smoothness (described in Remark 4 of Section II in [2]). We consider a connected mean-square-error (MSE) network of $N=50$ nodes and $M_{k}=L=5$. The $N$ nodes are placed randomly in the $[0,1]\times[0,1]$ square, and the weighted graph is then constructed according to a thresholded Gaussian kernel weighting function based on the distance between nodes. Particularly, the weight $c_{k\ell}$ of edge $(k,\ell)$ connecting nodes $k$ and $\ell$ that are $d_{k\ell}$ apart is:

[TABLE]

with $\sigma=0.12$ and $\kappa=0.33$. We assume real data case. Each agent is subjected to streaming data $\{\boldsymbol{d}_{k}(i),\boldsymbol{u}_{k,i}\}$ assumed to satisfy a linear regression model [4]:

[TABLE]

for some unknown $L\times 1$ vector $w^{\star}_{k}$ to be estimated with $\boldsymbol{v}_{k}(i)$ denoting a zero-mean measurement noise. For these networks, the risk functions take the form of mean-square-error costs:

[TABLE]

The processes $\{\boldsymbol{u}_{k,i},\boldsymbol{v}_{k}(i)\}$ are assumed to be zero-mean Gaussian with: i) $\mathbb{E}\boldsymbol{u}_{k,i}\boldsymbol{u}_{\ell,i}^{\top}=R_{u,k}=\sigma^{2}_{u,k}I_{L}$ if $k=\ell$ and zero otherwise; ii) $\mathbb{E}\boldsymbol{v}_{k}(i)\boldsymbol{v}_{\ell}(i)=\sigma^{2}_{v,k}$ if $k=\ell$ and zero otherwise; and iii) $\boldsymbol{u}_{k,i}$ and $\boldsymbol{v}_{k}(i)$ are independent of each other. The variances $\sigma^{2}_{u,k}$ and $\sigma^{2}_{v,k}$ are generated from the uniform distributions $\text{unif}(0.5,2)$ and $\text{unif}(0.2,0.8)$, respectively. Let ${\scriptstyle\mathcal{W}}^{\star}=\text{col}\{w_{1}^{\star},\ldots,w^{\star}_{N}\}$. The signal ${\scriptstyle\mathcal{W}}^{\star}$ is generated by smoothing a signal ${\scriptstyle\mathcal{W}}_{o}$ by a diffusion kernel. Particularly, we generate ${\scriptstyle\mathcal{W}}^{\star}$ according to ${\scriptstyle\mathcal{W}}^{\star}=[(Ve^{-\tau\Lambda}V^{\top})\otimes I_{L}]{\scriptstyle\mathcal{W}}_{o}$ with $\tau=30$, ${\scriptstyle\mathcal{W}}_{o}$ a randomly generated vector from the Gaussian distribution $\mathcal{N}(0.1\times\mathds{1}_{NL},I_{NL})$, and $\{V=[v_{1},\ldots,v_{N}],\Lambda=\text{diag}\{\lambda_{1},\ldots,\lambda_{N}\}\}$ are the matrices of eigenvectors and eigenvalues of $L_{c}=\text{diag}\{C\mathds{1}_{N}\}-C$ with $[C]_{k\ell}=c_{k\ell}$ given by (83). Figure 1 (right) illustrates the normalized squared $\ell_{2}$-norm of the spectral component $\overline{w}^{\star}_{m}=(v_{m}^{\top}\otimes I_{L}){\scriptstyle\mathcal{W}}^{\star}$. It can be observed that the signal is mainly localized in $[0,0.1]$. Note that, for MSE networks, it holds that $H_{k}(w_{k})=R_{u,k}$ $\forall w_{k}$. Furthermore, the gradient noise process (4) is given by:

[TABLE]

with covariance $R^{o}_{k}$ given by:

[TABLE]

where $W_{k}=(w^{\star}_{k}-w^{o}_{k})(w^{\star}_{k}-w^{o}_{k})^{\top}$, and where we used the fact that $\mathbb{E}[\boldsymbol{u}_{k,i}^{\top}\boldsymbol{u}_{k,i}W_{k}\boldsymbol{u}_{k,i}^{\top}\boldsymbol{u}_{k,i}]=2R_{u,k}W_{k}R_{u,k}+R_{u,k}\text{Tr}(R_{u,k}W_{k})$ since the regressors are zero-mean real Gaussian [27].

We run algorithm (2) for 5 different choices of matrix $\mathcal{U}$ in (1) with $\mathcal{U}=U\otimes I_{L}$: i) matrix $U$ chosen as the first eigenvector of the Laplacian $U=[v_{1}]=\frac{1}{\sqrt{N}}\mathds{1}_{N}$; ii) matrix $U$ chosen as the first two eigenvectors of the Laplacian $U=[v_{1}~{}v_{2}]$; iii) $U=[v_{1}~{}v_{2}~{}v_{3}]$; iv) $U=[v_{1}~{}\ldots~{}v_{4}]$; v) $U=[v_{1}~{}\ldots~{}v_{5}]$. Since $\mathcal{U}=U\otimes I_{L}$, the matrix $\mathcal{A}$ is of the form $A\otimes I_{L}$ with $A=[a_{k\ell}]$ an $N\times N$ matrix. In each case, the combination matrix $A$ is set as the solution of the optimization problem (68) ($\epsilon=0.01,\gamma=0$, ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A_{k\ell}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{1}=|a_{k\ell}|$), which is solved by the Douglas-Rachford algorithm (72) with $\eta=0.003$. Note that, for the $5$ different choices of $U$, the distributed implementation is feasible and the steady-state value of the cost in (68) is zero. We set $\mu=0.001$. We report the network $\text{MSD}^{\star}$ learning curves $\frac{1}{N}\mathbb{E}\|{\scriptstyle\mathcal{W}}^{\star}-\boldsymbol{{\scriptstyle\mathcal{W}}}_{i}\|^{2}$ in Fig. 2 (left). The results are averaged over $200$ Monte-Carlo runs. The learning curve of the non-cooperative solution, obtained from (2) by setting $\mathcal{A}=I_{LN}$, is also reported. The results show that the best performance is obtained when $\mathcal{U}=[v_{1}~{}v_{2}~{}v_{3}~{}v_{4}]\otimes I_{L}$. This is due to the fact that the columns of $\mathcal{U}$ constitute a basis spanning the useful signal subspace (see Fig. 1 (right)). As a consequence, a strong noise reduction may be obtained by projecting onto this subspace compared with the non-cooperative strategy where each agent estimates $w^{\star}_{k}$ without any cooperation. By forcing consensus (i.e., by choosing $U=[v_{1}]$), the resulting estimate $\boldsymbol{w}_{k,i}$ will be biased with respect to $w^{\star}_{k}$, which is not common across agents. The performance obtained when $U=[v_{1}~{}\ldots~{}v_{5}]$ is worse than the case where $U=[v_{1}~{}\ldots~{}v_{4}]$ due to a smaller noise reduction.

Finally, we illustrate Theorem 2 in Table II by reporting the steady-state $\text{MSD}=\limsup_{i\rightarrow\infty}\frac{1}{N}\mathbb{E}\|{\scriptstyle\mathcal{W}}^{o}-\boldsymbol{{\scriptstyle\mathcal{W}}}_{i}\|^{2}$ when $\mathcal{U}=[v_{1}~{}\ldots~{}v_{4}]\otimes I_{L}$ for $3$ different values of the step-size $\mu=\{10^{-2},10^{-3},10^{-4}\}$. A closed form solution for ${\scriptstyle\mathcal{W}}^{o}$ in (1) exists and is given by:

[TABLE]

where $\mathcal{H}=\text{diag}\{R_{u,k}\}_{k=1}^{N}$. We observe that, in the small adaptation regime, i.e., when $\mu\rightarrow 0$, the network MSD increases approximately $10$dB per decade (when $\mu$ goes from $\mu_{1}$ to $10\mu_{1}$). This means that the steady-state MSD is on the order of $\mu$. We also observe that, in the small adaptation regime, the distributed solution is able to attain the same performance as the centralized one. Finally, note that, for relatively large step-size ($\mu=10^{-2}$), expression (64) provides better results than (67) in the distributed case. This is due to neglecting the $O(1)$ term in (60) which is multiplied by $O(\mu^{2})$ (since $\mathcal{Y}=O(\mu^{2})$) when replaced in (64).

V Conclusion

In this paper, and its accompanying Part I [2], we considered inference problems over networks where agents have individual parameter vectors to estimate subject to subspace constraints that require the parameters across the network to lie in low-dimensional subspaces. Based on the gradient projection algorithm, we proposed an iterative and distributed implementation of the projection step, which runs in parallel with the stochastic gradient descent update. We showed that, for small step-size parameter, the network is able to approach the minimizer of the constrained problem to arbitrarily good accuracy levels. Furthermore, we derived a closed-form expressions for the steady-state mean-square-error (MSE) performance. These expressions revealed explicitly the influence of the gradient noise, data characteristics, and subspace constraints, on the network performance. Finally, among many possible convex formulations, we considered the design of feasible distributed solution that minimizes the number of edges to be added to the original graph.

Appendix A Proof of Theorem 1

To simplify the notation, we introduce the difference vector $\boldsymbol{z}_{i}\triangleq\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}_{i}^{e}-\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}.$ Using (44) in the expression for $\boldsymbol{\cal{B}}_{i-1}$ in (35), we can write:

[TABLE]

in terms of the constant coefficient matrix $\mathcal{B}$ in (42). Using (89) and (43), and subtracting (32) and (47), we then get:

[TABLE]

If we multiply both sides of (90) from the left by $(\mathcal{V}^{e}_{\epsilon})^{-1}$ we obtain:

[TABLE]

where we partitioned the vectors $(\mathcal{V}^{e}_{\epsilon})^{-1}\boldsymbol{z}_{i}$ and $\mu\Lambda_{\epsilon}^{e}(\mathcal{V}^{e}_{\epsilon})^{-1}\boldsymbol{c}_{i-1}$ into:

[TABLE]

with the leading vectors, $\{\overline{\boldsymbol{z}}_{i},\overline{\boldsymbol{c}}_{i-1}\}$, having dimensions $hP\times 1$ each. The matrix $\overline{\mathcal{B}}$ is given by:

[TABLE]

with the blocks $\{\mathcal{D}_{mn}\}$ given by:

[TABLE]

Recursion (91) has a form similar to the earlier recursion (66) in Part I [2] with three differences. First, the matrices $\{\mathcal{D}_{mn}\}$ in (92) are constant matrices; nevertheless, they satisfy the same bounds as the matrices $\{\boldsymbol{\cal{D}}_{mn,i-1}\}$ in eq. (66) in Part I [2]. In particular, from (115), (116), and (122) in Part I [2], it continues to hold that:

[TABLE]

for some positive constants $\sigma_{11},\sigma_{12},\sigma_{21},\sigma_{22}$ that are independent of $\mu$. Second, Third, the bias term $\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{b}]{\kern-0.3pt\bigwedge\kern-0.3pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{b}{\scalebox{-0.8}{\tmpbox}}^{e}$ in (66) in Part I [2] is absent from (92). Third, the gradient noise terms that appeared in recursion (66) in Part I [2] are now replaced by the perturbation sequences $\{\overline{\boldsymbol{c}}_{i-1},\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\boldsymbol{c}}]{\kern-0.3pt\bigwedge\kern-0.3pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\boldsymbol{c}}{\scalebox{-0.8}{\tmpbox}}_{i-1}\}$. However, these sequences can be bounded as follows:

[TABLE]

for some constant $r$ that is independent of $\mu$ since:

[TABLE]

To establish the above inequality, we start by noting that any cost $J_{k}(\cdot)$ satisfying (19) and (31) will also satisfy [4, Lemmas E.4, E.8]:

[TABLE]

for any $\Delta w_{k}$ and where \kappa^{\prime}_{d}\triangleq\max\{\kappa_{d},(\delta_{k}-\nu_{k})/(h{\color[rgb]{0,0,0}\epsilon})\}. Then, for each agent $k$ we have:

[TABLE]

Therefore,

[TABLE]

Now, replacing $\boldsymbol{c}_{i-1}$ in (101) by (43) and using (104) we conclude (101).

Repeating the argument that led to inequalities (129) and (130) in Part I [2] we obtain:

[TABLE]

and

[TABLE]

We can combine (105) and (106) into a single inequality recursion as follows:

[TABLE]

where $a=1-O(\mu)$, $b=O(\mu)$, $c=O(\mu^{2})$, $d=\rho(\mathcal{J}_{\epsilon})+\epsilon+O(\mu^{2})$, $e=O(\mu)$, and $f=O(\mu^{2})$. Using (30) and eq. (134) in Part I [2] we conclude that:

[TABLE]

and, hence, $\limsup_{i\rightarrow\infty}\mathbb{E}\|\boldsymbol{z}_{i}\|^{2}=O(\mu^{2}).$ It follows that $\limsup_{i\rightarrow\infty}\mathbb{E}\|\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}_{i}^{e}-\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}\|^{2}=O(\mu^{2}),$ which establishes (48). Finally, note that:

[TABLE]

where we used $|\mathbb{E}\boldsymbol{x}|\leq\mathbb{E}|\boldsymbol{x}|$ from Jensen’s inequality and where we applied Holder’s inequality:

[TABLE]

Hence, from (14) and (48) we get:

[TABLE]

since $\mu^{2}<\mu^{3/2}$ for small $\mu\ll 1$, which establishes (49).

From (14) and (111), it follows that:

[TABLE]

and, therefore, the long-term approximate model (47) is also mean-square stable.

Appendix B Low-rank approximation

From (95), we obtain:

[TABLE]

[TABLE]

where the block matrices $\{\mathcal{D}^{\top}_{mn},\mathcal{D}^{*}_{mn}\}$ are all on the order of $\mu$ with:

[TABLE]

of dimensions $hP\times hP$. Substituting (113) and (114) into (55) and using property $(A\otimes_{b}B)(C\otimes_{b}D)=AC\otimes_{b}BD,$ for block Kronecker products, we obtain:

[TABLE]

where we introduced:

[TABLE]

We partition $\mathcal{X}$ into the following block structure:

[TABLE]

where, for example, $\mathcal{X}_{11}$ is $(hP)^{2}\times(hP)^{2}$ and is given by:

[TABLE]

Since

[TABLE]

we proceed to evaluate $I-\mathcal{F}$. It follows that:

[TABLE]

and, in a manner similar to the way we assessed the size of the block matrices $\{\boldsymbol{\cal{D}}_{mn,i-1}\}$ in the proof of Theorem 1 in Part I [2], we can verify that:

[TABLE]

The matrix $I-\mathcal{X}$ is invertible since $I-\mathcal{F}$ is invertible; this is because $\rho(\mathcal{F})=[\rho(\mathcal{B})]^{2}<1$. Therefore, applying the block matrix inversion formula to $I-\mathcal{X}$ we get:

[TABLE]

where $\mathcal{Y}_{11}=(I-\mathcal{X}_{11})^{-1}\mathcal{X}_{12}\Delta^{-1}\mathcal{X}_{21}(I-\mathcal{X}_{11})^{-1}$, $\mathcal{Y}_{12}=(I-\mathcal{X}_{11})^{-1}\mathcal{X}_{12}\Delta^{-1}$, $\mathcal{Y}_{21}=\Delta^{-1}\mathcal{X}_{21}(I-\mathcal{X}_{11})^{-1}$, and $\Delta=(I-\mathcal{X}_{22})-\mathcal{X}_{21}(I-\mathcal{X}_{11})^{-1}\mathcal{X}_{12}$. The entries of $(I_{(hP)^{2}}-\mathcal{X}_{11})^{-1}$ are $O(1/\mu)$, while the entries in the second matrix on the right-hand side of the above equation are $O(1)$ when the step-size is small. That is, we can write:

[TABLE]

Moreover, since $O(1/\mu)$ dominates $O(1)$ for sufficiently small $\mu$, we can also write:

[TABLE]

Substituting (128) into (121) we arrive at (61). Since $Z=O(\mu)$, we conclude that (57) holds. We also conclude that (60) holds since:

[TABLE]

Appendix C Proof of Theorem 2

Consider the long-term model (47). Conditioning both sides of (47) on $\boldsymbol{\cal{F}}_{i-1}$, invoking the conditions on the gradient noise process from Assumption 2, and computing the conditional expectation, we obtain:

[TABLE]

where the term involving $\boldsymbol{s}_{i}^{e}(\boldsymbol{{\scriptstyle\mathcal{W}}}^{e}_{i-1})$ is eliminated because $\mathbb{E}[\boldsymbol{s}_{i}^{e}|\boldsymbol{\cal{F}}_{i-1}]=0$. Taking expectations again we arrive at:

[TABLE]

Since recursion (47) includes a constant driving term $\mu\mathcal{A}^{e}b^{e}$, we introduce the centered variable $\boldsymbol{y}_{i}\triangleq\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}-\mathbb{E}\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}.$ Subtracting (138) from (47), we find that $\boldsymbol{y}_{i}$ satisfies the following recursion:

[TABLE]

Although we are interested in evaluating $\limsup_{i\rightarrow\infty}\mathbb{E}\|\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}\|^{2}$, we can still rely on $\boldsymbol{y}_{i}$ since it holds for $i\gg 1$:

[TABLE]

where we used the fact that $\limsup_{i\rightarrow\infty}\|\mathbb{E}\widetilde{\boldsymbol{{\scriptstyle\mathcal{W}}}}^{e^{\prime}}_{i}\|=O(\mu)$ (see Appendix H). Therefore, from (49) and (140), we obtain:

[TABLE]

Let $\Sigma$ denote an arbitrary Hermitian positive semi-definite matrix that we are free to choose. Equating the squared weighted values of both sides of (139) and taking expectations conditioned on the past history gives:

[TABLE]

Taking expectations again, we get:

[TABLE]

From (54) and using same arguments as in [4, pp. 586], we can rewrite the second term on the R.H.S. of (143) as:

[TABLE]

for $i\gg 1$. Therefore, we obtain:

[TABLE]

In order to reduce the weighting matrix on the mean-square value of $\boldsymbol{z}_{i}$ in (145) to the identity, we need to select $\Sigma$ as the solution to the following Lyapunov equation:

[TABLE]

This equation has a unique Hermitian non-negative definite solution $\Sigma$ [4, pp. 772] since the matrix $\mathcal{B}$ is stable for sufficiently small step-size. Now, by applying the block vectorization operation to both sides of (146) and by using the property that:

[TABLE]

we find that:

[TABLE]

in terms of the matrix $\mathcal{F}$ defined in (55).

Now, substituting $\Sigma$ in (148) into (145), we obtain $\mathbb{E}\|\boldsymbol{y}_{i}\|^{2}$ on the left-hand side while the term $\text{Tr}(\Sigma\mathcal{Y})$ on the right-hand side becomes:

[TABLE]

where we used the property that:

[TABLE]

Using the fact that $(I-\mathcal{F})^{-1}=O(1/\mu)$ and following similar arguments as in [4, pp. 590], we can show that:

[TABLE]

Replacing (149) and (151) into (145) gives (63). Observe that the first term on the R.H.S. of (63) is $O(\mu)$ since $\|\mathcal{Y}\|=O(\mu^{2})$ and $\|(I-\mathcal{F})^{-1}\|=O(1/\mu)$. Therefore, this term dominates the factor $O(\mu^{1+\gamma_{m}})$.

Since $\mathcal{F}$ is a stable matrix for sufficiently small step-sizes, we can employ the expansion $(I-\mathcal{F})^{-1}=I+\mathcal{F}+\mathcal{F}^{2}+\mathcal{F}^{3}+\ldots,$ replace $\mathcal{F}$ by (55), and use properties (150) and (147) to write the first term on the R.H.S. of (63) as $\sum_{n=0}^{\infty}\text{Tr}(\mathcal{B}^{n}\mathcal{Y}(\mathcal{B}^{*})^{n})$.

Now, in order to establish (67), we shall use the low-rank approximation (61). Using definition (12) and (63), we obtain:

[TABLE]

From (61) we get:

[TABLE]

Using property (147), it is straightforward to verify that the last three terms combine into the following result:

[TABLE]

Let us therefore evaluate the matrix vector product $x\triangleq Z^{-1}\text{vec}(I_{hP}).$ Using the definition (62) for $Z$, the vector $x$ is therefore the unique solution to the linear system of equations:

[TABLE]

Let $X=\text{unvec}(x)$ denote the $hP\times hP$ matrix whose vector representation is $x$. Applying to each of the terms appearing on the left-hand side of the above expression the Kronecker product property (147), albeit using vec instead of bvec operation, we find that $(I_{hP}\otimes\mathcal{D}_{11}^{*})x=\text{vec}(\mathcal{D}_{11}^{*}X),$ and $(\mathcal{D}_{11}^{\top}\otimes I_{hP})x=\text{vec}(X\mathcal{D}_{11}).$ We conclude from these equalities and from (155) that $X$ is the unique solution to the (continuous-time) Lyapunov equation $\mathcal{D}_{11}^{*}X+X\mathcal{D}_{11}=I_{hP}.$ Since $\mathcal{D}_{11}$ in (96) is Hermitian, we obtain:

[TABLE]

Therefore, substituting into (153) gives:

[TABLE]

where we used the fact that $(\mathcal{U}^{e})^{*}\mathcal{A}^{e}=(\mathcal{U}^{e})^{*}$. Now substituting the above expression into the right-hand side of (152) and computing the limit as $\mu\rightarrow 0$, we arrive at expression (67).

Appendix D Projection onto $\Omega_{1}$ in (77) and $\Omega_{2}$ in (78)

The closed convex set $\Omega_{1}$ in (77) can be rewritten alternatively as:

[TABLE]

and the projection onto it is given by:

[TABLE]

The Lagrangian of the convex optimization problem in (159) is defined as:

[TABLE]

where $X\in\mathbb{R}^{M\times P}$, $Y\in\mathbb{R}^{P\times M}$, and $Z\in\mathbb{R}^{M\times M}$ are the matrices of Lagrange multipliers. From the Karush-Kuhn-Tucker (KKT) conditions, we obtain at the optimum $(\mathcal{A}^{o};X^{o},Y^{o},Z^{o})$:

[TABLE]

From (164), we obtain:

[TABLE]

Multiplying both sides of (165) by $\mathcal{U}$ and using the fact that $\mathcal{U}^{\top}\mathcal{U}=I$ from Assumption 3, we obtain:

[TABLE]

Combining the previous expression with (161), we get:

[TABLE]

Replacing (167) into (165) and using the fact that $\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}=\mathcal{U}\mathcal{U}^{\top}$, we arrive at:

[TABLE]

Pre-multiplying both sides of the previous equation by $\mathcal{U}^{\top}$ and using the fact that $\mathcal{U}^{\top}\mathcal{U}=I$, we obtain:

[TABLE]

Combining the previous expression with (162), we arrive at:

[TABLE]

Pre-multiplying both sides of the previous equation by $\mathcal{U}$ and using the fact that $\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}=\mathcal{U}\mathcal{U}^{\top}$, we obtain:

[TABLE]

Replacing (171) into (168), we arrive at:

[TABLE]

and thus,

[TABLE]

Combining (163) and the previous two equations, we obtain:

[TABLE]

Replacing the previous equation into (172), we arrive at:

[TABLE]

Now, projecting a symmetric matrix $\mathcal{C}$ onto $\Omega_{2}$ in (78) is given by:

[TABLE]

where $\Omega_{3}\triangleq\{\mathcal{A}|\|\mathcal{A}\|\leq 1-\epsilon\}$. In order to project the symmetric matrix $\mathcal{C}-\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}$ onto $\Omega_{3}$, we need to compute its eigenvalue decomposition $\mathcal{C}-\mathcal{P}_{{\scriptscriptstyle\mathcal{U}}}=\sum_{m=1}^{M}\lambda_{m}v_{m}v_{m}^{\top}$ and then threshold the eigenvalues $\{\lambda_{m}\}_{m=1}^{M}$ to have absolute magnitude at most $1-\epsilon$ [26, pp. 191–194]. Thus we obtain:

[TABLE]

where:

[TABLE]

Now, replacing the matrix $\mathcal{C}$ by (175), we obtain (81).

Appendix E Proof of Lemma 2

In order to establish Lemma 2, we first need to introduce Lemmas 3 and 4.

Lemma 3**.**

Let $\Omega$ denote a closed convex set. For any $\mathcal{C}\notin\Omega$, $\mathcal{A}^{o}=\Pi_{\Omega}(\mathcal{C})$ if and only if $\langle\mathcal{C}-\mathcal{A}^{o},\mathcal{A}-\mathcal{A}^{o}\rangle\leq 0$, $\forall\mathcal{A}\in\Omega$ where $\langle X,Y\rangle=\emph{\text{Tr}}(X^{\top}Y)$.

Proof.

($\Rightarrow$) Let $\mathcal{A}^{o}=\Pi_{\Omega}(\mathcal{C})$ for any given $\mathcal{C}\notin\Omega$, that is, suppose that $\mathcal{A}^{o}$ is the unique solution to the optimization problem. Let $\mathcal{A}\in\Omega$ be such that $\mathcal{A}\neq\mathcal{A}^{o}$. Let $\alpha\in(0,1)$. Since $\Omega$ is convex, $(1-\alpha)\mathcal{A}^{o}+\alpha\mathcal{A}=\mathcal{A}^{o}+\alpha(\mathcal{A}-\mathcal{A}^{o})\in\Omega$. By the assumed optimality of $\mathcal{A}^{o}$, we must have:

[TABLE]

and we obtain:

[TABLE]

Now, note that (180) holds for any $\alpha\in(0,1)$. Since the RHS of (180) can be made arbitrarily small for a given $\mathcal{A}$, the LHS can not be strictly positive. Thus, we conclude as desired:

[TABLE]

($\Leftarrow$) Let $\mathcal{A}^{o}\in\Omega$ be such that $\langle\mathcal{C}-\mathcal{A}^{o},\mathcal{A}-\mathcal{A}^{o}\rangle\leq 0,\forall\mathcal{A}\in\Omega$. We shall show that it must be the optimal solution. Let $\mathcal{A}\in\Omega$ and $\mathcal{A}\neq\mathcal{A}^{o}$. We have:

[TABLE]

Hence, $\mathcal{A}^{o}$ is the optimal solution to the optimization problem, and thus $\mathcal{A}^{o}=\Pi_{\Omega}(\mathcal{C})$ by definition. ∎

Lemma 4**.**

If $\Omega$ is further affine, then, for any $\mathcal{C}\notin\Omega$, $\mathcal{A}^{o}=\Pi_{\Omega}(\mathcal{C})$ if and only if $\langle\mathcal{C}-\mathcal{A}^{o},\mathcal{A}-\mathcal{A}^{o}\rangle=0$, $\forall\mathcal{A}\in\Omega$.

Proof.

($\Rightarrow$) Let $\mathcal{A}^{o}=\Pi_{\Omega}(\mathcal{C})$ for any given $\mathcal{C}\notin\Omega$, that is, suppose that $\mathcal{A}^{o}$ is the unique solution to the optimization problem. Let $\mathcal{A}\in\Omega$ be such that $\mathcal{A}\neq\mathcal{A}^{o}$. Let $\alpha\in\mathbb{R}$. Since $\Omega$ is affine, $(1-\alpha)\mathcal{A}^{o}+\alpha\mathcal{A}=\mathcal{A}^{o}+\alpha(\mathcal{A}-\mathcal{A}^{o})\in\Omega$. By the assumed optimality of $\mathcal{A}^{o}$, we must have:

[TABLE]

and we obtain:

[TABLE]

If $\alpha\geq 0$, we obtain:

[TABLE]

Now, note that (185) holds for any $\alpha\geq 0$. Since the RHS of (185) can be made arbitrarily small for a given $\mathcal{A}$, the LHS can not be strictly positive. Thus, we conclude:

[TABLE]

If $\alpha\leq 0$, we obtain:

[TABLE]

Now, note that (187) holds for any $\alpha\leq 0$. Since the RHS of (187) can be made arbitrarily large for a given $\mathcal{A}$, the LHS can not be strictly negative. Thus, we conclude:

[TABLE]

Combining (186) and (188), we conclude as desired:

[TABLE]

($\Leftarrow$) Let $\mathcal{A}^{o}\in\Omega$ be such that $\langle\mathcal{C}-\mathcal{A}^{o},\mathcal{A}-\mathcal{A}^{o}\rangle=0,\forall\mathcal{A}\in\Omega$. We shall show that it must be the optimal solution. Let $\mathcal{A}\in\Omega$ and $\mathcal{A}\neq\mathcal{A}^{o}$. We have:

[TABLE]

Hence, $\mathcal{A}^{o}$ is the optimal solution to the optimization problem, and thus $\mathcal{A}^{o}=\Pi_{\Omega}(\mathcal{C})$ by definition. ∎

Now we prove Lemma 2. Let $Y=\Pi_{\Omega_{1}}(\mathcal{C})$. From Lemma 4, we have:

[TABLE]

Let $Z=\Pi_{\Omega_{2}}(Y)$. From Lemma 3, we have:

[TABLE]

For $Z=\Pi_{\Omega_{2}}(\Pi_{\Omega_{1}}(\mathcal{C}))$ to be the projection of $\mathcal{C}$ onto $\Omega_{1}\cap\Omega_{2}$, we need to show from Lemma 3 that:

[TABLE]

under the conditions in Lemma 2. For any $\mathcal{A}\in\Omega_{1}\cap\Omega_{2}$, we have:

[TABLE]

which concludes the proof.

Appendix F Stability of fourth-order error moment

In this Appendix, we show that, under the same settings of Theorem 1 in [2] with the second-order moment condition (24) replaced by the fourth-order moment condition (29), the fourth-order moment of the network error vector is stable for sufficiently small $\mu$, namely, (30) holds for small enough $\mu$. We start by recalling that for any two complex column vectors $x$ and $y$, it holds that $\|x+y\|^{4}\leq\|x\|^{4}+3\|y\|^{4}+8\|x\|^{2}\|y\|^{2}+4\|x\|^{2}\text{Re}(x^{*}y)$ [4, pp. 523]. Applying this inequality to eq. (60) in [2], conditioning on $\boldsymbol{\cal{F}}_{i-1}$, computing the expectations of both sides, using Assumption 2, taking expectations again, and exploiting the convexity of $\|x\|^{4}$, we conclude that:

[TABLE]

for any arbitrary positive number $t\in(0,1)$. In the last inequality we used the bounds (115) and (116) in [2]. By selecting $t=\mu\sigma_{11}$, we arrive at:

[TABLE]

Applying similar arguments for relation (61) in [2] and using the relation $\|a+b+c\|^{4}\leq 27\|a\|^{4}+27\|b\|^{4}+27\|c\|^{4}$, we obtain:

[TABLE]

for any arbitrary positive number $t\in(0,1)$. In the last inequality we used relation (122) in [2]. Selecting $t=\rho(\mathcal{J}_{\epsilon})+\epsilon<1$, we arrive at:

[TABLE]

where $1-t=1-\rho(\mathcal{J}_{\epsilon})-\epsilon$.

In order to bound the fourth-order noise terms $\mathbb{E}\|\overline{\boldsymbol{s}}_{i}^{e}\|^{4}$ and $\mathbb{E}\|\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\boldsymbol{s}}]{\kern-0.3pt\bigwedge\kern-0.3pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\boldsymbol{s}}{\scalebox{-0.8}{\tmpbox}}_{i}^{e}\|^{4}$ appearing in (F) and (F), we first note from eq. (58) in [2] that:

[TABLE]

Now, applying Jensen’s inequality to the convex function $f(x)=x^{2}$, we can write:

[TABLE]

in terms of the individual gradient noise processes, $\mathbb{E}\|\boldsymbol{s}_{k,i}\|^{4}$. For each term $\boldsymbol{s}_{k,i}$, we have from (29) and from the Jensen’s inequality applied to the convex norm $\|x\|^{4}$:

[TABLE]

where $\bar{\beta}_{4,k}^{4}=8(\beta_{4,k}/h)^{4}$ and $\bar{\sigma}^{4}_{s4,k}=8(\beta_{4,k}/h)^{4}\|w^{o}_{k}\|^{4}+\sigma^{4}_{s4,k}$. Using the relations,

[TABLE]

[TABLE]

the term $\mathbb{E}\|\boldsymbol{s}_{i}\|^{4}$ in (195) can be bounded as follows:

[TABLE]

where $\beta^{4}_{4,\max}\triangleq N\max_{1\leq k\leq N}\bar{\beta}_{4,k}^{4}$, and $\sigma^{4}_{s4}\triangleq 4N\sum_{k=1}^{N}\bar{\sigma}^{4}_{s4,k}$. Substituting into (LABEL:eq:_fourth_order_intermediate_3), we get:

[TABLE]

Returning to (F), and using the bounds (128) in [2] and (199), we find that:

[TABLE]

Using the properties that, for any two random variables $\boldsymbol{a}$ and $\boldsymbol{c}$, it holds that [4, pp. 528]:

[TABLE]

we can write:

[TABLE]

so that:

[TABLE]

where:

[TABLE]

Returning to (F) and using similar arguments, we can verify that:

[TABLE]

so that,

[TABLE]

where the coefficients $\{c,d,c^{\prime},d^{\prime},f\}$ have the following form:

[TABLE]

Therefore, we can write

[TABLE]

in terms of the $2\times 2$ coefficient matrix $\Gamma$ of the form (132) in [2] which is stable matrix for sufficiently small $\mu$ and $\epsilon$. Moreover, using relation (135) in [2], we have:

[TABLE]

In this case, we can iterate (217) and use relation (134) in [2] to conclude that:

[TABLE]

and, therefore,

[TABLE]

Appendix G Stability of the coefficient matrix $\mathcal{B}$

In this Appendix, we show that, under the same settings of Theorem 1, the constant matrix $\mathcal{B}$ defined by (42) is stable for sufficiently small step-sizes. To establish this, we use similar argument as in [6, 4]. We first note that the matrix $\mathcal{B}$ in (42) is similar to the matrix $\overline{\mathcal{B}}$ in (95), and therefore has the same eigenvalues as the block matrix $\overline{\mathcal{B}}$ written as:

[TABLE]

where the blocks entries $\{\mathcal{D}_{mn}\}$ are given by (96)–(97). In a manner similar to the arguments used in the proof of Theorem 1 in [2], we can verify that:

[TABLE]

where $\sigma_{11}$ is a positive scalar independent of $\mu$. Thus, we obtain:

[TABLE]

Now recall that the matrix $\mathcal{J}_{\epsilon}^{e}$ defined in Table I is $h(M-P)\times h(M-P)$ and has a Jordan structure. We consider here the complex data case since the real data case can be easily deduced from the complex case by removing the block $(\mathcal{J}_{\epsilon}^{*})^{\top}$. It can be expressed in the following upper-triangular form:

[TABLE]

with scalars $\{\lambda_{a,\ell},\lambda_{a,\ell}^{*}\}$ on the diagonal, all of which have norms strictly less than one, and where the entries of the strictly upper-triangular matrix $\mathcal{K}$ are either $\epsilon$ or zero. It follows that:

[TABLE]

We introduce the eigen-decomposition of the Hermitian positive-definite matrix $\mathcal{D}_{11}$ and denote it by [6, 4]:

[TABLE]

where $U_{d}$ is unitary and $\Lambda_{d}$ has positive diagonal entries $\{\lambda_{k}\}$; the matrices $U_{d}$ and $\Lambda_{d}$ are $hP\times hP$. Using $U_{d}$, we further introduce the following block-diagonal similarity transformation:

[TABLE]

We now use (95) to get:

[TABLE]

where we introduced the $hP\times hP$ diagonal matrix:

[TABLE]

It follows that all off-diagonal entries of the above transformed matrix are at most $O(\mu^{1/hM})$. Although the factor $\mu^{1/hM}$ decays slower than $\mu$, it nevertheless becomes small for sufficiently small $\mu$. Calling upon the Gershgorin’s theorem222Consider an $N\times N$ matrix $A$ with scalar entries $\{a_{k\ell}\}$. With each diagonal entry $a_{kk}$ we associate a disc in the complex plane centered at $a_{kk}$ and with $r_{k}=\sum_{\ell=1,\ell\neq k}^{N}|a_{k\ell}|$. That is, $r_{k}$ is equal to the sum of the magnitudes of the non-diagonal entries on the same row as $a_{kk}$. We denote the disc by $D_{k}$; it consists of all points that satisfy $D_{k}=\{z\in\mathbb{C}\text{ such that }|z-a_{kk}|\leq r_{k}\}$. Gershgorin’s theorem states that the spectrum of $A$ (i.e., the set of all its eigenvalues, denoted by $\lambda(A)$) is contained in the union of all $N$ Gershgorin discs

$\lambda(A)\subset\cup_{k=1}^{N}D_{k}.$

A stronger statement of the Gershgorin theorem covers the situation in which some of the Gershgorin discs happen to be disjoint. Specifically, if the union of the $L$ discs is disjoint from the union of the remaining $N-L$ discs, then the theorem further asserts that $L$ eigenvalues of $A$ will lie in the first union of $L$ discs and the remaining $N-L$ eigenvalues of $A$ will lie in the second union of $N-L$ discs., we conclude that the eigenvalues of $\mathcal{B}$ are either located in the Gershgorin circles that are centered at the eigenvalues of $B$ with radii $O(\mu^{(hM+1)/hM})$ or in the Gershgorin circles that are centered at the $\{\lambda_{a,\ell},\lambda_{a,\ell}^{*}\}$ with radii $O(\mu^{1/M})$, namely,

[TABLE]

where $\lambda(\mathcal{B})$ and $\lambda(B)$ denote any of the eigenvalues of $\mathcal{B}$ and $B$, and $\ell=1,\ldots,L$. It follows that:

[TABLE]

Now since $\mathcal{J}_{\epsilon}$ is a stable matrix, we know that $\rho(\mathcal{J}_{\epsilon})<1$. We express this spectral radius as:

[TABLE]

where $\delta_{J}$ is positive and independent of $\mu$. We also know from (232) that:

[TABLE]

since $B=U^{*}_{d}(I_{hP}-\mathcal{D}_{11})U_{d}$. We conclude from (259) that:

[TABLE]

If we now select $\mu\ll 1$ small enough such that:

[TABLE]

then we would be able to conclude that $\rho(\mathcal{B})<1$ so that $\mathcal{B}$ is stable for sufficiently small step-sizes, as claimed.

If we exploit the structure of $\overline{\mathcal{B}}$ in (95) we can further show, for sufficiently small step-sizes, that:

[TABLE]

where the leading $(1,1)$ block in $(I-\overline{\mathcal{B}})^{-1}$ has dimensions $hP\times hP$.

To establish this we first note that, by similarity, the matrix $\overline{\mathcal{B}}$ is stable. Let

[TABLE]

where from (230)–(231), we have:

[TABLE]

The matrix $\mathcal{X}$ is invertible since $I-\overline{\mathcal{B}}$ is invertible. Moreover, $\mathcal{X}_{11}$ is invertible since $\mathcal{D}_{11}$ is Hermitian positive definite. Using the block matrix inversion formula, we can write:

[TABLE]

where $\Delta$ denotes the Schur complement of $\mathcal{X}$ relative to $\mathcal{X}_{11}$:

[TABLE]

We then use (273)–(274) to conclude (267).

Appendix H Stability of first-order error moment of (47)

In this Appendix, we show that, under the same settings of Theorem 1, the first-order moment of the long-term model (47) is stable for sufficiently small step-sizes, namely, it holds that:

[TABLE]

We first multiply both sides of recursion (138) from the left by $(\mathcal{V}_{\epsilon}^{e})^{-1}$ and use relation (59) in [2] to get:

[TABLE]

where the matrix $\overline{\mathcal{B}}$ in (95) is stable as shown in Appendix G. Recursion (288) can be written more compactly as:

[TABLE]

Since $\overline{\mathcal{B}}$ is stable and $\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{b}]{\kern-0.3pt\bigwedge\kern-0.3pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{b}{\scalebox{-0.8}{\tmpbox}}^{e}=O(\mu)$, we conclude from (291) and (267) that:

[TABLE]

It follows that

[TABLE]

and, hence,

[TABLE]