Resource-Efficient Distributed Recursive Gaussian Processes

# Resource-Efficient Distributed Recursive Gaussian Processes
Source: [https://arxiv.org/html/2609.26979](https://arxiv.org/html/2609.26979)
\{IEEEkeywords\}

Gaussian processes, multi\-agent systems, distributed inference, kernel methods, sensor fusion

Student MemberIEEE Ali Emre Balcı1\{\}^\{\\textbf\{1\}\}Student MemberIEEE Raj Thilak Rajan1\{\}^\{\\textbf\{1\}\}Senior MemberIEEEAffiliation:Signal Processing Systems, Faculty of EEMCS, Delft University of Technology, Mekelweg 4, 2628CD, Delft, The Netherlands

###### Abstract

Gaussian processes \(GPs\) provide a flexible framework for learning unknown functions from noisy measurements while quantifying predictive uncertainty, making them well suited for estimation in multi\-agent systems\. However, when measurements are collected by multiple agents, maintaining a unified GP model without centralized processing requires efficient distributed algorithms that can operate using local measurements and communication with neighboring agents\. In this work, we develop two distributed recursive GP \(RGP\) algorithms for multi\-output GP regression: ADMM\-RGP and PDMM\-RGP\. We analyze the stability and convergence of both algorithms and develop parameter selection strategies to accelerate convergence, thus reducing the communication burden\. The proposed methods are validated on a real\-world multi\-output wind dataset, and their convergence behavior is examined across communication graphs with varying connectivity\. Numerical experiments demonstrate that ADMM\-RGP and PDMM\-RGP can significantly reduce communication relative to the state of the art, while maintaining comparable estimation accuracy and network\-wide consensus\.

††corresponding:Corresponding author: Ali Emre Balcı \(a\.e\.balci@tudelft\.nl\)††note:This work is partially funded by the Sensor AI Lab through the AI Labs Program of the Delft University of Technology, and EU\-HORIZON\-KDT\-JU\-2023\-2\-RIA under grant agreement No 101139996, the ShapeFuture project\.## 1Introduction

\\IEEEPARstart

Multi agent systems have gained significant attention in recent years across a variety of applications, including search\-and\-rescue\[[1](https://arxiv.org/html/2609.26979#bib.bib1)\], agricultural monitoring\[[2](https://arxiv.org/html/2609.26979#bib.bib2)\], smart grids\[[3](https://arxiv.org/html/2609.26979#bib.bib3)\], and logistics\[[4](https://arxiv.org/html/2609.26979#bib.bib4)\]\. A key advantage of these systems is their ability to simultaneously collect data from multiple locations\[[5](https://arxiv.org/html/2609.26979#bib.bib5)\]\. Exploiting these measurements requires agents to combine locally acquired information to estimate an underlying quantity of interest, such as an environmental field or an unknown system function\. In many applications, measurements arrive sequentially and are spatially distributed across the network, motivating learning methods that can update their estimates online while operating with local computation and sparse communication\[[6](https://arxiv.org/html/2609.26979#bib.bib6),[7](https://arxiv.org/html/2609.26979#bib.bib7)\]\.

Gaussian processes \(GPs\) are a Bayesian supervised machine learning method for reconstructing unknown functions from noisy measurements\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\]\. As a nonparametric method, GPs do not assume a fixed functional form, making them flexible and advantageous in applications where dynamics are unknown or unpredictable\[[9](https://arxiv.org/html/2609.26979#bib.bib9)\]\. GPs are particularly attractive when measurements are limited or noisy due to their ability to quantify predictive uncertainty and incorporate prior knowledge\[[10](https://arxiv.org/html/2609.26979#bib.bib10)\]\. These properties make GPs well suited for multi\-agent systems, which often rely on sensing, estimation, and uncertainty quantification\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\]\. GPs have proven useful in a wide variety of multi\-agent tasks, including wireless traffic prediction\[[11](https://arxiv.org/html/2609.26979#bib.bib11)\], safe multi\-agent path planning\[[12](https://arxiv.org/html/2609.26979#bib.bib12),[7](https://arxiv.org/html/2609.26979#bib.bib7)\], and cooperative mapping\[[6](https://arxiv.org/html/2609.26979#bib.bib6)\]\.

Applying GPs online in multi\-agent systems nevertheless presents computational and communication challenges\. Centralized approaches require agents to transmit their measurements to a common processing unit, creating a single point of failure, while standard GP regression scales cubically with the number of observations\[[5](https://arxiv.org/html/2609.26979#bib.bib5)\]\. Communication constraints present an additional challenge, which is particularly important in unmanned aerial vehicles \(UAVs\), multi\-robot systems \(MRSs\), and wireless sensor networks \(WSNs\), where communication can constitute a significant fraction of the available energy budget\[[13](https://arxiv.org/html/2609.26979#bib.bib13),[14](https://arxiv.org/html/2609.26979#bib.bib14),[15](https://arxiv.org/html/2609.26979#bib.bib15)\]\.

Prior works on distributed GPs have largely focused on kernel hyperparameter training e\.g\.,\[[16](https://arxiv.org/html/2609.26979#bib.bib16),[17](https://arxiv.org/html/2609.26979#bib.bib17)\]propose proximal\-ADMM approaches for distributed hyperparameter training, while\[[5](https://arxiv.org/html/2609.26979#bib.bib5)\]introduce dataset augmentation to promote hyperparameter consensus\. Although these methods enable distributed GP training, they do not address the cubic computational complexity of GPs or support online regression in distributed settings\. The recursive GP \(RGP\) algorithm\[[18](https://arxiv.org/html/2609.26979#bib.bib18)\]was proposed to resolve the complexity problem by enabling GP regression from streaming data using a sparse, inducing point\-based approximation and Kalman filter\-like updates\. For multi\-agent networks, Consensus\-RGP\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\]extended RGP to model and predict multi\-output functions using distributed consensus algorithms over a network\. In a similar setting, D\-RF\-GP\[[19](https://arxiv.org/html/2609.26979#bib.bib19)\]adopts random Fourier features and uses ensembles of models to adapt to changing kernel length scales\. The more recent ROAD\-GP\[[20](https://arxiv.org/html/2609.26979#bib.bib20)\]further extends D\-RF\-GP to time\-varying functions and improved its robustness to outliers\.

Despite this progress, the communication and consensus efficiency of distributed online GP methods remain comparatively unexplored\. Motivated by this gap, we develop two communication\-efficient distributed recursive GP algorithms for multi\-output regression that we call ADMM\-RGP and PDMM\-RGP\. ADMM\-RGP adapts a communication\-efficient alternating direction method of multipliers \(ADMM\) scheme originally developed for distributed Kalman filtering\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\], while PDMM\-RGP applies the primal\-dual method of multipliers \(PDMM\) algorithm\[[22](https://arxiv.org/html/2609.26979#bib.bib22)\]to the distributed RGP problem\. Compared with the state\-of\-the\-art Consensus\-RGP\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\], we show that ADMM\-RGP significantly reduces communication while maintaining comparable accuracy and network\-wide consensus, whereas PDMM\-RGP further reduces communication and is particularly effective in sparse communication networks\. The main contributions of this work are as follows:

- •We develop ADMM\-RGP and PDMM\-RGP for distributed, recursive, multi\-output GP regression\.
- •We analyze the stability of both algorithms and provide parameter selection strategies for fast convergence\. With the proposed parameter selection, ADMM\-RGP is shown to converge faster than the state\-of\-the\-art\.
- •We evaluate the proposed methods on a real\-world, multi\-output wind dataset\[[23](https://arxiv.org/html/2609.26979#bib.bib23)\]to demonstrate their effectiveness in learning an unknown function\.
- •We compare the communication and computational complexity of the proposed methods with the state of the art across graphs of varying connectivity\.

In this study, our key focus lies in the distributed fusion layer rather than in the local GP approximation\. Consequently, ADMM\-RGP and PDMM\-RGP are not tied to a specific basis construction or kernel family\. Building on the inducing point formulation of Consensus\-RGP\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\], we isolate the communication efficiency aspect of distributed inference, and thus enable a controlled comparison with the state\-of\-the\-art methods\. The proposed framework can be readily extended to alternative GP representations such as random Fourier features\. Moreover, the inducing point formulation adopted here is kernel agnostic, and thus can accommodate a broader class of covariance functions\.

Notation:𝐈N∈ℝN×N\\mathbf\{I\}\_\{N\}\\in\\mathbb\{R\}^\{N\\times N\}denotes the identity matrix, and𝟏N∈ℝN\\mathbf\{1\}\_\{N\}\\in\\mathbb\{R\}^\{N\}denotes a vector of ones\. Theiith entry of a vector𝐱\\mathbf\{x\}is denoted by\[𝐱\]i\[\\mathbf\{x\}\]\_\{i\}, and the\(i,j\)\(i,j\)th entry of a matrix𝐗\\mathbf\{X\}is denoted by\[𝐗\]i​j\[\\mathbf\{X\}\]\_\{ij\}\. Probability density functions are denoted byp⁡\(⋅\)p\(\\cdot\), and the multivariate normal distribution with mean𝝁\\boldsymbol\{\\mu\}and covariance𝚺\\boldsymbol\{\\Sigma\}is denoted by𝒩⁡\(𝝁,𝚺\)\\mathcal\{N\}\(\\boldsymbol\{\\mu\},\\boldsymbol\{\\Sigma\}\)\. The Kronecker product is denoted by⊗\\otimes\. The set of eigenvalues of a matrix𝐗\\mathbf\{X\}is denoted byλ⁡\(𝐗\)\\lambda\(\\mathbf\{X\}\), whileλi\\lambda\_\{i\}denotes itsiith smallest eigenvalue, such thatλ1≤⋯≤λN\\lambda\_\{1\}\\leq\\cdots\\leq\\lambda\_\{N\}\. Finally,∥⋅∥\\\|\\cdot\\\|denotes the Euclidean norm for vectors and the spectral norm for matrices, while∥⋅∥F\\\|\\cdot\\\|\_\{F\}denotes the Frobenius norm\.

## 2Preliminaries

### 2\.1Multi\-Output Gaussian Processes

Consider a function𝐟:ℝD→ℝD′\\mathbf\{f\}:\\mathbb\{R\}^\{D\}\\rightarrow\\mathbb\{R\}^\{D^\{\\prime\}\}that may be observed with some noise

𝐲j=𝐟⁡\(𝐱j\)\+𝐯j,𝐯j∼𝒩⁡\(𝟎,𝐑v\),\\mathbf\{y\}\_\{j\}=\\mathbf\{f\}\(\\mathbf\{x\}\_\{j\}\)\+\\mathbf\{v\}\_\{j\},\\quad\\mathbf\{v\}\_\{j\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{R\}\_\{v\}\),\(1\)where𝐲j\\mathbf\{y\}\_\{j\}is the noisy measurement at input location𝐱j\\mathbf\{x\}\_\{j\}and𝐑v\\mathbf\{R\}\_\{v\}is the measurement noise covariance matrix\. With a collection of measurements we form a dataset𝒟≜\{𝐗,𝐲\}\\mathcal\{D\}\\triangleq\\\{\\mathbf\{X\},\\mathbf\{y\}\\\}, where𝐗≜\[𝐱1⋯𝐱J\]∈ℝD×J\\mathbf\{X\}\\triangleq\[\\mathbf\{x\}\_\{1\}\\ \\cdots\\ \\mathbf\{x\}\_\{J\}\]\\in\\mathbb\{R\}^\{D\\times J\}is the matrix of measurement locations and𝐲≜\[𝐲1⊤⋯𝐲J⊤\]⊤∈ℝD′​J\\mathbf\{y\}\\triangleq\[\\mathbf\{y\}\_\{1\}^\{\\top\}\\ \\cdots\\ \\mathbf\{y\}\_\{J\}^\{\\top\}\]^\{\\top\}\\in\\mathbb\{R\}^\{D^\{\\prime\}J\}is the corresponding vector of measurements\.

In particular, forD′=1D^\{\\prime\}=1, the unknown function may be modeled using a single\-output Gaussian process denoted as𝒢​𝒫​\(m⁡\(𝐱\),κ⁡\(𝐱,𝐱′\)\)\\mathcal\{GP\}\(m\(\\mathbf\{x\}\),\\kappa\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\)which is characterized by a mean function,m⁡\(𝐱\)m\(\\mathbf\{x\}\), and a kernel function,κ⁡\(𝐱,𝐱′\):ℝD×ℝD→ℝ\\kappa\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\):\\mathbb\{R\}^\{D\}\\times\\mathbb\{R\}^\{D\}\\rightarrow\\mathbb\{R\}, where the kernel captures the smoothness and range of values that the function is expected to exhibit\. The kernel function can be used to construct a kernel matrix\[𝐊⁡\(𝐗,𝐗′\)\]i​j≜κ⁡\(𝐱i,𝐱j′\)\[\\mathbf\{K\}\(\\mathbf\{X\},\\mathbf\{X\}^\{\\prime\}\)\]\_\{ij\}\\triangleq\\kappa\(\\mathbf\{x\}\_\{i\},\\mathbf\{x\}\_\{j\}^\{\\prime\}\)for two sets of input locations \(OPEN𝐗,𝐗′\)\\mathbf\{X\},\\mathbf\{X\}^\{\\prime\}\)\. The predicted distribution𝐟⁡\(𝐗∗\)∼𝒩⁡\(𝝁∗,𝚺∗\)\\mathbf\{f\}\(\\mathbf\{X\}\_\{\*\}\)\\sim\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{\*\},\\boldsymbol\{\\Sigma\}\_\{\*\}\)at a set of input locations𝐗∗\\mathbf\{X\}\_\{\*\}is then

𝝁∗\\displaystyle\\boldsymbol\{\\mu\}\_\{\*\}=𝐊⁡\(𝐗∗,𝐗\)​𝐊𝐲−1​𝐲,\\displaystyle=\\mathbf\{K\}\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)\\mathbf\{K\}\_\{\\mathbf\{y\}\}^\{\-1\}\\mathbf\{y\},\(2a\)𝚺∗\\displaystyle\\boldsymbol\{\\Sigma\}\_\{\*\}=𝐊⁡\(𝐗∗,𝐗∗\)−𝐊⁡\(𝐗∗,𝐗\)​𝐊𝐲−1​𝐊​\(𝐗,𝐗∗\),\\displaystyle=\\mathbf\{K\}\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{\*\}\)\-\\mathbf\{K\}\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)\\mathbf\{K\}\_\{\\mathbf\{y\}\}^\{\-1\}\\mathbf\{K\}\(\\mathbf\{X\},\\mathbf\{X\}\_\{\*\}\),\(2b\)where𝐊𝐲≜𝐊⁡\(𝐗,𝐗\)\+𝐈J⊗𝐑v\\mathbf\{K\}\_\{\\mathbf\{y\}\}\\triangleq\\mathbf\{K\}\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\mathbf\{I\}\_\{J\}\\otimes\\mathbf\{R\}\_\{v\}\. For multi\-output functions in general, the linear model of coregionalization \(LMC\) may be used to capture correlations between different output dimensions\[[24](https://arxiv.org/html/2609.26979#bib.bib24)\]\. In the LMC approach, each output dimension is represented as a linear combination ofQ≤D′Q\\leq D^\{\\prime\}latent single\-output functions,\{uq:ℝD→ℝ\}q=1Q\\\{u\_\{q\}:\\mathbb\{R\}^\{D\}\\rightarrow\\mathbb\{R\}\\\}\_\{q=1\}^\{Q\}\. Each latent function has a zero\-mean GP prior,uq​\(𝐱\)∼𝒢​𝒫​\(0,κq​\(𝐱,𝐱′\)\)u\_\{q\}\(\\mathbf\{x\}\)\\sim\\mathcal\{GP\}\(0,\\kappa\_\{q\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\), and a vector of weights𝐛q≜\[b1​q⋯bD′​q\]⊤\\mathbf\{b\}\_\{q\}\\triangleq\[b\_\{1q\}\\ \\cdots\\ b\_\{D^\{\\prime\}q\}\]^\{\\top\}, where𝐛i​q\\mathbf\{b\}\_\{iq\}is associated with output dimensionii\. The covariance of the outputs at two input locations is given by the matrix\-valued kernel function,𝓚⁡\(𝐱,𝐱′\)≜cov​\(𝐟⁡\(𝐱\),𝐟⁡\(𝐱′\)\)=∑q=1Qκq​\(𝐱,𝐱′\)​𝐁q\\boldsymbol\{\\mathcal\{K\}\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\\triangleq\\text\{cov\}\(\\mathbf\{f\}\(\\mathbf\{x\}\),\\mathbf\{f\}\(\\mathbf\{x\}^\{\\prime\}\)\)=\\sum\_\{q=1\}^\{Q\}\\kappa\_\{q\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\\mathbf\{B\}\_\{q\}, where𝐁q≜𝐛q​𝐛q⊤\\mathbf\{B\}\_\{q\}\\triangleq\\mathbf\{b\}\_\{q\}\\mathbf\{b\}\_\{q\}^\{\\top\}is the coregionalization matrix for latent functionqq\. The LMC block kernel matrix for two sets of input locations \(OPEN𝐗,𝐗′\)\\mathbf\{X\},\\mathbf\{X\}^\{\\prime\}\)is thus defined as\[𝐊⁡\(𝐗,𝐗′\)\]i​j≜𝓚⁡\(𝐱i,𝐱j′\)\.\[\\mathbf\{K\}\(\\mathbf\{X\},\\mathbf\{X\}^\{\\prime\}\)\]\_\{ij\}\\triangleq\\boldsymbol\{\\mathcal\{K\}\}\(\\mathbf\{x\}\_\{i\},\\mathbf\{x\}\_\{j\}^\{\\prime\}\)\.

### 2\.2Recursive Gaussian Processes

The recursive Gaussian process \(RGP\) algorithm provides an online estimate of latent function values𝐳≜𝐟⁡\(𝐗p\)\\mathbf\{z\}\\triangleq\\mathbf\{f\}\(\\mathbf\{X\}\_\{p\}\)at a fixed set of basis locations𝐗p∈ℝD×P\\mathbf\{X\}\_\{p\}\\in\\mathbb\{R\}^\{D\\times P\}\[[18](https://arxiv.org/html/2609.26979#bib.bib18)\]\. At each time steptt, a batch ofBt<JB\_\{t\}<Jnoisy measurements\{𝐗t,𝐲t\}\\\{\\mathbf\{X\}\_\{t\},\\mathbf\{y\}\_\{t\}\\\}is collected according to the measurement model in \([1](https://arxiv.org/html/2609.26979#S2.E1)\) and, using this batch and the latent prior𝐳∣𝐲1:t−1∼𝒩\(𝝁t−1,𝚺t−1\),\\mathbf\{z\}\\mid\\mathbf\{y\}\_\{1:t\-1\}\\sim\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{t\-1\},\\mathbf\{\\Sigma\}\_\{t\-1\}\),RGP obtains the posterior𝐳∣𝐲1:t∼𝒩\(𝝁t,𝚺t\)\\mathbf\{z\}\\mid\\mathbf\{y\}\_\{1:t\}\\sim\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{t\},\\mathbf\{\\Sigma\}\_\{t\}\), where𝐲1:t≜\{𝐲i\}i=1t\\mathbf\{y\}\_\{1:t\}\\triangleq\\\{\\mathbf\{y\}\_\{i\}\\\}\_\{i=1\}^\{t\}\. Given the prior distribution over the basis function values, the predictive distribution at the new measurement locations𝐟t≜𝐟⁡\(𝐗t\)∼𝒩⁡\(𝝁𝐟t,𝚺𝐟t\)\\mathbf\{f\}\_\{t\}\\triangleq\\mathbf\{f\}\(\\mathbf\{X\}\_\{t\}\)\\sim\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{\\mathbf\{f\}\_\{t\}\},\\mathbf\{\\Sigma\}\_\{\\mathbf\{f\}\_\{t\}\}\)is a Gaussian with

𝝁𝐟t\\displaystyle\\boldsymbol\{\\mu\}\_\{\\mathbf\{f\}\_\{t\}\}=𝐇t​𝝁t−1,\\displaystyle=\\mathbf\{H\}\_\{t\}\\boldsymbol\{\\mu\}\_\{t\-1\},\(3a\)𝚺𝐟t\\displaystyle\\boldsymbol\{\\Sigma\}\_\{\\mathbf\{f\}\_\{t\}\}=𝐊⁡\(𝐗t,𝐗t\)−𝐇t​𝐊​\(𝐗p,𝐗t\)\+𝐇t​𝚺t−1​𝐇t⊤,\\displaystyle=\\mathbf\{K\}\(\\mathbf\{X\}\_\{t\},\\mathbf\{X\}\_\{t\}\)\-\\mathbf\{H\}\_\{t\}\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{t\}\)\+\\mathbf\{H\}\_\{t\}\\mathbf\{\\Sigma\}\_\{t\-1\}\\mathbf\{H\}\_\{t\}^\{\\top\},\(3b\)𝐇t\\displaystyle\\mathbf\{H\}\_\{t\}≜𝐊⁡\(𝐗t,𝐗p\)​𝐊𝐗p−1,\\displaystyle\\triangleq\\mathbf\{K\}\(\\mathbf\{X\}\_\{t\},\\mathbf\{X\}\_\{p\}\)\\mathbf\{K\}\_\{\\mathbf\{X\}\_\{p\}\}^\{\-1\},\(3c\)𝐊𝐗p\\displaystyle\\mathbf\{K\}\_\{\\mathbf\{X\}\_\{p\}\}≜𝐊⁡\(𝐗p,𝐗p\)\+ϵ​𝐈,\\displaystyle\\triangleq\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{p\}\)\+\\epsilon\\mathbf\{I\},\(3d\)whereϵ\>0\\epsilon\>0is a small constant to ensure stability of the inversion\. The posterior mean and covariance can be calculated using a Kalman\-like correction update as

𝝁t\\displaystyle\\boldsymbol\{\\mu\}\_\{t\}=𝝁t−1\+𝐆t​\(𝐲t−𝐇t​𝝁t−1\),\\displaystyle=\\boldsymbol\{\\mu\}\_\{t\-1\}\+\\mathbf\{G\}\_\{t\}\(\\mathbf\{y\}\_\{t\}\-\\mathbf\{H\}\_\{t\}\\boldsymbol\{\\mu\}\_\{t\-1\}\),\(4a\)𝚺t\\displaystyle\\mathbf\{\\Sigma\}\_\{t\}=\(𝐈−𝐆t​𝐇t\)​𝚺t−1,\\displaystyle=\(\\mathbf\{I\}\-\\mathbf\{G\}\_\{t\}\\mathbf\{H\}\_\{t\}\)\\mathbf\{\\Sigma\}\_\{t\-1\},\(4b\)𝐆t\\displaystyle\\mathbf\{G\}\_\{t\}≜𝚺t−1​𝐇t⊤​\(𝐑t\+𝐇t​𝚺t−1​𝐇t⊤\)−1,\\displaystyle\\triangleq\\mathbf\{\\Sigma\}\_\{t\-1\}\\mathbf\{H\}\_\{t\}^\{\\top\}\(\\mathbf\{R\}\_\{t\}\+\\mathbf\{H\}\_\{t\}\\mathbf\{\\Sigma\}\_\{t\-1\}\\mathbf\{H\}\_\{t\}^\{\\top\}\)^\{\-1\},\(4c\)𝐑t\\displaystyle\\mathbf\{R\}\_\{t\}≜𝐊⁡\(𝐗t,𝐗t\)−𝐇t​𝐊​\(𝐗p,𝐗t\)\+𝐈Bt⊗𝐑v\.\\displaystyle\\triangleq\\mathbf\{K\}\(\\mathbf\{X\}\_\{t\},\\mathbf\{X\}\_\{t\}\)\-\\mathbf\{H\}\_\{t\}\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{t\}\)\+\\mathbf\{I\}\_\{B\_\{t\}\}\\otimes\\mathbf\{R\}\_\{v\}\.\(4d\)

### 2\.3Distributed Recursive Gaussian Processes

The distributed RGP problem considers a multi\-agent network ofN≥2N\\geq 2agents, whose communication topology is represented by a connected graph𝒢=\(𝒱,ℰ\)\\mathcal\{G\}=\(\\mathcal\{V\},\\mathcal\{E\}\), with vertices𝒱=\{1,⋯,N\}\\mathcal\{V\}=\\\{1,\\cdots,N\\\}representing the agents and edgesℰ⊆𝒱×𝒱\\mathcal\{E\}\\subseteq\\mathcal\{V\}\\times\\mathcal\{V\}representing communication links\. The set𝒩n≜\{m∈𝒱:\(n,m\)∈ℰ\}\\mathcal\{N\}\_\{n\}\\triangleq\\\{m\\in\\mathcal\{V\}:\(n,m\)\\in\\mathcal\{E\}\\\}represents the neighbors of nodenn\. At each time steptt, each agentnncollects its own batch ofBn,tB\_\{n,t\}independent noisy measurements\{𝐗n,t,𝐲n,t\}\\\{\\mathbf\{X\}\_\{n,t\},\\mathbf\{y\}\_\{n,t\}\\\}according to \([1](https://arxiv.org/html/2609.26979#S2.E1)\)\. A set ofPPbasis locations,𝐗p\\mathbf\{X\}\_\{p\}, are fixed and known to each agent\.

The objective is for each agent to obtain the RGP estimate that would result from jointly processing all measurements collected across the network\. Following the RGP formulation, the Gaussian distribution can be equivalently represented using the global information vector𝝃g,t\\boldsymbol\{\\xi\}\_\{g,t\}and global information matrix𝛀g,t\\mathbf\{\\Omega\}\_\{g,t\}, given by

𝝃g,t\\displaystyle\\boldsymbol\{\\xi\}\_\{g,t\}=𝚺t−1​𝝁t,𝛀g,t=𝚺t−1\.\\displaystyle=\\boldsymbol\{\\Sigma\}\_\{t\}^\{\-1\}\\boldsymbol\{\\mu\}\_\{t\},\\quad\\mathbf\{\\Omega\}\_\{g,t\}=\\boldsymbol\{\\Sigma\}\_\{t\}^\{\-1\}\.\(5\)Assuming the measurements at different agents are independent, the information\-form global RGP updates are

𝝃g,t\\displaystyle\\boldsymbol\{\\xi\}\_\{g,t\}=𝝃g,t−1\+∑n=1N𝐇n,t⊤​𝐑n,t−1​𝐲n,t,\\displaystyle=\\boldsymbol\{\\xi\}\_\{g,t\-1\}\+\\sum\_\{n=1\}^\{N\}\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{y\}\_\{n,t\},\(6\)𝛀g,t\\displaystyle\\mathbf\{\\Omega\}\_\{g,t\}=𝛀g,t−1\+∑n=1N𝐇n,t⊤​𝐑n,t−1​𝐇n,t,\\displaystyle=\\mathbf\{\\Omega\}\_\{g,t\-1\}\+\\sum\_\{n=1\}^\{N\}\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{H\}\_\{n,t\},\(7\)where𝐇n,t\\mathbf\{H\}\_\{n,t\}and𝐑n,t\\mathbf\{R\}\_\{n,t\}are defined by \([3c](https://arxiv.org/html/2609.26979#S2.E3.3)\) and \([4d](https://arxiv.org/html/2609.26979#S2.E4.4)\) with the measurement batch collected by agentnn\. Although separable, \([6](https://arxiv.org/html/2609.26979#S2.E6)\)\-\([7](https://arxiv.org/html/2609.26979#S2.E7)\) require measurement contributions from all agents, which is not feasible in a distributed system with no central processing station\. Each agentnntherefore maintains a local information vector𝝃n,t\\boldsymbol\{\\xi\}\_\{n,t\}and an information matrix𝛀n,t\\mathbf\{\\Omega\}\_\{n,t\}and exchanges information with neighboring agents to approximate the global solution\.

### 2\.4State of the Art: Consensus\-RGP

The Consensus\-RGP algorithm extends the RGP framework to model multi\-output functions in distributed, multi\-agent systems\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\]\. At each time steptt, each agentnnapplies updates \([4a](https://arxiv.org/html/2609.26979#S2.E4.1)\) and \([4b](https://arxiv.org/html/2609.26979#S2.E4.2)\) in the information form with their local set of measurements

𝝃n,t\\displaystyle\\boldsymbol\{\\xi\}\_\{n,t\}=𝝃n,t−1\+𝐇n,t⊤​𝐑n,t−1​𝐲n,t,\\displaystyle=\\boldsymbol\{\\xi\}\_\{n,t\-1\}\+\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{y\}\_\{n,t\},\(8a\)𝛀n,t\\displaystyle\\mathbf\{\\Omega\}\_\{n,t\}=𝛀n,t−1\+𝐇n,t⊤​𝐑n,t−1​𝐇n,t,\\displaystyle=\\mathbf\{\\Omega\}\_\{n,t\-1\}\+\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{H\}\_\{n,t\},\(8b\)where𝐇n,t≜𝐊⁡\(𝐗n,t,𝐗p\)​𝐊𝐗p−1\\mathbf\{H\}\_\{n,t\}\\triangleq\\mathbf\{K\}\(\\mathbf\{X\}\_\{n,t\},\\mathbf\{X\}\_\{p\}\)\\mathbf\{K\}\_\{\\mathbf\{X\}\_\{p\}\}^\{\-1\}and𝐑n,t≜𝐊⁡\(𝐗n,t,𝐗n,t\)−𝐇n,t​𝐊​\(𝐗p,𝐗n,t\)\+𝐈Bn,t⊗𝐑v\\mathbf\{R\}\_\{n,t\}\\triangleq\\mathbf\{K\}\(\\mathbf\{X\}\_\{n,t\},\\mathbf\{X\}\_\{n,t\}\)\-\\mathbf\{H\}\_\{n,t\}\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{n,t\}\)\+\\mathbf\{I\}\_\{B\_\{n,t\}\}\\otimes\\mathbf\{R\}\_\{v\}\. After the local updates, the agents performKKrounds of distributed average consensus\. In thekt​hk^\{th\}round, agentnnupdates its information parameters as

𝝃n,t\[k\]=∑m=1Nwn​m​𝝃m,t\[k−1\],𝛀n,t\[k\]=∑m=1Nwn​m​𝛀m,t\[k−1\],\\displaystyle\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\}=\\sum\_\{m=1\}^\{N\}w\_\{nm\}\\boldsymbol\{\\xi\}\_\{m,t\}^\{\[k\-1\]\},\\quad\\boldsymbol\{\\Omega\}\_\{n,t\}^\{\[k\]\}=\\sum\_\{m=1\}^\{N\}w\_\{nm\}\\boldsymbol\{\\Omega\}\_\{m,t\}^\{\[k\-1\]\},\(9\)wherewn​mw\_\{nm\}are the row\-stochastic consensus weights\. AfterKKrounds of average consensus, the global information form parameters𝝃g,t\\boldsymbol\{\\xi\}\_\{g,t\}and𝛀g,t\\boldsymbol\{\\Omega\}\_\{g,t\}can be approximated as

𝝃g,t\\displaystyle\\boldsymbol\{\\xi\}\_\{g,t\}=N​𝝃n,t\[K\],\\displaystyle=N\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[K\]\},\(10a\)𝛀g,t\\displaystyle\\mathbf\{\\Omega\}\_\{g,t\}=𝛀n,0\[K\]\+N⁡\(𝛀n,t\[K\]−𝛀n,0\)\.\\displaystyle=\\mathbf\{\\Omega\}\_\{n,0\}^\{\[K\]\}\+N\(\\boldsymbol\{\\Omega\}\_\{n,t\}^\{\[K\]\}\-\\boldsymbol\{\\Omega\}\_\{n,0\}\)\.\(10b\)

## 3Proposed Approaches

### 3\.1ADMM\-RGP

Inspired by a communication\-efficient distributed ADMM variant of the Distributed Kalman Filter \(DKF\)\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\], we propose a novel distributed ADMM\-RGP algorithm that presents the updates directly in information form, simplifying the convergence analysis\. Let𝐋⪰𝟎\\mathbf\{L\}\\succeq\\mathbf\{0\}be the Laplacian matrix of the communication graph\. To derive distributed information vector updates using the ADMM\-based approach of\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\], we formulate the separable average consensus optimization problem

###### Problem 3\.1\.

ξtminimize

1 2∑\_n=1^N∥χ\_n,t \-ξ\_n,t∥^2, subject toL\_ξξ\_t =0,

where𝝃t≜\[𝝃1,t⊤⋯𝝃N,t⊤\]⊤\\boldsymbol\{\\xi\}\_\{t\}\\triangleq\[\\boldsymbol\{\\xi\}\_\{1,t\}^\{\\top\}\\ \\cdots\\ \\boldsymbol\{\\xi\}\_\{N,t\}^\{\\top\}\]^\{\\top\},𝝌n,t≜N​𝐇n,t⊤​𝐑n,t−1​𝐲n,t\+𝝃n,t−1\\boldsymbol\{\\chi\}\_\{n,t\}\\triangleq N\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{y\}\_\{n,t\}\+\\boldsymbol\{\\xi\}\_\{n,t\-1\},𝐋ξ≜\(𝐋⊗𝐈Dξ\)\\mathbf\{L\}\_\{\\xi\}\\triangleq\(\\mathbf\{L\}\\otimes\\mathbf\{I\}\_\{D\_\{\\xi\}\}\), andDξ≜P​D′D\_\{\\xi\}\\triangleq PD^\{\\prime\}\. The equality constraint ensures consensus among all agents, as the nullspace of𝐋\\mathbf\{L\}is the consensus subspace, spanned by𝟏N\\mathbf\{1\}\_\{N\}\. This problem has solution

𝝃t∗=𝟏N⊗\(∑n=1N𝐇n,t⊤​𝐑n,t−1​𝐲n,t\+1N​∑n=1N𝝃n,t−1\)\.\\boldsymbol\{\\xi\}\_\{t\}^\{\*\}=\\mathbf\{1\}\_\{N\}\\otimes\\Big\(\\sum\_\{n=1\}^\{N\}\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{y\}\_\{n,t\}\+\\frac\{1\}\{N\}\\sum\_\{n=1\}^\{N\}\\boldsymbol\{\\xi\}\_\{n,t\-1\}\\Big\)\.\(11\)If the agents reach consensus on the information vector at each time step, then𝝃n,t=𝝃g,t\\boldsymbol\{\\xi\}\_\{n,t\}=\\boldsymbol\{\\xi\}\_\{g,t\}for alln∈𝒱n\\in\\mathcal\{V\}, as desired\. To solve[3\.1](https://arxiv.org/html/2609.26979#Thmtheorem1), we define the dual variable at nodenn, timettas𝝀n,t∈ℝDξ\\boldsymbol\{\\lambda\}\_\{n,t\}\\in\\mathbb\{R\}^\{D\_\{\\xi\}\}and let𝝀t≜\[𝝀1,t⊤⋯𝝀N,t⊤\]⊤\\boldsymbol\{\\lambda\}\_\{t\}\\triangleq\[\\boldsymbol\{\\lambda\}\_\{1,t\}^\{\\top\}\\ \\cdots\\ \\boldsymbol\{\\lambda\}\_\{N,t\}^\{\\top\}\]^\{\\top\}\. The augmented Lagrangian with penalty parameterτ\\tauis then

ℒτ​\(𝝃t,𝝀t\)=12​‖𝝌t−𝝃t‖22\+𝝀t⊤​𝐁ξ​𝝃t\+τ2​‖𝐁ξ​𝝃t‖22,\\displaystyle\\mathcal\{L\}\_\{\\tau\}\(\\boldsymbol\{\\xi\}\_\{t\},\\boldsymbol\{\\lambda\}\_\{t\}\)=\\frac\{1\}\{2\}\\\|\\boldsymbol\{\\chi\}\_\{t\}\-\\boldsymbol\{\\xi\}\_\{t\}\\\|\_\{2\}^\{2\}\+\\boldsymbol\{\\lambda\}\_\{t\}^\{\\top\}\\mathbf\{B\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}\+\\frac\{\\tau\}\{2\}\|\|\\mathbf\{B\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}\|\|\_\{2\}^\{2\},\(12\)where𝐁ξ\\mathbf\{B\}\_\{\\xi\}is the symmetric matrix square root of𝐋ξ\\mathbf\{L\}\_\{\\xi\}and𝝌t≜\[𝝌1,t⊤⋯𝝌N,t⊤\]⊤\\boldsymbol\{\\chi\}\_\{t\}\\triangleq\[\\boldsymbol\{\\chi\}\_\{1,t\}^\{\\top\}\\ \\cdots\\ \\boldsymbol\{\\chi\}\_\{N,t\}^\{\\top\}\]^\{\\top\}\. Differentiating \([12](https://arxiv.org/html/2609.26979#S3.E12)\) with respect to𝝃t\\boldsymbol\{\\xi\}\_\{t\}and𝝀t\\boldsymbol\{\\lambda\}\_\{t\}gives the following gradients\.

∇𝝃tℒτ\\displaystyle\\nabla\_\{\\boldsymbol\{\\xi\}\_\{t\}\}\\mathcal\{L\}\_\{\\tau\}=𝝃t−𝝌t\+𝐁ξ​𝝀t\+τ​𝐋ξ​𝝃t,\\displaystyle=\\boldsymbol\{\\xi\}\_\{t\}\-\\boldsymbol\{\\chi\}\_\{t\}\+\\mathbf\{B\}\_\{\\xi\}\\boldsymbol\{\\lambda\}\_\{t\}\+\\tau\\mathbf\{L\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\},\(13a\)∇𝝀tℒτ\\displaystyle\\nabla\_\{\\boldsymbol\{\\lambda\}\_\{t\}\}\\mathcal\{L\}\_\{\\tau\}=𝐁ξ​𝝃t\.\\displaystyle=\\mathbf\{B\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}\.\(13b\)Definingα\>0\\alpha\>0as the step size of the dual variable ascent, the primal and dual variable updates at loop iterationkkare

𝝀~t\[k\]\\displaystyle\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{t\}^\{\[k\]\}=𝝀~t\[k−1\]\+α​𝐋ξ​𝝃t\[k−1\],\\displaystyle=\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{t\}^\{\[k\-1\]\}\+\\alpha\\mathbf\{L\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\-1\]\},\(14a\)𝝃t\[k\]\\displaystyle\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}=𝝌t−𝝀~t\[k\]−τ​𝐋ξ​𝝃t\[k−1\],\\displaystyle=\\boldsymbol\{\\chi\}\_\{t\}\-\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{t\}^\{\[k\]\}\-\\tau\\mathbf\{L\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\-1\]\},\(14b\)where𝝀~t≜𝐁ξ​𝝀t\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{t\}\\triangleq\\mathbf\{B\}\_\{\\xi\}\\boldsymbol\{\\lambda\}\_\{t\}is an auxiliary dual variable\. Letln​m≜\[𝐋\]n​ml\_\{nm\}\\triangleq\[\\mathbf\{L\}\]\_\{nm\}, then the local update at nodennis therefore

𝝀~n,t\[k\]\\displaystyle\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{n,t\}^\{\[k\]\}=𝝀~n,t\[k−1\]\+α​∑m=1Nln​m​𝝃m,t\[k−1\],\\displaystyle=\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{n,t\}^\{\[k\-1\]\}\+\\alpha\\sum\_\{m=1\}^\{N\}l\_\{nm\}\\boldsymbol\{\\xi\}\_\{m,t\}^\{\[k\-1\]\},\(15a\)𝝃n,t\[k\]\\displaystyle\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\}=𝝌n,t−𝝀~n,t\[k\]−τ​∑m=1Nln​m​𝝃m,t\[k−1\]\.\\displaystyle=\\boldsymbol\{\\chi\}\_\{n,t\}\-\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{n,t\}^\{\[k\]\}\-\\tau\\sum\_\{m=1\}^\{N\}l\_\{nm\}\\boldsymbol\{\\xi\}\_\{m,t\}^\{\[k\-1\]\}\.\(15b\)
Now, to derive the distributed ADMM\-based information matrix updates, we define the half\-vectorization of the information matrix,𝝎≜vech​\(𝛀\)\\boldsymbol\{\\omega\}\\triangleq\\text\{vech\}\(\\boldsymbol\{\\Omega\}\), which has dimensionDω≜12​Dξ​\(Dξ\+1\)D\_\{\\omega\}\\triangleq\\tfrac\{1\}\{2\}D\_\{\\xi\}\(D\_\{\\xi\}\+1\)\. The distributed optimization problem is then

###### Problem 3\.2\.

ωtminimize

1 2∑\_n=1^N∥ϕ\_n,t \-ω\_n,t∥^2, subject toL\_ωω\_t =0,

where𝝎t≜\[𝝎1,t⊤⋯𝝎N,t⊤\]⊤\\boldsymbol\{\\omega\}\_\{t\}\\triangleq\[\\boldsymbol\{\\omega\}\_\{1,t\}^\{\\top\}\\ \\cdots\\ \\boldsymbol\{\\omega\}\_\{N,t\}^\{\\top\}\]^\{\\top\},ϕn,t≜vech​\(N​𝐇n,t⊤​𝐑n,t−1​𝐇n,t\)\+𝝎n,t−1\\boldsymbol\{\\phi\}\_\{n,t\}\\triangleq\\text\{vech\}\(N\\mathbf\{H\}\_\{n,t\}^\{\\top\}\\mathbf\{R\}\_\{n,t\}^\{\-1\}\\mathbf\{H\}\_\{n,t\}\)\+\\boldsymbol\{\\omega\}\_\{n,t\-1\}, and𝐋ω≜𝐋⊗𝐈Dω\\mathbf\{L\}\_\{\\omega\}\\triangleq\\mathbf\{L\}\\otimes\\mathbf\{I\}\_\{D\_\{\\omega\}\}\. Following a similar procedure as for the information vector, the information matrix updates are

𝝂~t\[k\]\\displaystyle\\tilde\{\\boldsymbol\{\\nu\}\}\_\{t\}^\{\[k\]\}=𝝂~t\[k−1\]\+α​𝐋ω​𝝎t\[k−1\],\\displaystyle=\\tilde\{\\boldsymbol\{\\nu\}\}\_\{t\}^\{\[k\-1\]\}\+\\alpha\\mathbf\{L\}\_\{\\omega\}\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\-1\]\},\(16a\)𝝎t\[k\]\\displaystyle\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\]\}=ϕt−𝝂~t\[k\]−τ​𝐋ω​𝝎t\[k−1\],\\displaystyle=\\boldsymbol\{\\phi\}\_\{t\}\-\\tilde\{\\boldsymbol\{\\nu\}\}\_\{t\}^\{\[k\]\}\-\\tau\\mathbf\{L\}\_\{\\omega\}\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\-1\]\},\(16b\)where𝝂~t∈ℝN​Dω\\tilde\{\\boldsymbol\{\\nu\}\}\_\{t\}\\in\\mathbb\{R\}^\{ND\_\{\\omega\}\}is the global auxiliary dual variable\. The fully distributed updates at nodennare

