Distributed Active State Estimation with User-Specified Accuracy

TL;DR

This paper presents a distributed method for mobile sensor networks to estimate hidden states with user-defined accuracy, using an innovative optimization approach and an information consensus filter to guide sensor positioning.

Contribution

It introduces a novel distributed random approximate projections method for LMI-constrained optimization in active state estimation tasks.

Findings

Robust convergence to desired estimation accuracy.

Effective localization of sparse landmarks.

Emergent behaviors in robotic sensor networks.

Abstract

In this paper, we address the problem of controlling a network of mobile sensors so that a set of hidden states are estimated up to a user-specified accuracy. The sensors take measurements and fuse them online using an Information Consensus Filter (ICF). At the same time, the local estimates guide the sensors to their next best configuration. This leads to an LMI-constrained optimization problem that we solve by means of a new distributed random approximate projections method. The new method is robust to the state disagreement errors that exist among the robots as the ICF fuses the collected measurements. Assuming that the noise corrupting the measurements is zero-mean and Gaussian and that the robots are self localized in the environment, the integrated system converges to the next best positions from where new observations will be taken. This process is repeated with the robots taking a sequence of observations until the hidden states are estimated up to the desired user-specified accuracy. We present simulations of sparse landmark localization, where the robotic team achieves the desired estimation tolerances while exhibiting interesting emergent behavior.

Equations

$${\mathbf{y}}_{i,s}(t)$$

$$\left(\mathcal{D}(t-1),\mathcal{O}(t-1)\right)\stackrel{\mathrm{ICF}}{⟼}\mathcal{D}(t).$$

$$S_i(t+1)=\stackrel{\text{from D(t)}}{\overbrace{S_i(t)}}+\underset{\text{from O(t)}}{\underbrace{{\sum }_{s\in \mathcal{S}}Q({\mathbf{r}}_s(t-1)+{\mathbf{u}}_s,\stackrel{\text{from D(t)}}{\overbrace{{\hat{\mathbf{x}}}_i(t)}})}}.$$

$$\underset{(\mathbf{\gamma },\mathbf{u})\in \Gamma \times \mathcal{U}}{\max }{\sum }_{i\in \mathcal{I}}{\gamma }_i$$

$$\text{s.t}{\gamma }_{i}I⪯S_i+{\sum }_{s\in \mathcal{S}}Q({\mathbf{r}}_s+{\mathbf{u}}_s,{\hat{\mathbf{x}}}_i)\forall i\in \mathcal{I},$$

$$h:{\mathbb{R}}^{m+qn}\times \left({\mathbb{R}}^{pm}\times {\left({\mathrm{Sym}}_{+}(p,\mathbb{R})\right)}^m\right)\times \mathcal{I}\to \mathrm{Sym}(p,\mathbb{R}),$$

$$h\left((\mathbf{\gamma },\mathbf{u});\mathcal{D},i\right)$$

$$-$$

$$\underset{{\left\{{\mathbf{z}}_s\right\}}_{s\in \mathcal{S}}\subset {\mathcal{X}}_0}{\min }$$

$$\text{s.t}h\left({\mathbf{z}}_s;\mathcal{D},i\right)$$

$$\displaystyle\!\!\!{\mathcal{D}}_{s,k}\leftarrow\Big{(}\big{(}(n\Xi_{1,s,k})^{-1}\boldsymbol{\xi}_{1,s,k},\dots,(n\Xi_{m,s,k})^{-1}\boldsymbol{\xi}_{m,s,k}\big{)},$$

$$\displaystyle\left(n\Xi_{1,s,k},\dots,n\Xi_{m,s,k}\right)\Big{)}$$

$$\Vert {\mathcal{D}}_{s,k}\Vert \le \underset{v\in \mathcal{S}}{\max }\Vert \left({\mathcal{D}}_{v,0}\right)\Vert ,$$

$${\mathsf{\Pi }}_{+}X=E\text{diag}\left[({\lambda }_1{)}_{+},\dots ,({\lambda }_p{)}_{+}\right]E^{\top },$$

$$h_{+}\left(\mathbf{z};\mathcal{D},i\right)$$

$$\nabla_{j}h_{+}({\mathbf{z}};{\mathcal{D}},i)=\psi_{ij}\Big{(}\mathsf{\Pi}_{+}\big{(}h\left({\mathbf{z}};{\mathcal{D}},i\right)\big{)},{\mathbf{z}},{\mathcal{D}}\Big{)},$$

$${\psi }_{ij}(X,\mathbf{z},\mathcal{D})\triangleq \{\begin{array}{ll}\frac{\mathrm{Tr}\left({\nabla }_{j}h(\mathbf{z},\mathcal{D},i)X\right)}{{\Vert X\Vert }_F}& \mathrm{if}{\Vert X\Vert }_F>0\\ d& \mathrm{else}\end{array}$$

$${\beta }_{s,k}=\frac{h_{+}({\mathbf{v}}_{s,k};{\mathcal{D}}_{s,k},{\omega }_{s,k})}{{\Vert h_{+}^{\prime }({\mathbf{v}}_{s,k};{\mathcal{D}}_{s,k},{\omega }_{s,k})\Vert }^2}$$

$${\nu }_{s,k}=h_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},{\omega }_{s,k}\right)-h_{+}\left({\mathbf{z}}_{s,k};\mathcal{D},{\omega }_{s,k}\right).$$

$${\mathbf{\delta }}_{s,k}=h_{+}^{\prime }\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},{\omega }_{s,k}\right)-h_{+}^{\prime }\left({\mathbf{z}}_{s,k};\mathcal{D},{\omega }_{s,k}\right).$$

$$\Vert {\mathbf{z}}_1-{\mathbf{z}}_2\Vert \le C_{\mathbf{z}}.$$

$$\Vert f^{\prime }(\mathbf{z})\Vert \le L_f,$$

$$|f({\mathbf{z}}_1)-f({\mathbf{z}}_2)|\le L_f\Vert {\mathbf{z}}_1-{\mathbf{z}}_2\Vert .$$

$$\Vert h_{+}^{\prime }(\mathbf{z};\mathcal{D},i)\Vert \le L_h.$$

$$\Vert h_{+}^{\prime }(\mathbf{z};{\mathcal{D}}_{s,k},i)\Vert =\Vert h_{+}^{\prime }(\mathbf{z};\mathcal{D},i)+{\mathbf{\delta }}_{s,k}\Vert \le L_h+D,$$

$$|h_{+}({\mathbf{z}}_1;{\mathcal{D}}_{s,k},i)-h_{+}({\mathbf{z}}_2;{\mathcal{D}}_{s,k},i)|\le (L_h+D)\Vert {\mathbf{z}}_1-{\mathbf{z}}_2\Vert .$$

$$|h_{+}(\mathbf{z};\mathcal{D},i)|\le C_h,$$

$$|h_{+}(\mathbf{z};{\mathcal{D}}_{s,k},i)|=|h_{+}(\mathbf{z};\mathcal{D},i)+{\nu }_{s,k}|\le C_h+N,$$

$$\Vert h_{+}^{\prime }(\mathbf{z};{\mathcal{D}}_{s,k},i)\Vert \ge c_h$$

$${\beta }_{s,k}=\frac{h_{+}(\mathbf{z};{\mathcal{D}}_{s,k},{\omega }_{s,k})}{\Vert h_{+}^{\prime }(\mathbf{z};{\mathcal{D}}_{s,k},{\omega }_{s,k}){\Vert }^2}\le \frac{C_h+N}{c_h^2}.$$

$${\sum }_{s\in \mathcal{S}}g({\mathbf{p}}_{s,k})\le {\sum }_{s\in \mathcal{S}}g({\mathbf{z}}_{s,k-1})$$

$$\displaystyle\sum_{s\in{\mathcal{S}}}g(\mathbf{p}_{s,k})\leq\leavevmode\nobreak\$$

$$\displaystyle=\leavevmode\nobreak\$$

$${\mathbf{\theta }}_{s,k}=\sum _{v\in \mathcal{S}}[W_k{]}_{sv}{\mathbf{\theta }}_{v,k-1}+{\mathbf{e}}_{s,k},\text{leavevmode}\forall s\in \mathcal{S},$$

Full text

I Introduction

In this paper, we control a robotic sensor network to estimate a collection of hidden states so that desired accuracy thresholds and confidence levels are satisfied. We require that estimation and control are completely distributed, i.e., belief states and control actions are decided upon locally through pairwise communication among robots in the network.

A common assumption in problems like the one discussed herein is that observations depend linearly on the state and are corrupted by Gaussian noise as well as that the sensors are self localized [1, 2, 3, 4, 5, 6, 7, 8]. Under these assumptions, a typical approach to estimate a set of hidden states is to use an Information Filter (IF). The role of the IF is to fuse new observations of the hidden states with a prior Gaussian distribution to create a minimum-variance posterior Gaussian distribution. Typically, the observations are subject to error covariance that is a function of, e.g., viewing position, angle of incidence, and the sensing modality. Several such models have been proposed in the distributed active sensing literature, see, e.g., [1, 7, 2, 3, 5, 6, 8], but regardless of the specific model used for the measurement covariance, the IF combines the data among the nodes in the network in an additive way, making the IF algorithm readily distributed; see, e.g., Information Consensus Filtering (ICF) [9].

The dependence of the measurement model on the sensor state has been widely used to obtain planning algorithms that, given a history of measurements, determine the next best set of observations of a hidden state. The objective is that this new set of observations optimizes an information-theoretic criterion of interest. In this paper, we choose to maximize the minimum eigenvalue of the posterior information matrix, effectively minimizing the directional variance of every hidden state. We express this requirement by reformulating the problem as a Linear Matrix Inequality (LMI) constrained optimization problem. Specifically, we introduce as many LMI constraints as the number of hidden states, that are coupled with respect to the sensors via an ICF that estimates the uncertainty of each hidden state using the sensor measurements. To solve this problem, we propose a new distributed optimization algorithm, which we call random approximate projections, that is robust to the state disagreement errors that exist among the robots as an ICF fuses the collected measurements. The ability to handle such errors is critical in combining distributed estimation and planning in one process.

The proposed distributed random approximate projections method falls in the class of so called consensus-based optimization algorithms, where the goal is to interleave a decentralized averaging step between local optimization steps so that the local estimates of all agents simultaneously obtain a network-wide consensus and optimality of the global problem. Existing algorithms require expensive optimization steps or projections on a complicated constraint set, e.g., sets of LMI constraints as in this paper, at every iteration. Instead, our method employs a subgradient step in the direction that minimizes the violation of randomly selected local constraints, significantly reducing computational cost compared to relevant methods that rely on projection operations. Assuming that each LMI constraint is selected with non-zero probability, we show that our proposed algorithm converges almost surely to the next best sensor positions from where a new set of observations minimizes the worst-case uncertainty of the hidden states. This process is repeated with the sensors taking a sequence of observations until all the hidden states are estimated up to the desired user-specified accuracy. For details on relevant distributed optimization methods, we refer the interested reader to [10] and references therein.

Related to our random projection technique is the Markov incremental random subgradient method [11]. The key difference is in the mechanism that is used to distribute the computations. In the Markov incremental algorithm, the agents sequentially update a single iterate sequence by passing the iterate to each other following a time inhomogeneous Markov chain. In our algorithm, each agent maintains a local iterate sequence and updates this by communicating with its neighbors.

In two recent papers [12, 13], a dual-decomposition method is proposed that can be used for problems similar to the distributed control problem considered here, with lower communication overhead and increased privacy. Generally, dual methods are well suited for network optimization problems, they come at increased local cost because of the need to compute the dual function at each node. The approaches described in [12, 13] additionally require a special problem structure; specifically, a decomposable objective and/or a star network topology.

Applications in the area of distributed sensor planning and estimation that are directly relevant to our work include SLAM [14], localization [4, 5, 6], coverage [7, 8], mobile target tracking [15, 3, 2, 1], classification [16], and inverse problems for PDE-constrained systems [17, 18]. Compared to state estimation, inverse problems for PDE-constrained systems address both state and source identification and are particularly difficult to solve.

A common characteristic of the applications discussed above is that the sensors must collectively reason over the possible posterior distributions resulting from

1 multiple hidden states being observed by a single sensor or

2 individual hidden states being observed by several sensors at once.

Without assuming Gaussian distributions, optimization over the possible posteriors can only be done approximately [16]. In this paper, and in much of the relevant literature [1, 2, 3, 4, 5, 6], it is assumed that the posteriors for the hidden states are Gaussian. Under this assumption, typical choices of an objective function to minimize in order to obtain the next best set of observations include the trace, the determinant, or the maximum eigenvalue of the posterior covariance matrix. Intuitively, the trace minimizes the sum of the uncertainties of the hidden states, the determinant the volume of the confidence ellipsoids of the hidden states, and the maximum eigenvalue the worst-case error. To our knowledge, in the context of the distributed active localization problem, existing works that control the maximum eigenvalue of the error covariance matrix either focus on sensor selection from a discrete set [19, 20, 21, 12, 13], are limited to two sensors [6], or do not factor in simultaneous measurements by other sensors when deciding where to go [22, 2, 4].

