# Distributed Estimation of Sparse Inverse Covariances

Tong Yao and Shreyas Sundaram Tong Yao and Shreyas Sundaram are with the
Elmore Family School of Electrical and Computer Engineering at Purdue
University. Email: {yao127<EMAIL_ADDRESS>This research was supported
by National Science Foundation (NSF) grant CMMI 1638311.

###### Abstract

Learning the relationships between various entities from time-series data is
essential in many applications. Gaussian graphical models have been studied to
infer these relationships. However, existing algorithms process data in a
batch at a central location, limiting their applications in scenarios where
data is gathered by different agents. In this paper, we propose a distributed
sparse inverse covariance algorithm to learn the network structure (i.e.,
dependencies among observed entities) in real-time from data collected by
distributed agents. Our approach is built on an online graphical alternating
minimization algorithm, augmented with a consensus term that allows agents to
learn the desired structure cooperatively. We allow the system designer to
select the number of communication rounds and optimization steps per data
point. We characterize the rate of convergence of our algorithm and provide
simulations on synthetic datasets.

##  I INTRODUCTION

Several applications involve analyzing many interacting entities, each
generating a large quantity of multi-variate time-series data. Learning the
relationships among these entities, especially in real-time, is essential for
pattern discovery, prediction, and correlation-based clustering. In
distributed scenarios, each agent (observer) can only collect data from a
subset of the entities and communicate with its neighbors within a specific
communication range to learn the relationships of the entire network in real-
time (i.e., online). Such distributed online correlation inference has many
applications, for example, in covariance-based and distributed clustering [1,
2], traffic network prediction and tracking [3, 4], social network inference
[5], and biomedical sensor data inference [6].

To learn the relationships between entities, one approach is to represent each
entity as a node in a graph, and each edge defines an interaction between the
nodes. Learning the network topology from data is studied extensively in the
literature [7, 8, 9, 10]. One method focuses on estimating the non-zero
elements of the sparse inverse covariance matrix of these random variables
(i.e., nodes) [7, 11] through inferring the edges of the graph by using an
$l1$-regularized Gaussian maximum likelihood estimator, assuming the random
variables are jointly Gaussian. In this case, an edge between nodes indicates
that the corresponding random variables are conditionally dependent, given all
the other variables. This is known as the graphical lasso problem [12, 13,
14]. Various optimization algorithms have been proposed to solve this problem,
including coordinate descent [13], proximal methods [15, 16], and alternating
minimization methods [17, 18]. However, these algorithms have not been studied
extensively for distributed and online inference. In [19], an online sparse
inverse covariance algorithm was proposed, but it is not directly applicable
to the distributed setting.

In this paper, given multi-variate time-series data gathered by each agent in
a network, we propose a peer-to-peer distributed algorithm for each agent to
estimate the underlying relationships between all the network variables in
real-time via distributed and online inverse covariance matrix estimation. We
provide theoretical guarantees on the convergence of the estimates and
characterize the asymptotic rate of convergence. We also demonstrate the
performance of our proposed algorithm through simulations.

##  II PROBLEM FORMULATION

Consider a set of $p$ random variables $X=\\{x_{1},x_{2},\ldots,x_{p}\\}$ that
are jointly Gaussian with zero mean and covariance $S^{*}$. These variables
are represented by a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\\{v_{1},\ldots,v_{p}\\}$ is the set of nodes, with each node
$v_{i}$ representing a random variable $x_{i}$. An edge
$(v_{i},v_{j})\in\mathcal{E}$ indicates that variable $x_{j}$ is conditionally
dependent on $x_{i}$, given all the other random variables. Conversely, if
$(v_{i},v_{j})\notin\mathcal{E}$, $v_{j}$ is conditionally independent of
$v_{i}$, given all the other variables. This lack of an edge corresponds to a
zero-entry in the inverse covariance matrix $(S^{*})^{-1}$ (e.g., see [13,
17]).

These relationships (i.e., graph structure) between the variables are unknown
a priori; the goal is to infer the edges of the graph based on samples of the
random variables. Specifically, at each time step $t\in\\{1,2,\ldots\\}$, the
network generates data
$X_{t}=\left[\begin{matrix}x_{1t}&x_{2t}&\cdots&x_{pt}\end{matrix}\right]^{T}\in\mathbb{R}^{p}$.
We assume each $X_{t}$ is independently and identically sampled from the
underlying Gaussian distribution, i.e., $X_{t}\sim\mathcal{N}(0,S^{*})$.

Consider a group of $n$ agents, $n\in\mathbb{N}$, with a known, unweighted,
undirected, and connected communication graph
$\mathcal{G}_{a}=(\mathcal{V}_{a},\mathcal{E}_{a})$, where
$\mathcal{V}_{a}=\\{1,2,\ldots,n\\}$ is the set of vertices representing the
agents and $\mathcal{E}_{a}\subseteq\mathcal{V}_{a}\times\mathcal{V}_{a}$ is
the set of edges. Note that the communication graph $\mathcal{G}_{a}$ is not
to be confused with the correlation graph $\mathcal{G}$. If an edge
$({i},{j})\in\mathcal{E}_{a}$, agent $i$ and $j$ can communicate with each
other. The neighbors of agent $i\in\mathcal{V}_{a}$ in graph $\mathcal{G}_{a}$
are represented by the set
$N_{i}=\\{j\in\mathcal{V}_{a}:(i,j)\in\mathcal{E}_{a}\\}$. We let
$\mathcal{N}_{i}$ denote the set of neighbors and agent $i$ itself, i.e.,
$\mathcal{N}_{i}=N_{i}\cup\\{i\\}$. Each agent $i\in\mathcal{V}_{a}$ observes
a subset of the random variables $X_{i}\subseteq X$ and collects the
corresponding time-series data $X_{i,t}\subseteq X_{t}$. The time-series data
collected by the neighbors of agent $i$ is denoted by $X_{N_{i},t}\subseteq
X_{t}.$

Given a set of data $\\{X_{1},X_{2},\ldots,X_{t}\\}$ up to time $t$, the
centralized maximum likelihood estimation problem is given by

$\operatorname*{minimize\,}_{\Omega_{t}\in\mathcal{S}_{++}^{p}}-\log\det\Omega_{t}+\operatorname{tr}(S_{t}\Omega_{t})+\lambda|\Omega_{t}|_{l1},$
(1)

where the set of $p\times p$ positive definite matrices is denoted by set
$\mathcal{S}_{++}^{p}$, and $S_{t}=\frac{1}{t}\sum_{i=1}^{t}X_{i}X_{i}^{T}$ is
the sample covariance matrix constructed from all of the data up to time $t$.
The terms $-\log\det\Omega_{t}+\operatorname{tr}(S_{t}\Omega_{t})$ are derived
from the Gaussian log-likelihood function [12], where $\operatorname{tr}$
denotes trace, and the term
$|\Omega_{t}|_{l1}=\sum_{i,j=1}^{p}|\Omega_{t}(i,j)|$ is the element-wise
$l_{1}$ norm, encouraging sparsity of the solution regulated by the penalty
parameter $\lambda\geq 0$.

For the distributed setting, no agent will have access to the entire sample
covariance matrix, since each agent only observes a subset of the variables.
Instead, we will have each agent $i\in\mathcal{V}_{a}$ maintain an estimate
$S_{i,t}$ of the sample covariance matrix, which it updates over time based on
its own measurements and information received from its neighbors. At each time
step $t\in\\{1,2,\ldots,\\}$, the optimization objective for agent
$i\in\mathcal{V}_{a}$ would then be to find an estimate $\Omega_{i,t}$ of
$(S^{*})^{-1}$ by solving

$\operatorname*{minimize\,}_{\Omega_{i,t}\in\mathcal{S}_{++}^{p}}-\log\det\Omega_{i,t}+\operatorname{tr}(S_{i,t}\Omega_{i,t})+\lambda|\Omega_{i,t}|_{l1},$
(2)

where we have simply updated (1) to show the explicit dependence of the
estimated sample covariance matrix $S_{i,t}$ of agent $i$ on the current time
$t$.

This paper aims to formulate a distributed online algorithm for each agent $i$
to calculate and update an estimate of the inverse covariance matrix
$\Omega^{*}=(S^{*})^{-1}$ in real-time. Given a sequence of observations of
agent $i$, denoted by $\\{X_{i,1},X_{i,2},X_{i,3},\ldots\\}$, the objective of
$i$ is to perform an online inference of the entire edge set $\mathcal{E}$,
through estimating the inverse covariance matrix $\Omega_{i,t}$ based on its
local data and by incorporating information from its neighbors. In particular,
we want all agents to reach consensus on their estimates asymptotically, i.e.,
as $t\to\infty$, the estimate of each agent converges asymptotically to the
agreement:
$\lim_{t\to\infty}\Omega_{1,t}=\Omega_{2,t}=\ldots=\Omega_{n,t}=\Omega^{*}.$

In the rest of the paper, we describe a distributed algorithm to solve this
problem, allowing the system designer to specify the number of communications
and optimization steps of the algorithm between the arrival of data points. We
start by describing the online alternating minimization algorithm for solving
problem (1) and subsequently describe the extension of the algorithm to the
distributed setting.

##  III Background

We build our approach on an alternating minimization algorithm proposed in
[18] (batch) and [19] (online) for solving problem (1); the batch algorithm
was shown in [18] to be significantly faster than other proximal methods such
as GISTA [15] and QUIC[16]; the online algorithm was shown in [19] to achieve
a similar result with fewer iterations in real-time settings. We first
describe the details of the centralized algorithm from [18, 19], and
subsequently discuss our modification to account for the distributed
estimation.

The approach in [18] formulates the primal and dual objective functions for
problem (1). The primal of (1) is:

$\displaystyle\begin{split}\operatorname*{minimize\,}_{\Omega_{t}\in\mathcal{S}^{p}_{++},\Phi_{t}\in\mathcal{S}^{p}_{++}}&-\log\det\Omega_{t}+\operatorname{tr}(S_{t}\Phi_{t})+\lambda|\Phi_{t}|_{l1}\\\
\operatorname*{subject\,to\,}&\Phi_{t}=\Omega_{t}.\end{split}$ (3)

The dual of (1) is given by

$\displaystyle\begin{split}\operatorname*{minimize\,}_{\Gamma_{t}\in\mathcal{S}_{++}^{p}}&-\log\det\Gamma_{t}-p\\\
\operatorname*{subject\,to\,}&|\Gamma_{t}(i,j)-S_{t}(i,j)|\leq\lambda\quad\forall
i,j,\end{split}$ (4)

where the symmetric positive definite matrix $\Gamma_{t}$ is the dual
variable, and $A(i,j)$ denotes the $(i,j)$-th element of matrix $A$.

Given the sample covariance matrix $S_{t}$, the alternating minimization
follows the iterative sequence of updates, where each iteration is indexed by
the variable $k\in\mathbb{N}$:

$\displaystyle\Omega^{k+1}_{t}$
$\displaystyle=\operatorname*{argmin\,}_{\Omega\in\mathcal{S}^{p}_{++}}-\log\det\Omega+\operatorname{tr}(\Gamma_{t}^{k}\Omega),$
(5) $\displaystyle\Phi^{k+1}_{t}$
$\displaystyle=\operatorname*{argmin\,}_{\Phi\in\mathcal{S}^{p}_{++}}\operatorname{tr}(S_{t}\Phi)+\lambda|\Phi|_{l1}-\operatorname{tr}(\Gamma^{k}_{t}\Phi)$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{\zeta_{t}^{k}}{2}\|\Omega_{t}^{k+1}-\Phi\|_{F}^{2},$
(6) $\displaystyle\Gamma^{k+1}_{t}$
$\displaystyle=\Gamma_{t}^{k}+\zeta_{t}^{k}(\Omega_{t}^{k+1}-\Phi_{t}^{k+1}).$
(7)

In the above equations, $\zeta_{t}^{k}$ is a step size, and
$\|A\|_{F}=\sqrt{\operatorname{tr}(AA^{T})}$ denotes the Frobenius norm of a
given matrix $A$. Taking the derivatives of the expressions for $\Omega^{k+1}$
and $\Phi^{k+1}$ and equating them to 0, we obtain the closed-form updates
[18]:

$\displaystyle\Omega^{k+1}_{t}$ $\displaystyle=(\Gamma^{k}_{t})^{-1},$ (8)
$\displaystyle\Phi^{k+1}_{t}$
$\displaystyle=\frac{1}{\zeta_{t}^{k}}\mathcal{S}_{\lambda}(\zeta_{t}^{k}\Omega^{k+1}_{t}-S_{t}+\Gamma^{k}_{t}).$
(9)

Here, $\mathcal{S}_{\lambda}(x)=\operatorname*{sign}(x)(\max(|x|-\lambda,0))$
is the soft-thresholding operator (applied element-wise to a matrix argument).
Following these update rules, $\Omega_{t}$ is interpreted as an approximately
sparse inverse covariance matrix, and $\Phi_{t}$ is interpreted as the
estimate of the sparse inverse covariance matrix. Substituting (8) – (9) for
the variables in (7), and using the clip function
$\mathcal{C}_{\lambda}(x)=\min(\max(x,-\lambda),\lambda)$ with the property
$x=\mathcal{S}_{\lambda}(x)+\mathcal{C}_{\lambda}(x)$, the dual update (7) can
be written as:

$\displaystyle\Gamma^{k+1}_{t}=\mathcal{C}_{\lambda}(\Gamma^{k}_{t}-S_{t}+\zeta_{t}^{k}(\Gamma^{k}_{t})^{-1})+S_{t}.$
(10)

In [19], the authors proposed an online algorithm (OGAMA) that iteratively
updates the sample covariance matrix at each time step $t$ (at the arrival of
data points $X_{t}$) by

$S_{t}=\frac{1}{t}\left(\left(t-1\right)S_{t-1}+X_{t}X_{t}^{T}\right).$

The algorithm allows the system designer to select the total number of
optimization iterations $K\in\mathbb{N}$ per new data point, and initializes
$\Gamma_{t_{0}}^{K}=S_{t_{0}}+\lambda I_{p}$, where
$I_{p}\in\mathbb{R}^{p\times p}$ is an identity matrix, with data up to a
user-defined $t_{0}\in\\{1,2,\ldots\\}$. The algorithm updates the estimates
by iterating through (10) $K$ times at each $t>t_{0}$, and updates (8) and (9)
at the end of each $t$ when $k=K$.

The step size $\zeta_{t}^{k}$ at each iteration $k$ is chosen to guarantee
convergence of the estimates $\Gamma_{t}^{k}$ to their desired quantities. It
was shown in [19] that $\forall t\geq t_{0}$ and $\forall
k\in\\{1,2,\ldots,K\\}$, the step size can be set as a constant
$\zeta=\zeta_{t}^{k}<a^{2}$ for some constant $a$.

We provide the pseudo-code implementation of the OGAMA algorithm in Algorithm
1. Note that we modified the variables to demonstrate the computation at each
agent $i$. In the centralized OGAMA, the input $S_{i,t}^{w}=S_{t}$ and the
output $\Gamma_{i,t}^{k}=\Gamma_{t}^{k}$, $\forall i\in\mathcal{V}_{a}$ and
$\forall w\in\\{1,2,\ldots,W\\}$, where $w$ denotes the communication
iterations which we discuss in the next section.

Input: $t,k,S_{i,t}^{w},\Gamma_{i,t}^{k-1},\zeta_{i,t}^{k-1}$

Result: $\Gamma_{i,t}^{k},\zeta_{i,t}^{k}$

Parameter : $K$, $\lambda$, $t_{0}$

if  _$t=t_{0}$_ then

$\Gamma_{i,t}^{k}=\Gamma_{i,t_{0}}^{K}=S_{i,t_{0}}^{w}+\lambda I_{p}$

else if  _$t >t_{0}$_ then

$\Gamma^{k}_{i,t}=\mathcal{C}_{\lambda}(\Gamma^{k-1}_{i,t}-S_{i,t}^{w}+\zeta_{i,t}^{k-1}(\Gamma_{i,t}^{k-1})^{-1})+S_{i,t}^{w}$

Choose $\zeta_{i,t}^{k}\in(0,(\lambda_{\min}(\Gamma_{i,t}^{k}))^{2})$

Algorithm 1 Online Graphical Alternating Minimization Algorithm (OGAMA)

##  IV Communication Protocols and Distributed Algorithm

In this section, we first discuss our assumptions on the problem. We then
propose the communication protocols in Algorithm 2 and the high level flow of
function executions for DGAMA in Algorithm 3.

###### Assumption 1

Each agent $i\in\mathcal{V}_{a}$ knows its set of measured variables $X_{i}$,
represented by a diagonal matrix $V_{i}\in\mathbb{R}^{p\times p}$. The $j$-th
diagonal of $V_{i}(j,j)=1$ if $x_{j}\in X_{i}$ and $V_{i}(j,j)=0$ otherwise.
Also, each agent $i$ knows its neighbors ${N_{i}}$ and the nodes measured by
its neighbors $V_{N_{i}}$.

###### Assumption 2

The agents have a synchronized clock such that all the agents have the same
value of $t$.

To tailor the algorithm for the real-time multi-agent setting, we allow the
system designer to select the maximum number of communication iterations $W$
and the maximum number of optimization iterations $K$ per time-step. We
discuss the range of $W$ and $K$ in Sec. V. The selection of $W$ and $K$ is
also based on the number of variables $p$, the hardware capabilities, the data
arrival rate, etc. When communication (or computation) is expensive or slow,
$W$ (or $K$) is set to be small; on the other hand, if communication (or
computation) is cheap or fast, $W$ (or $K$) can be set to be large. We let $W$
and $K$ remain the same for each time step for simplicity. We index each
communication iteration by the variable $w\in\\{1,2,\ldots,W\\}$ and index the
optimization iteration by the variable $k\in\\{1,2,\ldots,K\\}$. In
particular, we denote agent $i$’s estimate of $S_{t}$ after $w$ rounds of
communication at time-step $t$ by $S_{i,t}^{w}$.

We initialize each agent’s estimated sample covariance matrix as
$S_{i,0}=\textbf{0}\in\mathbb{R}^{p\times p}$. At each time step $t$, the
agents communicate with their neighbors to update $S_{i,t}^{w}$. The pseudo-
code implementation to communicate and update $S_{i,t}^{w}$ is in Algorithm 2.

At the beginning of time $t$, when $w=1$, upon the arrival of the observation
data $X_{i,t}$, each agent $i$ performs a round of communication by sending
its $X_{i,t}$ and receiving $X_{N_{i},t}$ from its neighbors. Based on this
communication data, we introduce the concept of an observable node in the
following definition.

###### Definition 1

We say node $x_{j}$ is observable by agent $i$ if $x_{j}\in
X_{\mathcal{N}_{i}}=X_{i}\cup X_{N_{i}}$, i.e., the node can be directly
observed by agent $i$ or indirectly observed through data from ${N_{i}}$.
Similarly, we say the $(l,m)$-th entry of sample covariance matrix
$S_{i}(l,m)$ is observable by agent $i$ if $x_{l}$ and $x_{m}$ are both
observable by agent $i$.

Observations of the pairwise relationships are required to solve our inverse
covariance estimation problem; thus, each pair of nodes must be observed by at
least one agent. With that in mind, we proposed the following definition.

###### Definition 2

The pairwise relationships of all nodes (i.e., random variables) are jointly
observable if

$\cup_{i}\\{X_{\mathcal{N}_{i}}\times X_{\mathcal{N}_{i}}\\}=X\times X.$

Joint observability ensures that for each pair of nodes, there is at least one
agent that can observe and compute the sample covariance between these nodes.
A demonstration of Def. 2 is presented in Fig. 1.

Figure 1: Four agents $(a_{1},\ldots,a_{4})$ observing five nodes
$(x_{1},\ldots,x_{5})$. The solid lines connecting the agents indicate
communication; the solid lines connecting the nodes indicate correlations; the
dashed lines between the agents and the nodes indicate that the agents
directly observe nodes. In this example, the sample covariance matrix is
jointly observable, i.e., for each pair of nodes, there is an agent that can
observe the pairwise relationship. The entries of the sample covariance matrix
observed by agent $1$ are omitted in the presentation.

After the communications among the agents, each agent computes the sample
covariance matrix $S_{i,t}^{w}$. If $w=1$, with data $X_{i,t}$ and
$X_{N_{i},t}$, each agent creates a data vector $\chi_{i,t}\in\mathbb{R}^{p}$.
If a variable $x_{j}$ is observable by agent $i$ at time $t$, the
corresponding entry $j$ of $\chi_{i,t}(j)=x_{jt}$, and $\chi_{i,t}(j)=0$
otherwise. With $\chi_{i,t}$, each agent computes $S_{i,t}^{1}$ with the
following updates:

$\displaystyle S_{il,t}$
$\displaystyle=\frac{1}{t}\left(\left(t-1\right)V_{i}S_{i,t-1}^{W}V_{i}+\chi_{i,t}\chi_{i,t}^{T}\right)$
(11) $\displaystyle S_{i,t}^{1}$
$\displaystyle=S_{il,t}+\frac{\sum_{j\in\mathcal{N}_{i}}S_{j,t-1}^{W}-V_{i}\sum_{j\in\mathcal{N}_{i}}S_{j,t-1}^{W}V_{i}}{|\mathcal{N}_{i}|}.$
(12)

Each agent updates its local sample covariance matrix, denoted by $S_{il,t}$,
computed recursively following (11), consisting of data up to time step $t$.
The outer product $\chi_{i,t}\chi_{i,t}^{T}$ creates a block matrix, such
that, if $x_{l}$ and $x_{m}$ are both observable by agent $i$, the $(l,m)$-th
entry $\chi_{i,t}\chi_{i,t}^{T}(l,m)=x_{lt}x_{mt}$, and 0 otherwise.
Similarly, the term $V_{i}S_{i,t-1}^{W}V_{i}$ creates a block matrix such that
$V_{i}S_{i,t-1}^{W}V_{i}(l,m)=S_{i,t-1}^{W}(l,m)$ if $x_{l}$ and $x_{m}$ are
both observable by $i$, and 0 otherwise. In (12), the sample covariance matrix
$S_{i,t}^{1}$ of $i$ is created by combining the observable entries from the
block matrix $S_{il,t}$, with the unobservable entries from
$\sum_{j\in\mathcal{N}_{i}}S_{j,t-1}^{W}/|\mathcal{N}_{i}|$. If the entries of
$S_{i,t}^{w}$ are observable by $i$, agent $i$ uses its information
exclusively; if the entries of $S_{i,t}^{w}$ are not observable by $i$, agent
$i$ uses the information from its neighbors. Agent $i$ then sends its
neighbors the new $S_{i,t}^{1}$ and receives $S_{N_{i},t}^{1}$ from its
neighbors.

When $w>1$, the update of $S_{i,t}^{w}$ follows

$S_{i,t}^{w}=V_{i}S_{i,t}^{w-1}V_{i}+\frac{\sum_{j\in\mathcal{N}_{i}}S_{j,t}^{w-1}-V_{i}\sum_{j\in\mathcal{N}_{i}}S_{j,t}^{w-1}V_{i}}{|\mathcal{N}_{i}|},$
(13)

where $V_{i}S_{i,t}^{w-1}V_{i}$ contains the observable entries of agent $i$
from the previous step (0 otherwise). Then agent $i$ computes $S_{i,t}^{w}$ by
integrating the observable entries from the previous iteration $S_{i,t}^{w-1}$
and unobservable entries given by
$\sum_{j\in\mathcal{N}_{i}}S_{j,t}^{w-1}/|\mathcal{N}_{i}|$. Agent $i$ then
sends its neighbors the new $S_{i,t}^{w}$ and receives $S_{N_{i},t}^{w}$ from
its neighbors.

Input: $t,w,S_{\mathcal{N}_{i},t}^{w-1}$

Result: $S_{\mathcal{N}_{i},t}^{w}=S_{i,t}^{w}\cup S_{N_{i},t}^{w}$

Parameter : $V_{i},N_{i}$

if  _$w=1$_ then

Send data $X_{i,t}$ and receive neighbor data $X_{N_{i},t}$

$X_{\mathcal{N}_{i},t}=X_{i,t}\cup X_{N_{i},t}$

Create $\chi_{i,t}$ with $X_{\mathcal{N}_{i},t}$ as described above (11)

Update $S_{il,t}$ as per (11)

Update $S_{i,t}^{1}$ as per (12)

else

Update $S_{i,t}^{w}$ as per (13)

Send $S_{i,t}^{w}$ and receive $S_{N_{i},t}^{w}$

Algorithm 2 Communication for Agent $i\in\mathcal{V}_{a}$ (COM)

In Algorithm 3, we present the flow of function executions. Each agent $i$
performs $W$ rounds of communication and obtains $S_{i,t}^{W}$ at each time
step $t$; subsequently, $i$ updates its local estimates $\Gamma_{i,t}^{K}$
after performing $K$ optimization iterations.

Input: $W,K,\lambda,t_{0}$

Initialize $S_{j,1}^{0}=S_{j,0}^{W}=\textbf{0}\in\mathbb{R}^{p\times
p}\;\forall j\in\mathcal{N}_{i}$

for  _$t\in\\{1,2,3,\ldots\\}$_ do

for  _$w=1:W$_ do

$S_{\mathcal{N}_{i},t}^{w}=$ COM($t,w,S_{\mathcal{N}_{i},t}^{w-1}$)
(consensus)

Extract $S_{i,t}^{W}$ from $S_{\mathcal{N}_{i},t}^{W}$

if  _$t=t_{0}$_ then

$\Gamma_{i,t}^{K},\zeta_{i,t}^{K}$ = OGAMA($t,K,S_{i,t}^{W},null,null$)

if  _$t >t_{0}$_ then

for  _$k=1:K$_ do

if  _$k=1$_ then

$\Gamma_{i,t}^{0}=\Gamma_{i,t-1}^{K}$

$\Gamma_{i,t}^{k},\zeta_{i,t}^{k}$ =
OGAMA($t,k,S_{i,t}^{W},\Gamma_{i,t}^{k-1},\zeta_{i,t}^{k-1}$)

Update the variables $\Omega_{t}^{K}$, $\Phi_{t}^{K}$ as per (8), (9)

Algorithm 3 Distributed Graphical Alternating Minimization Algorithm for Agent
$i\in\mathcal{V}_{a}$ (DGAMA)

##  V Convergence Analysis

In this section, we show that the dual variable estimate $\Gamma_{i,t}^{K}$ of
each agent $i$ converges to $\Gamma^{*}$, where $\Gamma^{*}=(\Omega^{*})^{-1}$
is the centralized fixed point solution of (10) and $\Omega^{*}$ is the
centralized optimal solution of (1), given the ground truth covariance matrix
$S^{*}$ [18]. The proofs of our results can be found in the Appendix.

###  V-A Error Bounds of Dual Variables

To simplify the analysis, we assume that $\forall i\in\mathcal{V}_{a},t\geq
t_{0}$, and $k\in\\{1,2,\ldots,K\\}$, we can select a constant step size
$\zeta=\zeta_{i,t}^{k}<a^{2}$, such that $0\prec
aI_{p}\preceq\Gamma_{i,t}^{k}\preceq b$ and $aI_{p}\preceq\Gamma^{*}\preceq
bI_{p}$. In Appendix -A, we show that constants $a$ and $b$ exist.

We provide the following results bounding the difference between the estimate
$\Gamma_{i,t}^{K}$ of each agent and the centralized optimal solution
$\Gamma^{*}$, after $K$ iterations at each time step $t$ (given the new data
point $X_{i,t}$), denoted by $\|\Gamma_{i,t}^{K}-\Gamma^{*}\|_{F}.$

###### Theorem 1

Assume $\forall i\in\mathcal{V}_{a}$, $\forall t\geq t_{0}$, and $\forall
k\in\\{1,2,\ldots,K\\}$, the iterates $\Gamma_{i,t}^{k}$ satisfy
$aI_{p}\preceq\Gamma_{i,t}^{k}\preceq bI_{p}$, and
$aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$, for some fixed constants $0<a<b$, and
$\zeta_{i,t}^{k}=\zeta$. After $W$ communication and $K$ optimization
iterations of Algorithm 3, at time step $t$, we have the following bound for
the dual variable:

$\|\Gamma_{i,t}^{K}-\Gamma^{*}\|_{F}\leq\beta^{K(t-t_{0})}\|\Gamma_{i,t_{0}}^{K}-\Gamma^{*}\|_{F}\\\
+2\sum_{m=1}^{K}\beta^{K-m}\sum_{l=t_{0}+1}^{t}\beta^{K(t-l)}\|S_{i,l}^{W}-S^{*}\|_{F},$
(14)

where
$\beta=\max\left\\{|1-\frac{\zeta}{a^{2}}|,|1-\frac{\zeta}{b^{2}}|\right\\}.$

###  V-B Convergence of Sample Covariance Matrix

To analyze the convergence of the estimates $\Gamma_{i,t}^{k}$ of the
algorithms, we show that the sample covariance matrix $S_{i,t}^{w}$ of each
agent converges to the true sample covariance matrix $S^{*}$. We consider each
entry $(l,m)$ of the sample covariance matrix $S_{t}$. We let
$L\subset\mathcal{V}_{a}$ denote the subset of agents that are observing
$S_{t}(l,m)$ and call these agents leaders. The other agents who cannot
observe $S_{t}(l,m)$ are called followers and denoted by
$F=\mathcal{V}_{a}\setminus L$. For our discussion below, we focus on a
specific pair of $(l,m)$, and without loss of generality, we let the first
$|L|$ agents $\\{1,2,\ldots,{|L|}\\}$ be leaders and let the last $|F|$ agents
$\\{{|L|+1},{|L|+2},\ldots,{n}\\}$ be followers. The same analysis holds for
all $l,m$.

For entry $(l,m)$, if $i\in L$, then $S_{i,t}^{w}(l,m)=S_{t}(l,m),\forall t$
and $\forall w=\\{1,2,\ldots,W\\}$, where $S_{t}$ is the sample covariance
matrix computed given data $\\{X_{1},X_{2},\ldots,X_{t}\\}$ up to time $t$. If
$S_{t}(l,m)$ is not observed by agent $i$, the update of unobserved entries in
(12) and (13) can be written as

$S_{i,t}^{w}(l,m)=\frac{1}{|\mathcal{N}_{i}|}\sum_{j\in\mathcal{N}_{i}}S_{j,t}^{w-1}(l,m)\;\forall
i\in F.$ (15)

We can write this update for the follower agents in the matrix form,

$\displaystyle S_{F,t}^{w}(l,m)$
$\displaystyle=\tilde{D}^{-1}(\tilde{I}+\tilde{A})S_{t}^{w-1}(l,m)$
$\displaystyle\triangleq P_{F,lm}S_{t}^{w-1}(l,m),$

where matrix
$\tilde{D}=\operatorname{diag}(|\mathcal{N}_{|L|+1}|,|\mathcal{N}_{|L|+2}|,\ldots,|\mathcal{N}_{n}|)$,
$\tilde{I}$ is the last $|F|$ rows of identity matrix $I_{n}$, $\tilde{A}$ is
the last $|F|$ rows of an unweighted adjacency matrix, and
$S_{t}^{w}(l,m)=[S_{1,t}^{w}(l,m),S_{2,t}^{w}(l,m),\ldots,S_{n,t}^{w}(l,m)]^{T}$.

Thus, for each $(l,m)$ entry of $S_{t}$, the updates of Algorithm 2 can be
written in the form

$\begin{bmatrix}S_{L,t}^{w}(l,m)\\\
S_{F,t}^{w}(l,m)\end{bmatrix}=\begin{bmatrix}P_{LL,lm}&P_{LF,lm}\\\
P_{FL,lm}&P_{FF,lm}\end{bmatrix}\begin{bmatrix}S_{L,t}^{w-1}(l,m)\\\
S_{F,t}^{w-1}(l,m)\end{bmatrix}.$

Since the leader agents keep their estimations constant, $P_{LL,lm}=I_{|L|}$
and $P_{LF,lm}=\textbf{0}$. The matrix $P_{FL,lm}$ is the first $|L|$ columns
and $P_{FF,lm}$ is the last $|F|$ columns of $P_{F,lm}$. It was shown that the
convergence rate of the communication network is dictated by the largest
eigenvalue of $P_{FF,lm}$ [20]. We provide the following lemma using
$P_{FF,lm}$.

###### Lemma 1

If all variables are jointly observable, the sample covariance matrix of each
agent $S_{i,t}^{w}$ converges to $S_{t}$ as $w\to\infty$, with a rate given by

$\|S_{i,t}^{w}-S_{t}\|_{F}\leq c\sigma^{w}\|S_{i,t}^{0}-S_{t}\|_{F},$ (16)

where $\sigma=\max_{lm}\lambda_{\max}(P_{FF,lm})$ denotes the largest
eigenvalue of $P_{FF,lm}$ for all $(l,m)$, $c\geq 1$ is some constant related
to the degrees of $\mathcal{G}_{a}$ [21], and $S_{i,t}^{0}=S_{i,t-1}^{w}$
denotes $i$’s initial estimate of the sample covariance matrix at time $t$.

Lemma 1 shows that the states of all agents converge at a rate given by the
largest eigenvalue of $P_{FF,lm}$. In this work, we will omit the discussion
on the bound of $\sigma$ and refer the readers to [22, 21, 20]. Since
$0\leq\sigma<1$ ($\mathcal{G}_{a}$ is connected), all agents reach consensus
when the rounds of communication $W\to\infty$, i.e.,
$\lim_{W\to\infty}S_{i,t}^{W}=S_{t},\forall i\in\mathcal{V}_{a}$ . To
guarantee the convergence of DGAMA, we let $W$ be large enough such that each
agent’s estimate of $S_{t}$ after $W$ rounds of communication $S_{i,t}^{W}$ is
better than the initial estimate $S_{i,t}^{0}$, i.e., $c\sigma^{W}<1$. Next,
we show that all agents reach consensus for such $W$ as the number of time
steps $t\to\infty$.

###### Theorem 2

If all variables are jointly observable, the estimates $S_{i,t}^{W}$ of sample
covariance matrices from Algorithm 2 converge to the ground truth sample
covariance matrix $S^{*}$ almost surely as $t\to\infty$, for all
$i\in\mathcal{V}_{a}$ and for all $W\in\mathbb{N}$ such that $c\sigma^{W}<1$:

$\lim_{t\to\infty}\|S_{i,t}^{W}-S^{*}\|_{F}=0\;a.s.$

###  V-C Convergence of Dual Variables

In this subsection, we show the convergence and the rate of convergence of the
dual variables $\Gamma_{i,t}^{K}$ as $t\to\infty$.

###### Corollary 1

Assume there exist constants $0<a<b$ such that, $\forall t\geq t_{0}$,
$\forall k\in\\{1,2,\ldots,K\\}$, and $\forall i\in\mathcal{V}_{a}$, the
quantities $\Gamma_{i,t}^{k}$ satisfy $aI\preceq\Gamma_{i,t}^{k}\preceq bI$,
$aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$, and $\zeta_{i,t}^{k}=\zeta<a^{2}$.
Then, for all $W\in\mathbb{N}$ such that $c\sigma^{W}<1$, as the number of
data points $t\to\infty$, the result $\Gamma_{i,t}^{K}$ converges to the
optimal solution $\Gamma^{*}$ almost surely:

$\lim_{t\to\infty}\|\Gamma_{i,t}^{K}-\Gamma^{*}\|=0\;a.s.$

Having established that the iterates provided by DGAMA converge, we
characterize the asymptotic rate of convergence for the dual variable of each
agent $\Gamma_{i,t}^{K}$ to $\Gamma^{*}$.

###### Corollary 2

Assume $\forall i\in\mathcal{V}_{a}$, $\forall t\geq t_{0}$ and $\forall
k\in\\{1,2,\ldots,K\\}$, the iterates $\Gamma_{i,t}^{k}$ satisfy $0\prec
aI_{p}\preceq\Gamma_{i,t}^{k}\preceq bI_{p}$, $\Gamma^{*}$ satisfies
$aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$, and $\zeta_{i,t}^{k}=\zeta$. Then for
all $\Delta\in(0,\frac{1}{2})$, for all sample paths in a set of measure 1,
there exists a $\bar{t}$, such that $\forall t\geq\bar{t}+1$, the update of
dual variable $\Gamma_{i,t}^{K}$ satisfies the following condition:

$\|\Gamma_{i,t}^{K}-\Gamma^{*}\|_{F}\leq\beta^{K(t+1-\bar{t})}\|\Gamma_{\bar{t}}^{K}-\Gamma^{*}\|_{F}\\\
+2\sum_{m=1}^{K}\beta^{K-m}\sum_{l=\bar{t}+1}^{t}\beta^{K(t-l)}\Big{[}(c\sigma^{W})^{l-\bar{t}}\|S_{i,\bar{t}}^{W}-S_{\bar{t}}^{*}\|_{F}\\\
+40p\max_{j}\left(S^{*}\left(j,j\right)\right)\\\
\cdot\big{(}(\frac{1}{l})^{\frac{1}{2}-\Delta}+2\sum_{j=\bar{t}+1}^{l}(c\sigma^{W})^{l+1-j}(\frac{1}{j-1})^{\frac{1}{2}-\Delta}\big{)}\Big{]},$
(17)

where
$\beta=\max\left\\{|1-\frac{\zeta}{a^{2}}|,|1-\frac{\zeta}{b^{2}}|\right\\}$.

Note from the above result that when $0<\beta<1$, the asymptotic rate at which
$\Gamma_{i,t}^{K}$ decreases to zero is dominated by the second term in (17).
Specifically, the rate is dominated by the convergence of the sample
covariance matrix. Also, increasing $K$ will speed up the convergence of the
estimates $\Gamma_{i,t}^{K}$. Comparing to the rate of convergence in the
centralized [19], increasing $W$ will decrease the deviation of each agent’s
estimates to the centralized estimates. We illustrate the dependence of this
bound on $t$, $K$, and $W$ in the next section with experiments.

(a) $K=1,W=1$

(b) $K=1,W=2$

(c) $K=2,W=2$

Figure 2: Simulations with varying parameters.

##  VI Experiments

We build the agent and node networks shown in Fig. 1 ($p=5,n=4$). We build the
node network by generating a sparse Erdos-Renyi network representing
$({S^{*}})^{-1}$, following the steps in [23]. We generate $40$ independent
and identically distributed samples from a Gaussian distribution
$\mathcal{N}(0,S^{*})$ and run DGAMA Algorithm with $\lambda=0.15$ and
$t_{0}=10$.

In Fig. 2, we show the convergence of our algorithm with different values of
$K$ (the number of iterations per data point) and $W$ (the number of
communication rounds per data point) in comparison with the centralized online
algorithm from [19]. The deviation of the iterates, measured by
$\|\Gamma_{i,t+1}^{K}-\Gamma^{*}\|_{F}$, converges to zero as the number of
data points $t$ increases. We see that the distributed algorithm provides
estimates comparable to that of the centralized online algorithm, with small
$K$ and $W$. Moreover, as $W$ increases, the estimates are closer to that of
the centralized algorithm; as $K$ increases, the estimates of the algorithm
converges to $\Gamma^{*}$ faster.

##  VII Conclusion

In this paper, we proposed a distributed and online sparse inverse covariance
matrix estimation algorithm. We provided convergence guarantees of each
agent’s estimates to what would have been obtained at a centralized agent if
all of the data were available simultaneously. Finally, we demonstrated the
empirical performance of our algorithm.

## References

  * [1] D. Hallac, S. Vare, S. Boyd, and J. Leskovec.  Toeplitz inverse covariance-based clustering of multivariate time series data.  In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’17, page 215–223, New York, NY, USA, 2017. Association for Computing Machinery.
  * [2] P. A. Forero, A. Cano, and G. B. Giannakis.  Distributed clustering using wireless sensor networks.  IEEE Journal of Selected Topics in Signal Processing, 5(4):707–724, 2011.
  * [3] S. Sun, R. Huang, and Y. Gao.  Network-scale traffic modeling and forecasting with graphical lasso and neural networks.  Journal of Transportation Engineering, 138(11):1358–1367, 2012\.
  * [4] C. Yang, M. Bakich, and E. Blasch.  Nonlinear constrained tracking of targets on roads.  In 2005 7th International Conference on Information Fusion, volume 1, pages 8 pp.–, 2005.
  * [5] M. Doostmohammadian and U. A. Khan.  Graph-theoretic distributed inference in social networks.  IEEE Journal of Selected Topics in Signal Processing, 8(4):613–623, 2014.
  * [6] L. Schwiebert, S. KS. Gupta, and J. Weinmann.  Research challenges in wireless networks of biomedical sensors.  In Proceedings of the 7th annual international conference on Mobile computing and networking, pages 151–165, 2001.
  * [7] E. Pavez, H. E. Egilmez, and A. Ortega.  Learning graphs with monotone topology properties and multiple connected components.  IEEE Transactions on Signal Processing, 66(9):2399–2413, 2018.
  * [8] S. Hassan-Moghaddam, N. K. Dhingra, and M. R. Jovanović.  Topology identification of undirected consensus networks via sparse inverse covariance estimation.  In IEEE 55th Conference on Decision and Control (CDC), pages 4624–4629, 2016.
  * [9] D. Materassi and M. V. Salapaka.  On the problem of reconstructing an unknown topology via locality properties of the Wiener filter.  IEEE Transactions on Automatic Control, 57(7):1765–1777, 2012.
  * [10] M. Ayazoglu, M. Sznaier, and N. Ozay.  Blind identification of sparse dynamic networks and applications.  In 50th IEEE Conference on Decision and Control and European Control Conference, pages 2944–2950, 2011.
  * [11] S. Sojoudi.  Equivalence of graphical lasso and thresholding for sparse graphs.  The Journal of Machine Learning Research, 17(1):3943–3963, 2016\.
  * [12] M. Yuan and Y. Lin.  Model selection and estimation in the Gaussian graphical model.  Biometrika, 94(1):19–35, 2007.
  * [13] J. Friedman, T. Hastie, and R. Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9(3):432–441, 2008.
  * [14] O. Banerjee, L. E. Ghaoui, and A. d’Aspremont.  Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data.  Journal of Machine Learning Research, 9(Mar):485–516, 2008.
  * [15] B. Rolfs, B. Rajaratnam, D. Guillot, I. Wong, and A. Maleki.  Iterative thresholding algorithm for sparse inverse covariance estimation.  In Advances in Neural Information Processing Systems, pages 1574–1582, 2012.
  * [16] C.-J. Hsieh, I. S. Dhillon, P. K. Ravikumar, and M. A. Sustik.  Sparse inverse covariance matrix estimation using quadratic approximation.  In Advances in Neural Information Processing Systems, pages 2330–2338, 2011.
  * [17] K. Scheinberg, S. Ma, and D. Goldfarb.  Sparse inverse covariance selection via alternating linearization methods.  In Advances in Neural Information Processing Systems, pages 2101–2109, 2010.
  * [18] O. Dalal and B. Rajaratnam.  Sparse Gaussian graphical model estimation via alternating minimization.  Biometrika, 104(2):379–395, 2017.
  * [19] T. Yao and S. Sundaram.  Online estimation of sparse inverse covariances.  In American Control Conference (ACC), pages 1935–1940, 2021.
  * [20] M. Pirani and S. Sundaram.  On the smallest eigenvalue of grounded Laplacian matrices.  IEEE Transactions on Automatic Control, 61(2):509–514, 2016.
  * [21] M. Pirani and S. Sundaram.  Spectral properties of the grounded Laplacian matrix with applications to consensus in the presence of stubborn agents.  In American Control Conference, pages 2160–2165, 2014.
  * [22] J. Ghaderi and R. Srikant.  Opinion dynamics in social networks: A local interaction game with stubborn agents.  In American control conference, pages 1982–1987. IEEE, 2013.
  * [23] K. Mohan, P. London, M. Fazel, D. Witten, and S. Lee.  Node-based learning of multiple Gaussian graphical models.  The Journal of Machine Learning Research, 15(1):445–488, 2014.
  * [24] A. Nedic, A. Ozdaglar, and P. A. Parrilo.  Constrained consensus and optimization in multi-agent networks.  IEEE Transactions on Automatic Control, 55(4):922–938, 2010.

###  -A Proof of Theorem 1

To prove Theorem 1, we first provide the following lemma showing that the
eigenvalues of the estimates are bounded.

###### Lemma 2

Let $t_{0}$ be a time such that the sample covariance matrix $S_{i,t_{0}}^{w}$
is nonsingular $\forall i\in\mathcal{V}_{a}$ and $\forall
w\in\\{1,2,\ldots,W\\}$. Define
$a=\min\\{\lambda_{\min}(S^{*}),\inf_{i,w,t\geq
t_{0}}\lambda_{\min}(S_{i,t}^{w})\\}-p\lambda$ and
$b=\max\\{\lambda_{\max}(S^{*}),\sup_{i,w,t\geq
t_{0}}\lambda_{\max}(S_{i,t}^{w})\\}+p\lambda$. Then $\forall t\geq
t_{0},\forall k\in\\{1,2,\ldots,K\\}$, and
$\lambda\in(0,\frac{1}{p}\min\\{\lambda_{\min}(S^{*}),\inf_{i,w,t\geq
t_{0}}\lambda_{\min}(S_{i,t}^{w})\\})$, the iterates $\Gamma_{i,t}^{k}$ of
Algorithm 1 satisfy $0\prec aI_{p}\preceq\Gamma_{i,t}^{k}\preceq bI_{p}$ and
$aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$.

###### Proof:

The proof is similar to [19], where we simply consider the infimum and
supremum of the eigenvalues of $S_{i,t}^{w}$, $\forall
i\in\mathcal{V}_{a},\forall w\in\\{1,2,\ldots,W\\}$, and $\forall t\geq
t_{0}$, for $a$ and $b$ respectively. ∎

The following lemma provides a bound of $\|\Gamma_{i,t}^{k}-\Gamma^{*}\|$
using the error from the previous iteration.

###### Lemma 3

Assume $\forall i\in\mathcal{V}_{a}$, $\forall t\geq t_{0}$, and $\forall
k\in\\{1,2,\ldots,K\\}$, iterates $\Gamma_{i,t}^{k}$ satisfy
$aI_{p}\preceq\Gamma_{i,t}^{k}\preceq bI_{p}$ and $\Gamma^{*}$ satisfies
$aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$ for some fixed constants $0<a<b$. At
the $k$-th iteration of time-step $t$, we have the following bound for the
dual variable:

$\|\Gamma_{i,t}^{k}-\Gamma^{*}\|_{F}\leq\beta_{i,t}^{k-1}\|\Gamma_{i,t}^{k-1}-\Gamma^{*}\|_{F}+2\|S_{i,t}^{W}-S^{*}\|_{F},$

where
$\beta_{i,t}^{k-1}=\max\left\\{|1-\frac{\zeta_{i,t}^{k-1}}{a^{2}}|,|1-\frac{\zeta_{i,t}^{k-1}}{b^{2}}|\right\\}$.

###### Proof:

Assume $aI_{p}\preceq\Gamma_{t}^{k}\preceq bI_{p},\forall t\geq t_{0}$,
$\forall k\in\\{1,2,\dots,K\\}$ and $aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$.
Using the dual update (10), matrix norm triangle inequality, and the non-
expansive property of the clip function $\mathcal{C}_{\lambda}$[18], we have:

$\displaystyle\|\Gamma_{i,t}^{k}-\Gamma^{*}\|_{F}=\|\mathcal{C}_{\lambda}(\Gamma_{i,t}^{k-1}+\zeta_{i,t}^{k-1}(\Gamma_{i,t}^{k-1})^{-1}-S_{i,t}^{W})+S_{i,t}^{W}$
$\displaystyle\qquad-\mathcal{C}_{\lambda}(\Gamma^{*}+\zeta_{i,t}^{k-1}(\Gamma^{*})^{-1}-S^{*})-S^{*}\|_{F}$
$\displaystyle\leq\|\mathcal{C}_{\lambda}(\Gamma_{i,t}^{k-1}+\zeta_{i,t}^{k-1}(\Gamma_{i,t}^{k-1})^{-1}-S_{i,t}^{W})$
$\displaystyle\qquad-\mathcal{C}_{\lambda}(\Gamma^{*}+\zeta_{i,t}^{k-1}(\Gamma^{*})^{-1}-S^{*})\|_{F}+\|S_{i,t}^{W}-S^{*}\|_{F}.$
$\displaystyle\leq\|(\Gamma_{i,t}^{k-1}+\zeta_{i,t}^{k-1}(\Gamma_{i,t}^{k-1})^{-1}-S_{i,t}^{W})$
$\displaystyle\qquad-(\Gamma^{*}+\zeta_{i,t}^{k-1}(\Gamma^{*})^{-1}-S^{*})\|_{F}+\|S_{i,t}^{W}-S^{*}\|_{F}$
$\displaystyle\leq\|(\Gamma_{i,t}^{k-1}+\zeta_{i,t}^{k-1}(\Gamma_{i,t}^{k-1})^{-1})-(\Gamma^{*}+\zeta_{i,t}^{k-1}(\Gamma^{*})^{-1})\|_{F}$
$\displaystyle\qquad+2\|S_{i,t}^{W}-S^{*}\|_{F}$
$\displaystyle\leq\max\left\\{\left|1-\frac{\zeta_{i,t}^{k-1}}{b^{2}}\right|,\left|1-\frac{\zeta_{i,t}^{k-1}}{a^{2}}\right|\right\\}\|\Gamma_{i,t}^{k-1}-\Gamma^{*}\|_{F}$
$\displaystyle\qquad+2\|S_{i,t}^{W}-S^{*}\|_{F},$

where the last inequality follows from Lemma in [15]. ∎

#### Proof of Theorem 1

Assume $\forall i\in\mathcal{V}_{a}$, $\forall t\geq t_{0}$, and $\forall
k\in\\{1,2,\dots,K\\}$, we can select a constant step size
$\zeta=\zeta_{i,t}^{k}$, where $0\prec aI_{p}\preceq\Gamma_{t}^{k}\preceq
bI_{p}$ and $aI_{p}\preceq\Gamma^{*}\preceq bI_{p}$. For
$t\in\\{t_{0}+1,t_{0}+2,\ldots\\}$, using the results from Lemma 3,

$\displaystyle\|\Gamma_{i,t_{0}+1}^{1}-\Gamma^{*}\|_{F}\leq\beta\|\Gamma_{i,t_{0}}^{K}-\Gamma^{*}\|_{F}+2\|S_{i,t_{0}+1}^{W}-S^{*}\|_{F}$
$\displaystyle\|\Gamma_{i,t_{0}+1}^{2}-\Gamma^{*}\|_{F}\leq\beta^{2}\|\Gamma_{i,t_{0}}^{K}-\Gamma^{*}\|_{F}+2\beta\|S_{i,t_{0}+1}^{W}-S^{*}\|_{F}$
$\displaystyle\qquad\qquad\qquad+2\|S_{i,t_{0}+1}^{W}-S^{*}\|_{F}.$

After $K$ iterations,

$\|\Gamma_{i,t_{0}+1}^{K}-\Gamma^{*}\|_{F}\leq\beta^{K}\|\Gamma_{i,t_{0}}^{K}-\Gamma^{*}\|_{F}\\\
+2\sum_{m=1}^{K}\beta^{K-m}\|S_{i,t_{0}+1}^{W}-S^{*}\|_{F}.$

At time $t$, at the end of $K$ iterations, we arrive at

$\displaystyle\|\Gamma_{i,t}^{K}-\Gamma^{*}\|_{F}\leq\beta^{K(t-t_{0})}\|\Gamma_{i,t_{0}}^{K}-\Gamma^{*}\|_{F}$
$\displaystyle\qquad\qquad+2\sum_{l=t_{0}+1}^{t}\beta^{K(t-l)}\sum_{m=1}^{K}\beta^{K-m}\|S_{i,l}^{W}-S^{*}\|_{F}.$

###  -B Proof of Lemma 1

Consider any entry $(l,m)$ of $S_{t}$. Under Def. 2, there is at least one
leader agent observing $S_{t}(l,m)$. If $i$ observes $S_{t}(l,m)$, then
$S_{i,t}^{w}(l,m)=S_{t}(l,m)$ $\forall i\in L$ and $\forall
w\in\\{1,2,\ldots,W\\}$. Since graph $\mathcal{G}_{a}$ is connected, the
values of every follower agent will converge to the same value as the leader
agents, i.e., $\lim_{w\to\infty}S_{i,t}^{w}(l,m)=S_{t}(l,m),\forall i\in F,$
geometrically with a rate at least equal to the largest eigenvalue of
$P_{FF,lm}$ [21, 22].

###  -C Proof of Theorem 2

To prove Theorem 2, we introduce the following lemma.

######  Lemma 4 ([24])

Let $0<\alpha<1$ and let $\\{\gamma_{t}\\}$ be a positive scalar sequence.
Assume that $\lim_{t\to\infty}\gamma_{t}=0$. Then

$\lim_{t\to\infty}\sum_{l=0}^{t}\alpha^{t-l}\gamma_{l}=0.$

#### Proof of Theorem 2

First, we characterize an error bound for $S_{i,t}^{W}$. Using triangle
inequality and applying Lemma 1, we have

$\displaystyle\|S_{i,1}^{W}-S^{*}\|_{F}$
$\displaystyle\leq\|S_{i,1}^{W}-S_{1}\|_{F}+\|S_{1}-S^{*}\|_{F}$
$\displaystyle\leq c\sigma^{W}\|S_{i,1}^{0}-S_{1}\|_{F}+\|S_{1}-S^{*}\|_{F}.$

Define $e_{t}=S_{t}-S_{t-1},\forall t\geq 1$. Recall that for initialization
$S_{i,0}^{w}=\textbf{0},\forall i,w$ and $S_{0}=\textbf{0}$. As a result,

$\|S_{i,1}^{W}-S_{1}\|_{F}\leq
c\sigma^{W}\|S_{0}+e_{1}\|_{F}=c\sigma^{W}\|e_{1}\|_{F},$

and

$\|S_{i,1}^{W}-S^{*}\|_{F}\leq c\sigma^{W}\|e_{1}\|_{F}+\|S_{1}-S^{*}\|_{F}.$

Similarly,

$\displaystyle\|S_{i,2}^{W}$ $\displaystyle-S^{*}\|_{F}\leq
c\sigma^{W}\|S_{i,2}^{0}-S_{2}\|_{F}+\|S_{2}-S^{*}\|_{F}$
$\displaystyle=c\sigma^{W}\|S_{i,1}^{W}-(S_{1}+e_{2})\|_{F}+\|S_{2}-S^{*}\|_{F}$
$\displaystyle\leq
c\sigma^{W}\|S_{i,1}^{W}-S_{1}\|_{F}+c\sigma^{W}\|e_{2}\|_{F}+\|S_{2}-S^{*}\|_{F}$
$\displaystyle\leq(c\sigma^{W})^{2}\|e_{1}\|_{F}+c\sigma^{W}\|e_{2}\|_{F}+\|S_{2}-S^{*}\|_{F}.$

At time step $t$, after $W$ rounds of communication, the difference between
the sample covariance estimation of agent $i\in\mathcal{V}_{a}$ and the ground
truth sample covariance $S^{*}$ is bounded by

$\|S_{i,t}^{W}-S^{*}\|_{F}\leq\|S_{t}-S^{*}\|_{F}+\sum_{j=1}^{t}(c\sigma^{W})^{t+1-j}\|e_{j}\|_{F}.$
(18)

From the law of large numbers, as $t\to\infty$, $S_{t}\to S^{*}$ almost surely
and $e_{t}\to 0$ almost surely. Thus, as $t\to\infty$, the first term of (18)
goes to zero almost surely. The second term is an instance of Lemma 4, where
$\alpha=c\sigma^{W}<1$ and $\gamma_{t}=\|e_{t}\|_{F}\to 0$ almost surely.

Figure 3: Demonstration of sample covariance convergence.

In Fig. 3, we use the experiment setup as in Sec.VI, let $W=1$ , and
demonstrate the convergence of sample covariance matrix of agents, evaluated
by $\|S_{i,t}^{W}-S^{*}\|_{F}$.

###  -D Proof of Corollary 1

From (14) in Theorem 1, if $\forall t\geq t_{0}$ and $\forall
k\in\\{1,2,\ldots,K\\}$, $\zeta_{i,t}^{k}=\zeta<a^{2}$, then $0<\beta<1$ . The
first term $\beta^{K(t-t_{0})}\|\Gamma_{i,t_{0}}^{K}-\Gamma^{*}\|_{F}$
converges to zero exponentially.

The sum $\sum_{l=t_{0}+1}^{t}\beta^{K(t-l)}\|S_{i,l}^{W}-S^{*}\|_{F}$ of the
second term is another instance of Lemma 4, where $\alpha=\beta^{K}<1$ and
$\gamma_{t}=\|S_{i,t}^{W}-S^{*}\|_{F}\to 0$ almost surely as $t\to\infty$ by
Theorem 2.

###  -E Proof of Corollary 2

We start with the following result on the convergence of the sample covariance
matrix.

######  Lemma 5 ([19])

Consider a sequence of independent and identically distributed $p$-dimensional
random vectors $\\{X_{1},X_{2},\ldots\\}$, where each $X_{t}$,
$t\in\mathbb{N}$, is drawn from $N(0,S^{*})$. For all $t\in\mathbb{N}$, let
$S_{t}=\frac{1}{t}\sum_{j=1}^{t}X_{j}X_{j}^{T}$ be the sample covariance
matrix for the data received up to time $t$. Then, for all
$\Delta\in(0,\frac{1}{2})$, there exists a set of sample paths of measure 1,
such that for each sample path in that set, there exists a finite time
$\bar{t}$, such that for all $t\geq\bar{t}$, the sample covariance matrix
satisfies the concentration inequality:

$\|S_{t}-S^{*}\|_{F}\leq
40p\max_{j}(S^{*}(j,j))\left(\frac{1}{t}\right)^{\frac{1}{2}-\Delta}.$

The above result, together with (18), leads directly to the following
characterization of the convergence rate of the sample covariance matrices
$S_{i,t}^{w}$ of each agent $i\in\mathcal{V}_{a}$.

###### Lemma 6

Consider a sequence of independent and identically distributed $p$-dimensional
random vectors $\\{X_{1},X_{2},\ldots\\}$, where each $X_{t}$,
$t\in\mathbb{N}$, is drawn from $N(0,S^{*})$. For all $t\in\mathbb{N}$, let
$S_{t}=\frac{1}{t}\sum_{j=1}^{t}X_{j}X_{j}^{T}$ be the sample covariance
matrix for the data received up to time $t$. Then, for all
$\Delta\in(0,\frac{1}{2})$, there exists a set of sample paths of measure 1,
such that for each sample path in that set, there exists a finite time
$\bar{t}$, such that for all $t\geq\bar{t}+1$, the sample covariance matrix
satisfies the concentration inequality:

$\|S_{i,t}^{W}-S^{*}\|_{F}\leq(c\sigma^{W})^{t-\bar{t}}\|S_{i,\bar{t}}^{W}-S_{\bar{t}}^{*}\|_{F}+40p\max_{j}(S^{*}(j,j))\\\
\cdot\Bigg{(}\left(\frac{1}{t}\right)^{\frac{1}{2}-\Delta}+2\sum_{j=\bar{t}+1}^{t}(c\sigma^{W})^{t+1-j}\left(\frac{1}{j-1}\right)^{\frac{1}{2}-\Delta}\Bigg{)}.$
(19)

###### Proof:

Using triangle inequality and applying Lemma 5, $\forall t\geq\bar{t}+1$, we
have

$\displaystyle\|e_{t}\|_{F}$
$\displaystyle=\|S_{t}-S_{t-1}\|_{F}\leq\|S_{t}-S^{*}\|_{F}+\|S_{t-1}-S^{*}\|_{F}$
$\displaystyle\leq
40p\max_{j}(S^{*}(j,j))\left(\left(\frac{1}{t}\right)^{\frac{1}{2}-\Delta}+\left(\frac{1}{t-1}\right)^{\frac{1}{2}-\Delta}\right)$
$\displaystyle\leq
80p\max_{j}(S^{*}(j,j))\left(\frac{1}{t-1}\right)^{\frac{1}{2}-\Delta}.$

To analyze the asymptotic behavior of the algorithm, we rewrite the update in
(18). Given the error $\|S_{i,\bar{t}}^{W}-S_{\bar{t}}^{*}\|_{F}$ at time
$\bar{t}$, we write out the error between $S_{i,\bar{t}+1}^{W}$ and $S^{*}$:

$\displaystyle\|S_{i,\bar{t}+1}^{W}$ $\displaystyle-S^{*}\|_{F}\leq
c\sigma^{W}\|S_{i,\bar{t}+1}^{0}-S_{\bar{t}+1}\|_{F}+\|S_{\bar{t}+1}-S^{*}\|_{F}$
$\displaystyle=c\sigma^{W}\|S_{i,\bar{t}}^{W}-(S_{\bar{t}}+e_{\bar{t}+1})\|_{F}+\|S_{\bar{t}+1}-S^{*}\|_{F}$
$\displaystyle\leq
c\sigma^{W}\left(\|S_{i,\bar{t}}^{W}-S_{\bar{t}}\|_{F}+\|e_{\bar{t}+1}\|_{F}\right)+\|S_{\bar{t}+1}-S^{*}\|_{F}.$

For all $t\geq\bar{t}+1$, we obtain

$\|S_{i,t}^{W}-S^{*}\|_{F}\leq(c\sigma^{W})^{t-\bar{t}}\|S_{i,\bar{t}}^{W}-S_{\bar{t}}^{*}\|_{F}\\\
+\sum_{j=\bar{t}+1}^{t}(c\sigma^{W})^{t+1-j}\|e_{j}\|_{F}+\|S_{t}-S^{*}\|_{F}.$
(20)

Substituting the corresponding terms in (20) with inequalities for
$\|e_{t}\|_{F}$ and in Lemma 5 , we obtain (19). ∎

The concentration inequality of $\|S_{i,t}^{W}-S^{*}\|_{F}$ in Lemma 6 and the
error bound $\|\Gamma_{i,t}^{K}-\Gamma^{*}\|_{F}$ in Theorem 1 immediately
lead to the convergence rate of the dual variable.