𝝂~n,t\[k\]\\displaystyle\\tilde\{\\boldsymbol\{\\nu\}\}\_\{n,t\}^\{\[k\]\}=𝝂~n,t\[k−1\]\+α​∑m=1Nln​m​𝝎m,t\[k−1\],\\displaystyle=\\tilde\{\\boldsymbol\{\\nu\}\}\_\{n,t\}^\{\[k\-1\]\}\+\\alpha\\sum\_\{m=1\}^\{N\}l\_\{nm\}\\boldsymbol\{\\omega\}\_\{m,t\}^\{\[k\-1\]\},\(17a\)𝝎n,t\[k\]\\displaystyle\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\]\}=ϕn,t−𝝂~n,t\[k\]−τ​∑m=1Nln​m​𝝎m,t\[k−1\],\\displaystyle=\\boldsymbol\{\\phi\}\_\{n,t\}\-\\tilde\{\\boldsymbol\{\\nu\}\}\_\{n,t\}^\{\[k\]\}\-\\tau\\sum\_\{m=1\}^\{N\}l\_\{nm\}\\boldsymbol\{\\omega\}\_\{m,t\}^\{\[k\-1\]\},\(17b\)where𝝂~n,t∈ℝDω\\tilde\{\\boldsymbol\{\\nu\}\}\_\{n,t\}\\in\\mathbb\{R\}^\{D\_\{\\omega\}\}is the auxiliary dual variable at nodenn\. AfterKKiterations of ADMM, the estimates of the mean𝝁n,t\\boldsymbol\{\\mu\}\_\{n,t\}and covariance𝚺n,t\\boldsymbol\{\\Sigma\}\_\{n,t\}can be obtained at each agentnnby computing