Methods that account for nonlinear dynamics and arbitrary error distributions [23, 15] are also relevant to the problem addressed herein. In these cases, updates to the hidden state pdfs may no longer be available in closed form, and so numerical methods are needed to evaluate the expectation integrals. To account for uncertainty in the sensor location and actions, decentralized Partially Observable Markov Decision Process (dec-POMDP) can be used to solve these problems. Dec-POMDPs have coupled value functions among the agents, and typically require a centralized planner [24]. At the same time, comparing the possible posterior distributions can be computationally expensive, so these methods do not usually scale well with the number of targets and agents. Although many of these works assume that the robots are self localized with respect to a common global reference frame [15, 16, 6, 1, 2, 3, 4], this is certainly not always the case; see, e.g., work on partially observable systems [24, 23] or the SLAM problem [14].

Contributions: In this paper we propose a distributed estimation and sensor planning method that can ensure desired estimation tolerances for large numbers of hidden states. To our knowledge, even though the distributed active sensing literature is well-developed, the ability to control worst-case estimation uncertainty in a decentralized way is new. This is possible by combining an ICF for distributed data fusion with an LMI-constrained optimization problem for sensor planning, for which we propose a new computationally inexpensive distributed method that is robust to the inexact pre-consensus filter data. Coordination among the sensors requires only weak connectivity of the communication graph and can, thus, handle occasional disruption of communication links.

The distributed planning algorithm itself is an extension of the authors’ previous work [25] and [26]. The difference is that the method developed in this paper is robust to state disagreement errors that appear as a result of the information filter. The ability to handle such errors is critical in integrating distributed estimation and planning. To the best of our knowledge, there are no distributed optimization methods that can handle efficiently LMI constraints under inexact data. It is also worth mentioning that a wide variety of convex constraints arising in control, such as quadratic inequalities, inequalities involving matrix norms, Lyapunov inequalities and quadratic matrix inequalities, can be cast as LMIs. Therefore, our developed algorithm can be used for other, potentially unrelated, control problems involving inexact data.

The rest of the paper is organized as follows. Section II formulates the problem under consideration. Then, Section III develops the proposed distributed method to maximize the minimum eigenvalue of the information matrix. Section IV states assumptions and the main convergence result of our proposed algorithm. Section V shows that the proposed algorithm converges to the optimal point. Section VI provides simulations, and Section VII concludes the paper.

II Problem Formulation

Consider the problem of estimating $m$ stationary hidden states $\left\{{\mathbf{x}}_{i}\in{\mathbb{R}}^{p}\right\}_{i\in{\mathcal{I}}}$ using noisy observations from $n$ mobile sensors, where $p$ is the dimension of the hidden states and ${\mathcal{I}}=\left\{1,\dots,m\right\}$. We assume that the hidden states are sparse so that there are no correlations between them. Denote the locations of the mobile sensors at time $t$ by $\left\{{\mathbf{r}}_{s}(t)\in{\mathbb{R}}^{q}\right\}_{s\in{\mathcal{S}}},$ where $q$ is the dimension of the sensors’ configuration space and ${\mathcal{S}}=\left\{1,\dots,n\right\}.$ The observation of state $i$ by sensor $s$ at time $t$ is given by

[TABLE]

We assume that the measurement errors are Normally distributed so that \boldsymbol{\zeta}_{i,s}(t)\sim N\big{(}{\mathbf{0}},Q({\mathbf{r}}_{s}(t),{\mathbf{x}}_{i})\big{)}, where $Q\colon{\mathbb{R}}^{q+p}\to{\rm Sym}_{+}(p,{\mathbb{R}})$ denotes the measurement precision matrix, or information matrix (not the covariance), and ${\rm Sym}_{+}(p,{\mathbb{R}})$ is the set of $p\times p$ symmetric positive semidefinite real matrices. We also assume that if signal ${\mathbf{x}}_{i}$ is out of range of sensor $s$, then the function $Q$ returns the information matrix ${\mathbf{0}}_{p\times p},$ corresponding to infinite variance. We provide an example of the measurement information function $Q$ in Section VI. Denote the set of all observations at time $t$ by ${\mathcal{O}}(t)\triangleq\left\{\left({\mathbf{y}}_{i,s}(t),Q({\mathbf{r}}_{s}(t),{\mathbf{x}}_{i})\right)\right\}_{i\in{\mathcal{I}},s\in{\mathcal{S}}}$

Moreover, let ${\hat{\mathbf{x}}}_{i}(t)$ denote the estimate of ${\mathbf{x}}_{i}$ at time $t$ and $S_{i}(t)$ denote the information matrix corresponding to the estimate ${\hat{\mathbf{x}}}_{i}(t).$ We assume that initial values for the pair $\left({\hat{\mathbf{x}}}_{i}(0),S_{i}(0)\right)$ are available a priori.

Let also ${\hat{\mathbf{x}}}(t)\in{\mathbb{R}}^{pm}$ denote the stack vector of all estimates ${\hat{\mathbf{x}}}_{i}(t)$ and $S(t)\in{\rm Sym}_{+}(p,{\mathbb{R}})\times\dots\times{\rm Sym}_{+}(p,{\mathbb{R}}):=\left({\rm Sym}_{+}(p,{\mathbb{R}})\right)^{m}$ denote the collection of all information matrices $S_{i}(t)$ at time $t$. Then, the global data at time $t$ can be captured by the set ${\mathcal{D}}(t)\triangleq({\hat{\mathbf{x}}}(t),S(t))\in{\mathbb{R}}^{pm}\times\left({\rm Sym}_{+}(p,{\mathbb{R}})\right)^{m}.$

Given prior data ${\mathcal{D}}(t-1)$ known by all sensors, and observations ${\mathcal{O}}(t-1)$, the current set of data ${\mathcal{D}}(t)$ can be determined in a distributed way using an Information Consensus Filter (ICF) [9] as111Explicit details on the mapping (2) are provided in Section III.

[TABLE]

Then, the goal is to control the next set of observations ${\mathcal{O}}(t)$ in order to reduce the future uncertainty in ${\mathcal{D}}(t+1)$. We can achieve this goal by controlling the mobile sensors and, therefore, the expected information $\left\{Q({\mathbf{r}}_{s}(t),{\mathbf{x}}_{i}(t)\right\}_{i\in{\mathcal{I}},s\in{\mathcal{S}}}$ associated with the observations $\left\{{\mathbf{y}}_{i,s}(t)\right\}_{i\in{\mathcal{I}},s\in{\mathcal{S}}}$. Specifically, let ${\mathbf{r}}_{s}(t)={\mathbf{r}}_{s}(t-1)+{\mathbf{u}}_{s}$ denote the motion model of sensor $s$, where ${\mathbf{u}}_{s}\in{\mathbb{R}}^{q}$ is a candidate control input used to drive sensor $s$ to a new position ${\mathbf{r}}_{s}(t)$. Then, our goal is to maximize the minimum eigenvalues of the information matrices in ${\mathcal{D}}(t+1)$ that, for each hidden state $i,$ are given by

[TABLE]

We assume that ${\mathbf{u}}_{s}$ lies in a compact set containing the origin. Denote this set by ${\mathcal{U}}_{s}.$ Note that there is a $\delta>0$ such that the closed ball of radius $\delta$ contains ${\mathcal{U}}_{s}.$ Define ${\mathcal{U}}\triangleq\prod_{s\in{\mathcal{S}}}{\mathcal{U}}_{s}$ and let ${\mathbf{u}}\in{\mathbb{R}}^{qn}$ denote the stack of all controllers ${\mathbf{u}}_{s}.$ Additionally, define the set $\Gamma\triangleq\prod_{i\in{\mathcal{I}}}[0,\tau_{i}]$, where $\tau_{i}$ represent the given user-specified estimation thresholds for each hidden state $i$. Define also $\boldsymbol{\gamma}=(\gamma_{1},\dots,\gamma_{m})\in{\mathbb{R}}^{m}.$ With the new notation, the problem we consider in this paper can be defined as

[TABLE]

where we have dropped dependence on the time $t$ for notational convenience. For the explicit references to time in (4b) see (3). The symbol $\preceq$ denotes an ordering on the negative semidefinite cone, i.e., $A\preceq B$ if and only if $A-B$ is negative semidefinite. The great challenge in solving (4) is that the constraints (4b) require knowledge of ${\mathcal{D}}(t)$, however, only ${\mathcal{D}}(t-1)$ and those parts of ${\mathcal{O}}(t-1)$ corresponding to local sensor measurements are available to any given sensor. In what follows, we propose an integrated approach that combines the ICF to determine the global data ${\mathcal{D}}(t)$ and the solution of the optimization (4) in one process. We show that our proposed method can handle pre-consensus disagreement errors due to the ICF and converges to the next best set of observations ${\mathcal{O}}(t)$, as desired.

III Concurrent Estimation and Sensor Planning

Problem (4) is a nonlinear semidefinite program. To simplify this problem, we replace the nonlinear function $Q({\mathbf{r}}_{s}+{\mathbf{u}}_{s},{\hat{\mathbf{x}}}_{i})$ in (4b) with its Taylor expansion about the input ${\mathbf{u}}_{s}=0$. Specifically, define the function

[TABLE]

where ${\rm Sym}(p,{\mathbb{R}})$ denotes the set of symmetric $p\times p$ dimensional matrices, by

[TABLE]

The notation $\nabla_{j}Q({\mathbf{r}}_{s},{\hat{\mathbf{x}}}_{i})$ is the partial derivative of the matrix-valued function $Q$ with respect to the $j$-th coordinate of ${\mathbb{R}}^{q+p},$ evaluated at $({\mathbf{r}}_{s},{\hat{\mathbf{x}}}_{i})$. Also, $[{\mathbf{u}}_{s}]_{j}$ denotes the $j$-th coordinate of the vector ${\mathbf{u}}_{s}$, and the semicolon notation in the arguments of $h$ separates the decision variables $(\boldsymbol{\gamma},{\mathbf{u}})$ from the parameters ${\mathcal{D}}$ and $i$. Note that the value of $h$ is negative semidefinite at a point $\left((\boldsymbol{\gamma},{\mathbf{u}});{\mathcal{D}},i\right)$ if and only if the linearized version of the constraint (4b) is satisfied for $\gamma_{i}$ using the data ${\mathcal{D}}=({\hat{\mathbf{x}}},S).$ The effect of the linear approximation is investigated in Section VI, Fig. 4.

To simplify notation, collect all decision variables in the vector ${\mathbf{z}}=\left(\boldsymbol{\gamma},{\mathbf{u}}\right)\in{\mathbb{R}}^{m+qn}$. Then, to decompose the linearized version of (4) over the set of sensors so that the problem can be solved in a distributed way, let ${\mathbf{z}}_{s}$ denote a copy of ${\mathbf{z}}$ that is local to sensor $s\in{\mathcal{S}}.$ Additionally, let $f({\mathbf{z}}_{s})=-\sum_{i\in{\mathcal{I}}}\gamma_{i,s}$ so that we may express the global objective as $\sum_{s\in{\mathcal{S}}}f({\mathbf{z}}_{s}).$ Then, problem (4) with linearized constraints can be expressed as

[TABLE]

where ${\mathcal{X}}_{0}=\Gamma\times{\mathcal{U}}$. In problem (6), the objective function (6a) is separable among the sensors and the constraints (6b) are linear in ${\mathbf{z}}_{s}$ and local to every sensor. In this form, problem (6) can be decomposed and solved in a distributed way. A necessary requirement is that, at the solution, the local variables ${\mathbf{z}}_{s}$ need to agree for all sensors $s,$ i.e., ${\mathbf{z}}_{s}={\mathbf{z}}_{v}$ for all $s,v\in{\mathcal{S}}$. As discussed in Section II, the main challenge in solving problem (6) is that the global data in ${\mathcal{D}}$ are not known to the sensors and need to be computed using the ICF concurrently with the optimization of the decision variables ${\mathbf{z}}_{s}$. The ICF is an iterative process by which every sensor updates a local copy of the global data until all sensors agree on the true values in ${\mathcal{D}}$. Therefore, before the ICF has converged, disagreement errors on the local copies of the data set ${\mathcal{D}}$ exist, which means that the gradients and objective function evaluations at every sensor node will not depend on the same data. In other words, gradients and objective function evaluations will be based on inexact data. An additional challenge is that the number of LMI constraints can grow very large with the number of sensors and hidden states. We address these challenges in the remainder of this section.

III-A Distributed Information Filtering

We solve problem (6) using an iterative approach. Let ${\mathcal{D}}_{s,k}$ denote the copy of the data set ${\mathcal{D}}$ that is local to sensor $s$ at iteration $k\in{\mathbb{N}}$, where ${\mathbb{N}}=\left\{1,2,3,\dots\right\}$ denotes the natural numbers. Additionally, let ${\mathbf{z}}_{s,k}$ denote the decision variables of the sensor $s$ at iteration $k$. The first step of the algorithm is for sensor $s$ to broadcast ${\mathbf{z}}_{s,k-1}$ and a summary of ${\mathcal{D}}_{s,k-1}.$ In Algorithm 1, the data summary for state $i$ is captured by the intermediate information matrix $\boldsymbol{\Xi}_{i,s,k}$ and information vector $\boldsymbol{\xi}_{i,s,k},$ defined in lines 2 and 3. We refer the reader to [9] for a detailed discussion of Information Consensus Filtering, but we note briefly here that if the sensors each compute $\boldsymbol{\Xi}_{i,s,0}$ and $\boldsymbol{\xi}_{i,s,0}$, and carry out the weighted averaging that is standard in all consensus algorithms, given in lines 7 and 8, then $n\boldsymbol{\Xi}_{i,s,k}$ and $\boldsymbol{\xi}_{i,s,k}$ converge geometrically to the consensus posterior information. In other words, under the ICF, we have that, for all $s\in{\mathcal{S}},$ $\left\{\left(n\Xi_{1,s,k},\dots,n\Xi_{m,s,k}\right)\right\}_{k\in{\mathbb{N}}}\to S$ for all $i\in{\mathcal{I}}$ and \left\{\big{(}(n\Xi_{1,s,k})^{-1}\boldsymbol{\xi}_{1,s,k},\dots,(n\Xi_{m,s,k})^{-1}\boldsymbol{\xi}_{m,s,k}\big{)}\right\}_{k\in{\mathbb{N}}}\to{\hat{\mathbf{x}}}. In particular, let ${\mathcal{G}}_{k}=({\mathcal{S}},{\mathcal{E}}_{k})$ represent a communication graph such that $(s,v)\in{\mathcal{E}}_{k}$ if and only if sensors $s$ and $v$ can communicate at iteration $k$. This defines the neighbor set, ${\mathcal{N}}_{s,k}\triangleq\left\{s\right\}\cup\left\{v\in{\mathcal{S}}\mid(v,s)\in{\mathcal{E}}_{k}\right\}.$ Sensor $s$ receives $\left\{{\mathbf{z}}_{v,k-1},\left\{\boldsymbol{\xi}_{i,v,k},\boldsymbol{\Xi}_{i,v,k}\right\}_{i\in{\mathcal{I}}}\mid v\in{\mathcal{N}}_{s,k}\right\}.$ With this new information, the sensor updates its own optimization variable and data summary via using weighted averaging. In particular, let $W_{k}$ be an $n\times n$ row stochastic matrix such that $[W_{k}]_{sv}>0\iff v\in{\mathcal{N}}_{s}.$ The exact assumptions for $W_{k}$ will be introduced later, when they are necessary for the proof of our main result. The consensus steps are depicted in lines 6, 7, and 8 of Algorithm 1. Then, the agent can compute ${\mathcal{D}}_{s,k}$ as shown in line 9, concluding the IF portion of iteration $k.$ Note that, given initial conditions $\left\{{\mathcal{D}}_{s,0}\right\}_{s\in{\mathcal{S}}},$ the data for any sensor $s$ at time $k$ is bounded by the initial conditions as

[TABLE]

where ${\left\|\cdot\right\|}$ is any norm in the space of data ${\mathbb{R}}^{pm}\times\left({\rm Sym}_{++}(p,{\mathbb{R}})\right)^{m}.$ This implies that the set of possible data is bounded with respect to any norm in the vector space ${\mathbb{R}}^{pm}\times\left({\rm Sym}(p,{\mathbb{R}})\right)^{m}.$ Note that, by the same arguments as can be found in [9], under any norm, all sensors’ data converge to the same consensus value, i.e., ${\left\|{\mathcal{D}}_{s,k}-{\mathcal{D}}_{v,k}\right\|}\to 0$ and ${\left\|{\mathcal{D}}_{s,k}-{\mathcal{D}}\right\|}\to 0$ for all $s,v\in{\mathcal{S}}.$

III-B Random Approximate Projections

The optimization part of the algorithm begins in line 10 by computing a negative subgradient of the $f$ at the local iterate, which we denote by $-f^{\prime}({\mathbf{z}}_{s,k})$, taking a gradient descent step, and projecting the iterate back to the simple constraint set ${\mathcal{X}}_{0}$. The positive step size $\alpha_{k}$ is diminishing, non-summable and square-summable, i.e, $\alpha_{k}\to 0$, $\sum_{k\in{\mathbb{N}}}\alpha_{k}=\infty$ and $\sum_{k\in{\mathbb{N}}}\alpha_{k}^{2}<\infty$. The last step of our method is not a projection to the feasible set, but instead a subgradient step in the direction that minimizes the violation of a single LMI constraint. Essentially, each node selects one of the LMI’s randomly, measures the violation of that constraint, and takes an additional gradient descent step with step size $\beta_{s,k}$ in order to minimize this violation. We call this process random approximate projection.

Specifically, line 11 randomly assigns $\omega_{s,k}\in{\mathcal{I}}$, effectively choosing to activate $h\left(\cdot;\cdot,\omega_{s,k}\right)$. The idea behind line 13 is to minimize a scalar metric that measures the amount of violation of $h\left(\cdot;\cdot,\omega_{s,k}\right)$, which we denote by $h_{+}\left(\cdot;\cdot,\omega_{s,k}\right),$ while remaining feasible. Let $\mathsf{\Pi}_{+}\colon{\rm Sym}(p,{\mathbb{R}})\to{\rm Sym}_{+}(p,{\mathbb{R}})$ denote the projection operator onto the set of positive semidefinite matrices, given by

[TABLE]

where $(\cdot)_{+}\triangleq\max\left\{\cdot,0\right\}$ and $\lambda_{1},\dots,\lambda_{p}$ are the eigenvalues and $E$ is a matrix with columns ${\mathbf{e}}_{1},\dots,{\mathbf{e}}_{p}$ that are the corresponding eigenvectors of $X$. A measure of the “positive definiteness” of $X$ is thus ${\left\|\mathsf{\Pi}_{+}X\right\|}_{F}\equiv\sqrt{(\lambda_{1})_{+}^{2}+\cdots+(\lambda_{p})_{+}^{2}}.$ This defines our measure of constraint violation,

[TABLE]

We use special notation to denote the vector of derivatives of $h_{+}$ with respect to the decision variables ${\mathbf{z}}$ evaluated at $({\mathbf{z}};{\mathcal{D}},i)$. In particular, define the gradient $h^{\prime}_{+}\left({\mathbf{z}};{\mathcal{D}},i\right)=\left(\nabla_{1}h_{+}({\mathbf{z}};{\mathcal{D}},i),\dots,\nabla_{m+qn}h_{+}({\mathbf{z}};{\mathcal{D}},i)\right)^{\top}=\nabla_{\mathbf{z}}h_{+}({\mathbf{z}};{\mathcal{D}},i).$ Defining the function $\psi_{ij}\colon{\rm Sym}(p,{\mathbb{R}})\times{\mathbb{R}}^{m+qn}\times\left({\mathbb{R}}^{pm}\times\left({\rm Sym}_{+}(p,{\mathbb{R}})\right)^{m}\right)\to{\mathbb{R}}$ for notational convenience, the $j$-th entry of $h^{\prime}_{+}\!\left({\mathbf{z}};{\mathcal{D}},i\right)$ is given by

[TABLE]

where

[TABLE]

and $d>0$ is some arbitrary constant. The function $\psi_{ij}$ is not continuous at $X={\mathbf{0}},$ and we set $d\neq 0$ for technical reasons discussed below. The step size

[TABLE]

given in line 12 is a variant of the Polyak step size [27]. Note that this is a well defined step size as we let $h_{+}^{\prime}({\mathbf{v}}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})=d{\bf 1}_{m+qn}$ (cf. Eq (11)) whenever $h_{+}({\mathbf{v}}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})=0$ and it is nonzero elsewhere by construction.

