1 Introduction↩︎

Networked AI systems, such as distributed sensing networks, edge intelligence, federated learning systems, and multi-agent autonomous platforms, require multiple agents to collaboratively learn and adapt models over communication networks [1]–[3]. A central challenge in these systems is how to enable each agent to update its local model efficiently while maintaining global coordination with other agents. This challenge becomes more pronounced when the learning objective has a hierarchical structure, where model parameters, hyperparameters, data weights, or task-specific adaptation variables must be optimized jointly. Such hierarchical learning problems can often be formulated as bilevel optimization problems. Bilevel optimization (BO) has attracted increasing attention in recent years. It has found important applications in meta-learning [4]–[6], hyperparameter optimization [4], [7], [8], reinforcement learning [9], and adversarial learning [10], [11]. To achieve state-of-the-art results in these domains, vast training datasets and distributed algorithms operating across multiple computing agents are often required.

This paper focuses on the following decentralized bilevel optimization (DBO) problem: \[\begin{align} \label{DBO} \mathop{\min}\limits_{x \in \mathbb{R}^d} & \;\Phi(x):= F(x,y^*(x)) := \frac{1}{n}\sum_{i=1}^{n}F_i(x,y^*(x))\quad\nonumber\\ \mathrm{s.t.} \; &y^*(x) \in \mathop{\textrm{argmin}}\limits_{y\in \mathbb{R}^q} f(x,y) := \frac{1}{n}\sum_{i=1}^{n}f_i(x,y), \end{align}\tag{1}\] where \(n\) collaborative nodes, connected via a communication graph \(\mathcal{G} = (\mathcal{V}, \mathcal{E})\), jointly solve this problem. Here, \(\mathcal{V}\) represents the set of agents, and \(\mathcal{E}\) represents the set of feasible communication links between the agents. Each node \(i\) privately holds an upper-level cost function \(F_i:\mathbb{R}^d\times\mathbb{R}^q\rightarrow\mathbb{R}\) and a lower-level cost function \(f_i:\mathbb{R}^d\times\mathbb{R}^q\rightarrow\mathbb{R}\). The decentralized framework eliminates the need for a central server and enables more efficient and scalable computation across distributed networks [12]. In this paper, we propose an efficient single-loop algorithm for solving 1 .

1.1 Related work↩︎

1.1.0.1 Bilevel optimization

When utilizing gradient methods for solving 1 , a central challenge is how to estimate the hypergradient \(\nabla \Phi(x)\), which is given by \[\label{hg} \small \begin{align} &\nabla\Phi(x)=\frac{1}{n}\sum_{i=1}^{n}\left[\nabla_1 F_i\left(x, y^*(x)\right) - \nabla^2_{12} f_i\left(x, y^*(x)\right) z(x)\right],\\ &z(x) = \left[\frac{1}{n}\sum_{i=1}^{n}\nabla_{22}^2 f_i\left(x, y^*(x)\right)\right]^{-1} \frac{1}{n}\sum_{i=1}^{n}\nabla_2 F_i\left(x, y^*(x)\right). \end{align}\tag{2}\] To address this challenge, various algorithms have been proposed in the literature, leveraging approaches such as approximate implicit differentiation (AID) [9], [13]–[18], iterative differentiation (ITD) [4], [6], [14]–[16], [19], [20] and Neumann series-based approach [13], [17]. Furthermore, both AID and ITD approaches involve heavy Hessian- and Jacobian-vector multiplications, which also consist of a complex double-loop structure. Recently, [18] introduced a single-loop framework, named SOBA, to solve BO with single agent, i.e., \(n=1\) in 1 . This framework approximates \(z(x)\) without requiring heavy matrix-vector multiplications, thus saving effort in solving the linear system. Specifically, the SOBA algorithm maintains three sequences, which are updated as follows: \[\begin{align} y^{k+1} = y^k - \beta_k D_y^k, v^{k+1} = v^k + \eta_k D_v^k, x^{k+1} = x^k - \alpha_k D_x^k, \end{align}\] where \(\alpha_k\), \(\beta_k\) and \(\eta_k\) are positive stepsizes, \(D_y^k, D_v^k\) and \(D_x^k\) are respectively unbiased stochastic estimators of \(\nabla_2f(x^k,y^k), \nabla_2F(x^k,y^k) - \nabla^2_{22}f(x^k,y^k)v^k\) and \(\nabla_1F(x^k,y^k) - \nabla^2_{12}f(x^k,y^k)v^k\). The SOBA framework was later extended by [21] to address the case where the lower-level problem is merely convex. [22] addressed the non-convexity of the lower-level function and proposed two techniques to tackle the resulting challenges.

1.1.0.2 Decentralized bilevel optimization

In the past few years, a series of decentralized optimization algorithms have been proposed to solve 1 . The first algorithm for solving 1 was proposed in [23]. They introduced the DSBO algorithm, which integrates a decentralized approach to solve the linear system in 2 . Subsequently, the same authors improved the DSBO algorithm by incorporating the moving average technique [24]. Other works on DBO include [25], in which a gossip-based DBO algorithm was proposed and its sample complexity was established, and [26], in which momentum and variance-reduced techniques were employed to solve distributed stochastic BO problems. All of these algorithms rely on inner-loop updates to assess the lower-level solution as well as the Hessian inverse, which demands extensive computational resources and leads to significant communication overhead, thus greatly reducing their practical applicability. [27] proposed a stochastic linearized augmented Lagrangian method for solving 1 , which, however, also requires heavy matrix-vector multiplications when approximating the hypergradient. Recently, [28] extended the SOBA algorithm [18] to the decentralized setting and proposed the first single-loop decentralized algorithm, named SLDBO, for solving 1 . Similar to SOBA, SLDBO does not require heavy matrix-vector multiplications when approximating the hypergradient. Additionally, SLDBO imposes no assumptions regarding gradient heterogeneity or the boundedness of \(\{v_k\}\). [29] introduced a decentralized SOBA with the moving average technique and established its transient iteration complexity. Similarly, [30] employed variance reduction and gradient tracking techniques to address the stochastic DBO problem. However, the consensus error was not taken into account in the convergence analysis. We point out that these single-loop algorithms still need to calculate gradients, Jacobian and Hessian matrices in each iteration, which can be time consuming for large-scale problems.

In this paper, we incorporate the snapshot gradient tracking (SSGT) technique [31] – a recent technique for decentralized optimization – to SLDBO and demonstrate that the new algorithm   is much more efficient than SLDBO. Unlike traditional methods that track the average gradient of decision variables, the SSGT tracks the average gradient of a snapshot point incorporating historical information. This technique employs intermittent gradient, Jacobian, and Hessian computations, offering a more effective approach to utilizing computational resources.

1.2 Main results↩︎

Contributions of this paper are summarized as follows.

• We propose a novel single-loop algorithm named Snapshot-SLDBO (), which integrates the snapshot idea from SSGT [31] into SLDBO [28] and intermittently skips the computations of gradients, Jacobians and Hessians, thus significantly reducing the overall computational cost.

• We analyze the iteration complexity of , demonstrating that it achieves \(\mathcal{O}(1/K)\) convergence rate, where \(K\) is the iteration counter. Moreover, we establish a high probability non-ergodic convergence result for .

• We validate the proposed method on hyperparameter optimization, data hyper-cleaning, and decentralized meta-learning tasks. The results demonstrate improved time efficiency and competitive learning performance in networked AI scenarios. Additionally, we also study the numerical impact of network topology and data heterogeneity on the performance of .