𝝁n,t=𝚺n,t​𝝃n,t\[K\],𝚺n,t=𝛀n,t−1,\\displaystyle\\boldsymbol\{\\mu\}\_\{n,t\}=\\mathbf\{\\Sigma\}\_\{n,t\}\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[K\]\},\\quad\\mathbf\{\\Sigma\}\_\{n,t\}=\\mathbf\{\\Omega\}\_\{n,t\}^\{\-1\},\(18\)where𝛀n,t=vech−1​\(𝝎n,t\[K\]\)\\mathbf\{\\Omega\}\_\{n,t\}=\\text\{vech\}^\{\-1\}\(\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[K\]\}\)\. The ADMM\-RGP algorithm is summarized in Algorithm[1](https://arxiv.org/html/2609.26979#alg1)\.

Algorithm 1ADMM\-RGP1:Inputs:

𝐋\\mathbf\{L\},

𝓚\\boldsymbol\{\\mathcal\{K\}\},

α\\alpha,

τ\\tau
2:

𝝎n,0←vech​\(\(𝐊⁡\(𝐗p,𝐗p\)\+ϵ​𝐈\)−1\)\\boldsymbol\{\\omega\}\_\{n,0\}\\leftarrow\\text\{vech\}\(\(\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{p\}\)\+\\epsilon\\mathbf\{I\}\)^\{\-1\}\)∀n∈𝒱\\forall n\\in\\mathcal\{V\}
3:

𝝃n,0←𝟎\\boldsymbol\{\\xi\}\_\{n,0\}\\leftarrow\\mathbf\{0\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}
4:for

t=1t=1to

TTdo

5:for

n∈𝒱n\\in\\mathcal\{V\}do

6:Get new batch of measurements

\{𝐗n,t,𝐲n,t\}\\\{\\mathbf\{X\}\_\{n,t\},\\mathbf\{y\}\_\{n,t\}\\\}
7:Calculate

𝝌n,t\\boldsymbol\{\\chi\}\_\{n,t\}and

ϕn,t\\boldsymbol\{\\phi\}\_\{n,t\}
8:

𝝃n,t\[0\]←𝝌n,t,𝝎n,t\[0\]←ϕn,t\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[0\]\}\\leftarrow\\boldsymbol\{\\chi\}\_\{n,t\},\\ \\boldsymbol\{\\omega\}\_\{n,t\}^\{\[0\]\}\\leftarrow\\boldsymbol\{\\phi\}\_\{n,t\}
9:

𝝀~n,t\[0\]←𝟎,𝝂~n,t\[0\]←𝟎\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{n,t\}^\{\[0\]\}\\leftarrow\\mathbf\{0\},\\tilde\{\\boldsymbol\{\\nu\}\}\_\{n,t\}^\{\[0\]\}\\leftarrow\\mathbf\{0\}
10:endfor

11:for

k=1k=1to

KKdo

12:for

n∈𝒱n\\in\\mathcal\{V\}do

13:Share

𝝃n,t\[k−1\],𝝎n,t\[k−1\]\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\-1\]\},\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\-1\]\}with all

m∈𝒩nm\\in\\mathcal\{N\}\_\{n\}
14:Receive

𝝃m,t\[k−1\],𝝎m,t\[k−1\]\\boldsymbol\{\\xi\}\_\{m,t\}^\{\[k\-1\]\},\\boldsymbol\{\\omega\}\_\{m,t\}^\{\[k\-1\]\}from all

m∈𝒩nm\\in\\mathcal\{N\}\_\{n\}
15:Update

𝝀~n,t\[k\]\\tilde\{\\boldsymbol\{\\lambda\}\}\_\{n,t\}^\{\[k\]\}and

𝝃n,t\[k\]\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\}following \([15a](https://arxiv.org/html/2609.26979#S3.E15.1)\)–\([15b](https://arxiv.org/html/2609.26979#S3.E15.2)\)

16:Update

𝝂~n,t\[k\]\\tilde\{\\boldsymbol\{\\nu\}\}\_\{n,t\}^\{\[k\]\}and

𝝎n,t\[k\]\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\]\}following \([17a](https://arxiv.org/html/2609.26979#S3.E17.1)\)–\([17b](https://arxiv.org/html/2609.26979#S3.E17.2)\)

17:endfor

18:endfor

19:

𝝃n,t←𝝃n,t\[K\]\\boldsymbol\{\\xi\}\_\{n,t\}\\leftarrow\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[K\]\},

𝝎n,t←𝝎n,t\[K\]\\boldsymbol\{\\omega\}\_\{n,t\}\\leftarrow\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[K\]\}for

∀n∈𝒱\\forall n\\in\\mathcal\{V\}
20:Calculate

𝝁n,t\\boldsymbol\{\\mu\}\_\{n,t\}and

𝚺n,t\\mathbf\{\\Sigma\}\_\{n,t\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}using \([18](https://arxiv.org/html/2609.26979#S3.E18)\)

21:endfor

22:Outputs:

𝝁n,T\\boldsymbol\{\\mu\}\_\{n,T\}and

𝚺n,T\\mathbf\{\\Sigma\}\_\{n,T\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}

### 3\.2PDMM\-RGP

The optimization problems in the ADMM\-RGP derivation are node\-separable with strongly convex, differentiable objectives and convex constraints, making them well suited for the primal–dual method of multipliers \(PDMM\) framework\[[25](https://arxiv.org/html/2609.26979#bib.bib25)\]\. PDMM has been empirically shown to converge faster than ADMM in certain settings, including average consensus problems\[[22](https://arxiv.org/html/2609.26979#bib.bib22)\]\.

To reformulate[3\.1](https://arxiv.org/html/2609.26979#Thmtheorem1)for the application of PDMM, we define two directed edges\(n→m\)\(n\\rightarrow m\)and\(m→n\)\(m\\rightarrow n\)for each\(n,m\)∈ℰ\(n,m\)\\in\\mathcal\{E\}\. Furthermore let𝐀n→m≜an→m​𝐈Dξ\\mathbf\{A\}\_\{n\\rightarrow m\}\\triangleq a\_\{n\\rightarrow m\}\\mathbf\{I\}\_\{D\_\{\\xi\}\}, where

an→m≜\{−\[𝐋\]n​mif​n<m,\[𝐋\]n​mif​n\>m,\\displaystyle a\_\{n\\rightarrow m\}\\triangleq\\begin\{cases\}\-\[\\mathbf\{L\}\]\_\{nm\}&\\text\{if \}n<m,\\\\ \[\\mathbf\{L\}\]\_\{nm\}&\\text\{if \}n\>m,\\end\{cases\}\(19\)which leads to the following reformulated problem

###### Problem 3\.3\.

ξtminimize

1 2∑\_n=1^N∥χ\_n,t \-ξ\_n,t∥^2, subject toA\_n→mξ\_n,t \+A\_m→nξ\_m,t=0∀\(n,m\)∈E\.

Now, we construct a lifted constraint matrix, by assigning each directed edge\(n→m\)\(n\\rightarrow m\)to a unique indexin→m∈\{1,…,2​M\}i\_\{n\\rightarrow m\}\\in\\\{1,\\dots,2M\\\}, whereM=\|ℰ\|M=\|\\mathcal\{E\}\|denotes the number of undirected edges in the communication graph\. Let the lifted constraint matrix be denoted by𝐂∈ℝ2​M×N\\mathbf\{C\}\\in\\mathbb\{R\}^\{2M\\times N\}then

\[𝐂\]r​n=\{an→mif​r=in→m,0otherwise\.\[\\mathbf\{C\}\]\_\{rn\}=\\begin\{cases\}a\_\{n\\rightarrow m\}&\\text\{if \}r=i\_\{n\\rightarrow m\},\\\\ 0&\\text\{otherwise\}\.\\end\{cases\}\(20\)Additionally, consider a2​M×2​M2M\\times 2Mpermutation matrix𝐏\\mathbf\{P\}that exchanges the components at indicesin→mi\_\{n\\rightarrow m\}andim→ni\_\{m\\rightarrow n\}for each undirected edge\(n,m\)∈ℰ\(n,m\)\\in\\mathcal\{E\}\. Lastly, let the auxiliary dual variables be denoted by𝝊n\|m,t∈ℝDξ\\boldsymbol\{\\upsilon\}\_\{n\|m,t\}\\in\\mathbb\{R\}^\{D\_\{\\xi\}\}at each directed edge\(n→m\)\(n\\rightarrow m\)of the communication graph, which are aggregated into𝝊t∈ℝ2​M​Dξ\\boldsymbol\{\\upsilon\}\_\{t\}\\in\\mathbb\{R\}^\{2MD\_\{\\xi\}\}following the same edge order used to build𝐂\\mathbf\{C\}and𝐏\\mathbf\{P\}\.

The PDMM updates\[[26](https://arxiv.org/html/2609.26979#bib.bib26)\]to solve[3\.3](https://arxiv.org/html/2609.26979#Thmtheorem3)are then given by

𝝃t\[k\]\\displaystyle\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}=arg⁡min𝝃​\(‖𝝌t−𝝃‖2\+\(𝐏ξ​𝝊t\[k−1\]\)⊤​𝐂ξ​𝝃\+c2​‖𝐂ξ​𝝃‖2\),\\displaystyle=\\arg\\underset\{\\boldsymbol\{\\xi\}\}\{\\min\}\(\\\|\\boldsymbol\{\\chi\}\_\{t\}\-\\boldsymbol\{\\xi\}\\\|^\{2\}\+\(\\mathbf\{P\}\_\{\\xi\}\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\-1\]\}\)^\{\\top\}\\mathbf\{C\}\_\{\\xi\}\\boldsymbol\{\\xi\}\+\\frac\{c\}\{2\}\\\|\\mathbf\{C\}\_\{\\xi\}\\boldsymbol\{\\xi\}\\\|^\{2\}\),𝝊t\[k\]\\displaystyle\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\]\}=𝐏ξ​𝝊t\[k−1\]\+2​c​𝐂ξ​𝝃t\[k\],\\displaystyle=\\mathbf\{P\}\_\{\\xi\}\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\-1\]\}\+2c\\mathbf\{C\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\},wherec\>0c\>0is a constant parameter,𝐂ξ≜𝐂⊗𝐈Dξ\\mathbf\{C\}\_\{\\xi\}\\triangleq\\mathbf\{C\}\\otimes\\mathbf\{I\}\_\{D\_\{\\xi\}\}, and𝐏ξ≜𝐏⊗𝐈Dξ\\mathbf\{P\}\_\{\\xi\}\\triangleq\\mathbf\{P\}\\otimes\\mathbf\{I\}\_\{D\_\{\\xi\}\}\. To obtain the𝝃t\\boldsymbol\{\\xi\}\_\{t\}update, take the gradient with respect to the primal variable and set it to zero\. Because𝐂⊤​𝐂\\mathbf\{C\}^\{\\top\}\\mathbf\{C\}is a diagonal matrix with entries that sum the squared weights corresponding to each node andc\>0c\>0,\(𝐈\+c​𝐂ξ⊤​𝐂ξ\)\(\\mathbf\{I\}\+c\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\\mathbf\{C\}\_\{\\xi\}\)is invertible, leading us to the global PDMM iterates for the information vector

𝝃t\[k\]\\displaystyle\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}=\(𝐈\+c​𝐂ξ⊤​𝐂ξ\)−1​\(𝝌t−𝐂ξ⊤​𝐏ξ​𝝊t\[k−1\]\),\\displaystyle=\(\\mathbf\{I\}\+c\\mathbf\{C\_\{\\xi\}^\{\\top\}C\_\{\\xi\}\}\)^\{\-1\}\(\\boldsymbol\{\\chi\}\_\{t\}\-\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\\mathbf\{P\}\_\{\\xi\}\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\-1\]\}\),\(22a\)𝝊t\[k\]\\displaystyle\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\]\}=𝐏ξ​𝝊t\[k−1\]\+2​c​𝐂ξ​𝝃t\[k\]\.\\displaystyle=\\mathbf\{P\}\_\{\\xi\}\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\-1\]\}\+2c\\mathbf\{C\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}\.\(22b\)Subsequently, the distributed iterates at each agentnnare