IV Main Results

In this section, we discuss our assumptions and state the main result. The detailed proofs are deferred to Section V. The feasible set is ${\mathcal{X}}={\mathcal{X}}_{0}\cap\{{\mathbf{z}}\in{\mathbb{R}}^{m+qn}\mid h\left({\mathbf{z}};{\mathcal{D}},i\right)\preceq 0,\leavevmode\nobreak\ \forall i\in{\mathcal{I}}\}.$ The optimal value and solution set of the problem (6) is $f^{*}=\min_{{\mathcal{X}}}f({\mathbf{z}})$ and ${\mathcal{X}}^{*}=\{{\mathbf{z}}\in{\mathcal{X}}\mid f({\mathbf{z}})=f^{*}\}.$ Note that ${\mathbf{0}}\in\Gamma\times{\mathcal{U}}$ and $h({\mathbf{0}};{\mathcal{D}},i)\preceq 0$ for all possible ${\mathcal{D}}$ and $i,$ thus there is always a feasible point.

At step $k\in{\mathbb{N}}$, ${\mathcal{D}}_{s,k}$ in line 9 of Algorithm 1 may not yet have converged to ${\mathcal{D}},$ so we denote the error in the constraint violation by

[TABLE]

Similar to the definition in (13), we define the error in the gradient of the constraint violation to be

[TABLE]

The vectors $\nu_{s,k}$ and $\boldsymbol{\delta}_{s,k}$ will enable us to quantify the joint effect of the function value error $\nu_{s,k}$ and the subgradient error $\boldsymbol{\delta}_{s,k}$ in studying convergence properties of Algorithm 1.

IV-A Bounds

The norms of $\nu_{s,k}$ and $\boldsymbol{\delta}_{s,k}$ are bounded, i.e., for some scalars $D,N>0$, there holds for all $k\in{\mathbb{N}}$ and $s\in{\mathcal{S}}$ a.s.: $|\nu_{s,k}|\leq N$ and $\|\boldsymbol{\delta}_{s,k}\|\leq D.$ We calculate the exact values of $D$ and $N$ in Corollaries A.2 and A.7, respectively.

The set ${\mathcal{X}}_{0}=\Gamma\times{\mathcal{U}}$ is convex and compact. Therefore, there must exist a constant $C_{\mathbf{z}}>0$ such that for any $\mathbf{z}_{1},\mathbf{z}_{2}\in\mathcal{X}_{0}$

[TABLE]

Also, note that the functions $f(\cdot)$ and $h(\mathbf{z};{\mathcal{D}},i)$ for $i\in{\mathcal{I}}$ are convex (not necessarily differentiable) over some open set that contains $\mathcal{X}_{0}$. A direct consequence of this and the compactness of $\mathcal{X}_{0}$ is that the subdifferentials $\partial f(\mathbf{z})$ and $\partial h_{+}(\mathbf{z};{\mathcal{D}},i)$ are nonempty over $\mathcal{X}_{0}$, where note that here we use the symbol $\partial$ to refer to the subdifferential set and not a partial derivative. It also implies that the subgradients $f^{\prime}(\mathbf{z})\in\partial f(\mathbf{z})$ and $h_{+}^{\prime}(\mathbf{z};{\mathcal{D}},i)\in\partial h_{+}(\mathbf{z};{\mathcal{D}},i)$ are uniformly bounded over the set $\mathcal{X}_{0}$. That is, there is a scalar $L_{f}$ such that for all $f^{\prime}(\mathbf{z})\in\partial f(\mathbf{z})$ and $\mathbf{z}\in\mathcal{X}_{0}$,

[TABLE]

Note that for the objective function $f$ in (6), we have $L_{f}=1.$ Also, there is a scalar $L_{h}$ such that for all $h_{+}^{\prime}(\mathbf{z};{\mathcal{D}},i)\in\partial h_{+}(\mathbf{z};{\mathcal{D}},i)$, $\mathbf{z}\in\mathcal{X}_{0}$, $i\in{\mathcal{I}}$,

[TABLE]

We calculate the exact value of $L_{h}$ in Corollary A.7. From relation (14) we have for any $i\in{\mathcal{I}}$, $s\in{\mathcal{S}}$, $k\in{\mathbb{N}}$, and $\mathbf{z}\in\mathcal{X}_{0}$

[TABLE]

The compactness of $\mathcal{X}_{0}$, the boundedness of the data sequences ${\mathcal{D}}_{s,k}$, and the continuity of the functions $h_{+}(\mathbf{z};{\mathcal{D}},i)$ for $i\in{\mathcal{I}}$ imply that there exist constants $C_{h}>0$ such that for any $s\in{\mathcal{S}}$, $k\in{\mathbb{N}}$, $i\in{\mathcal{I}}$ and $\mathbf{z}\in\mathcal{X}_{0}$

[TABLE]

where the second relation follows from (13). Furthermore, when $h_{+}(\mathbf{z};{\mathcal{D}}_{s,k},i)\neq 0$, we have $h_{+}^{\prime}(\mathbf{z};{\mathcal{D}}_{s,k},i)\neq{\mathbf{0}}$. Therefore, there must exist a constant $c_{h}>0$ such that

[TABLE]

for all $\mathbf{z}\in\mathcal{X}_{0}$, $s\in{\mathcal{S}}$, $k\in{\mathbb{N}}$, and $i\in{\mathcal{I}}$. This and relation (19b) imply that

[TABLE]

Note that when $h_{+}(\mathbf{z};{\mathcal{D}}_{s,k},\omega_{s,k})=0$, the above bound holds trivially.

IV-B Assumptions

In the preceding sections, we have made extensive use of the function $Q\colon{\mathbb{R}}^{q}\times{\mathbb{R}}^{p}\to{\rm Sym}_{+}(p,{\mathbb{R}}).$ For our algorithm to converge, we require the following to hold true regarding this information model.

Assumption IV.1**.**

We assume that the information function $Q$

1. (a)

is bounded, 2. (b)

is twice subdifferentiable 3. (c)