While   builds upon the foundational ideas of SLDBO and SSGT, their integration into a coherent and provably convergent algorithm presents non-trivial challenges. The primary difficulty lies in adapting the snapshot mechanism—originally conceived for single-level optimization—to the more intricate, three-sequence dynamic (\(x, y, v\)) inherent in bilevel algorithms like SLDBO. Our core contribution is therefore not merely the algorithmic proposal, but the rigorous convergence analysis that validates this complex integration. We demonstrate that despite intermittent computations across coupled sequences, the algorithm maintains a solid theoretical convergence guarantee, a key step toward making DBO more practical for large-scale problems. We also note that our current analysis is focused on the deterministic setting, where full gradients are accessible; extending this framework to the stochastic realm remains an important direction for future investigation.

1.3 Organization↩︎

The remainder of this paper is organized as follows. In Section 2, we introduce the notation, assumptions, and essential preliminaries which will be utilized in subsequent analysis. In Section 3, we propose the algorithm and analyze its convergence. Numerical results on hyperparameter optimization, data hyper-cleaning, and decentralized meta-learning are presented in Section 4 to demonstrate the effectiveness of S\(^3\)LDBO. Finally, concluding remarks are drawn in Section 5.

2 sec:Preliminaries↩︎

We denote the optimal value of 1 as \(F^*\). The gradients of \(f\) with respect to \(x\) and \(y\) are represented by \(\nabla_1f(x,y)\) and \(\nabla_2 f(x,y)\) respectively, and the Jacobian matrix of \(\nabla_1f\) and Hessian matrix of \(f\) with respect to \(y\) are denoted as \(\nabla^2_{12}f(x,y)\) and \(\nabla^2_{22}f(x,y)\) respectively. Unless otherwise specified, \(\|\cdot\|\) refers to the \(\ell_2\) norm for vectors and the Frobenius norm for matrices. The operator norm of a matrix \(Z\) is denoted by \(\|Z\|_{\textrm{op}}\). We also define the projection operator \(\mathcal{P}_{r}\), which projects a given point onto a Euclidean ball with radius \(r\geq 0\), i.e., \[\mathcal{P}_{r}[z]:= \arg\min\nolimits_{\|z'\|\leq r}\|z'-z\|=\min \left\{1, r/\|z\|\right\} z.\]

Throughout this paper, we adhere to the following standard assumptions in the literature of bilevel optimization and decentralized optimization problems.

The following results are instrumental in our convergence analysis. Lemma 1 is a well-known result in decentralized optimization, as noted in [32]. Lemma 2, adopted from Lemma 3.2 in [33], elucidates the relationship between the consensus error before and after the projection.

3 Our  Algorithm↩︎

In this section, we propose the Snapshot-SLDBO algorithm (), which improves the SLDBO algorithm [28] with the idea of SSGT [31]. We start with discussing related algorithms and explaining the derivation of .

3.1 Algorithm development↩︎

The SSGT algorithm [31] is designed for solving the following single-level decentralized problem: \[\min\nolimits_{x} h(x) = \frac{1}{n}\sum\nolimits_{i=1}^n h_i(x),\] where \(x\) is the global decision variable and each \(h_i\) is a smooth and strongly convex function. Let \(\{(\xi^k,\zeta^k): k=0,1,2,\ldots\}\) be a sequence of two-point random variables with \(\xi^k \sim Bernoulli(p)\) and \(\zeta^k \sim Bernoulli(l)/l\), where \(p,l\in(0,1)\). With the initial values \((y^0_i,z^0_i,u^0_i,s_i^0)\) and \(g_i^0 = \nabla h_i(s_i^0)\), \(i=1,\ldots,n\), SSGT iterates for \(k=0,1,2,\ldots\) as \[\begin{align} x_i^{k} & = (1-\alpha-\tau) y_i^{k} + \alpha z_i^{k} + \tau u_i^{k}, \\ z_i^{k+1} & = \frac{1}{1+\beta} \sum\nolimits_{j=1}^n w_{ij} \biggl\{z_j^{k} + \beta x_j^{k} \\ &\quad \;- \eta \Bigl[g_j^{k} + \zeta^{k} \left(\nabla h_j(x_j^{k}) - \nabla h_j(s_j^{k})\right)\Bigl]\biggl\},\\ y_i^{k+1} & = x_i^{k} + \gamma \left(z_i^{k+1} - z_i^{k}\right), \\ s_i^{k+1} & = \xi^{k} x_i^{k} + (1-\xi^{k}) s_i^{k}, \\ u_i^{k+1} & = \sum\nolimits_{j=1}^n w_{ij} \left(\xi^{k} x_j^{k} + (1-\xi^{k}) u_j^{k}\right), \\ g_i^{k+1} & = \sum\nolimits_{j=1}^n w_{ij} g_j^{k} + \xi^{k} \left(\nabla h_i(x_i^{k}) - \nabla h_i(s_i^{k})\right), \end{align}\] where \(\alpha\), \(\beta\), \(\gamma\), \(\eta\) and \(\tau\) are positive constants. In other gradient-tracking based algorithms, when \(\nabla h_i(x_i^{k+1}) - \nabla h_i(x_i^k)\) is large, a small stepsize is taken to control the consensus error, ultimately resulting in suboptimal convergence rates. The introduction of the snapshot point \(s_i^k\) is designed to overcome this limitation and mitigate its impact on convergence rates. Meanwhile, the update of \(x_i^k\) is expressed as a linear combination of three other variables \(y_i^k, z_i^k\), and \(u_i^k\), ensuring that the distance between \(x_i^k\) and \(s_i^k\) is not too large, and \(\nabla h_i(x_i^k) - \nabla h_i(s_i^k)\) is not too large. Most importantly, the SSGT algorithm skips the gradient computation when \(\xi^k = \zeta^k = 0\), which greatly reduces the computational burden.

The SLDBO algorithm [28] for solving 1 updates in the \(k\)th iteration as follows: \[\begin{align} &t_{y, i}^{k+1}=\sum\nolimits_{j=1}^n w_{i j} t_{y, j}^{k}+d_{y,i}^{k+1}-d_{y, i}^{k}, \\ &y_i^{k+1}=\sum\nolimits_{j=1}^n w_{i j} (y_j^k-\beta t_{y, j}^k), \tag{3}\\ &t_{v, i}^{k+1}=\sum\nolimits_{j=1}^n w_{i j} t_{v, j}^{k}+d_{v,i}^{k+1}-d_{v, i}^{k}, \\ &v_i^{k+1}=\mathcal{P}_{r_v}\left[\sum\nolimits_{j=1}^n w_{i j} (v_j^k+\eta t_{v, j}^k)\right], \tag{4}\\ &t_{x, i}^{k+1}=\sum\nolimits_{j=1}^n w_{i j}t_{x, j}^{k}+d_{x,i}^{k+1}-d_{x, i}^{k}, \\ &x_i^{k+1}=\sum\nolimits_{j=1}^n w_{i j} (x_j^k-\alpha t_{x, j}^k), \tag{5} \end{align}\] where \(r_v:= L_{F,0}/\sigma\), and \(d_{y,i}^k\), \(d_{v,i}^k\) and \(d_{x,i}^k\) are defined in [yd]-[xd], respectively. In each iteration, the SLDBO algorithm needs to evaluate one gradient, one Jacobian, and one Hessian of the lower-level function, and two gradients of the upper-level function. This brings heavy computational burden for large-scale problems. To alleviate the burden caused by computing gradients, Jacobians and Hessians in each iteration, we propose the  algorithm, which skips these computations intermittently using the SSGT idea. The details of our  algorithm are described in Algorithm 1. Inspired by SSGT, we also introduce a stochastic variable \(\xi^k\) in the \(k\)th iteration. It should be noted that whenever the random variable \(\xi^k=0\), the algorithm does not compute any gradients, Jacobians, and Hessians. This leads to substantial savings in computational time when dealing with large-scale problems. In particular, as shown in [28], the SLDBO algorithm computes three gradients, one Jacobian, and one Hessian in each iteration. Assuming that the amount of CPU time required for these computations is \(\mathcal{T}\) in each iteration, this totals \(\mathcal{T}K\) for \(K\) iterations. In contrast, only needs an expected computational cost of \(p\mathcal{T}K\) for the same number of iterations, where \(p\in(0,1)\). This significantly reduces computational time.

3.2 Convergence results↩︎

Next, we present convergence results of . First, we recall the definition of a stationary point as presented in [34].

Our primary convergence rate results for Algorithm 1 are summarized in Theorem 1 and Theorem 2, and their proofs are postponed to the Appendix 6.

4 Numerical experiments↩︎

In this section, we evaluate the proposed S\(^3\)LDBO algorithm on three representative decentralized bilevel learning tasks, including hyperparameter optimization, data hyper-cleaning, and decentralized meta-learning. The first two tasks are conducted on synthetic data, MNIST, and corrupted Fashion-MNIST, while the meta-learning experiment is conducted on miniImageNet. These tasks cover both classical bilevel learning problems and neural-network-based multi-agent learning scenarios, and are used to assess the computational efficiency and learning performance of the proposed snapshot mechanism.

For hyperparameter optimization, we compare S\(^3\)LDBO with SLDBO [28], MA-DSBO [24], and SLAM [27]. SLDBO is the most directly related single-loop DBO baseline, MA-DSBO is a representative double-loop method, and SLAM is a stochastic linearized augmented Lagrangian method for DBO. For data hyper-cleaning, we mainly compare S\(^3\)LDBO with SLDBO in order to isolate the effect of the snapshot mechanism against its closest single-loop counterpart. For decentralized meta-learning, we further include MAML [35] and ANIL [36] as representative meta-learning baselines. Other DBO algorithms are not included because they are designed for different stochastic settings, relying on different assumptions, or focusing on communication complexity rather than the local derivative-evaluation cost considered here. A major goal of the experiments is to isolate the computational savings achieved by the snapshot mechanism relative to the most direct single-loop counterpart, while also examining whether such savings translate into improved time efficiency in neural-network-based meta-learning.

Throughout all experiments, we set the number of nodes to \(n=8\) and adopt a ring topology to model the decentralized communication network. The corresponding weight matrix \(W=(w_{ij})\in\mathbb{R}^{n\times n}\) is defined as follows: for \(i,j=1,\ldots,n\), \(w_{ij}=w\) if \(i=j\); \(w_{ij}=(1-w)/2\) if \(i=j\pm 1\) or \((i,j)\in\{(1,n),(n,1)\}\); and \(w_{ij}=0\) otherwise. We set \(w=0.4\), under which each node is connected to two neighboring nodes and the associated spectral constant is \(\rho\approx 0.724\). Unless otherwise specified, all results are reported with respect to CPU time or the number of expensive local derivative evaluations, so as to directly reflect the computational savings brought by the snapshot mechanism.

All the experiments were performed within Python 3.8 running on a laptop with AMD Ryzen 5 3550H CPU at 2.1GHz with 16GB memory. We employed mpi4py [37] for parallel computing.

4.1 Hyperparameter optimization↩︎

4.1.1 Synthetic data↩︎

We first carry out numerical experiments on the logistic regression problem with \(\ell_2\) regularization. Let \(\psi(t) = \log (1+e^{-t})\) for \(t\in\mathbb{R}\) and \(d\) be the dimension of the data. Following [24] and [28], for each node \(i\), we have \[\begin{align} F_i(\lambda, \omega)&=\sum_{\left(x_e, y_e\right) \in \mathcal{D}_i^{\prime}} \psi\left(y_e x_e^{\top} \omega\right) ,\\ f_i(\lambda, \omega)&=\sum_{\left(x_e, y_e\right) \in \mathcal{D}_i} \psi\left(y_e x_e^{\top} \omega\right)+\frac{1}{2} \sum\nolimits_{j=1}^d e^{\lambda_j} \omega_j^2, \end{align}\] where \(\mathcal{D}_i\) and \(\mathcal{D}_i^{\prime}\) denote the training and testing datasets on node \(i\), respectively. We aim to determine the optimal hyperparameter \(\lambda\) such that \(\omega^*(\lambda)\) minimizes the logistic loss evaluated on the testing dataset.

We utilize synthetic heterogeneous data generated in the same manner as in [24] and [28]. Specifically, the data distribution of \(x_e\) on node \(i\) follows a normal distribution with mean \(0\) and variance \(i^2 \cdot r^2\), where \(r\) is the heterogeneity rate. In our experiments, we set \(r\) to \(1\). For the corresponding response variable, we let \(y_e=x_e^{\top} w+0.1 z\), where \(z\) is sampled from the standard normal distribution. The training and testing datasets used in all algorithms each consist of 20,000 samples. In our experiments, we compare   with SLDBO [28], MA-DSBO [24] and SLAM[27]. In SLDBO and , the parameter \(r_v\) is set to 20, which can be chosen via grid search, and the values of both \({\alpha}\) and \({\eta}\) are assigned as 0.5. Furthermore, we set \({\beta}\) to 1.2 and the maximum number of iterations is set to \(10^3\), which are the same as in [28]. MA-DSBO utilizes two key parameters, namely \(T\) and \(N\). Here, \(T\) stands for the number of iterations carried out in the inner loop, while \(N\) represents the number of Hessian-inverse-gradient product iterations. MA-DSBO requires sufficient inner-loop iterations to estimate the LL solution and the hypergradient accurately. SLAM has one key parameter \(b\), which is associated with matrix-vector products used for estimating the inverse of the Hessian. In our experiments, we set \(T=5\), \(N=5\) and \(b=5\).

To visualize the performance of the three algorithms being compared, we present the curves of train loss, test loss, classification accuracy, and hypergradient norm against CPU time (in seconds) in Fig. 2 and Fig. 3, corresponding to \(d = 60\) and \(d = 300\) respectively. As the dimension of the data increases, the computational demands for calculating the gradients, Jacobian, and Hessian matrices rise significantly. The algorithm attains higher numerical efficiency due to its ability to skip certain computationally intensive steps, thereby enhancing the overall performance. In this experiment, we tested two values of \(p\), namely \(p=0.3\) and \(p=0.75\). These values were chosen to represent a low-frequency (\(p=0.3\)) versus a high-frequency (\(p=0.75\)) update regime. This allows us to empirically study the trade-off between per-iteration cost and convergence speed discussed. It can be seen that \(p = 0.3\) yields a better performance compared to \(p = 0.75\). Furthermore, by comparing Fig. 2 and Fig. 3, it can be seen that the advantage of setting \(p = 0.3\) becomes clearer for problems with larger dimensions compared to when \(p = 0.75\).

Fig. 4 and Fig. 5 display the comparison results between two single-loop algorithms,  and SLDBO, with respect to train loss, test loss, accuracy, and hypergradient magnitude versus the number of lower-level function’s gradient evaluations, where Fig. 4 corresponds to \(d=60\) and Fig. 5 corresponds to \(d=300\). It is noteworthy that within each iteration of the SLDBO algorithm, the iteration requires one computation of the lower-level function’s gradient, two computations of the upper-level function’s gradient, one computation of the lower-level function’s Jacobian matrix, and one computation of the lower-level function’s Hessian matrix. In Fig. 4 and Fig. 5, the horizontal axis represents the number of lower-level function gradient computations in one node, based on which we present comparative curve results. It can be seen from both figures that  with \(p=0.3\) performs the best, followed by  with \(p=0.75\) and SLDBO.

Fig. 6 showcases the performance of across four graphs, which is showed in Fig. 7. It can be observed that the performance of improves as the connectivity of the graph increases. Specifically, performs best when using the complete graph, followed by the random, the grid and the line graphs. This indicates that a network with better connectivity can generally accelerate the convergence process.

4.1.2 Real-world data↩︎

Similar to the case of synthetic data, we define \(\mathcal{D}_i\) and \(\mathcal{D}_i^{\prime}\) as the training and testing datasets for node \(i\), respectively. Next, we apply , SLDBO and MA-DSBO to solve hyperparameter optimization problem using the MNIST dataset [38]. In particular, the objective and constraint functions in 1 are given by \[\begin{align} F_i(\lambda, \omega)&= \frac{1}{|\mathcal{D}_i^{\prime}|}\sum_{(x_e,y_e)\in\mathcal{D}_i^{\prime}} L(x_e^{\top}\omega, y_e), \\ f_i(\lambda, \omega)&= \frac{1}{|\mathcal{D}_i|}\sum_{(x_e,y_e)\in \mathcal{D}_i}L(x_e^{\top}\omega, y_e) + \frac{1}{cd}\sum_{{m=1}}^{c}\sum_{j=1}^{d}e^{\lambda_j}\omega^2_{mj}, \end{align}\] where \(\omega \in \mathbb{R}^{c\times d}\) is the model parameter, \(L\) represents the cross entropy loss, and \(|S|\) denotes the cardinality of a set \(S\). In this experiment, \(c\) represents the number of classes and \(d\) represents the number of features. The values of \(c\) and \(d\) are set to \(10\) and \(784\) respectively. The training and testing datasets each consist of 60,000 samples, with balanced representation across all classes. For both   and SLDBO, the hyperparameters were set as follows: \(r_v=100\), \({\alpha}={\eta}=6\), and \({\beta}=12\). For MA-DSBO, we set \(T=10\) and \(N=20\). For SLAM, we set \(b=5\). In the experiment, we only tested \(p=0.3\) for . The comparison results of test loss, train loss, and classification accuracy for the three algorithms with respect to CPU time are presented in Fig. 8. These results clearly demonstrate that all the compared algorithms can efficiently solve this problem. Moreover, our proposed algorithm   performs the best in terms of improved convergence rate and classification accuracy, followed by SLDBO, SLAM and MA-DSBO.

4.2 Data hyper-cleaning↩︎

In this subsection, we compare the proposed   with SLDBO on a data hyper-cleaning problem [19] for the Fashion-MNIST dataset [39]. The dataset consists of 60,000 images for training and 10,000 images for testing. Within the training phase, the 60,000 training images are partitioned into two distinct subsets: a training set of 50,000 images (denoted as \(\mathcal{S}_{\mathcal{T}}\)) and a validation set of 10,000 images (denoted as \(\mathcal{S}_{\mathcal{V}}\)). In the data hyper-cleaning problem, we aim to train a classifier in a corrupted setting. Here, each training data’s label is replaced by a random class number with probability prob (i.e. the corruption rate). We again use a decentralized ring network with \(n\) = 8 clients. This problem can be modeled as the bilevel optimization problem 1 with \[\begin{align} F_i(\lambda, w) &= \frac{1}{|\mathcal{S}^{i}_{\mathcal{V}}|}\sum_{(x_e, y_e)\in \mathcal{S}^{i}_{\mathcal{V}}} L(w^{\top}x_e, y_e),\\ f_i(\lambda, w) &= \frac{1}{|\mathcal{S}^{i}_{\mathcal{T}}|}\sum_{(x_e, y_e)\in \mathcal{S}^{i}_{\mathcal{T}}} \sigma(\lambda_e) L(w^{\top}x_e, y_e ) +\tilde{C}\|w\|^2. \end{align}\] Here, \(L\) represents the cross-entropy loss, \(\tilde{C}\) is a regularization parameter, \(\sigma(\cdot)\) is the sigmoid function, \(\mathcal{S}^{i}_{\mathcal{T}}\) and \(\mathcal{S}^{i}_{\mathcal{V}}\) respectively denote the training and validation sets of the \(i\)th client. Following the description in [40], we assume that on each client, training examples are independently drawn with class labels following a categorical distribution over \(M\) classes parameterized by a vector \(\mathbf{q} = (q_1,\ldots,q_M)\) where \(q_i \geq 0\) for all \(i\) and \(\|\mathbf{q}\|_1 = 1\). To synthesize a population of non-identical clients, we draw \(\mathbf{q}\) from a Dirichlet distribution \(\mathbf{q} \sim Dir (\tilde{\alpha} {\mathbf{p}_{cls}})\). Here, \({\mathbf{p}_{cls}}\) represents a prior class distribution over \(M\) classes, and \(\tilde{\alpha} > 0\) is a concentration parameter that controls the identicalness among clients. In the experiment, we set \(\tilde{C} = 0.005\).

In Fig. 9, we display the data class distributions for both the i.i.d. case and the Dirichlet case with \(\tilde{\alpha} = 0.1\), \(\mathbf{p}_{cls} = (1,\ldots,1)\) and \(M\) = 10. For both cases, our rule for corrupting the \(10\) labels is to contaminate them using any one of the other 9 types. Fig. 10 shows the results of upper-level train loss, test accuracy and F1-score of the two algorithms for the two types of data distributions with corruption rate prob\(~= 0.4\), from which it is clear that data heterogeneity does not influence the performance of the two algorithms. This is because neither of the two compared algorithms depends on the assumption of data heterogeneity.

To clearly illustrate the performance of the two compared algorithms, Fig. 11 shows the results of the upper-level train loss, classification accuracy and F1-score with corruption rates prob\(~= 0.4\) and prob\(~= 0.7\) as the CPU time progresses. Here, classification accuracy refers to the accuracy of classifying data whose labels have not been corrupted, and F1-score measures the quality of the data cleaner [41]. The detailed experimental results of the final test accuracy, F1-score, and the consumed CPU time (in seconds) are summarized in Table 1, with an additional corruption rate prob\(~= 0.1\). Compared to SLDBO, our algorithm achieves a 35% reduction in CPU time, enhances accuracy by approximately 4%, and increases the F1-score by 2 to 4 points. From these results, it is evident that outperforms SLDBO in terms of both speed and accuracy.

4.3 Decentralized Meta-Learning↩︎

To further evaluate the applicability of the proposed algorithm to networked AI systems, we conduct experiments on a decentralized meta-learning task [35], [36] using the miniImageNet dataset[42]. Meta-learning aims to learn a shared representation or initialization that can rapidly adapt to new tasks using only a small number of task-specific samples. This setting naturally leads to a bilevel structure, where the lower-level problem performs task-specific adaptation on the support set and the upper-level problem evaluates the adapted model on the query set.

We consider a collection of tasks \(\{\mathcal{T}_r:r=1,\ldots,R\}\). For each task \(\mathcal{T}_r\), let \(x\) denote the shared model parameter and \(y_r\) denote the task-specific parameter. The task loss is defined as \[\ell_r(x,y_r) = \mathbb{E}_{\xi\in\mathcal{D}_r} L(x,y_r;\xi),\] where \(L\) denotes the supervised learning loss. In the decentralized setting, the data associated with each task are distributed over \(n\) agents. For node \(i\) and task \(\mathcal{T}_r\), we denote the local dataset by \(\mathcal{D}_{r,i}\), which is further split into a support set \(\mathcal{D}^{\mathrm{tr}}_{r,i}\) and a query set \(\mathcal{D}^{\mathrm{val}}_{r,i}\).

The decentralized meta-learning problem [43] can be formulated as the DBO problem 1 . For each node \(i\), the upper-level and lower-level objectives are given by \[F_i(x,y) = \frac{1}{R} \sum_{r=1}^{R} \left[ \frac{1}{|\mathcal{D}^{\mathrm{val}}_{r,i}|} \sum_{\xi\in\mathcal{D}^{\mathrm{val}}_{r,i}} L(x,y_r;\xi) \right],\] and \[f_i(x,y) = \frac{1}{R} \sum_{r=1}^{R} \left[ \frac{1}{|\mathcal{D}^{\mathrm{tr}}_{r,i}|} \sum_{\xi\in\mathcal{D}^{\mathrm{tr}}_{r,i}} L(x,y_r;\xi) + \mathcal{R}(y_r) \right],\] where \(y=(y_1,\ldots,y_R)\) collects the task-specific parameters and \(\mathcal{R}(y_r)\) is a regularization term. The lower-level objective adapts the task-specific parameter using the support set, while the upper-level objective measures the query-set performance after adaptation. Compared with the previous logistic-regression-based experiments, this experiment is more representative of networked AI applications because it involves neural-network models, task-level adaptation, and distributed task data.

In our implementation, we use a 5-way 10-shot few-shot classification setting. The dataset is miniImageNet [42], which is generated from ImageNet [44] and consists of 100 classes with each class containing 600 images of size 84 × 84. Each task contains 5 classes, with 10 support samples and 10 query samples per class. The decentralized network consists of \(n=8\) nodes. We use 1000 training tasks, 200 validation tasks, and 600 testing tasks. The batch size is set to 32, and all algorithms are trained for 1500 iterations. The models are evaluated every 25 iterations, and each evaluation is averaged over 200 episodes. The same task splits, batch size, and evaluation protocol are used for all compared methods.

We compare S\(^3\)LDBO with SLDBO, MAML[35] and ANIL[36]. MAML[35] is a representative gradient-based meta-learning baseline, while ANIL[36] restricts task-specific adaptation mainly to the final layers. For a fair comparison, the stepsizes of all compared methods are tuned from the same candidate set on the validation tasks. For S\(^3\)LDBO, we set the upper-level stepsize to \(\alpha = 0.002\), the lower-level stepsize to \(\beta = 0.004\), and the auxiliary-variable stepsize to \(\eta = 0.001\). For SLDBO, we use the same stepsizes. For MAML and ANIL, the meta stepsize and the inner-loop adaptation stepsize are set to \(\alpha_{\mathrm{meta}} = 0.002\) and \(\alpha_{\mathrm{in}} = 0.002\), respectively. These values are selected to provide stable training behavior for all methods under the same miniImageNet task split. For the snapshot probability in S\(^3\)LDBO, we use \(p=0.6\) in this experiment. This value is chosen as a moderate snapshot frequency that balances computational saving and update accuracy: a smaller \(p\) further reduces derivative-evaluation cost but may introduce larger snapshot errors, whereas a larger \(p\) makes the method closer to SLDBO and weakens the computational advantage. Since the influence of \(p\) has already been investigated in the synthetic and data hyper-cleaning experiments, the decentralized meta-learning experiment focuses on validating whether the proposed snapshot mechanism remains effective in a more computationally demanding neural-network-based setting. Therefore, we report the representative result with \(p=0.6\).

Fig. 12 and Fig. 13 compare different algorithms on the meta-learning problem in terms of training and testing accuracy, where Fig. 12 reports the results with respect to the number of iterations and Fig. 13 reports the results with respect to CPU time. Under the iteration-based comparison, S\(^3\)LDBO achieves the highest training accuracy and obtains testing accuracy comparable to MAML, while outperforming SLDBO and ANIL. This indicates that the proposed method maintains strong task adaptation performance during the iterative optimization process. Under the CPU-time-based comparison, S\(^3\)LDBO achieves the best performance on both training and testing tasks within the same time budget. These results show that the snapshot mechanism not only preserves competitive learning performance in terms of iterations, but also improves practical time efficiency by reducing expensive derivative evaluations. Overall, the meta-learning experiment demonstrates that S\(^3\)LDBO is effective for decentralized task adaptation and is suitable for networked AI applications with limited computational resources.

5 Concluding remarks↩︎

In this paper, we propose a novel single-loop algorithm, , which integrates the SLDBO framework with the snapshot gradient tracking technique to efficiently solve decentralized bilevel optimization problems. Unlike other single-loop algorithms that require the computation of full gradients, Jacobian, and Hessian matrices in every iteration – a time-consuming task in practical settings,   allows agents to skip these calculations in some iterations, which significantly reduces computational costs in large-scale distributed optimization. We establish the ergodic iteration complexity and the high probability nonergodic iteration complexity of . Numerical experiments on hyperparameter optimization, data hyper-cleaning, and decentralized meta-learning confirm the efficacy of the proposed algorithm. These results suggest that snapshot-based DBO algorithm is a promising tool for scalable and adaptive learning in networked AI systems. A key avenue for future research is the extension of   to the stochastic setting. Developing an algorithm that can effectively handle variance from both stochastic data sampling and the snapshot mechanism would require new analytical tools and could significantly broaden the practical applicability of this work to large-scale machine learning problems.

6 Complete proof of the convergence results↩︎

In this appendix, we provide the proof of our convergence results, i.e., Theorem 1 and Theorem 2. Assumptions 1 and 2 are used throughout the proof.

Recall that \(r_v:= L_{F,0}/\sigma\) and \(\rho\) is defined in Assumption 2. Moreover, we define the following constants for the analysis: \[\label{def-constants-main-text} \left\{ \begin{array}{l} L_v:=\left(L_{F,1}+L_{f,2}r_v\right)\left(1+ L_{f,1}/\sigma \right), \quad L_1:=L_{F,1}+L_{f,2}r_v, \medskip \\ L_{\Phi}:= L_{F, 1}+\frac{2 L_{F, 1} L_{f, 1}+L_{f, 2} L_{F, 0}^2}{\sigma}+\frac{2 L_{f, 1} L_{F, 0} L_{f, 2}+L_{f, 1}^2 L_{F, 1}}{\sigma^2}+\frac{L_{f, 2} L_{f, 1}^2 L_{F, 0}}{\sigma^3}, \end{array} \right.\tag{6}\] where constants such as \(\sigma\), \(L_{F,0}\), \(L_{F,1}\), \(L_{f,1}\), \(L_{f,2}\) are defined in Assumptions 1 and 2.

6.0.1 Notation, constants, and roadmap↩︎

For convenience in the subsequent convergence proof, we rewrite the \(k\)th iteration of Algorithm 1 in a more concise form as follows, together with [yd], [vd] and [xd]: \[\tag{7} \begin{align} z_{y,i}^{k} &= t_{y,i}^k + \xi^k(d_{y,i}^k - d_{\tilde{y},i}^k)/p,& y_i^{k+1}&=\sum\nolimits_{j=1}^n w_{i j} (y_j^k-\beta z_{y, j}^k), \tag{8} \\ \tilde{y}_i^{k+1} &= \xi^k y_i^{k} + (1-\xi^k)\tilde{y}_i^{k},& t_{y, i}^{k+1}&=\sum\nolimits_{j=1}^n w_{i j} \big(t_{y, j}^{k}+ \xi^k (d_{y,j}^{k}-d_{\tilde{y}, j}^{k})\big), \tag{9} \\ z_{v,i}^{k} &= t_{v,i}^k + \xi^k (d_{v,i}^k - d_{\tilde{v},i}^k)/p,& v_i^{k+1}&=\mathcal{P}_{r_v}\left[\sum\nolimits_{j=1}^n w_{i j}( v_j^k+\eta z_{v, j}^k)\right], \tag{10} \\ \tilde{v}_i^{k+1} &= \xi^k v_i^{k} + (1-\xi^k)\tilde{v}_i^{k},& t_{v, i}^{k+1}&=\sum\nolimits_{j=1}^n w_{i j} \big(t_{v, j}^{k}+\xi^k (d_{v,j}^{k}-d_{\tilde{v}, j}^{k})\big), \tag{11} \\ z_{x,i}^{k} &= t_{x,i}^k + \xi^k (d_{x,i}^k - d_{\tilde{x},i}^k)/p,& x_i^{k+1}&=\sum\nolimits_{j=1}^n w_{i j} (x_j^k-\alpha z_{x,j}^k),\tag{12} \\ \tilde{x}_i^{k+1} &= \xi^k x_i^{k} + (1-\xi^k)\tilde{x}_i^{k},& t_{x, i}^{k+1}&=\sum\nolimits_{j=1}^n w_{i j}\big( t_{x, j}^{k}+\xi^k (d_{x,j}^k-d_{\tilde{x}, j}^{k})\big),\tag{13}\\ d_{\tilde{y},i}^{k+1} & = \xi^{k}d_{y,i}^{k} + (1-\xi^{k})d_{\tilde{y},i}^{k} ,& d_{\tilde{v},i}^{k+1} & = \xi^{k}d_{v,i}^{k} + (1-\xi^{k})d_{\tilde{v},i}^{k} ,\\ d_{\tilde{x},i}^{k+1} & = \xi^{k}d_{x,i}^{k} + (1-\xi^{k})d_{\tilde{x},i}^{k}. \end{align}\] From 7 , we know that \(d_{\tilde{y},i}^{k+1}, d_{\tilde{v},i}^{k+1}\) and \(d_{\tilde{x},i}^{k+1}\) have the following representations, which will be useful in the subsequent analysis: \[\begin{align} d_{\tilde{y},i}^{k+1} &= \nabla_2 f_i(\tilde{x}_i^{k+1},\tilde{y}_i^{k+1}), \tag{14}\\ d_{\tilde{v},i}^{k+1} &= \nabla_2 F_i(\tilde{x}_i^{k+1},\tilde{y}_i^{k+1}) - \nabla^2_{22} f_i(\tilde{x}_i^{k+1},\tilde{y}_i^{k+1}) \tilde{v}_i^{k+1},\tag{15}\\ d_{\tilde{x},i}^{k+1} &= \nabla_1 F_i(\tilde{x}_i^{k+1},\tilde{y}_i^{k+1}) - \nabla^2_{12} f_i(\tilde{x}_i^{k+1},\tilde{y}_i^{k+1}) \tilde{v}_i^{k+1}.\tag{16} \end{align}\] For the analysis, we further define \[\label{def:bar9} \left\{ \begin{array}{lll} {\bar{x}^{k}} := \frac{1}{n}\sum_{i=1}^{n}x_{i}^k, & {\bar{y}^{k}} := \frac{1}{n}\sum_{i=1}^{n}y_{i}^k, & {\bar{v}^{k}} := \frac{1}{n}\sum_{i=1}^{n}v_{i}^k, \smallskip\\ \bar{d}_x^k := \frac{1}{n}\sum_{i=1}^{n}d_{x,i}^k, & \bar{d}_y^k := \frac{1}{n}\sum_{i=1}^{n}d_{y,i}^k, & \bar{d}_v^k := \frac{1}{n}\sum_{i=1}^{n}d_{v,i}^k, \smallskip\\ \bar{d}_{\tilde{x}}^k := \frac{1}{n}\sum_{i=1}^{n}d_{\tilde{x},i}^k, & {\bar{d}_{\tilde{y}}^{k}} := \frac{1}{n}\sum_{i=1}^{n}d_{\tilde{y},i}^k, & \bar{d}_{\tilde{v}}^k := \frac{1}{n}\sum_{i=1}^{n}d_{\tilde{v},i}^k, \smallskip \\ \bar{t}_x^k := \frac{1}{n}\sum_{i=1}^{n}t_{x,i}^k, & \bar{t}_y^k := \frac{1}{n}\sum_{i=1}^{n}t_{y,i}^k, & \bar{t}_v^k := \frac{1}{n}\sum_{i=1}^{n}t_{v,i}^k, \smallskip \\ \bar{z}_x^k := \frac{1}{n}\sum_{i=1}^{n}z_{x,i}^k, & \bar{z}_y^k := \frac{1}{n}\sum_{i=1}^{n}z_{y,i}^k, & \bar{z}_v^k := \frac{1}{n}\sum_{i=1}^{n}z_{v,i}^k, \end{array} \right.\tag{17}\] and two additional constants \[\label{def-constants-proof} C_1 := 3\max\{(L_{F,1}+r_{v}L_{f,2})^2,L_{f,1}^2\} \text{~~and~~} C_2 := L_{f,1}^2.\tag{18}\]

6.0.1.1 Roadmap of the proof.

Here we briefly describe the roadmap of the proof of Theorem 1. First, we define the Lyapunov function \[\begin{align} V_k = \; & F(\bar{x}^k, y ^{*}(\bar{x}^k)) - F^{*} + a_1 \|\bar{y}^k - y^{*}(\bar{x}^k)\|^2 + a_2 \|\bar{v}^k - v^{*}(\bar{x}^k)\|^2 + \frac{a_3}{n}\sum_{i=1}^n \|x_i^k -\bar{x}^k\|^2 \nonumber\\ & + \frac{a_4}{n} \sum_{i=1}^n \|y_i^k -\bar{y}^k\|^2 +\frac{a_5}{n} \sum_{i=1}^n \|v_i^k - \bar{v}^k\|^2 + \frac{a_6 \alpha^2}{n} \sum_{i=1}^n\|z_{x,i}^k - \bar{z}_{x}^k\|^2 + \frac{a_7 \beta^2}{n}\sum_{i=1}^n \|z_{y,i}^k - \bar{z}_{y}^k\|^2 \nonumber \\ & + \frac{a_8 \eta^2}{n} \sum_{i=1}^n \|z_{v,i}^k - \bar{z}_{v}^k\|^2 + \frac{a_9}{n} \sum_{i=1}^n \|x_i^k - \tilde{x}_i^k\|^2 + \frac{a_{10}}{n} \sum_{i=1}^n \|y_i^k - \tilde{y}_i^k\|^2 + \frac{a_{11}}{n} \sum_{i=1}^n \|v_i^k -\tilde{v}_i^k\|^2,\label{Lyapunov} \end{align}\tag{19}\] where \(v^*(x) := \left[\nabla^2_{22} f(x,y^*(x))\right]^{-1}\nabla_2 F(x,y^*(x))\) and \(a_1, a_2, \ldots, a_{11}\) are constants. We will show that, for another set of positive constants \(A_1, A_2, \ldots, A_{12}\), the following inequality holds \[\begin{align} \mathbb{E}\Big[V_{k+1} - V_{k}\Big] \leq & -\frac{\alpha}{2} \mathbb{E}\Big[\|\nabla \Phi(\bar{x}_{k})\|^2\Big] -A_1\mathbb{E}\Big[\|\bar{d}_{x}^k\|^2\Big] -A_2\mathbb{E}\Big[\|\bar{y}^k - y^{*}(\bar{x}^k)\|^2\Big] - A_3 \mathbb{E}\Big[\|\bar{v}^k - v^{*}(\bar{x}^k)\|^2\Big] \\ & - \frac{A_4}{n}\mathbb{E}\Big[\sum_{i=1}^n\|x_i^k -\bar{x}^k\|^2\Big] - \frac{A_5}{n}\mathbb{E}\Big[\sum_{i=1}^n\|y_i^k -\bar{y}^k\|^2\Big] - \frac{A_6}{n}\mathbb{E}\Big[\sum_{i=1}^n\|v_i^k -\bar{v}^k\|^2\Big]\\ & - \frac{A_7}{n}\mathbb{E}\Big[\sum_{i=1}^n\|z_{x,i}^k - \bar{z}_{x}^k\|^2\Big] - \frac{A_8}{n}\mathbb{E}\Big[\sum_{i=1}^n\|z_{y,i}^k - \bar{z}_{y}^k\|^2\Big] - \frac{A_9}{n}\mathbb{E}\Big[\sum_{i=1}^n\|z_{v,i}^k - \bar{z}_{v}^k\|^2\Big]\\ & - \frac{A_{10}}{n}\mathbb{E}\Big[\sum_{i=1}^n\|x_i^k -\tilde{x}_{i}^k\|^2\Big] - \frac{A_{11}}{n}\mathbb{E}\Big[\sum_{i=1}^n\|y_i^k -\tilde{y}_{i}^k\|^2\Big] - \frac{A_{12}}{n} \mathbb{E}\Big[\sum_{i=1}^n \|v_i^k -\tilde{v}_{i}^k\|^2\Big]. \end{align}\] Then, the proof can be obtained by telescoping the above inequality.

6.0.2 Some lemmas for consensus error of Algorithm 1↩︎

Recall that \({\bar{x}^{k}}\), \({\bar{y}^{k}}\), \({\bar{v}^{k}}\), \(\bar{d}_x^k\), \(\bar{d}_y^k\), \(\bar{d}_v^k\), \(\bar{d}_{\tilde{x}}^k\), \(\bar{d}_{\tilde{y}}^k\), \(\bar{d}_{\tilde{v}}^k\), \(\bar{t}_x^k\), \(\bar{t}_y^k\), \(\bar{t}_v^k\), \(\bar{z}_x^k\), \(\bar{z}_y^k\) and \(\bar{z}_v^k\) are defined in 17 . Let \(\mathcal{F}_k\) be the sigma field generated by \(\{\xi^j\}_{0\leq j\leq k-1}\) and \(\mathbb{E}_k[\cdot]\) be the conditional expectation taken over \(\mathcal{F}_k\). Then, for any \(k \geq 0\), the variables \(\{x_i^k, y_i^k, v_i^k, d_{x,i}^k, d_{y,i}^k, d_{v,i}^k, \tilde{x}_i^k, \tilde{y}_i^k, \tilde{v}_i^k, d_{\tilde{x},i}^k, d_{\tilde{y},i}^k, d_{\tilde{v},i}^k, t_{x,i}^k, t_{y,i}^k, t_{v,i}^k\}\) are measurable with respect to \(\mathcal{F}_k\). We first prove some useful lemmas.

Then, we will bound the terms related to the consensus error in 19 .

Next, we bound the term \(\mathbb{E}_k\Big[\sum_{i=1}^n \|d_{x,i}^k-d_{\tilde{x},i}^k\|^2\Big]\) in 22 , as well as \(\mathbb{E}_k\Big[\sum_{i=1}^n \|d_{y,i}^k-d_{\tilde{y},i}^k\|^2\Big]\) and \(\mathbb{E}_k\Big[\sum_{i=1}^n \|d_{v,i}^k-d_{\tilde{v},i}^k\|^2\Big]\).

Next, we bound the term \(\mathbb{E}_k\Big[\sum_{i=1}^n \|v_i^{k+1} - v_i^k\|^2\Big]\) in 27 .

Next, we bound \(\mathbb{E}_k\Big[\sum_{i=1}^n (\|d_{\star,i}^k - d_{\tilde{\star},i}^k\|^2 + \|d_{\star,i}^{k+1} - d_{\tilde{\star},i}^{k+1}\|^2)\Big]\) for \(\star = x, y, v\) respectively. By combining 23 28 , we derive \[\begin{align} &\mathbb{E}_k\Big[\sum_{i=1}^n \big(\|d_{x,i}^k - d_{\tilde{x},i}^k\|^2 + \|d_{x,i}^{k+1} - d_{\tilde{x},i}^{k+1}\|^2\big)\Big] \nonumber \\ \leq &C_1 \mathbb{E}_k\Big[\sum_{i=1}^n \big(\|x_i^{k+1} - \tilde{x}_{i}^{k+1}\|^2 + \|y_i^{k+1} - \tilde{y}_{i}^{k+1}\|^2 + \|v_i^{k+1} - \tilde{v}_{i}^{k+1}\|^2 + \|x_i^{k} - \tilde{x}_{i}^{k}\|^2 + \|y_i^{k} - \tilde{y}_{i}^{k}\|^2 + \|v_i^{k} - \tilde{v}_{i}^{k}\|^2\big)\Big] \nonumber \\ \leq & C_1\sum_{i=1}^n\Big( \frac{1-p+\rho}{\rho}\big(\|x_i^k -\tilde{x}_{i}^k\|^2 + \|y_i^k -\tilde{y}_{i}^k\|^2 +\|v_i^k -\tilde{v}_{i}^k\|^2\big) \nonumber\\ & \quad + \frac{1}{1-\rho}\mathbb{E}_k\Big[\sum_{i=1}^n \|x_i^{k+1} - x_i^k\|^2+\sum_{i=1}^n \|y_i^{k+1} - y_i^k\|^2+\sum_{i=1}^n \|v_i^{k+1} - v_i^k\|^2\Big]\Big). \label{dxk43dxk431} \end{align}\tag{31}\] Similarly, we can deduce \[\begin{align} &\mathbb{E}_k\Big[\sum_{i=1}^n \big(\|d_{y,i}^k - d_{\tilde{y},i}^k\|^2 + \|d_{y,i}^{k+1} - d_{\tilde{y},i}^{k+1}\|^2\big)\Big] \nonumber \\ \leq \, & C_2{\sum_{i=1}^n\Big( \frac{1-p+\rho}{\rho}\big(\|x_i^k -\tilde{x}_{i}^k\|^2 + \|y_i^k -\tilde{y}_{i}^k\|^2\big) + \frac{1}{1-\rho}\mathbb{E}_k\Big[\sum_{i=1}^n \|x_i^{k+1} - x_i^k\|^2+\sum_{i=1}^n \|y_i^{k+1} - y_i^k\|^2\Big] \Big) }, \tag{32}\\ \text{and} \nonumber \\ &\mathbb{E}_k\Big[\sum_{i=1}^n \big(\|d_{v,i}^k - d_{\tilde{v},i}^k\|^2 + \|d_{v,i}^{k+1} - d_{\tilde{v},i}^{k+1}\|^2\big)\Big] \nonumber \\ \leq \, & C_1\sum_{i=1}^n\Big( \frac{1-p+\rho}{\rho}\big(\|x_i^k -\tilde{x}_{i}^k\|^2 + \|y_i^k -\tilde{y}_{i}^k\|^2 + \|v_i^k -\tilde{v}_{i}^k\|^2\big) \nonumber\\ & \quad + \frac{1}{1-\rho}\mathbb{E}_k\Big[\sum_{i=1}^n \|x_i^{k+1} - x_i^k\|^2+\sum_{i=1}^n \|y_i^{k+1} - y_i^k\|^2+\sum_{i=1}^n \|v_i^{k+1} - v_i^k\|^2\Big]\Big). \tag{33} \end{align}\] Next, we recall a Lemma from [28], which is useful for bounding \(\|\bar{d}_{y}^k\|^2\) and \(\|\bar{d}_{v}^k\|^2\).

Now, by combining 22 , the inequalities in Remark 3, Lemma 8, Remark 6 and 31 35 , we can establish the following results for \(\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{\star,i}^{k+1} - \bar{z}_{\star}^{k+1}\|^2\Big]\) for \(\star = x, y, v\) respectively. Again, details are omitted for simplicity. \[\begin{align} &\mathbb{E}_k\Big[ \frac{1}{n}\sum_{i=1}^n\|z_{x,i}^{k+1} - \bar{z}_{x}^{k+1}\|^2\Big] \nonumber\\ \leq & \Big(\rho + \frac{4\alpha^2 C}{(1-\rho)^3 \rho p}\Big)\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{x,i}^{k} - \bar{z}_{x}^{k}\|^2\Big] + \frac{4\beta^2 C}{(1-\rho)^3\rho p}\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{y,i}^{k} - \bar{z}_{y}^{k}\|^2\Big] + \frac{4\eta^2 C}{(1-\rho)^3\rho p}\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{v,i}^{k} - \bar{z}_{v}^{k}\|^2\Big] \nonumber\\ & + \Big(\frac{(1-p+\rho)C}{(1-\rho)^2 \rho^2 p} + \frac{4(\alpha^2+\beta^2+\eta^2)C^2}{(1-\rho^3)p^2}\Big){\frac{1}{n}\sum_{i=1}^n\big(\|x_i^k - \tilde{x}_i^k\|^2 + \|y_i^k - \tilde{y}_i^k\|^2 + \|v_i^k - \tilde{v}_i^k\|^2\big)} + \frac{4\alpha^2 C}{(1-\rho)^3 \rho p}{\|\bar{d}_{x}^k\|^2} \nonumber\\ & + \frac{(8+8L_{f,1}^2\beta^2 + 20\eta^2 L_1^2)C}{(1-\rho)^3 \rho p} \Big({\frac{1}{n}\sum_{i=1}^n\|x_i^k - \bar{x}^k\|^2} + {\frac{1}{n}\sum_{i=1}^n\|y_i^k - \bar{y}^k\|^2}\Big) + \frac{8C}{(1-\rho)^3 \rho p}{\frac{1}{n}\sum_{i=1}^n\|v_i^k - \bar{v}^k\|^2} \nonumber\\ & + \frac{(8L_{f,1}^2\beta^2 + 20\eta^2 L_1^2 )C}{(1-\rho)^3 \rho p} {\|\bar{y}^k - y^*(\bar{x}^k)\|^2} + \frac{20\eta^2 L_1^2 C}{(1-\rho)^3 \rho p} {\|\bar{v}^k - v^*(\bar{x}^k)\|^2}, \label{bound95for95zx-barz} \end{align}\tag{36}\] \[\begin{align} &\mathbb{E}_k\Big[ \frac{1}{n}\sum_{i=1}^n\|z_{y,i}^{k+1} - \bar{z}_{y}^{k+1}\|^2\Big] \nonumber\\ \leq & \Big(\rho + \frac{4\beta^2 C}{(1-\rho)^3 \rho p}\Big)\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{y,i}^{k} - \bar{z}_{y}^{k}\|^2\Big] + \frac{4\alpha^2 C}{(1-\rho)^3\rho p}\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{x,i}^{k} - \bar{z}_{x}^{k}\|^2\Big] \nonumber\\ & + \Big(\frac{(1-p+\rho)C}{(1-\rho)^2 \rho^2 p} + \frac{4(\alpha^2+\beta^2+\eta^2)C^2}{(1-\rho^3)p^2}\Big){\frac{1}{n}\sum_{i=1}^n\big(\|x_i^k - \tilde{x}_i^k\|^2 + \|y_i^k - \tilde{y}_i^k\|^2\big)} + \frac{8L_{f,1}^2\beta^2 C}{(1-\rho)^3 \rho p} {\|\bar{y}^k - y^*(\bar{x}^k)\|^2} \nonumber\\ & + \frac{(8+8L_{f,1}^2\beta^2)C}{(1-\rho)^3 \rho p} \Big({\frac{1}{n}\sum_{i=1}^n\|x_i^k - \bar{x}^k\|^2} + {\frac{1}{n}\sum_{i=1}^n\|y_i^k - \bar{y}^k\|^2}\Big) + \frac{4\alpha^2 C}{(1-\rho)^3 \rho p}{\|\bar{d}_{x}^k\|^2}, \tag{37}\\ &\mathbb{E}_k\Big[ \frac{1}{n}\sum_{i=1}^n\|z_{v,i}^{k+1} - \bar{z}_{v}^{k+1}\|^2\Big] \nonumber\\ \leq & \Big(\rho + \frac{4\eta^2 C}{(1-\rho)^3 \rho p}\Big)\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{v,i}^{k} - \bar{z}_{v}^{k}\|^2\Big] + \frac{4\alpha^2 C}{(1-\rho)^3\rho p}\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{x,i}^{k} - \bar{z}_{x}^{k}\|^2\Big] + \frac{4\beta^2 C}{(1-\rho)^3\rho p}\mathbb{E}_k\Big[\frac{1}{n}\sum_{i=1}^n\|z_{y,i}^{k} - \bar{z}_{y}^{k}\|^2\Big] \nonumber\\ & + \Big(\frac{(1-p+\rho)C}{(1-\rho)^2 \rho^2 p}+ \frac{4(\alpha^2+\beta^2+\eta^2)C^2}{(1-\rho^3)p^2}\Big) {\frac{1}{n}\sum_{i=1}^n\big(\|x_i^k - \tilde{x}_i^k\|^2 + \|y_i^k - \tilde{y}_i^k\|^2 + \|v_i^k - \tilde{v}_i^k\|^2\big)} + \frac{4\alpha^2 C}{(1-\rho)^3 \rho p}{\|\bar{d}_{x}^k\|^2} \nonumber\\ & + \frac{(8+8L_{f,1}^2\beta^2 + 20\eta^2 L_1^2)C}{(1-\rho)^3 \rho p} \Big({\frac{1}{n}\sum_{i=1}^n\|x_i^k - \bar{x}^k\|^2} + {\frac{1}{n}\sum_{i=1}^n\|y_i^k - \bar{y}^k\|^2}\Big) + \frac{8C}{(1-\rho)^3 \rho p}{\frac{1}{n}\sum_{i=1}^n\|v_i^k - \bar{v}^k\|^2} \nonumber\\ & + \frac{(8L_{f,1}^2\beta^2 + 20\eta^2 L_1^2 )C}{(1-\rho)^3 \rho p} {\|\bar{y}^k - y^*(\bar{x}^k)\|^2} + \frac{20\eta^2 L_1^2 C}{(1-\rho)^3 \rho p} {\|\bar{v}^k - v^*(\bar{x}^k)\|^2},\tag{38} \end{align}\] where \(C = \max \{C_1,C_2\}\) and \(L_1\) is defined in 6 .