𝝃n,t\[k\]\\displaystyle\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\}=𝝌n,t−∑m∈𝒩nan→m​𝝊m\|n,t\[k−1\]1\+c​∑m∈𝒩nan→m2,\\displaystyle=\\frac\{\\boldsymbol\{\\chi\}\_\{n,t\}\-\\sum\_\{m\\in\\mathcal\{N\}\_\{n\}\}a\_\{n\\rightarrow m\}\\boldsymbol\{\\upsilon\}\_\{m\|n,t\}^\{\[k\-1\]\}\}\{1\+c\\sum\_\{m\\in\\mathcal\{N\}\_\{n\}\}a\_\{n\\rightarrow m\}^\{2\}\},\(23a\)𝝊n\|m,t\[k\]\\displaystyle\\boldsymbol\{\\upsilon\}\_\{n\|m,t\}^\{\[k\]\}=𝝊m\|n,t\[k−1\]\+2​c​an→m​𝝃n,t\[k\]\.\\displaystyle=\\boldsymbol\{\\upsilon\}\_\{m\|n,t\}^\{\[k\-1\]\}\+2ca\_\{n\\rightarrow m\}\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\}\.\(23b\)
Note that the agentnnrequires the value of𝝊m\|n,t\\boldsymbol\{\\upsilon\}\_\{m\|n,t\}from agentmm\. Rather than transmitting these values individually over each edge using unicast communication, agentmmcan broadcast its primal variable𝝃m,t\\boldsymbol\{\\xi\}\_\{m,t\}to all of its neighbors, allowing each neighbor to locally reconstruct the value of𝝊m\|n,t\\boldsymbol\{\\upsilon\}\_\{m\|n,t\}\. This broadcast strategy reduces the communication overhead\.Along similar lines,[3\.2](https://arxiv.org/html/2609.26979#Thmtheorem2)can be solved with PDMM by defining auxiliary dual variables𝝍n\|m,t∈ℝDω\\boldsymbol\{\\psi\}\_\{n\|m,t\}\\in\\mathbb\{R\}^\{D\_\{\\omega\}\}for each directed edge\(n→m\)\(n\\rightarrow m\)in the communication graph\. The local updates at nodennare

𝝎n,t\[k\]\\displaystyle\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\]\}=ϕn,t−∑m∈𝒩nan→m​𝝍m\|n,t\[k−1\]1\+c​∑m∈𝒩nan→m2,\\displaystyle=\\frac\{\\boldsymbol\{\\phi\}\_\{n,t\}\-\\sum\_\{m\\in\\mathcal\{N\}\_\{n\}\}a\_\{n\\rightarrow m\}\\boldsymbol\{\\psi\}\_\{m\|n,t\}^\{\[k\-1\]\}\}\{1\+c\\sum\_\{m\\in\\mathcal\{N\}\_\{n\}\}a\_\{n\\rightarrow m\}^\{2\}\},\(24a\)𝝍n\|m,t\[k\]\\displaystyle\\boldsymbol\{\\psi\}\_\{n\|m,t\}^\{\[k\]\}=𝝍m\|n,t\[k−1\]\+2​c​an→m​𝝎n,t\[k\]\.\\displaystyle=\\boldsymbol\{\\psi\}\_\{m\|n,t\}^\{\[k\-1\]\}\+2ca\_\{n\\rightarrow m\}\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\]\}\.\(24b\)The PDMM\-RGP algorithm using a broadcast communication protocol is summarized in Algorithm[2](https://arxiv.org/html/2609.26979#alg2)\.

Algorithm 2PDMM\-RGP1:Inputs:

𝐋\\mathbf\{L\},

𝓚\\boldsymbol\{\\mathcal\{K\}\},

cc
2:

𝝎n,0←vech​\(\(𝐊⁡\(𝐗p,𝐗p\)\+ϵ​𝐈\)−1\)\\boldsymbol\{\\omega\}\_\{n,0\}\\leftarrow\\text\{vech\}\(\(\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{p\}\)\+\\epsilon\\mathbf\{I\}\)^\{\-1\}\),

∀n∈𝒱\\forall n\\in\\mathcal\{V\}
3:

𝝃n,0←𝟎\\boldsymbol\{\\xi\}\_\{n,0\}\\leftarrow\\mathbf\{0\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}
4:for

t=1t=1to

TTdo

5:for

n∈𝒱n\\in\\mathcal\{V\}do

6:Get new batch of measurements

\(𝐗n,t,𝐲n,t\)\(\\mathbf\{X\}\_\{n,t\},\\mathbf\{y\}\_\{n,t\}\)
7:Calculate

𝝌n,t\\boldsymbol\{\\chi\}\_\{n,t\}and

ϕn,t\\boldsymbol\{\\phi\}\_\{n,t\}
8:

𝝃n,t\[0\]←𝝌n,t,𝝎n,t\[0\]←ϕn,t\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[0\]\}\\leftarrow\\boldsymbol\{\\chi\}\_\{n,t\},\\ \\boldsymbol\{\\omega\}\_\{n,t\}^\{\[0\]\}\\leftarrow\\boldsymbol\{\\phi\}\_\{n,t\}
9:

𝝊n\|m,t\[0\]←𝟎,𝝍n\|m,t\[0\]←𝟎\\boldsymbol\{\\upsilon\}\_\{n\|m,t\}^\{\[0\]\}\\leftarrow\\mathbf\{0\},\\boldsymbol\{\\psi\}\_\{n\|m,t\}^\{\[0\]\}\\leftarrow\\mathbf\{0\}for

m∈𝒩nm\\in\\mathcal\{N\}\_\{n\}\(at

nn\)

10:

𝝊m\|n,t\[0\]←𝟎,𝝍m\|n,t\[0\]←𝟎\\boldsymbol\{\\upsilon\}\_\{m\|n,t\}^\{\[0\]\}\\leftarrow\\mathbf\{0\},\\boldsymbol\{\\psi\}\_\{m\|n,t\}^\{\[0\]\}\\leftarrow\\mathbf\{0\}for

m∈𝒩nm\\in\\mathcal\{N\}\_\{n\}\(at

nn\)

11:endfor

12:for

k=1k=1to

KKdo

13:for

n∈𝒱n\\in\\mathcal\{V\}do

14:Update

𝝃n,t\[k\]\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\}and

𝝊n\|m,t\[k\]\\boldsymbol\{\\upsilon\}\_\{n\|m,t\}^\{\[k\]\}using \([23a](https://arxiv.org/html/2609.26979#S3.E23.1)\)\-\([23b](https://arxiv.org/html/2609.26979#S3.E23.2)\)

15:Update

𝝎n,t\[k\]\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\]\}and

𝝍n\|m,t\[k\]\\boldsymbol\{\\psi\}\_\{n\|m,t\}^\{\[k\]\}using \([24a](https://arxiv.org/html/2609.26979#S3.E24.1)\)\-\([24b](https://arxiv.org/html/2609.26979#S3.E24.2)\)

16:Share

𝝃n,t\[k\],𝝎n,t\[k\]\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[k\]\},\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[k\]\}with all

m∈𝒩nm\\in\\mathcal\{N\}\_\{n\}
17:Receive

𝝃m,t\[k\],𝝎m,t\[k\]\\boldsymbol\{\\xi\}\_\{m,t\}^\{\[k\]\},\\boldsymbol\{\\omega\}\_\{m,t\}^\{\[k\]\}from all

m∈𝒩nm\\in\\mathcal\{N\}\_\{n\}
18:

𝝊m\|n,t\[k\]←𝝊n\|m,t\[k−1\]\+2​c​am→n​𝝃m,t\[k\]\\boldsymbol\{\\upsilon\}^\{\[k\]\}\_\{m\|n,t\}\\leftarrow\\boldsymbol\{\\upsilon\}\_\{n\|m,t\}^\{\[k\-1\]\}\+2ca\_\{m\\rightarrow n\}\\boldsymbol\{\\xi\}\_\{m,t\}^\{\[k\]\}\(at

nn\)

19:

𝝍m\|n,t\[k\]←𝝍n\|m,t\[k−1\]\+2​c​am→n​𝝎m,t\[k\]\\boldsymbol\{\\psi\}^\{\[k\]\}\_\{m\|n,t\}\\leftarrow\\boldsymbol\{\\psi\}\_\{n\|m,t\}^\{\[k\-1\]\}\+2ca\_\{m\\rightarrow n\}\\boldsymbol\{\\omega\}\_\{m,t\}^\{\[k\]\}\(at

nn\)

20:endfor

21:endfor

22:

𝝃n,t←𝝃n,t\[K\]\\boldsymbol\{\\xi\}\_\{n,t\}\\leftarrow\\boldsymbol\{\\xi\}\_\{n,t\}^\{\[K\]\},

𝝎n,t←𝝎n,t\[K\]\\boldsymbol\{\\omega\}\_\{n,t\}\\leftarrow\\boldsymbol\{\\omega\}\_\{n,t\}^\{\[K\]\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}
23:Calculate

𝝁n,t\\boldsymbol\{\\mu\}\_\{n,t\}and

𝚺n,t\\mathbf\{\\Sigma\}\_\{n,t\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}using \([18](https://arxiv.org/html/2609.26979#S3.E18)\)

24:endfor

25:Outputs:

𝝁n,T\\boldsymbol\{\\mu\}\_\{n,T\}and

𝚺n,T\\mathbf\{\\Sigma\}\_\{n,T\}∀n∈𝒱\\forall n\\in\\mathcal\{V\}

## 4Analysis

In this section, we establish the convergence properties of the proposed ADMM\-RGP and PDMM\-RGP algorithms\.

### 4\.1ADMM\-RGP Convergence

We begin with the analysis in\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\]\. Recall the ADMM\-RGP iterates from \([14b](https://arxiv.org/html/2609.26979#S3.E14.2)\) and \([16b](https://arxiv.org/html/2609.26979#S3.E16.2)\), which can be written as a discrete\-time, second\-order system of the form

𝝃t\[k\]\\displaystyle\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}=\(𝐈−\(α\+τ\)​𝐋ξ\)​𝝃t\[k−1\]\+τ​𝐋ξ​𝝃t\[k−2\],\\displaystyle=\(\\mathbf\{I\}\-\(\\alpha\+\\tau\)\\mathbf\{L\}\_\{\\xi\}\)\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\-1\]\}\+\\tau\\mathbf\{L\}\_\{\\xi\}\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\-2\]\},\(25a\)𝝎t\[k\]\\displaystyle\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\]\}=\(𝐈−\(α\+τ\)​𝐋ω\)​𝝎t\[k−1\]\+τ​𝐋ω​𝝎t\[k−2\]\.\\displaystyle=\(\\mathbf\{I\}\-\(\\alpha\+\\tau\)\\mathbf\{L\}\_\{\\omega\}\)\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\-1\]\}\+\\tau\\mathbf\{L\}\_\{\\omega\}\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\-2\]\}\.\(25b\)We further define error terms

𝐞t,ξ\[k\]≜𝝃t∗−𝝃t\[k\],𝐞t,ω\[k\]≜𝝎t∗−𝝎t\[k\],\\displaystyle\\mathbf\{e\}\_\{t,\\xi\}^\{\[k\]\}\\triangleq\\boldsymbol\{\\xi\}\_\{t\}^\{\*\}\-\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\},\\quad\\mathbf\{e\}\_\{t,\\omega\}^\{\[k\]\}\\triangleq\\boldsymbol\{\\omega\}\_\{t\}^\{\*\}\-\\boldsymbol\{\\omega\}\_\{t\}^\{\[k\]\},\(26\)𝐞t\[k\]≜\[𝐞t,ξ\[k\]⊤𝐞t,ω\[k\]⊤\]⊤\.\\displaystyle\\mathbf\{e\}\_\{t\}^\{\[k\]\}\\triangleq\\begin\{bmatrix\}\\mathbf\{e\}\_\{t,\\xi\}^\{\[k\]\\top\}&\\mathbf\{e\}\_\{t,\\omega\}^\{\[k\]\\top\}\\end\{bmatrix\}^\{\\top\}\.\(27\)Since the error𝐞t\[k\]\\mathbf\{e\}\_\{t\}^\{\[k\]\}vanishes in the consensus subspace associated withλ1​\(𝐋\)=0\\lambda\_\{1\}\(\\mathbf\{L\}\)=0, we remove this component and diagonalize the remaining dynamics to obtain the transformed error𝐞~t\[k\]∈ℝ\(N−1\)​\(Dξ\+Dω\)\\tilde\{\\mathbf\{e\}\}\_\{t\}^\{\[k\]\}\\in\\mathbb\{R\}^\{\(N\-1\)\(D\_\{\\xi\}\+D\_\{\\omega\}\)\}, with dynamics

\[𝐞~t\[k\]𝐞~t\[k−1\]\]=\[\(𝐈−\(α\+τ\)​𝚲¯\)τ​𝚲¯𝐈𝟎\]⏟𝐌​\[𝐞~t\[k−1\]𝐞~t\[k−2\]\],\\begin\{bmatrix\}\\tilde\{\\mathbf\{e\}\}\_\{t\}^\{\[k\]\}\\\\ \\tilde\{\\mathbf\{e\}\}\_\{t\}^\{\[k\-1\]\}\\end\{bmatrix\}=\\underbrace\{\\begin\{bmatrix\}\(\\mathbf\{I\}\-\(\\alpha\+\\tau\)\\bar\{\\mathbf\{\\Lambda\}\}\)&\\tau\\bar\{\\mathbf\{\\Lambda\}\}\\\\ \\mathbf\{I\}&\\mathbf\{0\}\\end\{bmatrix\}\}\_\{\\mathbf\{M\}\}\\begin\{bmatrix\}\\tilde\{\\mathbf\{e\}\}\_\{t\}^\{\[k\-1\]\}\\\\ \\tilde\{\\mathbf\{e\}\}\_\{t\}^\{\[k\-2\]\}\\end\{bmatrix\},\(28\)where𝚲¯≜diag​\(λ2​\(𝐋\),…,λN​\(𝐋\)\)⊗𝐈Dξ\+Dω\\bar\{\\boldsymbol\{\\Lambda\}\}\\triangleq\\text\{diag\}\(\\lambda\_\{2\}\(\\mathbf\{L\}\),\\dots,\\lambda\_\{N\}\(\\mathbf\{L\}\)\)\\otimes\\mathbf\{I\}\_\{D\_\{\\xi\}\+D\_\{\\omega\}\}\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\]\. Then, the eigenvalues of𝐌\\mathbf\{M\}are given by the eigenvalues of𝐌i\\mathbf\{M\}\_\{i\}fori=2,⋯,Ni=2,\\cdots,N, where