has bounded subdifferentials up to the second order, and 4. (d)

*has relatively few critical points, i.e., the sets of critical points of $Q$ and its partial derivatives up to the second order are measure zero. *

Denote the bound on the magnitude of $Q$ by $\eta_{0}$, the bound on its first derivatives by $\eta_{1}$, and the bound on its second derivatives by $\eta_{2}.$

The technical restrictions on the information model $Q$ are not overly restrictive in practice. In particular, items a, b, and c imply that one cannot obtain infinite information and that, by changing sensor configuration, the information rate (and the rate of change of the rate) cannot change infinitely quickly. Item d essentially allows us to distinguish between sensor configurations, i.e., the set of configurations that offer “optimal” information rates are relatively sparse with respect to the space of signal source positions and sensor configurations ${\mathbb{R}}^{q}\times{\mathbb{R}}^{p}.$

The next two assumptions are related to the random variables $\omega_{s,k}$. At each iteration $k$ of the inner-loop, recall that each sensor $s$ randomly generates $\omega_{s,k}\in{\mathcal{I}}$. We assume that they are i.i.d. samples from some probability distribution on ${\mathcal{I}}$.

Assumption IV.2**.**

*In the $k$-th iteration of the inner-loop, each element $i$ of ${\mathcal{I}}$ is generated with nonzero probability, i.e., for any $s\in{\mathcal{S}}$ and $k\in{\mathbb{N}}$ it holds that $\pi_{i}\triangleq\textrm{Pr}\{\omega_{s,k}=i\}>0.$ *

Assumption IV.2 is easy to satisfy because ${\mathcal{I}}$ is a finite set.

Assumption IV.3**.**

For all $s\in{\mathcal{S}}$ and $k\in{\mathbb{N}}$, there exists a constant $\kappa>0$ such that for all $\mathbf{z}\in\mathcal{X}_{0}$ $\mathrm{dist}^{2}(\mathbf{z},\mathcal{X})\leq\kappa\mathsf{E}\left[h_{+}^{2}(\mathbf{z};{\mathcal{D}}_{s,k},\omega_{s,k})\right].$

The upper bound in Assumption IV.3 is known as global error bound and is crucial for the convergence analysis of Algorithm 1. Sufficient conditions for this bound have been shown in [28, 29, 30, 31], which includes the case when each function $h(\cdot;\cdot,i)$ is linear in ${\mathbf{z}}$, or when the feasible set $\mathcal{X}$ has a nonempty interior.

Assumption IV.4**.**

For all $k\in{\mathbb{N}}$, the weighted graphs $\mathcal{G}_{k}=({\mathcal{S}},\mathcal{E}_{k},W_{k})$ satisfy:

1. (a)

There exists a scalar $\eta\in(0,1)$ such that $[W_{k}]_{sj}\geq\eta$ if $j\in\mathcal{N}_{s,k}$. Otherwise, $[W_{k}]_{sj}=0$. 2. (b)

The weight matrix $W_{k}$ is doubly stochastic, i.e., $\sum_{s\in{\mathcal{S}}}[W_{k}]_{sj}=1$ for all $j\in{\mathcal{S}}$ and $\sum_{j\in{\mathcal{S}}}[W_{k}]_{sj}=1$ for all $s\in{\mathcal{S}}$. 3. (c)

There exists a scalar $Q>0$ such that the graph $\left({\mathcal{S}},\cup_{\ell=0,\ldots,Q-1}\mathcal{E}_{t+\ell}\right)$ is strongly connected for any $t\geq 1$.

This assumption ensures a balanced communication between sensors. It also ensures that there exists a path from one sensor to every other sensor infinitely often even if the underlying graph topology is time-varying.

IV-C Main result

Our main result shows the almost sure convergence of Algorithm 1. Specifically, the result states that the sensors asymptotically reach an agreement to a random point ${\mathbf{z}}^{*}$ which is in the optimal set ${\mathcal{X}}^{*}$ almost surely (a.s.), as given in the following theorem.

Theorem IV.5** (Asymptotic Convergence Under Noise).**

Let Assumptions IV.1 - IV.4 hold. Let also the stepsizes be such that $\sum_{k=1}^{\infty}\alpha_{k}=\infty$ and $\sum_{k=1}^{\infty}\alpha_{k}^{2}<\infty$. Then, the iterates $\{{\mathbf{z}}_{s,k}\}$ generated by each agent $s\in{\mathcal{S}}$ via Algorithm 1 converge almost surely to the same point in the optimal set ${\mathcal{X}}^{*}$, i.e., for a random point ${\mathbf{z}}^{*}\in{\mathcal{X}}^{*}$ $\lim_{k\to\infty}{\mathbf{z}}_{s,k}={\mathbf{z}}^{*},\quad\forall s\in{\mathcal{S}}\quad a.s.$

V Convergence of Algorithm 1

V-A Preliminary Results

Lemma V.1** (Non-expansiveness [32]).**

Let $\mathcal{X}\subseteq\mathbb{R}^{d}$ be a nonempty closed convex set. The function $\mathsf{\Pi}_{\mathcal{X}}:\mathbb{R}^{d}\rightarrow\mathcal{X}$ is nonexpansive, i.e., $\|\mathsf{\Pi}_{\mathcal{X}}[{\mathbf{a}}]-\mathsf{\Pi}_{\mathcal{X}}[{\mathbf{b}}]\|\leq\|{\mathbf{a}}-{\mathbf{b}}\|$ for all ${\mathbf{a}},{\mathbf{b}}\in\mathbb{R}^{d}.$

Lemma V.2** (​​[33, Lemma 10-11, p. 49-50]).**

Let $\mathcal{F}_{k}\triangleq\left\{v_{\ell},u_{\ell},a_{\ell},b_{\ell}\right\}_{\ell=0}^{k}$ denote a collection of nonnegative real random variables for $k\in{\mathbb{N}}\cup\left\{\infty\right\}$ such that $\mathsf{E}[v_{k+1}|\mathcal{F}_{k}]\leq(1+a_{k})v_{k}-u_{k}+b_{k}$ for all $k\in\left\{0\right\}\cup{\mathbb{N}}\quad a.s.$ Assume further that $\left\{a_{k}\right\}$ and $\left\{b_{k}\right\}$ are a.s. summable. Then, we have a.s. that (i) $\left\{u_{k}\right\}$ is summable and (ii) there exists a nonnegative random variable $v$ such that $\left\{v_{k}\right\}\to v.$

In the next lemma, we relate the two iterates $\mathbf{p}_{s,k}$ and $\mathbf{z}_{s,k}$ in Line 6 and 13 of Algorithm 1. In particular, we show a relation of $\mathbf{p}_{s,k}$ and $\mathbf{z}_{s,k-1}$ associated with any convex function $g$ which will be often used in the analysis. For example, $g(\mathbf{z})=\|\mathbf{z}-\mathbf{a}\|^{2}$ for some $\mathbf{a}\in\mathbb{R}^{n}$ or $g(\mathbf{z})=\mathrm{dist}^{2}(\mathbf{z},\mathcal{X})$.

Lemma V.3** (Convexity and Double Stochasticity).**

For any convex function $g:\mathbb{R}^{n}\to\mathbb{R}$, we have

[TABLE]

Proof.

The double stochasticity of the weights plays a crucial role in this lemma. From the definition of $\mathbf{p}_{s,k}$ in Line 6 of Algorithm 1, we obtain

[TABLE]

∎

Lastly, for the convergence proof of our algorithm, we use a result from [34], which ensures successful consensus in the presence of a well behaved disturbance sequence.

Lemma V.4** (Perturbed Consensus).**

Let Assumption IV.4 hold. Consider the iterates generated by

[TABLE]

Suppose there exists a nonnegative nonincreasing scalar sequence $\{\alpha_{k}\}$ such that $\sum_{k=1}^{\infty}\alpha_{k}\|\mathbf{e}_{s,k}\|<\infty\text{ for all }s\in{\mathcal{S}}.$ Then, for all $s,v\in{\mathcal{S}}$, $\sum_{k=1}^{\infty}\alpha_{k}\|\boldsymbol{\theta}_{s,k}-\boldsymbol{\theta}_{v,k}\|<\infty.$

In addition to the well-known results of Lemmas V.1-V.4, we need the following intermediate results before the main result in Theorem IV.5 can be proven. The proofs of all of these Lemmas are deferred to Appendices A and B. In the first two of them, we posit that the two error sequences $\{\nu_{s,k}\}$ and $\{\boldsymbol{\delta}_{s,k}\}$ are summable.

Lemma V.5** (Summable Constraint Errors).**

For almost every bounded sequence $\left\{{\mathbf{z}}_{s,k}\right\},$ the error in the constraint $\{|\nu_{s,k}|\}_{k\in{\mathbb{N}}}$ is summable a.s., i.e., $\sum_{k=0}^{\infty}\left|\nu_{s,k}\right|<\infty.$

Lemma V.6** (Summable Constraint Violation Errors).**

For almost every bounded sequence $\left\{{\mathbf{z}}_{s,k}\right\},$ the error in the gradient of the constraint $\{\|\boldsymbol{\delta}_{s,k}\|\}_{k\in{\mathbb{N}}}$ is summable a.s., i.e., $\sum_{k=0}^{\infty}{\left\|\boldsymbol{\delta}_{s,k}\right\|}<\infty.$

In the next lemma, we relate the two iterates $\mathbf{p}_{s,k}$ and $\mathbf{z}_{s,k}$ in Line 6 and 13 of Algorithm 1.

Lemma V.7** (Basic Iterate Relation).**

Consider the sequences $\{\mathbf{z}_{s,k}\}_{k=0}^{\infty}$ and $\{\mathbf{p}_{s,k}\}_{k=0}^{\infty}$ for $s\in{\mathcal{S}}$ generated by Algorithm 1. Then, for any $\mathbf{z},\check{\mathbf{z}}\in\mathcal{X}$, $s\in{\mathcal{S}}$ and $k\in{\mathbb{N}}$, we have a.s.:

[TABLE]

where $A_{\eta,\tau}=1+4\eta+\tau$ and $\eta,\tau>0$ are arbitrary.

Since we use a random approximate projection, we cannot guarantee the feasibility of the sequences $\mathbf{p}_{s,k}$ and $\mathbf{z}_{s,k}$ for every $k\in{\mathbb{N}}$ and $s\in{\mathcal{S}}$. In the next lemma, we prove that $\mathbf{p}_{s,k}$ for all $s\in{\mathcal{S}}$ asymptotically achieve feasibility even under the effect of the disturbances $\{\nu_{s,k}\}$ and $\{\boldsymbol{\delta}_{s,k}\}$.

Lemma V.8** (Asymptotic Feasibility Under Noise).**

Let Assumption IV.2 and IV.3 hold. Let $\alpha_{k}$ be square summable Consider the sequences $\{\mathbf{p}_{s,k}\}$ for $s\in{\mathcal{S}}$ generated by Algorithm 1. Then, for any $s\in{\mathcal{S}}$, the sequence $\left\{\mathrm{dist}(\mathbf{p}_{s,k},\mathcal{X})\right\}_{k\in{\mathbb{N}}}$ is a.s. square summable.

In the final Lemma, we show that $\|\mathbf{z}_{s,k}-\mathbf{p}_{s,k}\|$ eventually converges to zero for all $s\in{\mathcal{S}}$. This result combined with Lemma V.8 implies that the sequences $\left\{\mathbf{z}_{s,k}\right\}_{k\in{\mathbb{N}}}$ also achieve asymptotic feasibility.

Lemma V.9** (Network Disagreement Under Noise).**

Let Assumptions IV.2-IV.4 hold. Let the sequence $\{\alpha_{k}\}$ is such that $\sum_{k=1}^{\infty}\alpha_{k}^{2}<\infty$. Define $\mathbf{e}_{s,k}=\mathbf{z}_{s,k}-\mathbf{p}_{s,k}$ for all $s\in{\mathcal{S}}$ and $k\geq 1$. Then, we have a.s.:

1. (a)

$\sum_{k=1}^{\infty}\|\mathbf{e}_{s,k}\|^{2}<\infty\quad\text{for all }s\in{\mathcal{S}}$, 2. (b)

\sum_{k=1}^{\infty}\alpha_{k}\|\check{\mathbf{z}}_{s,k}-\bar{\check{\mathbf{z}}}_{k}\|<\infty\quad\hbox{for all s\in{\mathcal{S}}},

where $\check{\mathbf{z}}_{s,k}=\mathsf{\Pi}_{\mathcal{X}}[\mathbf{p}_{s,k}]$ and $\bar{\check{\mathbf{z}}}_{k}=\frac{1}{n}\sum_{s\in{\mathcal{S}}}\check{\mathbf{z}}_{s,k}$.

Note that the sequences $\{\mathbf{z}_{s,k}\}_{k\in{\mathbb{N}}}$ for $s\in{\mathcal{S}}$ generated by Algorithm 1 can be represented as relation (22). That is, $\mathbf{z}_{s,k}=\sum_{s\in{\mathcal{S}}}[W_{k}]_{sj}\mathbf{z}_{j,k-1}+\mathbf{e}_{s,k}$ for all $s\in{\mathcal{S}}$ and $\mathbf{e}_{s,k}=\mathbf{z}_{s,k}-\mathbf{p}_{s,k}$. Therefore, from Lemma V.9(a), we have $\sum_{k=1}^{\infty}\alpha_{k}\|\mathbf{e}_{s,k}\|=\frac{1}{2}\sum_{k=1}^{\infty}\alpha_{k}^{2}+\frac{1}{2}\sum_{k=1}^{\infty}\|\mathbf{e}_{s,k}\|^{2}<\infty.$ Invoking Lemma V.4, we see that there is consensus among $\mathbf{z}_{s,k}$ for $s\in{\mathcal{S}}$.