6.0.3 Some lemmas for convergence rate of Algorithm 1↩︎

Recall that \(\Phi(x) = F(x,y^*(x))\) denotes the overall objective function. For any \(k>0\), \(\Phi(\bar{x}^k)\) and \(\nabla \Phi(\bar{x}^k)\) are also measurable with respect to \(\mathcal{F}_k\).

The next two lemmas bound \(\mathbb{E}_k\Big[\Vert\bar{y}^{k+1}-y^*(\bar{x}^k)\Vert^2\Big]\) and \(\mathbb{E}_k\Big[\Vert \bar{v}^{k+1} - v^* (\bar{x}^k)\Vert^2\Big]\), respectively.

6.0.4 Proof of Theorem 1.↩︎

6.0.5 Proof of Theorem 2.↩︎

References↩︎

---

1. Chao Yin is with the School of Mathematics, Hohai University, Nanjing, P. R. China (e-mail: yinchao@hhu.edu.cn).↩︎

2. Youran Dong, Bofan Wang, and Junfeng Yang are with the School of Mathematics, Nanjing University, Nanjing, P. R. China (e-mail: yrdong@smail.nju.edu.cn; wangbf@smail.nju.edu.cn; jfyang@nju.edu.cn).↩︎

3. Shiqian Ma is with the Department of Computational Applied Mathematics and Operations Research, Rice University, Houston, USA (e-mail: sqma@rice.edu).↩︎