𝐌i\\displaystyle\\mathbf\{M\}\_\{i\}≜\[1−\(α\+τ\)​λi​\(𝐋\)τ​λi​\(𝐋\)10\],\\displaystyle\\triangleq\\begin\{bmatrix\}1\-\(\\alpha\+\\tau\)\\mathbf\{\\lambda\}\_\{i\}\(\\mathbf\{L\}\)&\\tau\\mathbf\{\\lambda\}\_\{i\}\(\\mathbf\{L\}\)\\\\ 1&0\\end\{bmatrix\},\(29\)λ⁡\(𝐌i\)\\displaystyle\\lambda\(\\mathbf\{M\}\_\{i\}\)=12​\(α~i±α~i2\+4​τ​λi​\(𝐋\)\),\\displaystyle=\\frac\{1\}\{2\}\\Big\(\\tilde\{\\alpha\}\_\{i\}\\pm\\sqrt\{\\tilde\{\\alpha\}\_\{i\}^\{2\}\+4\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)\}\\ \\Big\),\(30\)whereα~i≜1−\(α\+τ\)​λi​\(𝐋\)\\tilde\{\\alpha\}\_\{i\}\\triangleq 1\-\(\\alpha\+\\tau\)\\lambda\_\{i\}\(\\mathbf\{L\}\)\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\]\. Since the transformed error dynamics constitute a discrete\-time linear system, convergence of the iterates is guaranteed if𝕄\\mathbb\{M\}is Schur stable i\.e\., if all eigenvalues of𝐌\\mathbf\{M\}lie strictly inside the unit circle\. In this case,𝐌\\mathbf\{M\}is shown to be Schur stable provided thatα\+2​τ<2​\(λN​\(𝐋\)\)−1\\alpha\+2\\tau<2\(\\lambda\_\{\\text\{N\}\}\(\\mathbf\{L\}\)\)^\{\-1\}andα,τ\>0\\alpha,\\tau\>0\[[21](https://arxiv.org/html/2609.26979#bib.bib21)\]\.

Note that the values ofα\\alphaandτ\\taumay be selected to minimize the spectral radius,ρ⁡\(𝐌\)\\rho\(\\mathbf\{M\}\), for fast convergence\. Assuming positiveα,τ\>0\\alpha,\\tau\>0, the eigenvalues of𝐌\\mathbf\{M\}are real, and the spectral radius can be expressed as

ρ⁡\(𝐌\)=maxi,j⁡\|λj​\(𝐌i\)\|=maxi⁡12​\(\|α~i\|\+α~i2\+4​τ​λi​\(𝐋\)\)\.\\rho\(\\mathbf\{M\}\)=\\max\_\{i,j\}\|\\lambda\_\{j\}\(\\mathbf\{M\}\_\{i\}\)\|=\\max\_\{i\}\\ \\frac\{1\}\{2\}\(\|\\tilde\{\\alpha\}\_\{i\}\|\+\\sqrt\{\\tilde\{\\alpha\}\_\{i\}^\{2\}\+4\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)\}\)\.The4​τ​λi​\(𝐋\)4\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)term increases the maximum eigenvalue, meaning thatτ\\taushould be chosen as close to00as possible to minimize the spectral radius\. However, whenτ=0\\tau=0, the ADMM iterates are identical to first\-order average consensus\. Therefore, we observe that forτ\>0\\tau\>0, ADMM\-RGP is guaranteed to be slower than Consensus\-RGP\. However, forτ<0\\tau<0, faster convergence could be achieved while preserving complexity\. The augmented Lagrangian in \([12](https://arxiv.org/html/2609.26979#S3.E12)\) may be written as

ℒτ​\(𝝃,𝝀\)\\displaystyle\\mathcal\{L\}\_\{\\tau\}\(\\boldsymbol\{\\xi\},\\boldsymbol\{\\lambda\}\)=12​𝝃⊤​\(𝐈\+τ​𝐋ξ\)​𝝃−\(𝝌⊤\+𝝀⊤​𝐁ξ\)​𝝃\+12​𝝌⊤​𝝌,\\displaystyle=\\frac\{1\}\{2\}\\boldsymbol\{\\xi\}^\{\\top\}\(\\mathbf\{I\}\+\\tau\\mathbf\{L\}\_\{\\xi\}\)\\boldsymbol\{\\xi\}\-\(\\boldsymbol\{\\chi\}^\{\\top\}\+\\boldsymbol\{\\lambda\}^\{\\top\}\\mathbf\{B\}\_\{\\xi\}\)\\boldsymbol\{\\xi\}\+\\frac\{1\}\{2\}\\boldsymbol\{\\chi\}^\{\\top\}\\boldsymbol\{\\chi\},which is convex in the primal variable if\(𝐈\+τ​𝐋ξ\)⪰𝟎\(\\mathbf\{I\}\+\\tau\\mathbf\{L\}\_\{\\xi\}\)\\succeq\\mathbf\{0\}\. Forτ<0\\tau<0, the smallest eigenvalue of\(𝐈\+τ​𝐋ξ\)\(\\mathbf\{I\}\+\\tau\\mathbf\{L\}\_\{\\xi\}\)is1\+τ​λN​\(𝐋\)1\+\\tau\\lambda\_\{N\}\(\\mathbf\{L\}\)\. Thus, the convexity of the augmented Lagrangian is preserved forτ\>−\(λN​\(𝐋\)\)−1\\tau\>\-\(\\lambda\_\{N\}\(\\mathbf\{L\}\)\)^\{\-1\}, leading to Theorem[4\.4](https://arxiv.org/html/2609.26979#Thmtheorem4)\.

###### Theorem 4\.4\.

Consider the error dynamics matrix𝐌\\mathbf\{M\}defined in \([28](https://arxiv.org/html/2609.26979#S4.E28)\) and its eigenvalues given by \([30](https://arxiv.org/html/2609.26979#S4.E30)\)\. Ifα\\alphaandτ\\tausatisfy

α\>0,−1λN​\(𝐋\)<τ<0,α\+2​τ<2λN​\(𝐋\),\\alpha\>0,\\quad\\ \\frac\{\-1\}\{\\lambda\_\{N\}\(\\mathbf\{L\}\)\}<\\tau<0,\\quad\\alpha\+2\\tau<\\frac\{2\}\{\\lambda\_\{N\}\(\\mathbf\{L\}\)\},\(31\)then𝐌\\mathbf\{M\}is Schur stable and ADMM\-RGP converges\.

###### Proof 4\.5\.

See Appendix[7](https://arxiv.org/html/2609.26979#S7)\.

The parametersα\\alphaandτ\\taumay be selected using a grid search over the range of values given in Theorem[4\.4](https://arxiv.org/html/2609.26979#Thmtheorem4)to find the values that jointly minimizeρ⁡\(𝐌\)\\rho\(\\mathbf\{M\}\)\. To reduce communication costs,KKmay be small enough that transient behavior dominates the error dynamics, in which case the parameters may be selected to minimize‖𝐌K‖F\\\|\\mathbf\{M\}^\{K\}\\\|\_\{F\}\.

### 4\.2PDMM\-RGP Convergence

We begin with the observation that the objective functions of[3\.1](https://arxiv.org/html/2609.26979#Thmtheorem1)and[3\.2](https://arxiv.org/html/2609.26979#Thmtheorem2)are closed, strongly convex, proper, and differentiable, so the PDMM iterations are guaranteed to converge for anyc\>0c\>0\[[25](https://arxiv.org/html/2609.26979#bib.bib25)\]\. Whileccdoes not affect asymptotic stability, note that it does influence convergence speed, a property of particular interest in communication\-constrained systems where the number ofkk\-iterations per time step may be limited\. To analyze the convergence of PDMM\-RGP, we define an auxiliary dual variable,𝜻t\[k\]≜𝐏ξ​𝝊t\[k\]\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\]\}\\triangleq\\mathbf\{P\}\_\{\\xi\}\\boldsymbol\{\\upsilon\}\_\{t\}^\{\[k\]\}, which has dynamics

𝜻t\[k\]\\displaystyle\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\]\}=𝐀ζ​𝜻t\[k−1\]\+2​c​𝐏ξ​𝐂ξ​𝐂¯c​𝝌t,\\displaystyle=\\mathbf\{A\}\_\{\\zeta\}\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\-1\]\}\+2c\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\\bar\{\\mathbf\{C\}\}\_\{c\}\\boldsymbol\{\\chi\}\_\{t\},\(32a\)𝐂¯c\\displaystyle\\bar\{\\mathbf\{C\}\}\_\{c\}≜\(𝐈\+c​𝐂ξ⊤​𝐂ξ\)−1,\\displaystyle\\triangleq\(\\mathbf\{I\}\+c\\mathbf\{C\_\{\\xi\}^\{\\top\}\\mathbf\{C\}\_\{\\xi\}\}\)^\{\-1\},\(32b\)𝐀ζ\\displaystyle\\mathbf\{A\}\_\{\\zeta\}≜𝐏ξ−2​c​𝐏ξ​𝐂ξ​𝐂¯c​𝐂ξ⊤\.\\displaystyle\\triangleq\\mathbf\{P\}\_\{\\xi\}\-2c\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\\bar\{\\mathbf\{C\}\}\_\{c\}\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\.\(32c\)Furthermore, we define subspacesΨ≜ran​\(𝐂ξ\)\+ran​\(𝐏ξ​𝐂ξ\)⊆ℝ2​M​Dξ\\Psi\\triangleq\\text\{ran\}\(\\mathbf\{C\}\_\{\\xi\}\)\+\\text\{ran\}\(\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\)\\subseteq\\mathbb\{R\}^\{2MD\_\{\\xi\}\}andΨ⟂≜ker​\(𝐂⊤\)∩ker​\(\(𝐏𝐂\)⊤\)⊆ℝ2​M​Dξ\\Psi^\{\\perp\}\\triangleq\\text\{ker\}\(\\mathbf\{C\}^\{\\top\}\)\\cap\\text\{ker\}\(\(\\mathbf\{PC\}\)^\{\\top\}\)\\subseteq\\mathbb\{R\}^\{2MD\_\{\\xi\}\}\. Let𝜻Ψ,t\[k\]\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\]\}and𝜻Ψ⟂,t\[k\]\\boldsymbol\{\\zeta\}\_\{\\Psi^\{\\perp\},t\}^\{\[k\]\}denote theΨ\\PsiandΨ⟂\\Psi^\{\\perp\}components of𝜻t\[k\]\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\]\}, respectively\. According to Lemma 5\.1 in\[[27](https://arxiv.org/html/2609.26979#bib.bib27)\],Ψ\\PsiandΨ⟂\\Psi^\{\\perp\}are invariant over𝐏ξ\\mathbf\{P\}\_\{\\xi\}\. In other words,𝐱∈Ψ⟹𝐏ξ​𝐱∈Ψ,𝐲∈Ψ⟂⟹𝐏ξ​𝐲∈Ψ⟂\\mathbf\{x\}\\in\\Psi\\implies\\mathbf\{P\}\_\{\\xi\}\\mathbf\{x\}\\in\\Psi,\\ \\mathbf\{y\}\\in\\Psi^\{\\perp\}\\implies\\mathbf\{P\}\_\{\\xi\}\\mathbf\{y\}\\in\\Psi^\{\\perp\}\. Furthermore,Ψ\\PsiandΨ⟂\\Psi^\{\\perp\}are invariant under𝐀ζ\\mathbf\{A\}\_\{\\zeta\}\.

𝐀ζ​𝜻Ψ,t\[k\]\\displaystyle\\mathbf\{A\}\_\{\\zeta\}\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\]\}=\(𝐏ξ−2​c​𝐏ξ​𝐂ξ​𝐂¯c​𝐂ξ⊤\)​𝜻Ψ,t\[k\]∈Ψ,\\displaystyle=\(\\mathbf\{P\}\_\{\\xi\}\-2c\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\\bar\{\\mathbf\{C\}\}\_\{c\}\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\)\\boldsymbol\{\\zeta\}^\{\[k\]\}\_\{\\Psi,t\}\\in\\Psi,\(33a\)𝐀ζ​𝜻Ψ⟂,t\[k\]\\displaystyle\\mathbf\{A\}\_\{\\zeta\}\\boldsymbol\{\\zeta\}\_\{\\Psi^\{\\perp\},t\}^\{\[k\]\}=\(𝐏ξ−2​c​𝐏ξ​𝐂ξ​𝐂¯c​𝐂ξ⊤\)​𝜻Ψ⟂,t\[k\]\\displaystyle=\(\\mathbf\{P\}\_\{\\xi\}\-2c\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\\bar\{\\mathbf\{C\}\}\_\{c\}\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\)\\boldsymbol\{\\zeta\}^\{\[k\]\}\_\{\\Psi^\{\\perp\},t\}\(33b\)=𝐏ξ​𝜻Ψ⟂,t\[k\]∈Ψ⟂\.\\displaystyle=\\mathbf\{P\}\_\{\\xi\}\\boldsymbol\{\\zeta\}^\{\[k\]\}\_\{\\Psi^\{\\perp\},t\}\\in\\Psi^\{\\perp\}\.\(33c\)Thus, the𝜻t\[k\]\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\]\}iterations decompose into,

𝜻Ψ,t\[k\]\\displaystyle\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\]\}=𝐀Ψ​𝜻Ψ,t\[k−1\]\+𝚷Ψ​\(2​c​𝐏ξ​𝐂ξ​𝐂¯c​𝝌t\),\\displaystyle=\\mathbf\{A\}\_\{\\Psi\}\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\-1\]\}\+\\mathbf\{\\Pi\}\_\{\\Psi\}\(2c\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\\bar\{\\mathbf\{C\}\}\_\{c\}\\boldsymbol\{\\chi\}\_\{t\}\),\(34a\)𝜻Ψ⟂,t\[k\]\\displaystyle\\boldsymbol\{\\zeta\}\_\{\\Psi^\{\\perp\},t\}^\{\[k\]\}=𝐀Ψ⟂​𝜻Ψ⟂,t\[k−1\]\+𝚷Ψ⟂​\(2​c​𝐏ξ​𝐂ξ​𝐂¯c​𝝌t\),\\displaystyle=\\mathbf\{A\}\_\{\\Psi^\{\\perp\}\}\\boldsymbol\{\\zeta\}\_\{\\Psi^\{\\perp\},t\}^\{\[k\-1\]\}\+\\mathbf\{\\Pi\}\_\{\\Psi^\{\\perp\}\}\(2c\\mathbf\{P\}\_\{\\xi\}\\mathbf\{C\}\_\{\\xi\}\\bar\{\\mathbf\{C\}\}\_\{c\}\\boldsymbol\{\\chi\}\_\{t\}\),\(34b\)where𝐀Ψ≜𝚷Ψ​𝐀ζ​𝚷Ψ\\mathbf\{A\}\_\{\\Psi\}\\triangleq\\mathbf\{\\Pi\}\_\{\\Psi\}\\mathbf\{A\}\_\{\\zeta\}\\mathbf\{\\Pi\}\_\{\\Psi\}and𝐀Ψ⟂≜𝚷Ψ⟂​𝐀ζ​𝚷Ψ⟂\\mathbf\{A\}\_\{\\Psi^\{\\perp\}\}\\triangleq\\mathbf\{\\Pi\}\_\{\\Psi^\{\\perp\}\}\\mathbf\{A\}\_\{\\zeta\}\\mathbf\{\\Pi\}\_\{\\Psi^\{\\perp\}\}are the restrictions of𝐀ζ\\mathbf\{A\}\_\{\\zeta\}to the subspaces and𝚷Ψ\\mathbf\{\\Pi\}\_\{\\Psi\}and𝚷Ψ⟂\\mathbf\{\\Pi\}\_\{\\Psi^\{\\perp\}\}are the orthogonal projections ontoΨ\\PsiandΨ⟂\\Psi^\{\\perp\}\.

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/Original_Wind_Field_reconstructed_field.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/Centralized_GP_reconstructed_field.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/Centralized_RGP_reconstructed_field.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/Consensus-RGP_reconstructed_field.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/ADMM-RGP_reconstructed_field.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/PDMM-RGP_reconstructed_field.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/colorbar1.png)