V-B Proof of Theorem IV.5

We invoke Lemma V.7 with $\check{\mathbf{z}}=\check{\mathbf{z}}_{s,k}$ (Note that $\check{\mathbf{z}}_{s,k}$ is defined as $\mathsf{\Pi}_{\mathcal{X}}[\mathbf{p}_{s,k}]$ in Lemma V.9), $\tau=4$ and $\eta=\kappa(L_{h}+D)^{2}$. We also let $\mathbf{z}=\mathbf{z}^{*}$ for an arbitrary $\mathbf{z}^{*}\in\mathcal{X}^{*}$. Therefore, for any $\mathbf{z}^{*}\in\mathcal{X}^{*}$, $s\in{\mathcal{S}}$ and $k\in{\mathbb{N}}$, we have a.s.:

[TABLE]

where $A=5+4\kappa(L_{h}+D)^{2}$. From Assumption IV.3, we know that $\mathrm{dist}^{2}(\mathbf{p}_{s,k},\mathcal{X})\leq\kappa\mathsf{E}\left[h^{2}_{+}(\mathbf{p}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})\mid\mathcal{F}_{k-1}\right].$

Denote by $\mathcal{F}_{k}$ the sigma-field induced by the history of the algorithm up to time $k$, i.e., $\mathcal{F}_{k}=\{\mathbf{z}_{s,0},\leavevmode\nobreak\ (\omega_{s,\ell},\leavevmode\nobreak\ 1\leq\ell\leq k),\leavevmode\nobreak\ s\in{\mathcal{S}}\}\quad\text{for }k\in{\mathbb{N}},$ and $\mathcal{F}_{0}=\{\mathbf{z}_{s,0},\leavevmode\nobreak\ s\in{\mathcal{S}}\}$. Taking the expectation conditioned on $\mathcal{F}_{k-1}$ in relation (V-B), summing this over $s\in{\mathcal{S}}$, and using the above relation, we obtain

[TABLE]

where we used Lemma V.3 for the first term on the right-hand side. Recall that $f(\mathbf{z})=\sum_{s\in{\mathcal{S}}}f(\mathbf{z})$ and $\bar{\check{\mathbf{z}}}_{k}\triangleq\frac{1}{n}\sum_{\ell\in{\mathcal{S}}}\check{\mathbf{z}}_{\ell,k}$. Using $\bar{\check{\mathbf{z}}}_{k}$ and $f$, we can rewrite the term $\sum_{s\in{\mathcal{S}}}(f(\check{\mathbf{z}}_{s,k})-f(\mathbf{z}^{*}))$ as follows:

[TABLE]

where (25a) follows from adding and subtracting $f(\bar{\check{\mathbf{z}}}_{k})$; (25b) follows from the convexity of the function $f$; (25c) follows from the Schwarz inequality; and (25d) follows from relation (16a) and $\bar{\check{\mathbf{z}}}_{k}\in\mathcal{X}\subseteq\mathcal{X}_{0}$. Combining (25d) with (24), we obtain

[TABLE]

where we omitted the negative term.

Since $\bar{\check{\mathbf{z}}}_{k}\in\mathcal{X}$, we have $f(\bar{\check{\mathbf{z}}}_{k})-f^{*}\geq 0$. Thus, under the assumption $\sum_{k=0}^{\infty}\alpha_{k}^{2}<\infty$ and Lemma V.9(b), the above relation satisfies all the conditions of Lemma V.2. Using this, we have the following results.

Result 1: For some ${\mathbf{z}}^{*}\in{\mathcal{X}}^{*}$ and any $s\in{\mathcal{S}}$, the sequence $\{\|\mathbf{z}_{s,k}-\mathbf{z}^{*}\|\}$ is convergent a.s.

Result 2: $\sum_{k=1}^{\infty}\alpha_{k}(f(\bar{\check{\mathbf{z}}}_{k})-f^{*})<\infty\quad a.s.$

As a direct consequence of Result 1, we know that the sequences $\{\|\mathbf{p}_{s,k}-\mathbf{z}^{*}\|\}$ and $\{\|\check{\mathbf{z}}_{s,k}-\mathbf{z}^{*}\|\}$ are also convergent a.s. (This is straightforward from Line 6 of Algorithm 1 and Lemma V.8 with the knowledge that $\mathrm{dist}^{2}(\mathbf{p}_{s,k},\mathcal{X})=\|\mathbf{p}_{s,k}-\check{\mathbf{z}}_{s,k}\|^{2}$.) Since $\|\bar{\check{\mathbf{z}}}_{k}-\mathbf{z}^{*}\|\leq\frac{1}{N}\sum_{i\in{\mathcal{S}}}\|\check{\mathbf{z}}_{s,k}-\mathbf{z}^{*}\|$, we further know that the sequence $\{\|\bar{\check{\mathbf{z}}}_{k}-\mathbf{z}^{*}\|\}$ is also convergent a.s.

As a direct consequence of Result 2, since $\alpha_{k}$ is not summable, $\liminf_{k\rightarrow\infty}(f(\bar{\check{\mathbf{z}}}_{k})-f^{*})=0\;a.s.$ From this relation and the continuity of $f$, it follows that the sequence $\{\bar{\check{\mathbf{z}}}_{k}\}$ must have one accumulation point in the set $\mathcal{X}^{*}$ a.s. This and the fact that $\{\|\bar{\check{\mathbf{z}}}_{k}-\mathbf{z}^{*}\|\}$ is convergent a.s. for every $\mathbf{z}^{*}\in\mathcal{X}^{*}$ imply that

[TABLE]

Also, by Lemma V.9(b), we have

[TABLE]

The fact that $\{\|\check{\mathbf{z}}_{s,k}-\mathbf{z}^{*}\|\}$ is convergent a.s. for all $s\in{\mathcal{S}}$, together with (26) and (27) implies that

[TABLE]

Since $\|\mathbf{z}_{s,k}-\check{\mathbf{z}}_{s,k}\|\leq\|\mathbf{z}_{s,k}-\mathbf{p}_{s,k}\|+\|\mathbf{p}_{s,k}-\check{\mathbf{z}}_{s,k}\|$, by invoking Lemma V.8 and V.9(a), we have $\lim_{k\to\infty}\|\mathbf{z}_{s,k}-\check{\mathbf{z}}_{s,k}\|=0$ for all $s\in{\mathcal{S}}\quad a.s.$ This and relation (28) give us the desired result, which is $\lim_{k\rightarrow\infty}\mathbf{z}_{s,k}=\mathbf{z}^{\star}$ for all $s\in{\mathcal{S}}\quad a.s.$

VI Numerical Experiments

In this section, we illustrate Algorithm 1 on a network of $n$ mobile sensors that collaboratively localize a collection of $m$ landmarks. We use a distributed connectivity maintenance controller [35] to ensure that the sensor network can share information and collaborate for all time. The controller operates by signaling for a strong attractive force whenever a link is going to be deleted that can break connectivity. The controller also deletes or replaces weak links as long as connectivity is preserved.

VI-A Sensor Model

We assume that the landmarks live in ${\mathbb{R}}^{2}$ and denote the configuration of landmark $i$ by ${\mathbf{x}}_{i}.$ We also assume that the configuration space of the robot is ${\mathbb{R}}^{2}$.

Generic models for the measurement covariance matrix for sparse landmark localization have been proposed [1, 7, 2]. These models have two important characteristics that apply to a wide variety of localization sensors:

(i) Measurement quality is inversely proportional to viewing distance, and

(ii) The direction of maximum uncertainty is the viewing direction, i.e., the sensors are more proficient at sensing bearing than range.

The general idea is to use the vector ${\hat{\mathbf{x}}}_{i}(t)-{\mathbf{r}}_{s}(t)$ to construct a diagonal matrix in a coordinate frame relative to the sensor, then rotate the matrix to a global coordinate frame so that it can be compared to other observations.

Following the accepted generic sensor models, we express $Q({\mathbf{r}}_{s},{\hat{\mathbf{x}}}_{i})$ by its eigenvalue decomposition $Q=R_{\phi_{i}}\Lambda_{i}R_{\phi_{i}}^{\top},$ where the argument $({\mathbf{r}}_{s},{\hat{\mathbf{x}}}_{i})$ of $Q$, $\Lambda_{i}$, and $R_{\phi_{i}}$ is implicit. The angle $\phi_{i}\in\left(-\frac{\pi}{2},\frac{\pi}{2}\right]$ is defined so that the first column of $R_{\phi_{i}}$, which by convention is $[\cos\phi_{i}\;\sin\phi_{i}]^{\top}\!\!,\;$ is parallel to the subspace spanned by the vector ${\hat{\mathbf{x}}}_{i}-{\mathbf{r}}_{s}$. This angle is given by $\phi_{i}=\tan^{-1}\left(([{\hat{\mathbf{x}}}_{i}]_{2}-[{\mathbf{r}}_{s}]_{2})/([{\hat{\mathbf{x}}}_{i}]_{1}-[{\mathbf{r}}_{s}]_{1})\right),$ where the subscripts outside the brackets refer to the first and second coordinate of the vector representing the location of the landmark and sensor located at ${\hat{\mathbf{x}}}_{i}$ and ${\mathbf{r}}_{s}$. The eigenvalue matrix is given by $\Lambda_{i}=\textrm{diag}(\lambda_{1,i},\lambda_{2,i}),$ where $\lambda_{1,i}=c_{0}(/1+c_{1}{\left\|{\hat{\mathbf{x}}}_{i}-{\mathbf{r}}_{s}\right\|}_{2})\lambda_{2,i}=\rho\lambda_{1,i}.$ In the body frame of sensor $s$, the eigenvalues $\lambda_{1,i}$ and $\lambda_{2,i}$ control the shape of the confidence ellipses for individual measurements of the $i$-th target. The parameter $c_{0}>0$ represents the overall sensor quality and scales the whole region equally, $c_{1}>0$ controls the sensitivity to depth, and $\rho>1$ controls eccentricity of the confidence regions associated with the measurements. Figure 1 illustrates the measurement model for one sensor and two targets.

Since $\tan^{-1}$ is not continuous at zero, we need to impose the following restriction so that Assumption IV.1 holds. Essentially, we assume that, for the selected $\delta>0$, which recall has an affect on the robot’s ability to translate in ${\mathbb{R}}^{2},$ we choose a $\tau_{i}$ such that $\tau_{i}<c_{0}/(1+c_{1}\delta).$ Thus, if the robot is within $\delta$ of a target $i$, the constraint $h({\mathbf{z}},{\mathcal{D}},i)$ is trivially satisfied, and therefore we will never evaluate $Q$ or any of its derivatives in this region.

VI-B Results

Figure 2 shows an example of the resulting active sensing trajectories for a network of four sensors that cooperatively localize a uniformly random distribution of 100 landmarks in a square $100\,{{\rm m}}^{2}$ workspace to a desired accuracy of $1\,{{\rm m}}$. This corresponds to an eigenvalue tolerance of $\tau_{1}=\cdots=\tau_{m}=9\,{{\rm m}}^{-2}.$ The prior distribution on the landmarks was uninformative, i.e., the initial mean estimate and information matrix are given by the initial observations $\left\{{\mathbf{y}}_{i,s}(0)\mid i\in{\mathcal{I}},s\in{\mathcal{S}}\right\}$ and the corresponding information matrices were set to $\left\{Q({\mathbf{r}}_{s}(0),{\mathbf{y}}_{i,s}(0))\right\}_{i\in{\mathcal{I}},s\in{\mathcal{S}}}.$ The parameters of the sensing model $Q$ were set to $c_{0}=\frac{1}{2},c_{1}=10,\rho=30$. The feasible control set ${\mathcal{U}}$ was set to ${\mathcal{U}}=\left\{{\mathbf{u}}_{s}\in{\mathbb{R}}^{2}\mid{\left\|{\mathbf{u}}_{s}\right\|}\leq 1\,{{\rm m}}\right\}.$ In the figure, note that between $t=194$ and $t=242$, the agent in the bottom right moves north in order to maintain connectivity while the other agents move to finish the localization task. It can be seen from the four snapshots of the trajectories in Fig. 2, that the robots effectively “divide and conquer” the task based on proximity to unlocalized targets.

Fig. 3 shows the evolution of the minimum eigenvalues of the information matrices $S_{i}(t)$ from the simulation shown in Fig. 2. Note that, even after the minimum eigenvalue associated with a landmark has surpassed $\tau$, more information is still collected, but it does not affect the control signal. Note also that the maximal uncertainty in the bottom panel does not go below the desired threshold, even though the minimal eigenvalue does go above the threshold in the top panel. This is expected because we use the 95% confidence interval to generate the threshold $\tau$, thus we expect 5% of landmarks to have errors larger than the desired threshold.

In this case study, we selected the termination criterion for Algorithm 1 to be at most 30 iterations. When testing the algorithm with $\alpha_{k}=1/k$, we observed that the relative change in the iterates ${\mathbf{z}}_{s,k}$ and the data ${\mathcal{D}}_{s,k}$ would often reach levels at or below 10% inside of thirty iterations, so we allowed the message passing to terminate early if this criterion was met.

Figure 4 measures the amount of constraint violation of the original nonlinear constraints for the solution obtained from the linearized version that we observed in a simulation with twenty landmarks. For the ${\mathcal{U}}$ and the parameters for the function $Q$ used in these simulations, the constraint violation was very small, and driven to zero typically for $t>10.$ We tested larger feasible sets ${\mathcal{U}}$ and (predictably) observed a general increase in the constraint violation, although we still observed reasonable robot behaviors.

VII Conclusion

In this paper, we addressed the problem of controlling a network of mobile sensors so that a set of hidden states are estimated up to a user-specified accuracy. We proposed a new algorithm that combines the ICF with control. The first major contribution is the capability of controlling worst-case error in the hidden states. The second major contribution is the novel numerical technique, which is robust toward noisy gradients, computationally light, can handle large systems of coupled LMIs. We showed that our method converges to the optimal solution and presented simulations for sparse landmark localization, where the team was able to achieve the desired estimation tolerances while exhibiting interesting emergent behaviors, notably “dividing and conquering” the distributed estimation task, observing landmarks from diverse viewing angles.

Appendix A Proof of Lemmas V.5 and V.6

To prove Lemmas V.5 and V.6, consider the following.

Proposition A.1** (Summable Disagreement in Constraints).**

Let Assumptions IV.1 - IV.4 hold. For any bounded sequence $\left\{{\mathbf{z}}_{s,k}\right\}$ and any $i\in{\mathcal{I}},$

[TABLE]

Proof.

In this proof, we will essentially be comparing the function $h$ for the unconverged local data ${\mathcal{D}}_{s,k}$ and the consensus data ${\mathcal{D}}.$ Note that at some iteration $k\in{\mathbb{N}},$ the ICF is actually computing the information summary $(\boldsymbol{\xi}_{1,s,k},\dots,\boldsymbol{\xi}_{m,s,k},\Xi_{1,s,k},\dots,\Xi_{m,s,k})$, which is used to compute ${\mathcal{D}}_{s,k}$ in line 9 of Algorithm 1.

In particular, we will need to refer to the unconverged mean and information matrix, which is

[TABLE]

In the proof, we will drop the reference to the sensor $s$, as we only consider one sensor at a time. In other words, the unconverged data will be written as

[TABLE]

for $k\in{\mathbb{N}},$ i.e., this is implicitly the data local to sensor $s$ in the proof. For the consensus data ${\mathcal{D}}_{\infty}=\left(\left({\hat{\mathbf{x}}}_{1,\infty},\dots,{\hat{\mathbf{x}}}_{m,\infty}\right),\left(S_{1,\infty},\dots,S_{m,\infty}\right)\right)$, we also drop the $k$ subscript, i.e., we write the consensus data as ${\mathcal{D}}=({\hat{\mathbf{x}}},S)=\left(\left({\hat{\mathbf{x}}}_{1},\dots,{\hat{\mathbf{x}}}_{m}\right),\left(S_{1},\dots,S_{m}\right)\right)$.

Let us fix an arbitrary $i\in{\mathcal{I}}$, and prove the proposition for the sequence $\omega_{k}=\left\{i,i,i,\dots\right\}.$ If the proof holds for this sequence, it thus holds for the random sequences $\omega_{k}\subset{\mathcal{I}}$ because

[TABLE]

and if each sequence ${\left\|h\left({\mathbf{z}}_{k};{\mathcal{D}}_{k},i\right)-h\left({\mathbf{z}}_{k};{\mathcal{D}},i\right)\right\|}_{F}$ is summable, so is their sum.

Recall the definition of $h$ from (5). We have that

[TABLE]

where $S_{i}$ and ${\hat{\mathbf{x}}}_{i}$ are the eventual consensus variables that compose the consensus data ${\mathcal{D}},$ i.e., $S_{i}=S_{i,\infty}$ and ${\hat{\mathbf{x}}}_{i}={\hat{\mathbf{x}}}_{i,\infty}.$ The remainder of the proof will show that each of the constituent terms, specifically, (29a), (29b), and (29c), are summable, and the proof will follow from that immediately.

(29a) For this term, recall from line 9 of Algorithm 1 that $S_{i,k}\leftarrow n\Xi_{i,k}$, where $n$ is the number of agents in the network. It is well known that for the consensus sequence $\left\{\Xi_{i,k}\right\}_{k\in{\mathbb{N}}}$ the errors $\left\{\Xi_{i,k}-\Xi_{i,\infty}\right\}_{k\in{\mathbb{N}}}$ are summable; cf., e.g., [36]. Noting that $S_{i}=n\Xi_{i,\infty}$ and $S_{i,k}=n\Xi_{i,k},$ we have the result.

(29b) By Assumption IV.1, ${\left\|\nabla_{q+j}Q({\mathbf{r}}_{v},{\mathbf{x}}_{i})\right\|}_{F}\leq\eta_{1}$ for all ${\mathbf{r}}_{v}$, ${\mathbf{x}}_{i}$ and $j$. Then, for each $v\in{\mathcal{S}}$ and $k\in{\mathbb{N}},$ using the taylor expansion of $Q({\mathbf{r}}_{v},\cdot)$ about the point ${\hat{\mathbf{x}}}_{i}$, we have that

[TABLE]

which is summable because ${\left\|{\hat{\mathbf{x}}}_{i,k}-{\hat{\mathbf{x}}}\right\|}$ is the error sequence of a consensus filter, which is summable.

(29c) First, since ${\mathbf{u}}_{v}<\delta$ for all $v\in{\mathcal{S}},$ for any $j=1,\dots,p$, we have that

[TABLE]

Also due to Assumption IV.1, $\eta_{2}$ is the upper bound for the norm of the second derivative of $Q$. In particular, let $\nabla_{q+\ell}\nabla_{j}Q\left({\mathbf{r}}_{v},{\hat{\mathbf{x}}}_{i}\right)$ denote the mixed second partial derivative of $Q$ with respect to the $\ell$-th component of the ${\mathbf{x}}_{i}$ coordinates and the $j$-th component of the ${\mathbf{r}}_{v}$ coordinates of the function $Q$, evaluated at the point $({\mathbf{r}}_{v},{\hat{\mathbf{x}}}_{i}).$ This notation is useful for the following inequality, which uses the first order Taylor expansion with respect to the ${\mathbf{x}}_{i}$ coordinates of $\nabla_{j}Q$. In particular, the norm on the right hand side of (30) is bounded above by

[TABLE]

The sequence above is summable again because ${\left\|{\hat{\mathbf{x}}}_{i,k}-{\hat{\mathbf{x}}}_{i}\right\|}$ is the error sequence of a consensus filter, which is summable.

Since (29a), (29b), and (29c) are each summable, their sum is also a summable sequence, and the proof follows.∎

Proof of Lemma V.5.

Suppose $\omega_{s,k}=\left\{i,i,i,\dots\right\}$. If the proof holds for this arbitrary $i$, then it holds for arbitrary sequences $\omega_{s,k}\subset{\mathcal{I}}$ for the same reasons as in Proposition A.1. Let $X_{k}=h\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},i\right)$ and $Y_{k}=h\left({\mathbf{z}}_{s,k};{\mathcal{D}},i\right)$. From the definition of $\nu_{s,k},$

[TABLE]

which is summable by Proposition A.1. In the above manipulations, (31b) is less than or equal to (31c) by the reverse triangle inequality and (31c) is less than or equal to (31d) because ${\rm Sym}_{+}(p,{\mathbb{R}})$ is a convex set. Thus, $\left|\nu_{s,k}\right|$ is element-wise less than or equal to a summable sequence, and the proof follows. ∎

Corollary A.2** (Bound for $\left|\nu_{s,k}\right|$).**

The constant $N$, which is the bound for $\left|\nu_{s,k}\right|$, is given by $N=\frac{1}{n}\max_{v\in{\mathcal{S}}}{\left\|(\Xi_{i,0})_{v}\right\|}_{F}+m(\eta_{0}+2\eta_{1}q\delta).$

Proof.

From the proof of Lemma V.5, each term $\left|\nu_{s,k}\right|$ is bounded above by ${\left\|X_{k}-Y_{k}\right\|}_{F},$ which, in view of (29), is the sum of three terms. The corollary follows by bounding each term, then taking the sum of the bounds. We have that

[TABLE]

The first term (32a) is the error sequence of a consensus filter, which is bounded by the largest initial condition in the network $\frac{1}{n}\max_{v\in{\mathcal{S}}}{\left\|(\Xi_{i,0})_{v}\right\|}_{F}$. For the second term (32b), recall Assumption IV.1 bounded ${\left\|Q\right\|}_{F}$ by $\eta_{0}.$ Given that $Q$ is positive semidefinite, (32b) is bounded by $m\eta_{0}.$ Similarly, the third term is bounded by $2\eta_{1}q\delta$, and, by taking the sum, the proof follows. ∎

The proof of Lemma V.6 requires the help of several intermediate results. First, we present two basic results in classical analysis. Both are direct results of the Implicit Function Theorem; see, e.g., [37].

Corollary A.3**.**

Let $p\colon{\mathbb{R}}^{n}\to{\mathbb{R}}$ be a nonzero polynomial function. Then, the set of roots of $p$ has measure zero in ${\mathbb{R}}^{n}.$

Corollary A.4**.**

Let $f\colon{\mathbb{R}}^{n+m}\to{\mathbb{R}}^{m}$ be a smooth function, $f^{\prime}$ surjective, and $A\subset{\mathbb{R}}^{m}$ a measure zero set. Then, the preimage of $A$ under $f$ is measure zero in ${\mathbb{R}}^{n+m}.$

Proposition A.5** (Set of Reachable Information Matrices that Are Singular or Have Nonsimple Eigenvalues is Measure Zero).**

Let

[TABLE]

Now, recall the function $h$ defined in (5) and fix arbitrary $i\in{\mathcal{I}}.$ The set

[TABLE]

is measure zero.

Proof.

Both the determinant and the discriminant of the characteristic polynomial are nonvanishing polynomials on ${\rm Sym}(p,{\mathbb{R}}).$ Thus, ${\mathcal{Z}}_{1}$ and ${\mathcal{Z}}_{2}$ represent the zero locii of their respective polynomials, and they are measure zero sets by Corollary A.3.

To complete the proof, we need only to check surjectivity of $h^{\prime}$ and apply Theorem A.4. Recall the first definition of $h$ with respect to the constant ${\mathbf{r}}$ and the variables $\gamma,S,{\hat{\mathbf{x}}},$ and ${\mathbf{u}}$, which were later shortened to the inputs ${\mathbf{z}}$ and the data ${\mathcal{D}}.$ Omitting unnecessary subscripts from the original definition in (5), slightly abusing notation, we have

[TABLE]

Note that, in this notation, $h:{\mathbb{R}}\times{\rm Sym}_{++}(p,{\mathbb{R}})\times{\mathbb{R}}^{p+qn}\to{\rm Sym}(p,{\mathbb{R}}).$ Since ${\rm Sym}_{++}(p,{\mathbb{R}})$ is homeomorphic to ${\mathbb{R}}^{p(p+1)/2}$, without loss of generality, we may write $h$ as a map from ${\mathbb{R}}^{a+b}$ to ${\mathbb{R}}^{b}$ with $a=1+p(p+1)/2+p+qn$ and $b=p(p+1)/2.$ If we write the derivative $h^{\prime}:{\mathbb{R}}^{a+b}\to{\mathbb{R}}^{b}$ as a $b\times a$ matrix $\nabla h(\gamma,S,{\hat{\mathbf{x}}}_{i},{\mathbf{u}}),$ the partial derivatives corresponding to the vectorized version of $S$, i.e., $\nabla_{S}h(\gamma,S,{\hat{\mathbf{x}}}_{i},{\mathbf{u}})$ will be $I_{b}$, so that $\nabla h(\gamma,S,{\hat{\mathbf{x}}}_{i},{\mathbf{u}})$ will be easily be rank $b,$ which is the required condition for surjectivity of $h^{\prime}$. ∎

Proposition A.6** (Bounded Entries of $\nabla h_{+}$).**

For any input ${\mathbf{z}}$ and data ${\mathcal{D}}$, the function $\psi_{ij}(\cdot,{\mathbf{z}},{\mathcal{D}})$, defined in (11), representing the $j$-th entry of the vector $\nabla h_{+}({\mathbf{z}};{\mathcal{D}},i)$, is bounded for all $i\in{\mathcal{I}}$ and $j\in\left\{1,\dots,m+qn\right\}.$

Proof.

Recall the definition of $\psi_{ij}(\cdot,{\mathbf{z}},{\mathcal{D}})\colon{\rm Sym}(p,{\mathbb{R}})\to{\mathbb{R}}$, repeated here for convenience,

[TABLE]