Figure 2:Original and reconstructed wind fields\.If𝜻t\[0\]\\boldsymbol\{\\zeta\}\_\{t\}^\{\[0\]\}is initialized randomly, it has been shown that only the component𝜻Ψ,t\[k\]\\boldsymbol\{\\zeta\}^\{\[k\]\}\_\{\\Psi,t\}converges to a fixed point𝜻Ψ,t∗\\boldsymbol\{\\zeta\}^\{\*\}\_\{\\Psi,t\}ask→∞k\\rightarrow\\infty\[[27](https://arxiv.org/html/2609.26979#bib.bib27)\]\. Thus, define the error of the𝜻t\[k\]\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\]\}iterations as𝐞ζ,t\[k\]≜𝜻Ψ,t∗−𝜻Ψ,t\[k\]=𝐀Ψ​\(𝜻Ψ,t∗−𝜻Ψ,t\[k−1\]\)=𝐀Ψk​𝐞ζ,t\[0\],\\mathbf\{e\}\_\{\\zeta,t\}^\{\[k\]\}\\triangleq\\boldsymbol\{\\zeta\}^\{\*\}\_\{\\Psi,t\}\-\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\]\}=\\mathbf\{A\}\_\{\\Psi\}\(\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\*\}\-\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\-1\]\}\)=\\mathbf\{A\}\_\{\\Psi\}^\{k\}\\mathbf\{e\}\_\{\\zeta,t\}^\{\[0\]\},\. Thus the𝜻\\boldsymbol\{\\zeta\}\-iterations is governed by the spectral radius of𝐀Ψ\\mathbf\{A\}\_\{\\Psi\}\. The non\-convergence of theΨ⟂\\Psi^\{\\perp\}component does not affect the convergence of𝝃t\[k\]\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}, as it is annihilated by𝐂ξ⊤\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\. Therefore,𝝃t\[k\]\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}will converge to a fixed point,𝝃t∗\\boldsymbol\{\\xi\}\_\{t\}^\{\*\}\. Now, observe the error of the𝝃t\[k\]\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}iterate is given by

𝐞ξ,t\[k\]\\displaystyle\\mathbf\{e\}\_\{\\xi,t\}^\{\[k\]\}≜𝝃t∗−𝝃t\[k\]=−𝐂¯c​𝐂ξ⊤​𝜻Ψ,t∗\+𝐂¯c​𝐂ξ⊤​𝜻Ψ,t\[k−1\]\\displaystyle\\triangleq\\boldsymbol\{\\xi\}\_\{t\}^\{\*\}\-\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}=\-\\bar\{\\mathbf\{C\}\}\_\{c\}\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\*\}\+\\bar\{\\mathbf\{C\}\}\_\{c\}\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\\boldsymbol\{\\zeta\}\_\{\\Psi,t\}^\{\[k\-1\]\}\(35a\)=−𝐂¯c​𝐂ξ⊤​𝐀Ψk−1​𝐞ζ,t\[0\]\.\\displaystyle=\-\\bar\{\\mathbf\{C\}\}\_\{c\}\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\\mathbf\{A\}\_\{\\Psi\}^\{k\-1\}\\mathbf\{e\}\_\{\\zeta,t\}^\{\[0\]\}\.\(35b\)Observe that the𝝃t\[k\]\\boldsymbol\{\\xi\}\_\{t\}^\{\[k\]\}error is thus governed by the𝜻t\[k\]\\boldsymbol\{\\zeta\}\_\{t\}^\{\[k\]\}error, which converges asymptotically to𝟎\\mathbf\{0\}for anyc\>0c\>0\. To maximize the convergence rate,ccmay be chosen to minimizeρ⁡\(𝐀Ψ\)\\rho\(\\mathbf\{A\}\_\{\\Psi\}\); since this minimization is non\-convex incc, standard convex optimization methods cannot be used, andccmust instead be selected via a search over candidates\. While the spectral radius governs asymptotic convergence, transient behavior may still dominate for smallKK, in which caseccmay instead be selected to minimize‖𝐂ξ⊤​𝐀ΨK−1‖F\\\|\\mathbf\{C\}\_\{\\xi\}^\{\\top\}\\mathbf\{A\}\_\{\\Psi\}^\{K\-1\}\\\|\_\{F\}\.

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/Communication_Graph_1.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/Communication_Graph_2.png)

Figure 1:Communication graphs used in simulations\.

## 5Experiments

The dataset used in the experiments consists of 10\-meter winduu\- andvv\-components indexed by GPS location, obtained from the “ERA5 post\-processed daily statistics on single levels from 1940 to present” dataset provided by the Copernicus Climate Data Store\[[23](https://arxiv.org/html/2609.26979#bib.bib23)\]\. Simulations useN=10N=10agents, whose fixed positions in the 2D wind field are shown in Figure[1](https://arxiv.org/html/2609.26979#S4.F1)with two different communication graphs\. At each time step, each agentnndraws sample locations from𝒩⁡\(𝐩n,σs2​𝐈\)\\mathcal\{N\}\(\\mathbf\{p\}\_\{n\},\\sigma\_\{s\}^\{2\}\\mathbf\{I\}\), where𝐩n\\mathbf\{p\}\_\{n\}is the agent’s position andσs=0\.25\\sigma\_\{s\}=0\.25\. Given a measurement location𝐱\\mathbf\{x\}, the measurement is drawn from𝒩⁡\(𝐟⁡\(𝐱\),𝐑v\)\\mathcal\{N\}\(\\mathbf\{f\}\(\\mathbf\{x\}\),\\mathbf\{R\}\_\{v\}\), where𝐑v=diag​\(0\.01,0\.01\)\\mathbf\{R\}\_\{v\}=\\text\{diag\}\(0\.01,0\.01\)\. Each agent collects2020measurements per time step forT=50T=50time steps, giving10,00010,000training data points\. The set of 400 basis locations𝐗p\\mathbf\{X\}\_\{p\}form a 20\-by\-20 grid over the input space, andQ=2Q=2latent functions are used, with the first function’s kernel hyperparameters pre\-trained on theuuwind speeds and the second trained on thevvwind speeds\. 2,500 prediction locations𝐗∗\\mathbf\{X\}\_\{\*\}, which form a 50\-by\-50 grid over the input space, are used to evaluate the trained model\. Two Laplacian weightings are considered: the unweighted graph Laplacian and the symmetric “optimal” weights obtained by solving the spectral\-norm minimization problem in\[[28](https://arxiv.org/html/2609.26979#bib.bib28)\]\. The parametersα\\alpha,τ\\tau, andccare tuned to optimize convergence speed for each algorithm\.

### 5\.1Performance Metrics

Given the posterior covariances𝚺n,t=𝛀n,t−1\\boldsymbol\{\\Sigma\}\_\{n,t\}=\\boldsymbol\{\\Omega\}\_\{n,t\}^\{\-1\}and means𝝁n,t=𝚺n,t​𝝃n,t\\boldsymbol\{\\mu\}\_\{n,t\}=\\boldsymbol\{\\Sigma\}\_\{n,t\}\\boldsymbol\{\\xi\}\_\{n,t\}of agentnnat timettat the basis locations, the predictive distribution at a set of test locations𝐗∗∈ℝD×P∗\\mathbf\{X\}\_\{\*\}\\in\\mathbb\{R\}^\{D\\times P\_\{\*\}\}is𝒩⁡\(𝝁∗n,t,𝚺∗n,t\)\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{\*n,t\},\\boldsymbol\{\\Sigma\}\_\{\*n,t\}\), where𝝁∗n,t=𝐊⁡\(𝐗∗,𝐗p\)​𝐊𝐗p−1​𝝁n,t\\boldsymbol\{\\mu\}\_\{\*n,t\}=\\mathbf\{K\}\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{p\}\)\\mathbf\{K\}\_\{\\mathbf\{X\}\_\{p\}\}^\{\-1\}\\boldsymbol\{\\mu\}\_\{n,t\}, and𝚺∗n,t=𝐊⁡\(𝐗∗,𝐗∗\)\+𝐊⁡\(𝐗∗,𝐗p\)​\(𝐊𝐗p−1​𝚺n,t−𝐈\)​𝐊𝐗p−1​𝐊​\(𝐗p,𝐗∗\)\\mathbf\{\\Sigma\}\_\{\*n,t\}=\\mathbf\{K\}\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{\*\}\)\+\\mathbf\{K\}\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{p\}\)\(\\mathbf\{K\}\_\{\\mathbf\{X\}\_\{p\}\}^\{\-1\}\\boldsymbol\{\\Sigma\}\_\{n,t\}\-\\mathbf\{I\}\)\\mathbf\{K\}\_\{\\mathbf\{X\}\_\{p\}\}^\{\-1\}\\mathbf\{K\}\(\\mathbf\{X\}\_\{p\},\\mathbf\{X\}\_\{\*\}\)\. Now, letting𝐟∗∈ℝP∗​D′\\mathbf\{f\}\_\{\*\}\\in\\mathbb\{R\}^\{P\_\{\*\}D^\{\\prime\}\}denote the true function values at the test locations, we use the root mean square error \(RMSE\) defined as

RMSE​\(𝐟∗,𝝁∗n,t\)≜\(𝝁∗n,t−𝐟∗\)⊤​\(𝝁∗n,t−𝐟∗\)P∗​D′,\\text\{RMSE\}\(\\mathbf\{f\}\_\{\*\},\\boldsymbol\{\\mu\}\_\{\*n,t\}\)\\triangleq\\sqrt\{\\dfrac\{\(\\boldsymbol\{\\mu\}\_\{\*n,t\}\-\\mathbf\{f\}\_\{\*\}\)^\{\\top\}\(\\boldsymbol\{\\mu\}\_\{\*n,t\}\-\\mathbf\{f\}\_\{\*\}\)\}\{P\_\{\*\}D^\{\\prime\}\}\},\(36\)to evaluate the accuracy of the predictions\. To evaluate the level of consensus among agents, we use mean variance of predictions \(MVOP\) defined as

MVOP​\(𝝁∗1,t,…​𝝁∗N,t\)≜1P∗​D′​∑i=1P∗​D′σi,t2,\\displaystyle\\text\{MVOP\}\(\\boldsymbol\{\\mu\}\_\{\*1,t\},\\dots\\boldsymbol\{\\mu\}\_\{\*N,t\}\)\\triangleq\\frac\{1\}\{P\_\{\*\}D^\{\\prime\}\}\\sum\\limits\_\{i=1\}^\{P\_\{\*\}D^\{\\prime\}\}\\sigma^\{2\}\_\{i,t\},\(37\)whereσi,t2≜1N−1​∑n=1N\(\[𝝁∗n,t\]i−μ¯i,t\)2\\sigma\_\{i,t\}^\{2\}\\triangleq\\frac\{1\}\{N\-1\}\\sum\_\{n=1\}^\{N\}\(\[\\boldsymbol\{\\mu\}\_\{\*n,t\}\]\_\{i\}\-\\bar\{\\mu\}\_\{i,t\}\)^\{2\}is the unbiased sample variance of theit​hi^\{th\}prediction acrossNNagents at timettandμ¯i,t=1N​∑n=1N\[𝝁∗n,t\]i\\bar\{\\mu\}\_\{i,t\}=\\frac\{1\}\{N\}\\sum\_\{n=1\}^\{N\}\[\\boldsymbol\{\\mu\}\_\{\*n,t\}\]\_\{i\}is the sample mean\.

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/ksweep_graph1.png)

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/ksweep_graph2.png)

Figure 3:Plots of consensus metric MVOP vs\. number of communication roundsKKusing communication graph 1 \(top\) and 2 \(bottom\) andT=20T=20time steps\.![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/alg_conn_sweep.png)Figure 4:Plot of consensus metric MVOP vs\. algebraic connectivity of the communication graph withT=20T=20\.
### 5\.2Results and Discussion