Note that $\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i)$ is a constant matrix for any $j$ because $h$ is a linear function of ${\mathbf{z}}.$ To simplify notation, call this matrix $A_{ij}$. Let $\boldsymbol{\lambda}=[\lambda_{1}\;\cdots\;\lambda_{p}]^{\top}$ denote the vector of eigenvalues of $X$ and ${\mathbf{v}}_{j}$ denote the eigenvector corresponding to $\lambda_{j}.$ We can write $X=\sum_{\alpha=1}^{p}\lambda_{\alpha}{\mathbf{v}}_{\alpha}{\mathbf{v}}_{\alpha}^{\top}.$ If $X\neq{\mathbf{0}}$, then ${\left\|\boldsymbol{\lambda}\right\|}\neq 0$ and

[TABLE]

Let $B_{ij}\triangleq\textrm{diag}(v_{1}^{\top}A_{ij}v_{1},\dots,v_{p}^{\top}A_{ij}v_{p}),$ and note that $\psi_{ij}={{\bf 1}}^{\top}B_{ij}\boldsymbol{\lambda}/{\left\|\boldsymbol{\lambda}\right\|}.$ Note that the limit of $\psi_{ij}(X,{\mathbf{z}},{\mathcal{D}})$ therefore exists for sequences $\left\{X_{k}\right\}\to{\mathbf{0}}$. With respect to bounding the magnitude of $\psi_{ij}$, we have that

[TABLE]

Therefore, $\psi_{ij}$ is bounded. Additionally, $\left|\psi_{ij}\right|$ is maximized if $\left\{{\mathbf{v}}_{\alpha}\right\}$ represent a set of orthogonal eigenvectors of $A_{ij}$, that is, if $X$ has the same eigenvectors as $A_{ij}$. In this case, the bound is tight and equal to ${\left\|\boldsymbol{\mu}_{ij}\right\|},$ where $\boldsymbol{\mu}_{ij}$ is the vector of eigenvalues of the matrix $A_{ij}\triangleq\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i).$

It is left to show that $A_{ij}$ is bounded. This obvious when we write out the individual terms of $h$ from (5). It can be seen that if $j\leq m,$ then $A_{ij}=I.$ Otherwise, it is the matrix of partial derivatives of $Q$, which are bounded by Assumption IV.1. ∎

Corollary A.7** (Bound for ${\left\|\boldsymbol{\delta}_{s,k}\right\|}$).**

Let $A_{ij}\triangleq\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i).$ The constant $L_{h}$ and the constant $D$, which are bounds for ${\left\|h^{\prime}_{+}({\mathbf{z}},{\mathcal{D}},i)\right\|}$ and ${\left\|\boldsymbol{\delta}_{s,k}\right\|}$, respectively, are given by $L_{h}=\sqrt{\bar{\mu}^{2}_{1}+\cdots+\bar{\mu}^{2}_{m+qn}},\;\;D=2L_{h},$ where $\bar{\mu}_{j}\triangleq\max_{i\in{\mathcal{I}}}\left\{{\left\|\boldsymbol{\mu}_{ij}\right\|}\mid\boldsymbol{\mu}_{ij}\right\}$ is the vector of eigenvalues of $A_{ij}.$

Proof.

The $j$-th entry of $h^{\prime}_{+}$ is a restriction of $\psi_{ij}$ to a subset of its domain. Thus, bounds for $\psi_{ij}$ also hold for the $j$-th entry of $h^{\prime}_{+}.$ In the proof of Proposition A.6, we saw that $\psi_{ij}$ is bounded by the 2-norm of the eigenvalues of $A_{ij}\triangleq\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i).$ Therefore, $L_{h}={\left\|\boldsymbol{\mu}_{ij}\right\|}_{2}.$ Taking the norm of both sides of (14), we have that

[TABLE]

Thus, the $j$-th entry of $\boldsymbol{\delta}_{s,k}$ is bounded by $2\bar{\mu}_{j}$ and ${\left\|\boldsymbol{\delta}_{s,k}\right\|}$. Taking the norm, we retrieve the result of the Corollary. ∎

Proposition A.8** (Bounded Derivatives of $\psi_{ij}$).**

Assume $X\in{\rm Sym}_{+}(p,{\mathbb{R}}).$ Then, the magnitudes of the partial derivatives ${\left\|\nabla_{\mathcal{D}}\psi_{ij}(X,{\mathbf{z}},{\mathcal{D}})\right\|}$ and ${\left\|\nabla_{X}\psi_{ij}(X,{\mathbf{z}},{\mathcal{D}})\right\|}$ are point-wise bounded for all $i\in{\mathcal{I}}$ and $j\in\left\{1,\dots,m+qn\right\}.$

Proof.

If ${\left\|X\right\|}=0,$ then $\psi_{ij}$ is constant, and thus the magnitudes of its partial derivatives are zero. The rest of the proof deals with ${\left\|X\right\|}>0.$

To simplify notation, we refer to $\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i)$ by the shorthand $A_{ij}.$ Choose $\epsilon\in{\mathbb{R}}$ such that ${\left\|X\right\|}>\epsilon>0.$ Fix arbitrary $i$. Since $X$ and $A_{ij}$ are symmetric, for any $j$ we have that $\nabla_{X}\!\psi_{ij}(X,{\mathbf{z}},{\mathcal{D}})\!=\!\frac{1}{{\left\|X\right\|}}A_{ij}-\frac{{{\rm Tr}}A_{ij}X}{{\left\|X\right\|}^{3}}X$ and

[TABLE]

which is strictly less than 1 for all ${\left\|X\right\|}>{\left\|A_{ij}\right\|}.$ For ${\left\|X\right\|}\leq{\left\|A_{ij}\right\|},$ $({\left\|A_{ij}\right\|}/\epsilon)^{4}$ is an upper bound for (36).

For the term $\nabla_{\mathcal{D}}\psi_{ij}(X,{\mathbf{z}},{\mathcal{D}}),$ note that the only part of $\psi_{ij}$ that depends on ${\mathcal{D}}$ is $A_{ij}.$ Thus by the chain rule,

[TABLE]

The first term is given by $\nabla_{A_{ij}}\psi_{ij}(X,{\mathbf{z}},{\mathcal{D}})=\frac{1}{{\left\|X\right\|}_{F}}I\prec\frac{1}{\epsilon}I.$ For the second term, recall that ${\mathcal{D}}=({\hat{\mathbf{x}}},S).$ Split the partial derivative into two parts $\nabla_{\mathcal{D}}=(\nabla_{{\hat{\mathbf{x}}}},\nabla_{S})$ so that

[TABLE]

For any $j$, $\nabla_{S}\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i)=I_{p},$ and obviously ${\left\|I_{p}\right\|}\leq p.$ Similarly, $\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i)=I_{p}$ for $j\leq m.$ If $j>m,$ then $\nabla_{{\hat{\mathbf{x}}}}\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i)$ results in second mixed partial derivatives of $Q$, which we assumed to be bounded in Assumption IV.1.∎

Proposition A.9** (Local Differentiability of Constraint Gradient).**

Consider an arbitrary ${\mathbf{z}}\in{\mathbb{R}}^{m+qn}$ and ${\mathcal{D}}=({\hat{\mathbf{x}}},S)$ with ${\hat{\mathbf{x}}}_{i}\in{\mathbb{R}}^{p}$ and $S_{i}\in{\rm Sym}_{+}(p,{\mathbb{R}}).$ Suppose $h\left({\mathbf{z}};{\mathcal{D}},i\right)$ is nonsigular with distinct eigenvalues. Then, $\exists$ an open neighborhood $\Omega$ containing $h\left({\mathbf{z}};{\mathcal{D}},i\right)$ such that \psi_{ij}\left(\mathsf{\Pi}_{+}\left(\cdot\right),{\mathbf{z}},{\mathcal{D}}\right)\big{|}_{\Omega} is differentiable.

Proof.

Notation Note that $f\in C^{\infty}({\mathcal{A}},{\mathcal{B}})$ means that the $n$-th derivative of $f$ exists at all points in ${\mathcal{A}}$ for all $n\in{\mathbb{N}}.$ We also use the notation $f|_{\Omega}$ to denote the restriction of $f$ to $\Omega.$

Roughly speaking, the proof uses a well known result [38] to show that, if $h\left({\mathbf{z}};{\mathcal{D}},i\right)$ has distinct eigenvalues, then for $\alpha=1,\dots,p$ we can find an open neighborhood $\Omega$ in which there are differentiable maps $\lambda_{\alpha}\in C^{\infty}(\Omega_{\alpha},{\mathbb{R}})$ and ${\mathbf{v}}_{\alpha}\in C^{\infty}(\Omega_{\alpha},{\mathbb{R}}^{p})$ that represent the $\alpha$-th eigenvalue-eigenvector pair of each $X\in\Omega$ and agree with the eigenvalues and eigenvectors of $h\left({\mathbf{z}};{\mathcal{D}},i\right)$. Since $\psi_{ij}\left(\mathsf{\Pi}_{+}(\cdot),{\mathbf{z}},\Sigma\right)$ is well behaved in an open neighborhood, it is not difficult to arrive at the desired result. The details of the proof are as follows, and generalize the line of reasoning in [38].

Let $(\bar{\lambda}_{\alpha},\bar{{\mathbf{v}}}_{\alpha})$ be an eigenvalue-eigenvector pair of $h\left({\mathbf{z}};{\mathcal{D}},i\right).$ Since $h\left({\mathbf{z}};{\mathcal{D}},i\right)\in{\rm Sym}(p,{\mathbb{R}})$, we know that $\bar{\lambda}_{\alpha}$ is real valued and distinct. Consider the mapping $\boldsymbol{\phi}\colon{\mathbb{R}}\times{\mathbb{R}}^{p}\times{\rm Sym}(p,{\mathbb{R}})\to{\mathbb{R}}^{p+1}$ given by

[TABLE]

By definition, $\boldsymbol{\phi}(\lambda,{\mathbf{u}},X)={\mathbf{0}}$ if and only if $(\lambda,{\mathbf{v}})$ are a (normalized) eigenvalue-eigenvector pair of $X$. It is also a well known properties of eigenvalues and eigenvectors that

[TABLE]

is full rank if and only if $\bar{\lambda}_{\alpha}$ is simple, which is true by assumption. By the implicit function theorem, there is an open neighborhood $\Omega_{\alpha}$ of $h\left({\mathbf{z}};{\mathcal{D}},i\right)$ in which we can find differentiable functions $\lambda_{\alpha}\in C^{\infty}(\Omega_{\alpha},{\mathbb{R}})$ and ${\mathbf{v}}_{\alpha}\in C^{\infty}(\Omega_{\alpha},{\mathbb{R}}^{p})$ such that $\lambda_{\alpha}(h\left({\mathbf{z}};{\mathcal{D}},i\right))=\bar{\lambda}_{\alpha}$ and ${\mathbf{v}}_{\alpha}(h\left({\mathbf{z}};{\mathcal{D}},i\right))=\bar{{\mathbf{v}}}_{\alpha}$ Let

[TABLE]

If $\bar{\lambda}_{\alpha}$ is positive for all $\alpha=1,\dots,p,$ every $\Omega_{\alpha}^{+}$ is an open set that contains $h\left({\mathbf{z}};{\mathcal{D}},i\right).$ Let

[TABLE]

which is nonempty and open. When restricted to $\Omega$, the mapping $\mathsf{\Pi}_{+}$ is the identity map, and thus (at least for this case) the proof follows if \psi_{ij}\left(\cdot,{\mathbf{z}},{\mathcal{D}}\right)\big{|}_{\Omega} is differentiable. To verify that it is, consider arbitrary $X\in\Omega$. Recall again the definition of $\psi_{ij}(\cdot,{\mathbf{z}},{\mathcal{D}})\colon{\rm Sym}(p,{\mathbb{R}})\to{\mathbb{R}}$ from (11). Note that $\nabla_{j}h({\mathbf{z}},{\mathcal{D}},i)$ is a constant matrix because $h$ is a linear function of ${\mathbf{z}}.$ Call this matrix $A_{ij}$. We have that

[TABLE]

which uses the same algebraic manipulation as the derivation of (35b). By inspection, $\psi_{ij}(\cdot,{\mathbf{z}},{\mathcal{D}})$ is differentiable at all $X\in\Omega$, and the result follows.

If $\bar{\lambda}_{\alpha}$ is negative for all $\alpha=1,\dots,p,$ every $\Omega_{\alpha}^{-}$ is an open set that contains $h\left({\mathbf{z}};{\mathcal{D}},i\right).$ Let

[TABLE]

which is again open and nonempty. Note that $\mathsf{\Pi}_{+}$ projects every matrix $X\in\Omega$ to the matrix of zeros. The map $\psi_{ij}({\mathbf{0}},{\mathbf{z}},{\mathcal{D}})$ is thus constant w.r.t. $X\in\Omega$ and is thus differentiable.

If $h\left({\mathbf{z}};{\mathcal{D}},i\right)$ has some positive and negative eigenvalues, then we can split them into two groups. Let $\alpha$ index the positive eigenvalues and $\beta$ denote the negative eigenvalues. Similar to before, let

[TABLE]

which is again nonempty, open, and contains $h\left({\mathbf{z}};{\mathcal{D}},i\right).$ Note that $\mathsf{\Pi}_{+}$ does not change $\lambda_{\alpha}$ or ${\mathbf{v}}_{\alpha}$ for matrices $X\in\Omega$, i.e.,

[TABLE]

The eigenvalues corresponding to $\beta$ of $\mathsf{\Pi}_{+}(X)$, on the other hand, are always projected to zero. Thus we write that $\lambda_{\beta}\left(\mathsf{\Pi}_{+}(X)\right)\equiv 0,$ and we have, omitting some elementary algebra, that

[TABLE]

By inspection, $\psi_{ij}\left(\mathsf{\Pi}_{+}(\cdot),{\mathbf{z}},{\mathcal{D}}\right)$ is differentiable at $X$. Since $X$ was an arbitrary matrix from the set $\Omega$, the result follows.

Since $h\left({\mathbf{z}};{\mathcal{D}},i\right)$ was assumed to be nonsingular, all possibilities of combinations eigenvalues have been covered, and, for each case, we found an open set $\Omega$ in which the desired result holds. ∎

Proof of Lemma V.6.

Recall the definition

[TABLE]

where $\left\{\omega_{s,k}\right\}_{k\in{\mathbb{N}}}\subset{\mathcal{I}}.$ Consider sequences of the form $\left\{i,i,i,\dots\right\}.$ We will use the fact that

[TABLE]

to show that that $\boldsymbol{\delta}_{s,k}$ is summable, provided that each sequence ${\left\|h^{\prime}_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},i\right)-h^{\prime}_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}},i\right)\right\|}$ is summable.

To show that $h^{\prime}_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},i\right)-h^{\prime}_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}},i\right)$ is summable, we will show that the sequence of the entries for an arbitrary entry $j\in\left\{1,\dots,m+qn\right\}$ is summable. If this is true, the norm of the whole vector will be summable as well because ${\left\|\cdot\right\|}_{a}\geq{\left\|\cdot\right\|}_{b}.$ if $a\leq b.$ Fix $j\in\left\{1,\dots,m+qn\right\}$.

Case 1 If $j\leq m,$ then this entry corresponds to the variable from the original problem that we denoted as $\gamma_{j}$. Note that, for this choice of $j$, we have that

[TABLE]

for all possible variables ${\mathbf{z}}$ or data ${\mathcal{D}}.$ Therefore, the $j$-th entry of $h^{\prime}_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},i\right)-h^{\prime}_{+}\left({\mathbf{z}}_{s,k};{\mathcal{D}},i\right)$, which is the difference of two functions

[TABLE]

because these functions are actually the same constant.

Case 2 If $j>m$, then we are looking at an entry of the gradient with respect to the control variables ${\mathbf{u}}_{s},$ which are non trivial. In this case, recall from the definition in (11) that $\nabla_{j}h_{+}({\mathbf{z}},{\mathcal{D}},i)=\psi_{ij}\left(\mathsf{\Pi}_{+}\left(h\left({\mathbf{z}};{\mathcal{D}},i\right)\right),{\mathbf{z}},{\mathcal{D}}\right).$ In other words, we can look at $\nabla_{j}h_{+}$ as essentially a composition of three functions: $\psi_{ij},\mathsf{\Pi}_{+},\;h.$ In the preceding lemmas and propositions, we have established useful facts about each of these constituent functions that will enable us to complete the proof.

Let $X_{k}=h\left({\mathbf{z}}_{s,k};{\mathcal{D}}_{s,k},i\right)$ and $Y_{k}=h\left({\mathbf{z}}_{s,k};{\mathcal{D}},i\right),$ and note that $X_{k}\to Y_{k}$ because $h$ is continuous and ${\mathcal{D}}_{s,k}\to{\mathcal{D}}.$ By Proposition A.5, both sequences $\left\{{\mathbf{z}}_{s,k},{\mathcal{D}}_{s,k}\right\}_{k\in{\mathbb{N}}}$ and $\left\{{\mathbf{z}}_{s,k},{\mathcal{D}}\right\}_{k\in{\mathbb{N}}}$ almost surely lie in the preimage of the set of nonsingular matrices that have distinct eigenvalues with respect to the mapping $h(\cdot;\cdot,i).$ Let $\Omega_{k}$ be an open neighborhood of $Y_{k}$ such that \psi_{ij}\left(\mathsf{\Pi}_{+}\left(\cdot\right),{\mathbf{z}},{\mathcal{D}}\right)\big{|}_{\Omega_{k}} is differentiable. Such a neighborhood is guaranteed to be available by Proposition A.9 as long as $Y_{k}$ converges, and the limit point is nonsingular and does not have a simple spectrum. We know that $Y_{k}$ converges by construction of the decreasing step sizes in Algorithm 1, and w.p.1. this limit point will be nonsingular and have a nonsimple spectrum by Proposition A.5. Then, we can find $K\in{\mathbb{N}}$ such that $\forall\ell\geq K,$ $X_{\ell}\in\Omega_{\ell}.$

By Proposition A.8, the derivatives of $\psi_{ij}\left(\mathsf{\Pi}_{+}\left(X\right),{\mathbf{z}},{\mathcal{D}}\right)$ with respect to both $X$ and ${\mathcal{D}}$ are bounded. Denote the larger bound by $\mu_{\ell}$. Proposition A.6 bounds $\psi_{ij}$ globally almost surely. We have that

[TABLE]

where $\langle\cdot,\cdot\rangle_{{\rm Sym}(p,{\mathbb{R}})}$ denotes the matrix inner product. The final term is summable if and only if ${\left\|X_{\ell}-Y_{\ell}\right\|}$ is summable. To verify that it is, note that

[TABLE]

which is summable by Proposition A.1.∎

Appendix B Proof of Lemmas V.7-V.9

Proof of Lemma V.7.

From Line 13 of Algorithm 1, the following chain of relations hold for any feasible point $\mathbf{z}\in\mathcal{X}\subseteq\mathcal{X}_{0}$ a.s.:

[TABLE]

where (46a) follows from the nonexpansiveness of the projection operator; (46b) follows from the expansion of $\|\cdot\|^{2}$; (46c) follows from relation (14) and the Schwarz inequality; and (46d) follows from the convexity of the function $h_{+}(\cdot;\cdot,\omega_{s,k})$, and $h_{+}(\mathbf{z};{\mathcal{D}},\omega_{s,k})=0$ due to the feasibility of the point $\mathbf{z}$, i.e., $h_{+}(\mathbf{v}_{s,k};{\mathcal{D}},\omega_{s,k})\leq\langle h^{\prime}_{+}(\mathbf{v}_{s,k};{\mathcal{D}},\omega_{s,k}),\mathbf{v}_{s,k}-\mathbf{z}\rangle.$ From relation (13) and the definition of $\beta_{s,k}$ in (12), we further have

[TABLE]

where the last inequality follows from relation (18a). Combining this with inequality (46d), we obtain for any $\mathbf{z}\in\mathcal{X}$ a.s.:

[TABLE]

For the last term on the right-hand side, we can rewrite $h_{+}(\mathbf{v}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})=h_{+}(\mathbf{v}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})-h_{+}(\mathbf{p}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})+h_{+}(\mathbf{p}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k}).$ Therefore,

[TABLE]

Regarding the first term on the right-hand side of (B), the following chain of relations holds:

[TABLE]

where (49a) follows from the Schwarz inequality; (49b) follows from relation (18b); (49c) follows from Line 10 of Algorithm 1; and (49d) follows from using $2|a||b|\leq\tau a^{2}+\frac{1}{\tau}b^{2}$ with $a=\alpha_{k}(L_{h}+D)L_{f}$ and $b=h_{+}(\mathbf{p}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})$, and $\tau>0$ is arbitrary. Using relations (B)-(49) in (47), we obtain for all $\mathbf{z}\in\mathcal{X}$ a.s.:

[TABLE]

Similarly to (46a)-(46d), from Line 10 of Algorithm 1, the following chain of relations hold for any $\mathbf{z}\in\mathcal{X}\subseteq\mathcal{X}_{0}$ a.s.:

[TABLE]

where (51a) follows from the nonexpansive projection; (51b) follows from the expansion of $\|\cdot\|^{2}$; (51c) follows from the convexity of the function $f$ and relation (16a). For the second term on the right-hand side, we further have for any $\check{\mathbf{z}}\in\mathcal{X}$:

[TABLE]

where (52a) follows from adding and subtracting $f(\check{\mathbf{z}})$; (52b) follows from the convexity of the function $f$; (52c) follows from the Schwarz inequality; and (52d) follows from using $2|a||b|\leq 4\eta a^{2}+\frac{1}{4\eta}b^{2}$ with $a=\alpha_{k}L_{f}$ and $b=\|\mathbf{p}_{s,k}-\check{\mathbf{z}}\|$. Combining (52d) in (51c), we have for any $\mathbf{z},\check{\mathbf{z}}\in\mathcal{X}\subseteq\mathcal{X}_{0}$ a.s.:

[TABLE]

Substituting this inequality in relation (50), we obtain

[TABLE]

where $A_{\eta,\tau}=1+4\eta+\tau$. From relation (15), we have $\|\mathbf{v}_{s,k}-\mathbf{z}\|\leq C_{\mathbf{z}}.$ Using this and the upper estimate for $\beta_{s,k}$ in (21) for bounding the three error terms completes the proof. ∎

Proof of Lemma V.8.

We use Lemma V.7 with $\check{\mathbf{z}}=\mathbf{z}=\mathsf{\Pi}_{\mathcal{X}}[\mathbf{p}_{s,k}]$. Therefore, for any $s\in{\mathcal{S}}$, and $k\geq 1$, we obtain a.s.:

[TABLE]

where $A_{\eta,\tau}=1+4\eta+\tau$ and $\eta,\tau>0$ are arbitrary. By the definition of the projection, we have $\mathrm{dist}(\mathbf{p}_{s,k},\mathcal{X})=\|\mathbf{p}_{s,k}-\mathsf{\Pi}_{\mathcal{X}}[\mathbf{p}_{s,k}]\|$ and $\mathrm{dist}(\mathbf{z}_{s,k},\mathcal{X})=\|\mathbf{z}_{s,k}-\mathsf{\Pi}_{\mathcal{X}}[\mathbf{z}_{s,k}]\|\leq\|\mathbf{z}_{s,k}-\mathsf{\Pi}_{\mathcal{X}}[\mathbf{p}_{s,k}]\|.$ Upon substituting these estimates in relation (53), we obtain

[TABLE]

Taking the expectation conditioned on $\mathcal{F}_{k-1}$ and noting that $\mathbf{p}_{s,k}$ is fully determined by $\mathcal{F}_{k-1}$, we have for any $s\in{\mathcal{S}}$ and $k\geq 1$ a.s.:

[TABLE]

We now choose $\tau=4$, $\eta=\kappa(L_{h}+D)^{2}$ and use Assumption IV.3 to yield

[TABLE]

where $A=5+4\kappa(L_{h}+D)^{2}$. Finally, by summing over all $s\in{\mathcal{S}}$ and using Lemma V.3 with $h(x)=\mathrm{dist}^{2}(x,\mathcal{X})$, we arrive at the following relation:

[TABLE]

Therefore, from $\sum_{k\geq 0}\alpha_{k}^{2}<\infty$ and Lemmas V.5 and V.6, all the conditions in the convergence theorem (Theorem V.2) are satisfied and the desired result follows. ∎

Proof of Lemma V.9(a).

From Line 6-10 and 13 of Algorithm 1, we define $\mathbf{e}_{s,k}\triangleq\mathbf{z}_{s,k}-\mathbf{p}_{s,k}$ for $s\in{\mathcal{S}}$ and $k\geq 0$. Hence, $\mathbf{e}_{s,k}$ can be viewed as the perturbation that we make on $\mathbf{p}_{s,k}$ after the network consensus step in Line 6 of Algorithm 1. Consider $\|\mathbf{e}_{s,k}\|$, for which we can write

[TABLE]

where the second inequality follows from fact that $\mathbf{p}_{s,k},\mathbf{v}_{s,k}\in\mathcal{X}_{0}$, the nonexpansiveness of the projection operator, relations (16a) and (20).

Applying $(a+b)^{2}\leq 2a^{2}+2b^{2}$ in inequality (B), we have for all $s\in{\mathcal{S}}$ and $k\geq 1$,

[TABLE]

For the term $h_{+}^{2}(\mathbf{v}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})$ in (58), we have the following chain of relations:

[TABLE]

where (59) follows from adding and subtracting $h_{+}(\mathbf{p}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})$; (60) follows from relation (18b); (61) follows from Line 10 of Algorithm 1; (62) follows from using the relation $2|a||b|\leq a^{2}+b^{2}$ with $a=\alpha_{k}(L_{h}+D)L_{f}$ and $b=h_{+}(\mathbf{p}_{s,k};{\mathcal{D}}_{s,k},\omega_{s,k})$. By combining (58) and (62), we obtain

[TABLE]

where the last inequality is from (17). Therefore, from Lemma V.8 and $\sum_{k=1}^{\infty}\alpha_{k}^{2}<\infty$, we conclude that $\sum_{k=1}^{\infty}\|\mathbf{e}_{s,k}\|^{2}<\infty$ for all $s\in{\mathcal{S}}\quad a.s.$ ∎

Proof of Lemma V.9(b).

The proof of this part of Lemma will coincide with [26, Lemma 6(b)]. ∎