Consensus\-RGP, ADMM\-RGP, PDMM\-RGP, and centralized RGP were all simulated across 100 Monte Carlo \(MC\) simulations with different random seeds\. All four algorithms achieved an average RMSE of0\.08620\.0862with a 95% confidence interval of\[0\.0860,0\.0864\]\[0\.0860,0\.0864\]\. Figure[2](https://arxiv.org/html/2609.26979#S4.F2)shows the reconstructed fields from one simulation alongside that of a centralized, non\-recursive GP\. The high number of basis locations and large grid used here were chosen deliberately to focus on communication efficiency effects, which is the goal of this work\. Hence, a systematic sweep over basis location count is therefore outside the scope of this study\. The behavior of these methods under a reduced basis location regime — where distributed and centralized solutions are known to diverge — has been characterized in\[[8](https://arxiv.org/html/2609.26979#bib.bib8)\]\.

Figure[3](https://arxiv.org/html/2609.26979#S5.F3)shows the averagelog10⁡\(MVOP\)\\log\_\{10\}\(\\text\{MVOP\}\)versusKK, the number of communication rounds per time step, over 100 MC runs\. For communication graph 1 \(top\), ADMM\-RGP achieves lower MVOP values than the other algorithms for any givenKK\. In particular, ADMM\-RGP withK=7K=7achieves a consensus level comparable to Consensus\-RGP withK=10K=10, reducing communication by 30% without sacrificing performance\. Although PDMM\-RGP performs less favorably on graph 1, it outperforms the other algorithms on graph 2, achieving MVOP values up to three orders of magnitude lower than Consensus\-RGP\. On graph 2, both PDMM\-RGP and ADMM\-RGP achieve performance withK=5K=5comparable to Consensus\-RGP withK=10K=10, corresponding to a 50% reduction in communication\. Figure[4](https://arxiv.org/html/2609.26979#S5.F4)plots the averagelog10⁡\(MVOP\)\\log\_\{10\}\(\\text\{MVOP\}\)over 100 MC runs versus the algebraic connectivity of 20 different communication graphs\. PDMM\-RGP outperforms the other algorithms in sparse graphs, while ADMM\-RGP consistently outperforms Consensus\-RGP and performs particularly well in densely connected graphs\. Figure[5](https://arxiv.org/html/2609.26979#S5.F5)shows box plots of computation time per communication round over 250 rounds as a function of the agent’s number of neighbors\. All three algorithms exhibit linear scaling, although PDMM\-RGP has a steeper increase\. While Consensus\-RGP has the lowest per\-round computational complexity, ADMM\-RGP and PDMM\-RGP reduce communication and can reduce total computation by requiring fewer rounds\.

![Refer to caption](https://arxiv.org/html/2609.26979v1/Figures/comp_time.png)Figure 5:Box plots of computation time from 250 communication rounds using communication graph 1\. Markers denote the median times\.

## 6Conclusion

This paper presented two distributed, recursive, multi\-output GP algorithms with reduced communication and computational complexity\. We analyzed their convergence and provided practical guidance for parameter selection to ensure fast convergence\. Experiments show the proposed algorithms require fewer resources than the state of the art while achieving comparable reconstruction performance, with PDMM\-RGP more effective on sparse graphs and ADMM\-RGP on denser graphs\. Future work will consider asynchronous communication, time\-varying functions and communication graphs, and extensions to other GP representations such as random\-feature and multiple\-model formulations\. The proposed fusion framework may also be applied to distributed trajectory and dynamical state estimation, as well as active sensing and cooperative control problems\.\\appendices

## 7Proof of Theorem[4\.4](https://arxiv.org/html/2609.26979#Thmtheorem4)

Using the propertydet​\(𝐀\)=∏iλi​\(𝐀\)\\text\{det\}\(\\mathbf\{A\}\)=\\prod\_\{i\}\\lambda\_\{i\}\(\\mathbf\{A\}\)\[[29](https://arxiv.org/html/2609.26979#bib.bib29)\], we havedet\(𝐌i\)=−τ​λi​\(𝐋\)=λ1​\(𝐌i\)​λ2​\(𝐌i\)\\det\(\\mathbf\{M\}\_\{i\}\)=\-\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)=\\lambda\_\{1\}\(\\mathbf\{M\}\_\{i\}\)\\lambda\_\{2\}\(\\mathbf\{M\}\_\{i\}\)\. Becauseλi​\(𝐋\)\\lambda\_\{i\}\(\\mathbf\{L\}\)is positive,0<−τ​λi​\(𝐋\)<10<\-\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)<1\. If the eigenvalues of𝐌i\\mathbf\{M\}\_\{i\}are real, the characteristic equation,

pi​\(λ\)=λ2−\(1−\(α\+τ\)​λi​\(𝐋\)\)​λ−τ​λi​\(𝐋\)=0,p\_\{i\}\(\\lambda\)=\\lambda^\{2\}\-\(1\-\(\\alpha\+\\tau\)\\lambda\_\{i\}\(\\mathbf\{L\}\)\)\\lambda\-\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)=0,has two real roots, which satisfy\|λ1​λ2\|=\|−τ​λi​\(𝐋\)\|<1\|\\lambda\_\{1\}\\lambda\_\{2\}\|=\|\-\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\)\|<1\. Therefore, both roots have magnitude<1<1provided that neither one crosses the±1\\pm 1boundary\. Evaluating the characteristic equation at11and−1\-1, we find

pi​\(1\)=α​λi​\(𝐋\)\>0,pi​\(−1\)=2−\(α\+2​τ\)​λi​\(𝐋\)\>0,p\_\{i\}\(1\)=\\alpha\\lambda\_\{i\}\(\\mathbf\{L\}\)\>0,\\qquad p\_\{i\}\(\-1\)=2\-\(\\alpha\+2\\tau\)\\lambda\_\{i\}\(\\mathbf\{L\}\)\>0,under the assumed bounds onτ\\tau\. Therefore, both eigenvalues lie strictly within\(−1,1\)\(\-1,1\)\. In the case that the eigenvalues are complex conjugate pairs of the forma±b​ja\\pm bj, we havedet\(𝐌i\)=\(a\+b​j\)​\(a−b​j\)=a2\+b2\\det\(\\mathbf\{M\}\_\{i\}\)=\(a\+bj\)\(a\-bj\)=a^\{2\}\+b^\{2\}with squared magnitude\|λ1​\(𝐌i\)\|2=\|λ2​\(𝐌i\)\|2=a2\+b2\|\\lambda\_\{1\}\(\\mathbf\{M\}\_\{i\}\)\|^\{2\}=\|\\lambda\_\{2\}\(\\mathbf\{M\}\_\{i\}\)\|^\{2\}=a^\{2\}\+b^\{2\}\. Therefore, if the eigenvalues of𝐌i\\mathbf\{M\}\_\{i\}are complex, they have squared magnitudes equal todet\(𝐌i\)=−τ​λi​\(𝐋\)\\det\(\\mathbf\{M\}\_\{i\}\)=\-\\tau\\lambda\_\{i\}\(\\mathbf\{L\}\), which we have established is in the range\(0,1\)\(0,1\)\. Thus, we have shown that the eigenvalues of𝐌i\\mathbf\{M\}\_\{i\}lie strictly inside of the unit circle and thus𝐌\\mathbf\{M\}is Schur stable\.

## References

- \[1\]R\. N\. Darmanin and M\. K\. Bugeja, “A review on multi\-robot systems categorised by application domain,” in*2017 25th Mediterranean Conference on Control and Automation \(MED\)*, Jul\. 2017, pp\. 701–706, iSSN: 2473\-3504\.
- \[2\]M\. Davoodi, S\. Faryadi, and J\. M\. Velni, “A Graph Theoretic\-Based Approach for Deploying Heterogeneous Multi\-agent Systems with Application in Precision Agriculture,”*Journal of Intelligent & Robotic Systems*, vol\. 101, no\. 1, p\. 10, Dec\. 2020\.
- \[3\]O\. P\. Mahela, M\. Khosravy, N\. Gupta, B\. Khan, H\. H\. Alhelou, R\. Mahla, N\. Patel, and P\. Siano, “Comprehensive overview of multi\-agent systems for controlling smart grids,”*CSEE Journal of Power and Energy Systems*, vol\. 8, no\. 1, pp\. 115–131, Jan\. 2022\.
- \[4\]A\. Dorri, S\. S\. Kanhere, and R\. Jurdak, “Multi\-Agent Systems: A Survey,”*IEEE Access*, vol\. 6, pp\. 28 573–28 593, 2018\.
- \[5\]G\. P\. Kontoudis and D\. J\. Stilwell, “Scalable, Federated Gaussian Process Training for Decentralized Multi\-Agent Systems,”*IEEE Access*, vol\. 12, pp\. 77 800–77 815, 2024\.
- \[6\]T\. Ding, R\. Zheng, S\. Zhang, and M\. Liu, “Resource\-Efficient Cooperative Online Scalar Field Mapping via Distributed Sparse Gaussian Process Regression,”*IEEE Robotics and Automation Letters*, vol\. 9, no\. 3, pp\. 2295–2302, Mar\. 2024\.
- \[7\]K\. Jakkala and S\. Akella, “Multi\-Robot Informative Path Planning from Regression with Sparse Gaussian Processes,” in*2024 IEEE International Conference on Robotics and Automation \(ICRA\)*, May 2024, pp\. 12 382–12 388\.
- \[8\]Y\. P\. K\. Rao, T\. Keviczky, and R\. T\. Rajan, “Consensus\-based recursive multi\-output gaussian process,” in*2026 60th Asilomar Conference on Signals, Systems, and Computers*\. Pacific Grove, CA, USA: IEEE, 2026, accepted for publication\.
- \[9\]A\. E\. Balcı and R\. T\. Rajan, “Multiple Model Recursive Gaussian Process for Robust Target Tracking,”*IEEE Open Journal of Signal Processing*, pp\. 1–9, 2025\.
- \[10\]A\. Lederer, Z\. Yang, J\. Jiao, and S\. Hirche, “Cooperative Control of Uncertain Multiagent Systems via Distributed Gaussian Processes,”*IEEE Transactions on Automatic Control*, vol\. 68, no\. 5, pp\. 3091–3098, May 2023\.
- \[11\]Y\. Xu, F\. Yin, W\. Xu, J\. Lin, and S\. Cui, “Wireless Traffic Prediction With Scalable Gaussian Process: Framework, Algorithms, and Verification,”*IEEE Journal on Selected Areas in Communications*, vol\. 37, no\. 6, pp\. 1291–1306, Jun\. 2019\.
- \[12\]Z\. Zhu, E\. Biyik, and D\. Sadigh, “Multi\-Agent Safe Planning with Gaussian Processes,” in*2020 IEEE/RSJ International Conference on Intelligent Robots and Systems \(IROS\)*\. Las Vegas, NV, USA: IEEE, Oct\. 2020, pp\. 6260–6267\.
- \[13\]L\. Gupta, R\. Jain, and G\. Vaszkun, “Survey of Important Issues in UAV Communication Networks,”*IEEE Communications Surveys & Tutorials*, vol\. 18, pp\. 1–1, Nov\. 2015\.
- \[14\]J\. Gielis, A\. Shankar, and A\. Prorok, “A Critical Review of Communications in Multi\-robot Systems,”*Current Robotics Reports*, vol\. 3, no\. 4, pp\. 213–225, 2022\.
- \[15\]M\. A\. Razzaque and S\. Dobson, “Energy\-Efficient Sensing in Wireless Sensor Networks Using Compressed Sensing,”*Sensors \(Basel, Switzerland\)*, vol\. 14, no\. 2, pp\. 2822–2859, Feb\. 2014\.
- \[16\]A\. Xie, F\. Yin, Y\. Xu, B\. Ai, T\. Chen, and S\. Cui, “Distributed Gaussian Processes Hyperparameter Optimization for Big Data Using Proximal ADMM,”*IEEE Signal Processing Letters*, vol\. 26, no\. 8, pp\. 1197–1201, Aug\. 2019\.
- \[17\]P\. Zhai and R\. T\. Rajan, “Distributed Gaussian Process Hyperparameter Optimization for Multi\-Agent Systems,” in*ICASSP 2023 \- 2023 IEEE International Conference on Acoustics, Speech and Signal Processing \(ICASSP\)*, Jun\. 2023, pp\. 1–5\.
- \[18\]M\. F\. Huber, “Recursive Gaussian process: On\-line regression and learning,”*Pattern Recognition Letters*, vol\. 45, pp\. 85–91, Aug\. 2014\.
- \[19\]F\. Llorente, D\. Waxman, and P\. M\. Djurić, “Decentralized Online Ensembles of Gaussian Processes for Multi\-Agent Systems,” in*ICASSP 2025 \- 2025 IEEE International Conference on Acoustics, Speech and Signal Processing \(ICASSP\)*, Apr\. 2025, pp\. 1–5, iSSN: 2379\-190X\.
- \[20\]F\. Llorente, D\. Waxman, S\. Jantre, N\. M\. Urban, and S\. E\. Minkoff, “Robust, Online, and Adaptive Decentralized Gaussian Processes,” in*ICASSP 2026 \- 2026 IEEE International Conference on Acoustics, Speech and Signal Processing \(ICASSP\)*, May 2026, pp\. 22 487–22 491, iSSN: 2379\-190X\.
- \[21\]M\. Iqbal, K\. Kumar, and S\. Särkkä, “Communication\-Efficient Distributed Kalman Filtering Using ADMM,”*IEEE Transactions on Automatic Control*, vol\. 71, no\. 3, pp\. 1916–1923, Mar\. 2026\.
- \[22\]G\. Zhang and R\. Heusdens, “Distributed Optimization Using the Primal\-Dual Method of Multipliers,”*IEEE Transactions on Signal and Information Processing over Networks*, vol\. 4, no\. 1, pp\. 173–187, Mar\. 2018\.
- \[23\]Copernicus Climate Change Service, “ERA5 post\-processed daily statistics on single levels from 1940 to present,” Oct\. 2024\.
- \[24\]M\. A\. Álvarez, L\. Rosasco, and N\. D\. Lawrence, “Kernels for Vector\-Valued Functions: A Review,”*Foundations and Trends in Machine Learning*, vol\. 4, no\. 3, pp\. 195–266, Mar\. 2012\.
- \[25\]T\. W\. Sherson, R\. Heusdens, and W\. B\. Kleijn, “Derivation and Analysis of the Primal\-Dual Method of Multipliers Based on Monotone Operator Theory,”*IEEE Transactions on Signal and Information Processing over Networks*, vol\. 5, no\. 2, pp\. 334–347, Jun\. 2019\.
- \[26\]R\. Heusdens and G\. Zhang, “Distributed Optimisation Via The Generalised Primal\-Dual Method of Multipliers Under Unreliable and Quantised Communication,” in*ICASSP 2026 \- 2026 IEEE International Conference on Acoustics, Speech and Signal Processing \(ICASSP\)*\. Barcelona, Spain: IEEE, May 2026, pp\. 281–285\.
- \[27\]Q\. Li, R\. Heusdens, and M\. G\. Christensen, “Communication efficient privacy\-preserving distributed optimization using adaptive differential quantization,”*Signal Processing*, vol\. 194, p\. 108456, May 2022\.
- \[28\]L\. Xiao and S\. Boyd, “Fast linear iterations for distributed averaging,”*Systems & Control Letters*, vol\. 53, no\. 1, pp\. 65–78, Feb\. 2004\.
- \[29\]K\. B\. Petersen and M\. S\. Pedersen, “The Matrix Cookbook,” Nov\. 2012\.