Regret-Optimal Federated Transfer Learning for Kernel Regression – with Applications in American Option Pricing

Xuwei Yang<EMAIL_ADDRESS>
Anastasis Kratsios<EMAIL_ADDRESS>
Florian Krach<EMAIL_ADDRESS>
Matheus Grasselli<EMAIL_ADDRESS>
Aurelien Lucchi<EMAIL_ADDRESS>
We propose an optimal iterative scheme for federated transfer learning, where a central planner has access to datasets ${\cal D}_1,\dots,{\cal D}_N$ for the same learning model $f_{\theta}$.  Our objective is to minimize the cumulative deviation of the generated parameters $\{\theta_i(t)\}_{t=0}^T$ across all $T$ iterations from the specialized parameters
$\theta^\star_{1},\ldots,\theta^\star_N$ obtained for each dataset, while respecting the loss function for the model $f_{\theta(T)}$ produced by the algorithm upon halting.  We only allow for continual communication between each of the specialized models (nodes/agents) and the central planner (server), at each iteration (round).
For the case where the model $f_{\theta}$ is a finite-rank kernel regression, we derive explicit updates for the regret-optimal algorithm.  By leveraging symmetries within the regret-optimal algorithm, we further develop a nearly regret-optimal heuristic that runs with $\mathcal{O}(Np^2)$ fewer elementary operations, where $p$ is the dimension of the parameter space.
Additionally, we investigate the adversarial robustness of the regret-optimal algorithm showing that an adversary which perturbs $q$ training pairs by at-most $\varepsilon>0$, across all training sets, cannot reduce the regret-optimal algorithm's regret by more than $\mathcal{O}(\varepsilon q \bar{N}^{1/2})$, where $\bar{N}$ is the aggregate number of training pairs. To validate our theoretical findings, we conduct numerical experiments in the context of American option pricing, utilizing a randomly generated finite-rank kernel.

§ INTRODUCTION

In classical regression analysis, one typically selects a regression model $f:\Theta \times \mathbb{R}^d\rightarrow \mathbb{R}$ by optimizing a performance criterion on a (finite) training dataset $\mathcal{D}_1\subset \mathbb{R}^d\times \mathbb{R}$. The model $f$ is parameterized by $\Theta\subseteq \mathbb{R}^p$, with $p,d\in \mathbb{N}_+$, and for fixed $\theta \in \Theta$ we write $f_\theta \eqdef f(\theta, \cdot)$ for the function that maps the input $x\in \mathbb{R}^d$ to the model's output $y\in \mathbb{R}$.
Selecting a regression model then amounts to identifying a parameter vector $\theta_1 \in \Theta$ which is typically chosen to minimize the penalized sum-of-squares errors (SSE) on the given dataset

\begin{equation}
\label{eq:Regression_Klassische0}
\inf_{\theta \in \Theta}\,
\sum_{(x,y)\in \mathcal{D}_1}\,
\kappa \,
\|\theta\|_2.
\end{equation}

Here, $\kappa \geq 0$ is a hyperparameter controlling the regularity of the model, which is often implicitly reflected by the Euclidean norm of its parameter vector $\|\theta\|_2$, though other regularizers are available, for example in <cit.>.  Even if $y$ is obscured by “noise”, there is an underlying homogeneity assumption underpinning (<ref>), namely that the training data all belongs to the same dataset and the elements of that dataset represent the same task.  Typically (see e.g. <cit.>), one assumes that there is a suitable, yet unknown, “target” function $f^\star:\mathbb{R}^d\rightarrow \mathbb{R}$ and a “noise” generating probability distribution $\mu$ on $\mathbb{R}$ such that for every $(x,y)$ in the training dataset $\mathcal{D}_1$ one has $y=f^{\star}(x)+\varepsilon_y$ where $\{\varepsilon_y\}_{(x,y)\in \mathcal{D}_1}$ are i.i.d. random variables with law $\mu$.

A parameter vector optimizing (<ref>) need not be available in closed-form.  Instead, provided it exists, it is typically identified approximately via an iterative gradient descent-type procedure.  Even if a closed-form expression exists, it is often more computationally efficient to resort to iterative algorithms <cit.>.  A deterministic iterative algorithm initialized at $\theta_0\in\mathbb{R}^p$ and running for $T\in\mathbb{N}_+$ iterations is characterized by its iterates.  Thus, any such iterative deterministic algorithm is identifiable with a sequence in $\mathbb{R}^p$ of length $1+T$; i.e.

\begin{equation}
\label{eq:algorithms}
\theta^{0\cdots T}\eqdef \left\{\theta(t)\right\}_{t=0}^T,
\qquad
\theta(0)=\theta_0.
\end{equation}

With this identification and a measure of algorithmic efficiency, e.g. regret, one can characterize the most efficient algorithm as in <cit.>.

The situation becomes more complicated in modern machine learning, where one often has access to large amounts of data, typically drawn from several similar, yet distinct sources.
These types of comparable datasets arise in meta-learning <cit.>), in multi-task learning <cit.>, in quantitative finance due to natural market cycles <cit.>, or in robust finance approaches to model uncertainty resulting in several plausible stochastic time-series models wherewith data can be simulated <cit.>.

This paper focuses on the optimization aspect of the regression problem under multi-sourced data, with a focus on applications in mathematical finance, instead of the more often explored questions of learning-theoretic nature such as the risk of a federated learning procedure as in <cit.>.  Specifically, we construct a regret-optimal federated optimization algorithm, as defined shortly, under no communication restrictions between the server (central planner) and the local nodes (agents) at each round.

Optimizing with Multi-Sourced Data:
In situations where we have access to several distinct additional datasets $\mathcal{D}_2,\dots,\mathcal{D}_N$ whose information we wish to leverage, we can improve the predictive performance of our regression model past what can be achieved using the focal or primary dataset $\mathcal{D}_1$ in the classical approach (<ref>).
Here, the underpinning assumption is that the datasets $\mathcal{D}_i$, for $i\neq 1$, are generated by a different, but comparable, mechanism from which $\mathcal{D}_1$ was generated.

The introduction of multi-sourced data yields a spectrum of interpretations of optimality extending (<ref>), which can be naturally organized from a game-theoretic point of view.  In this perspective, we can consider several agents (or nodes), each representing the influence of a dataset $\mathcal{D}_i$ on the model selection problem, who propose parameter vectors $\theta_1,\dots,\theta_N$ optimized for their own datasets. The agents are then  assumed to be organized by a central planner (or server), representing the user. The objective of this central planner is to maximize the predictive power of the selected model on new data sampled from the same process which generated the focal dataset $\mathcal{D}_1$.

The Central Planner:
The purpose of the central planner is to avoid the situation whereby the structure of the different datasets is ignored in the model selection problem, by merging them into a single dataset $\mathcal{D}\eqdef \cup_{i=1}^N\,\mathcal{D}_i$ and thereby reducing the problem to (<ref>). Though intuitively simple, such merging  approach can be particularly disadvantageous in the presence of heterogeneity among the datasets, since a dataset that is potentially very different from the focal dataset $\mathcal{D}_1$ would be equally influential in the regression problem (<ref>).
Instead, the user, acting as a central planner, organizes the model selection problem through a cooperative objective

\begin{align}
\label{eq:terminal_time_objective}
\begin{aligned}
l(\Theta, w;\,\mathcal{D} )
\eqdef &
\sum_{i=1}^N
\sum_{j=1}^{|\mathcal{D}_i|}
- y^i_j
\\
\mbox{ with }
\,
\theta^w\eqdef & \sum_{i=1}^N
\end{aligned}
\end{align}

where $\Theta\eqdef (\theta_1^{\top}, \cdots, \theta_N^{\top})^\top\in \mathbb{R}^{Np \times 1}$. Here, the central planner ascribes the influence of each agent on the model selection problem through the weight vector $w=(w_1, \cdots, w_N)\in [0,1]^N$ with $\sum_{i=1}^N\,w_i=1$.
As can be seen in (<ref>), this influence is two-fold: on the tone hand, it impacts the cooperative objective $l$ in the aggregation of the SSE of each player; on the other hand, it affects the proportion of each player's preferred parameter choice entering into the collectively selected parameter $\theta^{w}$.

The choice of $w$ is approached in different ways in the literature, e.g. by bagging

reminiscent of mean-field games <cit.>, by data-driven weighting procedures <cit.>, Bayesian aggregation <cit.>, or with dictionaries <cit.>.

The Agents:
When acting in isolation, each agent's preferred model is determined by optimizing a version of  (<ref>) using only their dataset.  Under mild conditions[E.g. coercivity and lower semi-continuity in $\theta$.] on $f$, this corresponds to an individual parameter selection

\begin{align}
\theta^{\star}_i \in
\operatorname{argmin}_{\theta \in \Theta}
\,
\sum_{(x,y)\in \mathcal{D}_i}\,
\kappa \,
\|\theta\|_2,
\label{eq:Regression_Klassische1}
\end{align}

for $i=1,\dots,N$, which we assume to be fixed and known to all agents prior to their collective parameter selection. Taking into account the other agents, each agent acts by specifying an iterative deterministic algorithm of the form (<ref>), with the intent of maximizing their influence on the joint model selection problem (<ref>).  Thus, the $i^{th}$ player wants the jointly selected model, specified by $\theta^{w}(T)=\sum_{n=1}^N\, w_i\,\theta_i(T)$, to be as close to $\theta_i^{\star}$ as possible, thereby maximally encoding the characteristics of the $i^{th}$ dataset into the selected model $f_{\theta^w(T)}$. The accumulated regret which the $i^{th}$ player incurs by deploying an algorithm $\theta^{0\cdots T}_{i}$ that deviates from their preferred selection $\theta^*_i$ is measured by

\begin{equation}
\label{eq:penalized_regret_functional__REGRETFORM___Individual}
\begin{aligned}
%     \mathcal{R}_i(\theta^{0\cdots T}_{i})
% \eqdef &
\sum_{t=0}^{T-1}
\,
\underbrace{
\lambda \,
\left\| \theta_i(t+1) - \theta^\star_{i} \right\|_2^2
}_{\text{Preference Strength}}
\underbrace{
\beta \, \left\| \Delta \theta_i(t) \right\|_2^2
}_{\text{Algorithm Stability}}
\end{aligned}
\end{equation}

where $\lambda>0$ is a hyperparameter quantifying the player's attachment towards their preferred model choice $\theta^{\star}_i$ and the hyperparameter $\beta > 0$ quantifies the the stability of the algorithm during the iterations, where  $\Delta \theta_i(t)\eqdef \theta_i(t+1)-\theta_i(t)$.

The central planner then organizes the action $\Theta^{0\cdots T}\eqdef
\left\{\Theta(t)\right\}_{t=0}^{T}$ of the system of $N$ players to reach the objective (<ref>) while encoding their individual regrets (<ref>), where

\begin{equation}
\Theta (t) = (\theta_1(t)^\top,\ldots,\theta_N(t)^\top)^\top \in \mathbb{R}^{Np \times 1}. \notag
\end{equation}

This is achieved by coupling the individual agent's regret functionals (<ref>) through the systemic regret functional

\begin{equation}
\label{eq:penalized_regret_functional__REGRETFORM}
\begin{aligned}
\mathcal{R}(\Theta^{0\cdots T})
\eqdef &
\underbrace{
l (\Theta(T), w;\,\mathcal{D} )
}_{\text{Cooperative Objective}}
\,\,  l^{\star}
\sum_{t=0}^{T-1}
\,
\big(
\underbrace{
\lambda \left\| \Theta(t+1) - \Theta^\star \right\|_2^2
}_{\text{Preference Strength}}
\underbrace{
\beta  \left\| \Delta \Theta(t) \right\|_2^2
}_{\text{Algo. Stability}}
\big)
\end{aligned}
\end{equation}

\begin{equation}
\label{Theta_star}
\Theta^\star\eqdef (\theta_1^{\star \top}, \cdots, \theta_N^{\star\top})^\top
\in \mathbb{R}^{Np \times 1}
\end{equation}

encodes the individual preferences of the players,

\begin{equation}
\Delta \Theta(t)\eqdef \Theta(t+1)-\Theta(t)
\label{Theta-simple__nocontrol}
\end{equation}

quantifies the integrative updates of any candidate optimizing sequence $\Theta^{0\cdots T}$, and the ideal terminal loss is

\[
\eqdef
\underbrace{
\min_{\Xi \in \mathbb{R}^{Np\times 1}}\,
l(\Xi,w; \mathcal{D})
% + \lambda \|\Xi - \Theta^{\star} \|_2
}_{\text{Cooperative Sub-optimality}}
\]

The term $(l(\cdot,\cdot;\cdot)-l^{\star})$ measures the central planner's regret in failing to optimize (<ref>).  The algorithm stability term in (<ref>) appears in forward-backward proximal splitting algorithms with quadratic objectives (e.g. <cit.>) and several standard proximal algorithms (see <cit.>).

The systemic regret functional (<ref>) acts as a performance criterion for any optimization algorithm, played by $N$ players, initialized at an arbitrary $\Theta(0)$, and running for $T$ iterations.
In what follows, we take the initial condition of the algorithm as $\Theta(0)=\Theta^\star$ for simplicity, but other choices are possible. In other words, the regret of each player is initialized at zero, as they all start from their individually preferred parameter $\theta^\star_i$, and starts to accumulate as the system moves away from $\Theta^\star$ towards the optimizer of  (<ref>).

§.§ Contributions

We exhibit a regret optimal algorithm, by which we mean an optimizer of (<ref>) when $f_{\theta}:\mathbb{R}^d\rightarrow \mathbb{R}$ is linear in its trainable parameter
[We emphasize that all vectors are represented as column-vectors and not row vectors.]
$\theta\in \mathbb{R}^p$ but non-linear in its inputs $x\in \mathbb{R}^d$ through
a feature map $\phi:\mathbb{R}^d\rightarrow \mathbb{R}^p$, typically with $p\gg d$, namely

\begin{equation}
\label{eq:FRKR}
\end{equation}

In this case, $f_\theta$ is called a finite-rank kernel ridge regressor (fKRR) and a closed-form solution to the ridge-regression problem (<ref>) exists and is given by (<ref>) in Appendix <ref>.

A key insights is that (<ref>) has the same minimizers as the energy functional

\begin{align}
\label{eq:penalized_regret_functional}
{\cal L}(\Theta^{0\cdots T})
\eqdef &
\underbrace{
\sum_{t=0}^{T-1}
\,
\left[
\lambda \left\| \Theta(t+1) - \Theta^\star \right\|_2^2
\beta \left\| \Delta \Theta(t) \right\|_2^2
\right]
\underbrace{
l (\Theta(T), w;\,\mathcal{D} )
Terminal Cost
\end{align}

Every algorithm (<ref>) can be regarded as driven by a “control”, and the functional ${\cal L}$ has the form of a “cost function” of an optimal control problem,  decomposable into the sum of a “running/operational cost” and a “terminal/objective cost”.  This allows us to leverage tools from optimal control <cit.> to derive a  closed-form expression of the unique regret-optimal algorithm driven by the unique “optimal control” as defined in Algorithm <ref>.

We prove that this procedure is adversarially robust, meaning that the value of (<ref>) varies by $\mathcal{O}(\varepsilon \sqrt{\bar{N}}q)$ if a malicious adversary can perturb at most $1/q$ percent of all the data in each of the dataset $\mathcal{D}_i$ by at-most $\varepsilon\ge 0$, where $\bar{N}=\sum_{i=1}^N |\mathcal{D}_i|$.
Furthermore, we show that our regret-optimal algorithm has a computational complexity of $\mathcal{O}(N^2p^3 + T\,(Np)^{2.373})$.
Moreover, we characterize the optimal choice of the weights $w$

which the central planner must use to maximize the selected model's performance on the primary/focal dataset $\mathcal{D}_1$.

Finally, we show that if all the datasets $\mathcal{D}_1$, $\dots$, $\mathcal{D}_N$ are equally dissimilar (but not necessarily similar) then the regret-optimal algorithm can be accelerated.  This acceleration requires $\mathcal{O}(N\,p^2)$ fewer elementary operations.  We also show that if this “accelerated” algorithm is used in situations when the datasets differ by a little, then it does not deviate far from the regret-optimal algorithm.  Thus it can be used as an accelerated heuristic.

§.§ Outline

The paper is organized as follows. In Section <ref>, we introduce the regret-optimal algorithm and demonstrate the optimality, computational complexity, and adversarial robustness of the algorithm. We  further introduce an accelerated algorithm of lower computational complexity that can achieve near-regret optimality.
In Section <ref>, we conduct a convergence analysis of the regret-optimal algorithm and compare it to the standard gradient descent approach. Then we show the transfer learning capabilities of the algorithm for American option pricing.
Appendix <ref> is denoted to proving the optimality, adversarial robustness, and computational complexity of the regret-optimal algorithm. The near-regret optimality and computational complexity of the accelerated algorithm are also proven. Appendix <ref> briefly introduces the background of finite-rank kernel ridge regressors.

§.§ Notation

We briefly introduce some additional notation.
Each dataset has a finite number of samples which we write as $\mathcal{D}_i = \{ (x_j^i, y_j^i) \, | \, 1 \leq j \leq |\mathcal{D}_i| \}$. We assume that some feature map $\phi$ is fixed and denote the features of each sample of $\mathcal{D}_i$ as

$$u_j^i \eqdef \phi(x_j^i).$$

Moreover, we denote the joint inputs and outputs for each dataset as

\begin{align}
X^i \eqdef \left( x^i_1, \cdots , x^i_{|\mathcal{D}_i|} \right)^\top \in \mathbb{R}^{|\mathcal{D}_i| \times d} \quad \text{ and } \quad
Y^i \eqdef (y^i_1, \cdots, y^i_{|\mathcal{D}_i|})^\top \in \R^{|\mathcal{D}_i|}, \notag
\end{align}

and the joint features as

\begin{align}
U^i \eqdef \Phi(X^i) \eqdef
\left( \phi(x^i_1), \cdots , \
\phi(x^i_{|\mathcal{D}_i|}) \right)^\top \in \mathbb{R}^{|\mathcal{D}_i| \times p}. \notag
\end{align}

Finally, we use $I_p$ to denote a $p\times p$ identity matrix.

§.§ Related Work

Examples of $f_\theta$ of the form (<ref>) include extreme learning machine <cit.> or any finite-rank kernel ridge regressor (fKRR) <cit.>.  Random feature models, and finite-rank ridge regressors, are typical in contemporary quantitative finance <cit.> or reservoir computing <cit.>, where it is favourable to rapidly deploy a deep learning model due to time or computing-power limitations/constraints.  Finite-rank kernel ridge regressors are also typical in the setting of standard transfer <cit.> and multi-task learning <cit.> pipelines where one often has a pre-trained deep neural network model and, for any novel task, the user freezes the hidden layers in the deep neural network model and only fine-tunes/trains the deep neural network's final linear layer.
This is useful, since the hidden layers of most deep neural network architectures are known to process general, cruder features and the final layer of any such model captures finer, task-specific features <cit.>.

The optimization of a model trained from multiple data sources lies at the heart of multi-task learning.  In  <cit.>, the authors cast the task of optimizing (<ref>) as a multi-objective optimization problem.  They propose to solve it via a multi-gradient descent algorithm which, under mild conditions, converges to a Pareto stationary point; i.e. which is essentially a critical point of each of the $i$-th player's penalized SSE.
Though iterative schemes for multi-objective optimization can efficiently reach critical points <cit.>, this problem's critical point is not necessarily a minimum since the problem (<ref>) is not convex when $w$ is shared amongst all $N$ players.

This proposed algorithm reflects the ideas of the federated stochastic gradient descent (FedSGD) algorithm of <cit.> and of its online counterpart FedOGD which is an online federated version of online gradient descent[See <cit.> for details on online gradient descent] and more efficient online versions thereof such as <cit.>.  From the federated learning perspective, the problem of training individual fKRR on each dataset, and then centralizing them by optimizing the mixing parameter $w$ in (<ref>) and its origins date back to local gradient descent introduced in the optimization and control literature in <cit.>.  Though several authors (e.g. <cit.>) have proposed various progressive improvements of federated gradient-descent-type algorithms for training a learner from several data sources, the problem of studying the most efficient iterative federated learning algorithm for regret minimization has not yet been tackled.

We note that (stochastic) optimal control tools have revealed new insights into machine learning.  Examples include an optimal step-size control <cit.> and batch-size control <cit.>, developing online subgradient methods that adapt to the data's observed geometry <cit.>, developing new maximum principle-based training algorithms for deep neural networks <cit.>, or unifying (stochastic) optimization frameworks <cit.>.

§ MAIN RESULTS

This section contains the description of the unique regret-optimal algorithm and an analysis of its theoretical properties. Throughout our manuscript, we maintain the following two assumptions.

There are constants $K_x,K_y>0$ such that all
$(x^i_j, y^i_j) \in \mathcal{D}$ satisfy $\|u^i_j\|_2 \leq K_x$ and $|y^i_j|^2 \leq K_y$, where $u^i_j = \phi (x^i_j)$.

The data ${\cal D} = \cup_{i=1}^N\,\mathcal{D}_i$ and the weights $w^{\star}$ satisfy

\begin{align}
& \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|} u^i_j u^{i\top}_j \geq 0 ,
\label{Sum(uij)Tuij>=0}\\
& [w^\star_1 I_p ,\cdots, w^\star_N I_p]^{\top}
\left( \sum_{i=1}^N
w^\star_i \,
\sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j \right)
[ w^\star_1 I_p, \cdots, w^\star_N I_p ] \geq 0 .
\label{WeightedSum(uij)Tuij>=0}
\end{align}

With the notation introduced before, the
object $[w_1^\star I_p, \cdots, w_N^\star I_p ]$ above is a $p \times Np$ matrix.

§.§ The Algorithms

Step 1: Adaptive “Prior” Sharing: The $i^{th}$ dataset $\mathcal{D}_i$ is shared with the main dataset $\mathcal{D}_1$ if the best model $f_{\theta^\star_i}$, optimizing (<ref>) on dataset $\mathcal{D}_i$, can outperform the optimal model trained on the main dataset, up to a slack factor $\eta>0$.  The relative importance of each dataset is proportional to $f_{\theta^\star_i}$'s performance relative to the training performance of the best model trained on the main dataset.
When $\eta$ is small, sharing is unlikely; when $\eta$ is large, nearly all datasets are pooled together with approximately the same relative importance.

Step 2: Optimized “Posterior” Sharing: We then optimize the relative importance/weights given to the shared datasets from the naive adaptive sharing.  The optimization is a closed-form update of the adaptive “prior” weights which optimize a certain maximum expected performance criterion subject to a relative entropy penalization against the weights obtained from adaptive prior sharing.

Visual Summary of Optimal Information Sharing (Algorithm <ref>).

We begin by describing our procedure for determining the optimal mixture weights defining the cooperative objective function (<ref>) used in  (<ref>).  This procedure is illustrated graphically in Figure <ref> and detailed in Algorithm <ref>.

Optimal Information Sharing - Initialization to Algorithm <ref>

Datasets $\mathcal{D}_1,\dots,\mathcal{D}_N$, finite-rank kernel $\phi$, hyperparameters[Different ridge hyperparameters $\kappa$ can be used to define each regressor in Algorithm  <ref>] $\kappa,\eta>0$.
Initialize Locally-Optimal Learners and Record Scores on Main Dataset

ForParallelFor $i:1,\dots,N$ in parallelend
$\theta^\star_i \leftarrow
( \Phi^\top(X^i) \Phi(X^i)
+ \kappa I_p )^{-1}
\Phi^\top(X^i) Y^i$
*Optimize each fKRR
\leftarrow
\frac1{|\mathcal{D}_1|}\sum_{j=1}^{|\mathcal{D}_1|}\,
\big\|
-  y_j^1
\big\|_2^2

*Score predictive power of fKRR on main dataset

Get Adaptive Prior
ForParallelFor $i:1,\dots,N$ in parallelend
\leftarrow
\frac{
}{ \sum_{j=1}^N\, [s^1+\eta-s^j]_+ }
Get Optimized Posterior
ForParallelFor $i:1,\dots,N$ in parallelend
\leftarrow
\frac{
\,
\sum_{j=1}^N\,
Optimal Weights $w^\star$ and Locally-Optimal fKRR Parameters $\theta^\star_1,\dots,\theta^\star_N$.

Suppose now that we have used Algorithm <ref> to compute weights $w^\star=(w^\star_1, \cdots, w^\star_N)$ and the locally optimal fKRR parameters $\theta^\star_1$, ..., $\theta^\star_N$$w^{\star}$. The main novelty in the next algorithm follows a similar spirit to <cit.> and <cit.>, by casting the problem of constructing a regret-optimal algorithm as an optimal control problem.

As before, denote

\begin{align}
& \Theta(t)\eqdef [ \theta_1^{\top}(t), \cdots, \theta_N^{\top}(t) ]^{\top}, \quad
\Theta^\star \eqdef [ \theta^{\star \top}_1, \cdots, \theta^{\star \top}_N ]^{\top} ,
\notag \\
& \boldsymbol\alpha(t) \eqdef [\alpha_1^{\top}(t), \cdots, \alpha_N^{\top}(t)]^{\top}, \notag
\end{align}

which are all elements of $\R^{Np}$, and introduce the dynamics that $\Theta$ follows
\begin{align}
& \Theta(0)=\Theta^\star,
\quad
\Theta(t+1) = \Theta(t) + \boldsymbol\alpha(t) ,
\quad 0\leq t \leq T-1 ,
\label{Theta-simple}
\end{align}

which is equivalently written as
\begin{align}
\Theta(t+1) = \Theta^\star + \sum_{u=0}^{t} \boldsymbol\alpha(u) , \quad 0\leq t \leq
T-1 .
\label{Theta(t)Theta(0)}
\end{align}

Then, since the optimizers of the energy (<ref>) coincide with the optimizers of the systemic regret functional (<ref>), searching for a regret-optimal algorithm reduces to solving for an optimal control $\mathbf{\alpha}$ to minimize the cost
(<ref>) subject to (<ref>).
The optimal control problem
(<ref>) and (<ref>) is a discrete-time linear quadratic optimal control problem, and can be solved by dynamic programming methods that are standard in the literature. For instance, by <cit.>, we obtain the solution for the regret-optimal algorithm in closed-form as solutions of a system of Riccati equations (<ref>) and (<ref>), given in
Theorem <ref>.

Algorithm <ref> below implements the optimal control approach in two steps. In Step 1, we solve the equation system (<ref>) and (<ref>) backwards for $t=T$, $T-1$, ..., $1$. Then in Step 2, we obtain the optimal control $\boldsymbol\alpha(t)$ and update the parameter $\Theta(t)$ for $t=0$, $1$, ... $T$.

Regret-Optimal Optimization Algorithm

Datasets $\mathcal{D}_1,\dots,\mathcal{D}_N$, $N$. Iterations $T\in \mathbb{N}_+$, finite-rank kernel $\phi$ and hyperparameters $\lambda, \beta, \kappa, \eta >0$.
Get Initialize Weights and Locally-Optimal fKRR Parameters

$\theta^\star_1,\dots,\theta^\star_N,w^{\star}$ $\leftarrow$ Run: Algorithm <ref> with $\mathcal{D}_1,\dots,\mathcal{D}_N$, $\phi$, and $\kappa,\eta$.

Initialize Updates
= [w_1^\star I_p , \cdots, w_N^\star I_p ]^\top \left(\sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u_j^{i\top}\right) [w_1^\star I_p, \cdots, w_N^\star I_p ] $

$S(T) = -[w_1^\star I_p , \cdots, w_N^\star I_p ]^\top \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j y^i_j$

$\Theta(0) = \Theta^\star \eqdef (\theta_1^{\star \top}, \cdots, \theta_N^{\star \top} )^\top$

Generate Iterates
Update Driving Parameters

$P(t)= \beta I_{Np} - \beta^2 [(\lambda+\beta)I_{Np} + P(t+1)]^{-1}$

$S(t) = \beta [(\lambda + \beta )I_{Np} + P(t+1)]^{-1} (S(t+1) - \lambda \Theta^\star)$

$t=0, ..., T-1$
Update Control

$\boldsymbol{\alpha}(t) = -[(\lambda+\beta )I_{Np} + P(t+1)]^{-1}
[(\lambda I_{Np} + P(t+1))\Theta(t) - \lambda \Theta^\star + S(t+1) ]$

Return the Optimized fKRR $f_{\theta^{w^\star}}$

In the next subsections, we state our primary results, whose proofs are presented in Appendix <ref>. We first quantify the optimality of Algorithm <ref>. Next, we establish the optimality, complexity, and adversarial robustness of Algorithm <ref>.

§.§ Properties of Algorithm <ref>

By decoupling the optimization of the “dataset-relevant weights” $w$ in Algorithm <ref> from the optimization of the parameter $\theta$ in Algorithm  <ref>, we avoid open-loop forward-backward iterations, which iterate between optimizing $w$ and $\theta$.  We now study the optimality of the weight $w^{\star}$ generated by Algorithm <ref>.  The weights $w^{\star}$ represent the optimal probability distributions of randomly selecting amongst the $N$ optimized regressors $\{ f_{\theta^\star_i} \}_{i=1}^N$, each optimized locally on their dataset $\mathcal{D}_i$ by solving (<ref>).  Here, optimality is quantified in terms of the best-expected prediction on the focal dataset $\mathcal{D}_1$ subject to a penalty for deviating too far from a reference distribution.  This reference distribution serves to automatically delete models whose datasets are too dissimilar to $\mathcal{D}_1$.  Inspired by most meta-learning problems, e.g. <cit.>, similarity is quantified by extracting a “similarity score” from each learner.

Let $\mathcal{P}(\{\theta^\star_i\}_{i=1}^N)$ consist of all probability measures on the set $\{\theta^\star_i \}_{i=1}^N$.  Any probability measure $\mathbb{P}_w$ in $\mathcal{P}(\{\theta^\star_ i\}_{i=1}^N)$ is uniquely determined by a weight $w$ in the $N$-simplex; consisting of all $w\in [0,1]^N$ whose entries sum to $1$.  Let $\theta$ be a random variable with law $\mathbb{P}_{w}$, for some $\mathbb{P}_w\in \mathcal{P}(\{\theta^\star_i \}_{i=1}^N)$.  The average performance of the random learner $f_{\theta}$ on the dataset $\mathcal{D}_1$ is

\begin{equation}
\label{eq:optimizingQ}
\mathbb{E}_{\theta \sim \mathbb{P}_w}\,
\left[
\frac1{|\mathcal{D}_1|}\sum_{j=1}^{|\mathcal{D}_1|}\,
\big\|
\big\|_2^2
\right]
\sum_{i=1}^N\,
%%% Score
\underbrace{
\biggl[
\frac1{|\mathcal{D}_1|}\sum_{j=1}^{|\mathcal{D}_1|}\,
\big\|
\big\|_2^2
\biggl]
}_{i^{th} \mbox{ learner's score }(s_i)} .
\end{equation}

Fix an information sharing level $ \eta > 0$.  The information sharing level $\eta$ adaptively determines which of the datasets in $\mathcal{D}$ are implicitly shared with the principle dataset $\mathcal{D}_1$.  The information sharing is formalized through the probability measure

\[
\bar{\mathbb{P}}^{\eta}(\theta^\star_i) \eqdef
\frac{
[s^1+\eta -s^i]_+
\sum_{i=1}^N\, [s^1 + \eta - s^i]_+
\]

on the set of $N$ learners, or equivalently on the parameters $\theta^\star_1,\dots,\theta^\star_N$.
The measure $\bar{\mathbb{P}}^{\eta}$ thus allows us to minimize the regularized version of the loss in (<ref>),

\begin{equation}
\label{eq:optimizingQ__regularized}
\min_{
\mathbb{P}_w
\,
\,
\underbrace{
\mathbb{E}_{\theta \sim \mathbb{P}_w}\biggl[
% \Big[
\frac1{|\mathcal{D}_1|}\sum_{j=1}^{|\mathcal{D}_1|}\,
% w_i\,
\big\|
\big\|_2^2
\biggr]
}_{i^{th} \mbox{ learner's score }(s_i)}
\underbrace{
% (1-\eta)
\eta\,
\operatorname{D}_{\operatorname{KL}}\big(
\mathbb{P}_w
\|
\bar{\mathbb{P}}^{\eta}
\big)
\underset{\mbox{information sharing}}{\mbox{Inhomogeneity of}}
\end{equation}

where the infimum is taken over all $\mathbb{P}_w$ in $\mathcal{P}(\{\theta^\star_i\}_{i=1}^N)$ and  $\operatorname{D}_{\operatorname{KL}}$ denotes the Kulbeck-Leibler divergence with respect to $\bar{P}^{\eta}$ given for any probability measure $\mathbb{P}_w$ on $\{\theta^\star_i\}_{i=1}^N$, namely

\[
\operatorname{D}_{\operatorname{KL}}\big(\mathbb{P}_w\| \bar{\mathbb{P}}^{\eta})
\eqdef
% \begin{cases}
\sum_{i=1}^N\,
% \Biggl(
\log\biggl(
% \sum_{j=1}^N\,[s^1+\eta-s^j]_+
\frac{w_i}{
% \Big)
% +
% \log\big(
% \sum_{j=1}^N\,[s^1+\eta-s^j]_+
% \big)
\biggr)
% \Biggr)
%         & \mbox{ if } \mathbb{P}\ll \mathbb{Q} \\
%     \infty & \mbox{ otherwise }
% \end{cases}
\]

where $C_{\eta}\eqdef  \sum_{j=1}^N\,[s^1+\eta-s^j]_+$, if $\mathbb{P}_w$ is absolutely continuous with respect to $\bar{\mathbb{P}}^{\eta}$ and otherwise $\operatorname{D}_{\operatorname{KL}}\big(\mathbb{P}_w\| \bar{\mathbb{P}}^{\eta})=\infty$.
There is a unique probability measure optimizing the penalized functional (<ref>) determined by the weight $w^{\star}$ in Algorithm <ref>.

Fix $\eta>0$.  The unique weights $w^{\star}_i$ on the $N$-simplex minimizing (<ref>) are given by

\begin{equation}
\label{eq:prop:Optimal_Bayes_Weighting__OptimalWeights}
\begin{aligned}
\frac{
\,
\sum_{j=1}^N\,
, \qquad  i=1,\dots,N.
\end{aligned}
\end{equation}

§.§ Guarantees for Algorithm <ref>

Our next result shows that Algorithm <ref>, and its iterates, are regret-optimal.

Suppose that Assumptions <ref> and <ref> hold and fix $\Theta^{\star}\in \mathbb{R}^{N\,p}$ as in (<ref>).  Then, the systemic regret functional $\mathcal{R}$ defined in (<ref>) admits a unique minimizer $\hat\Theta_{\cdot}\eqdef (\hat\Theta(t))_{t=0}^T$ of the form (<ref>) with $\hat\Theta(0) = \Theta^{\star}$ aand increments

\begin{align}
\Delta \hat\Theta(t)
= & -[(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
\big[ (\lambda I_{Np} + P(t+1)) \hat\Theta(t)
- \lambda  \Theta^\star + S(t+1)  \big] ,   \label{optimala 0}
\end{align}

where, for $t=1,\dots,T$,  the matrices $P(t)$ and $S(t)$ in (<ref>) are determined by

\begin{align}
& \begin{cases}
P(T) =  [w^\star_1 I_p ,\cdots, w^\star_N I_p]^{\top}
\left( \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}   u^i_j u^{i\top}_j \right)
[ w^\star_1 I_p, \cdots, w^\star_N I_p] \\
P(t)  =  \beta I_{Np} - \beta^2 [(\lambda+\beta) I_{Np} + P(t+1) ]^{-1}  \label{eqnP(t) 0}
\end{cases},   \\
& \begin{cases}
S(T) =  - [ w^\star_1 I_p, \cdots, w^\star_N I_p]^{\top} \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  {u^i_j} y^i_j \\
S(t) =  \beta \big[ (\lambda+\beta)I_{Np} + P(t+1) \big]^{-1} (S(t+1) - \lambda \Theta^\star )
\label{eqnS(t) 0}
\end{cases} .
\end{align}

In particular, for each $t=0,\dots,T-1$

\[
\Delta \hat\Theta(t)
\boldsymbol{\alpha}(t)
\]

where $\boldsymbol{\alpha}$ is as in Algorithm <ref>.

Instead of the sums we can equivalently use our matrix notation to write $\sum_{j=1}^{|\mathcal{D}_i|} u^i_j u^{i\top}_j = {U^i}^\top U^i$ and $\sum_{j=1}^{|\mathcal{D}_i|}  {u^i_j} y^i_j = U^{i\top} Y^i$.

The explicit update rules for the regret-optimal optimization algorithm presented above permit us to compute its complexity.  Since we will be assuming that we work with floating-point arithmetic, all elementary arithmetic operations (e.g. addition, subtraction, multiplication) are of constant (i.e. $\mathcal{O}(1)$) complexity.  We are most interested in the regime where the number of datasets is much larger than the typical number of instances/data per dataset; i.e. when

\begin{equation}
\label{eq:moredatasets_than_dataperdataset}
\frac1{N}\sum_{i=1}^N\,|\mathcal{D}_i|
\le
\end{equation}

In the setting of Theorem <ref>, the computational complexity[We emphasize that we have assumed that computational complexity of all elementary arithmetic operations (e.g. addition, subtraction, multiplication) is $\mathcal{O}(1)$.] of computing the sequence $(\hat \Theta(t))_{t=0}^T$ defining the regret-optimal algorithm is $
\mathcal{O}\big( N^2 p^3 +
\,
%N^{2.373}\, p^2
\big)

Next, we quantify the sensitivity of Algorithm <ref> to adversarial attacks.  We study the robustness of Algorithm <ref> with respect to a family of perturbations that switch a given percentage $0\le q\le 1$ of all training data by selecting fake data points at a distance of at-most $\varepsilon\ge 0$ for any perturbed datapoint.  Thus, for a fixed set of datasets
$\{ \mathcal{D}_i \}_{i=1}^N$, and $q,\varepsilon$ as above,
we define the class $\mathbb{D}^{q,\varepsilon}$ of adversarially generated datasets of “magnitude” $(q,\varepsilon)$ as all sets of datasets
$\{ \widetilde{\mathcal{D}}_i \}_{i=1}^N$, $\widetilde{\mathcal{D}}_i=\{ \{(\tilde{x}^i_j,\tilde{y}^i_j) \, | \, 1 \leq j \leq |\mathcal{D}_i| \}\subseteq \mathbb{R}^{d}\times \mathbb{R}$ satisfying Assumptions <ref> and <ref>, as well as

\[
%%%% Number of Perturbed Instances
\underbrace{
\frac{
\#\{
(i,j):\, u^i_j\neq \tilde{u}^i_j \mbox{ or } y^i_j\neq \tilde{y}^i_j
\}
\sum_{i=1}^N\,|\mathcal{D}_i|
\le
}_{\mbox{Attack Persistence}}
\mbox{ and }
%%%% Maximum Magnitude of Adversarial Perturbations
\underbrace{
\max_{i,j} \left\{ \,
\max\left(
\|u^i_j-\tilde{u}^i_j\|_2
\right) \right\}
\le
\varepsilon
}_{\mbox{Attack Severity}}
\]

where the indices $(i,j)$ run over $\{(i,j) \, | \, 1 \leq i \leq N,\,1 \leq j \leq |\mathcal{D}_i|\}$ and $\tilde{u}^i_j\eqdef \phi(\tilde{x}^i_j)$.

In view of Theorem <ref>, the updates of the unique regret-optimal algorithm $\hat \Theta^{0\cdots T}$ are determined by the sequence $\boldsymbol{\alpha}_{\cdot}\eqdef (\boldsymbol{\alpha}(t))_{t=0}^T$ computed by Algorithm <ref>.  It is convenient to emphasize this connection on $\boldsymbol{\alpha}_{\cdot}$ and the connection with the dataset $\mathcal{D}$ when expressing the following theorem, with the energy (<ref>) determining the value of the systematic regret functional $\mathcal{R}$.  We therefore write $
{\cal L}(\boldsymbol{\alpha};\mathcal{D})
\eqdef
{\cal L}(\hat\Theta^{0\cdots T})

Fix a $\mathcal{D}$ satisfying Assumptions <ref> and <ref>.  For any level of “attack persistence” $0\le q\le 1$ and any degree of “attack severity” $\varepsilon\ge 0$ we have that

\[
\big| {\cal L} (  \boldsymbol\alpha ; \mathcal{D} )
{\cal L}(  \widetilde{\boldsymbol\alpha} ;  \widetilde{\mathcal{D}} )  \big|
\le
\mathcal{O}\Big(
\varepsilon
\sqrt{\bar{N}\,q}
\Big)
\]

$\widetilde{\boldsymbol\alpha}$ is the output of Algorithm <ref> for any given adversarially-generated set of datasets
$\{ \widetilde{\mathcal{D}}_i \}_{i=1}^N\in \mathbb{D}^{q,\varepsilon}$ and $\bar{N}=\sum_{i=1}^N |\mathcal{D}_i|$.  Here $\mathcal{O}$ hides a constant depending only on Assumption <ref>.

§.§ An Accelerated Heuristic

We now discuss a situation in which the computational complexity of the regret-optimal algorithm can be reduced.  This is achieved by exploiting symmetries in the optimal policy computed by Algorithm <ref>.

Examining the explicit update rules, we notice that if the weights $w$ defining the energy (<ref>) were to be $\overline{1/N}\eqdef (1/N,\dots,1/N)$ then the matrices defining the optimal policy/control's updates would become highly symmetric Toeplitz matrices, depending on exactly three $p$-dimensional block sub-matrices.  This allows us to parameterize the entire regret-optimal algorithm, for the sub-optimal weights $\overline{1/N}$, using only three low-dimensional systems whose dimension is independent of the number of datasets $N$.  In line with this, Algorithm <ref> implements the regret-optimal algorithm for the weight specification $\overline{1/N}$ by only keeping track of these three optimal low-dimensional systems, and then efficiently recombining them by leveraging the Toeplitz structure defining the optimal policy/control.  These low-dimensional systems are defined using helper functions, listed in equations (<ref>)-(<ref>).

The approximate near-optimality of Algorithm <ref> only depends on the stability of the optimal policy/control, defining the algorithm's updates as a function of the weights $w$.  If the optimal weights defining the systemic regret functional (<ref>) are close to $\overline{1/N}$, then Algorithm <ref> is nearly regret-optimal.  Otherwise, the sub-optimality it incurs in favor of computational efficiency is explicitly quantifiable and depends only on the difference between the optimal weights computed in Algorithm <ref> and the symmetric weights $\overline{1/N}$.

Note that, we record the matrix-valued functions $\gamma_1(\cdot)$ and $\gamma_2(\cdot)$ in Algorithm <ref>, directly below it.

Accelerated (Nearly) Regret-Optimal Algorithm

Datasets  $\mathcal{D}_1,\dots,\mathcal{D}_N$, $N$. Iterations $T\in \mathbb{N}_+$, finite-rank kernel $\phi$, hyperparameters $\lambda, \beta, \kappa, \eta >0$.
Get Initialize Weights and Locally-Optimal fRKR Parameters

$\theta^\star_1,\dots,\theta^\star_N,w^{\star}$ $\leftarrow$ Run: Algorithm <ref> with $\mathcal{D}_1,\dots,\mathcal{D}_N$, $\phi$, and $\kappa,\eta$.
Initialize Updates
= \pi_2(T)
= (1/N^3) \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j$
=  - (1/N^2) \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|}  {u^i_j} y^i_j $
$\theta_1(0) = \theta^\star_1$, $\theta_2(0) = \theta^\star_2$, $\dots$, $\theta_N(0) =\theta^\star_N$
$\theta^\star_{(N)} = \frac{1}{N} \sum_{i=1}^N \theta^\star_i .$
Generate Iterates
Update Low-Dimensional Driving Parameters

$\pi_1(t)= \beta I_p- \beta^2 \gamma_1(t+1)$
$\gamma_1(\cdot)$ is defined by (<ref>)

$\pi_2(t) = - \beta^2 \gamma_2(t+1)$
$\gamma_2(\cdot)$ is defined by (<ref>)

= \beta \big[ \gamma_1(t+1)
+ (N-1) \gamma_2(t+1)
\big] ( \pi_3(t+1) -\lambda \theta^\star_{(N)})$

$t=0, ..., T-1$
Update Low-Dimensional Control
$i=1, ..., N$
=  \widehat{\alpha}\big(\theta_i(t) ,\
(\theta_j(t))_{j = 1,j\neq i}^N ,\ (\gamma_i(t+1))_{i=1}^2,
\ (\pi_i(t+1))_{i=1}^3 \big)$

$\theta_i(t+1) =\theta_i(t) + \widehat{\alpha}_i(t)$

Synchronized fKRR

The matrix-valued functions $\gamma_1(\cdot)$ and $\gamma_2(\cdot)$ in Algorithm <ref> are defined as

\begin{align}
& \gamma_1(\cdot) \eqdef
\big\{ (\lambda+\beta )I_p
+\pi_1(\cdot) \notag \\
& \hspace{2cm} - (N-1) \pi_2(\cdot)\big[ (\lambda+ \beta )
I_p + \pi_1(\cdot)
+ (N-2)\pi_2(\cdot)]^{-1} \pi_2(\cdot) \big\}^{-1} ,
\label{gamma1} \\
& \gamma_2(\cdot)  \eqdef
- \gamma_1(\cdot) \pi_2(\cdot)
\big[ (\lambda+\beta )I_p + \pi_1(\cdot)
+ (N-2)\pi_2(\cdot) \big]^{-1},  \label{gamma2}
\end{align}

\begin{align*}
\pi_1(T) = \pi_2(T) = \frac{1}{N^3} \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j
\qquad
\mbox{ and }
\qquad
\pi_3(T) =  - \frac{1}{N^2} \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|}  {u^i_j} y^i_j
\end{align*}

and, for $t=T-1,\dots,0$,

\begin{align*}
\pi_1(t) &= \beta I_p- \beta^2 \gamma_1(t+1) \\
\pi_2(t) &= - \beta^2 \gamma_2(t+1) \\
\pi_3(t) &= \beta \big[ \gamma_1(t+1)
+ (N-1) \gamma_2(t+1)
\big] \left( \pi_3(t+1) -\lambda \theta^\star_{(N)}\right)
\end{align*}

where $\theta^\star_{(N)} = \frac{1}{N} \sum_{i=1}^N \theta^\star_i$.

The updates of the low-dimensional systems which Algorithm <ref> exploits are implemented via the helper function $\hat \alpha$, defined as

\begin{align}
\label{eq:RegretOptimizationHeuristic_Helper__1}
\widehat{\alpha}\big( \theta_i, ((\theta_j)_{j=1,j\neq i}^N) , (\gamma_i)_{i=1}^2,(\pi_i)_{i=1}^3 \big)  \eqdef &
- (\gamma_1-\gamma_2)(\lambda I_p + \pi_1 - \pi_2 )\theta_i \notag \\ &
- \big\{ \gamma_1 \pi_2 + \gamma_2 [ \lambda I_p + \pi_1 + (N-2)\pi_2 ] \big\} \sum_{j=1}^N \theta_j \notag \\
& - \big[ \gamma_1 + (N-1) \gamma_2 \big]
(\pi_3 - \lambda \theta^\star_{(N)}) .
\end{align}

The degree to which Algorithm <ref> is nearly regret-optimal depends on how far the weights $w$ determined by Algorithm <ref> are from the equal weighting $\overline{1/N}$.  This is precisely quantified by the following results.

Suppose that $\mathcal{D}$ satisfies Assumptions <ref> and <ref>.  and let $\hat\Theta^{0\cdots T}$ denote the unique regret-optimal algorithm[That is, the unique optimizer of (<ref>)] computed according to Algorithm <ref>. Then, the output $\vartheta^{0\cdots T}$ of Algorithm <ref> satisfies

(i) Near-Optimality: The difference $\big|
{\cal L}(\Theta^{0\cdots T})
{\cal L}(\vartheta^{0\cdots T})
\big|$ is of the order

\begin{align}
%    & \big|
%         {\cal L}(\Theta^{0\cdots T})
%         -
%         {\cal L}(\vartheta^{0\cdots T})
%     \big|
% \in
\mathcal{O}\Big(
\bar{N} N^3 \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i\Big|
+ \bar{N} N^{3/2} \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|
\notag \\
& \hspace{5cm} + \bar{N}  \max_{1\leq i \leq N} \Big|\frac{1}{N} - w_i^\star \Big|
+ N \big\|\Theta^\star - \Theta^\star_{(N)} \big\|_2
\Big)
\notag
\end{align}

(ii) Near-Optimal Increments: The difference $\|\Delta \Theta^{0\cdots T-1}-\Delta \vartheta^{0\cdots T-1} \|_{\ell^1}$ is of the order

\begin{align}
%    \|\Delta \Theta^{0\cdots T-1}-\Delta \vartheta^{0\cdots T-1} \|_{\ell^1}
% \in
\mathcal{O}\Big(
\bar{N} N^{3/2} \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big| \notag \\
& %\hspace{5cm}
+ \bar{N} N^{1/2} \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|
+  \big\|\Theta^\star - \Theta^\star_{(N)}\big\|_2
\Big),
\notag
\end{align}

where $\mathcal{O}$ hides a constant depending only on $\lambda$, $\beta$,
$\Theta^{\star}$, $(w_1^\star, ..., w_N^\star)$, $T$, $K_x$, and in $K_y$.

The next result shows that Algorithm <ref> accelerates the regret-optimal algorithm by a linear factor in $N$ and by a quadratic factor in $p$, when $N$ is large.

With the convention that floating-point arithmetic operations are $\mathcal{O}(1)$, we have that:

(i) Complexity: Algorithm <ref> has a complexity of $\mathcal{O}(TNp\max\{N,p\}))$

(ii) Acceleration over Regret-Optimal Algorithm: The regret-optimal algorithm has a complexity of $\mathcal{O}\big(TN^3 p^3\big)$.

When $N\ge p$, Algorithm <ref> requires $\mathcal{O}(N\,p^2)$ fewer operations than the regret-optimal algorithm.

Finally, we establish that Algorithm <ref> is the regret-optimal algorithm, optimizing the systemic regret functional (<ref>), in the “symmetric” case where the optimal weights $w^{\star}$ are all $(1/N,\dots,1/N)$.

Suppose that $w^{\star}_i=1/N$ for $i=1,\dots,N.$  Then, $\hat\Theta^{0\cdots T}=\vartheta^{0\cdots T}$.

§ EXPERIMENTS

We now experimentally verify our theoretical results, with a focus on quantitative finance.
In particular, we implemented the regret-optimal algorithm as well as its accelerated version and provide them together with the code for all experiments described in this section at <https://github.com/FlorianKrach/RegretOptimalFederatedTransferLearning>.

First, we conduct a convergence analysis of our regret-optimal algorithm, comparing it to standard gradient descent. We then demonstrate the transfer learning capabilities of our algorithm in the context of American option pricing.
We use the randomized neural networks approach of <cit.> as our fKRR model, since it best balances computational efficiency and expressiveness amongst the available non-linear versions of the <cit.> algorithm for American option pricing.

§.§ Convergence of the Regret-Optimal Algorithm

We first study the evolution of the energy functional  (<ref>) and loss (<ref>) achieved in each iteration by our regret-optimal Algorithm <ref> with weights initialized by Algorithm <ref>. That is to say, we compute ${\cal L}(\Theta^{0\cdots t})$ as defined in (<ref>) but for intermediate iterations $t$, as well as the loss $l(\Theta(t), w;\,\mathcal{D})$.

We compare the regret-optimal algorithm with its accelerated version, given by Algorithm <ref> with weights given by Algorithm <ref>, where we use $\Theta(0) = \Theta^\star_{(N)}$ and replace $\Theta^\star$ by $\Theta^\star_{(N)}$ in (<ref>).
Moreover, we compare to the results when using gradient descent to optimize the loss (<ref>) directly.
For all algorithms we use a ridge coefficient $\kappa=0$, regret coefficients $\lambda = 0$, $\beta=1$  and Algorithm <ref> with information sharing level $\eta=10$ to get the corresponding optimal weights $w^\star$ and to initialize $\Theta(0)$ with the locally optimal parameters.
We note that the comparison of our algorithms to gradient descent optimizing the loss (<ref>) only makes sense in the $\lambda=0$ regime, since otherwise the objectives and therefore the behaviour are different. On the other hand, using $\beta > 0$ doesn't change the objective, but only forces our algorithm to make equally sized steps $\Delta \Theta (t)$, where the step size depends on $\beta$ and $T$.

For our convergence analysis we use $N=5$ randomly generated datasets. In particular, each dataset $\mathcal{D}_i$ corresponds to a neural network $g_{\omega_i}$ with fixed randomly-generated hidden weights and biases.  Here $\omega_i$ denotes the vector of all these hidden weights and bias parameters. The networks have $d=5$ dimensional input, one hidden layer with $10$ hidden nodes, $\operatorname{ReLU}$ activation function and map to a $1$-dimensional output. The dataset is generated by sampling $S_i = 100$ inputs $x^i_j$ from a $d$-dimensional standard normal distribution and the corresponding target values are given by $g_{\omega_i}(x^i_j)$, i.e.,

\begin{equation*}
\mathcal{D}_i = \{ (x^i_j, g_{\omega_i}(x^i_j)) \, | \, 1 \leq i \leq S_i, x^i_j \sim N(0,I_d) \}.
\end{equation*}

Our feature map $\phi$ in (<ref>) is a randomized neural network
[Also called a random feature model or an extreme learning machine]
with one hidden layer, $500$ hidden nodes and $\operatorname{ReLU}$ activation as studied by <cit.>. In particular, the weights of the hidden layers are randomly initialized and fixed and only the weights of the last layer, corresponding to $\theta$ in (<ref>), are optimized. We note that by our choice for the number of hidden nodes to equal the total number of samples, the algorithms can in principle learn to have perfect replication on the training set, i.e., loss $l(\Theta(t), w;\,\mathcal{D}) = 0$.
For a fair comparison between the algorithms, we use the same randomly generated datasets and the same randomized neural networks for each of the methods.

Figure <ref> shows the averaged results for 10 runs each with a different set of randomly generated datasets and a different randomized neural network as fKRR. For the regret-optimal algorithm and for its accelerated version we use $T=10^3$ iterations, while we use $10^5$ steps for gradient descent with learning rate $7 \cdot 10^{-5}$. We note that for learning rates $\geq 8 \cdot 10^{-5}$ the training becomes unstable such that the loss explodes.
The gradient descent method gradually but slowly reduces the loss.
Since the gradient descent method makes so small parameter updates (and since $\lambda=0)$, the regret is nearly $0$ and the energy ${\cal L}$ nearly equal to the loss. Even though it runs $100$ times the number of iterations of our algorithms, it would still need much more training to achieve the same loss.
On the other hand, the regret-optimal algorithm and its accelerated version behave nearly indistinguishably, reducing the loss nearly to $0$ and piling up a little bit of regret. We note that by increasing the number of iterations, our algorithms could reduce the loss further by making smaller parameter updates and therefore piling up less regret, allowing them to focus even more on reducing the loss.

Loss and energy ${\cal L}$ for gradient descent (GD), the regret-optimal algorithm (RO) and the accelerated regret-optimal algorithm (ARO). Means and standard deviations over 10 runs are shown. Left: all iterations; middle: first 1000 iterations; right: iterations 500 to 1000 for our algorithms.

In the regime where $\lambda > 0$ the behaviour of the algorithm changes drastically. In particular, since deviating from $\Theta(0) = \Theta^\star$ is penalised, the algorithm remains at the starting point in the initial steps and only starts to converge when approaching $T$. The exact behaviour depends on the choice of $\lambda$ (and is additionally influenced by the choice of $\beta$). In Figure <ref> we compare the behaviour for $\lambda=10^{-4}$ and $\lambda = 2$ while otherwise using the same setup as before.  For the smaller value of $\lambda$, the algorithm starts the decent earlier and decreases the loss by more compared to the larger value for $\lambda$, where the algorithm deviates from $\Theta^\star$ only in the very last steps.

Loss and energy ${\cal L}$ for the regret-optimal algorithm (RO) and the accelerated regret-optimal algorithm (ARO). Means and standard deviations over 10 runs are shown. Left: $\lambda = 10^{-4}$; right: $\lambda=2$.

§.§ American Option Pricing

We test our algorithm in the context of the optimal stopping problem for American option pricing, a computationally highly challenging problem which has recently become tractable also in higher dimensional settings through deep learning <cit.>. This is a good example for our purpose for two reasons.
On the one hand, <cit.> did show that randomized neural networks, which are one class of fKRR, are well suited to (approximately) solve these optimal stopping problems. Moreover, it was shown that using randomized neural networks can heavily improve the computation time, which is also of importance to us since we need to solve several optimizations problems for a potentially large number of datasets.
On the other hand, to solve the optimal stopping problem one needs to generate synthetic datasets. Hence, this problem allows us to generate multiple synthetic datasets of different sizes and varying similarity on which we can test the transfer-learning capabilities of our algorithm. The performance can be measured by comparing the achieved prices with those of several baseline methods.

We first provide a short recap of the optimal stopping problem in American option pricing and then present our transfer-learning experiments.

§.§.§ Background: American Option Pricing and Optimal Stopping

An American option is a financial derivative that gives its owner the right to execute a certain predetermined trade at any time between the start date and the maturity of the option. For example, a call (put) option gives its owner the right to buy (sell) a certain stock from the writer of the option for a fixed price $K$, called the strike. If the current stock price $X_t$ is higher (lower) than $K$ then the option owner can make a profit of $(X_t - K)_+$ (or $(K - X_t)_+$ respectively) by directly selling (buying) the stock again in the open market after having exercised their right to buy (sell) from the option writer.
While American options allow for continuous-in-time execution, it is common practice to approximate them via Bermudan options, which can be executed on a predetermined time grid $t_0 < \dotsb t_M$. We briefly overview the American option pricing problem, details of which can be found in <cit.>.

In order to capture the computational challenges posed by the American option pricing problem in higher dimensions, we consider a stochastic process $(X_t)_{t=0}^M \in (\R^d)^{(M+1)}$ describing the price of $d$ stocks on the time grid

and let $(Z_t)_{t=0}^M \in \R^{(M+1)}$ be the corresponding discounted payoff process of an American option written on a one-dimensional process obtained from $X_t$ (e.g. the largest price among the $d$ stocks). We define the filtration $\mathbb{F} = (\mathcal{F}_t)_{t=0}^M$ generated by $X$ via $\mathcal{F}_t = \sigma(X_s | s\leq t)$. Then the price of the American option is given by the starting value $U_0$ of the Snell envelope <cit.>, given by

\begin{equation}\label{equ:snell envelope 2}
U_m = \operatorname{sup}_{\tau \in \mathcal{T}_m} \E[ Z_{\tau} \, | \, \mathcal{F}_m ],
\end{equation}

where $\mathcal{T}_m$ is the set of all stopping times $\tau \geq m$.
The smallest optimal stopping time is

\begin{align}\label{equ: optimal stopping time}
\tau_M &:= M, \\
\tau_m &:= \begin{cases}
m, & \text{if } Z_m \geq \E[U_{m+1}  \, | \, \mathcal{F}_m],\\
\tau_{m+1}, & \text{otherwise},
\end{cases} \nonumber
\end{align}

hence, to compute the price of the American option it suffices to compute the continuation values $\E[U_{m+1}  \, | \, \mathcal{F}_m]$ for all $m$.
In the following, we will use the randomized least squares Monte Carlo (RLSM) algorithm of <cit.>, which utilises randomized neural networks to approximate the continuation values via backward induction (note the backward recursive scheme that (<ref>) and (<ref>) imply).
In particular, this leads to $M$ backward-recursively defined regression problems that need to be solved sequentially. Since (<ref>) is a maximisation problem and since we only approximate the optimal stopping time, we compute lower bounds of the price $U_0$.
The standard procedure is to fix some law for the process $X$ and to sample a dataset of $\nu = \nu_1+\nu_2$ i.i.d. paths of $X$, where $\nu_1$ of them are used for training the randomized neural networks (approximating the continuation value functions) and $\nu_2$ of them are used to evaluate the resulting price $U_0$ (using these trained networks).
Through the independence of the two parts of the dataset it is clear that any over-fitting bias the networks might learn will not lead to prices that are too large. This is important, since the quality of a method can therefore be evaluated by the resulting price, using a higher price indicating a a better result.

§.§.§ Transfer Learning in American Option Pricing

In order to test the transfer learning capabilities of our algorithm, we do not generate a very large dataset of the most plausible distribtution of $X$ as is standard practice. Instead we generate many smaller, more or less similar datasets $\mathcal{D}_1, \dotsc, \mathcal{D}_N$ out of which  $\mathcal{D}_1$ is considered our main dataset on which we evaluate the performance.

We always use the RLSM algorithm to price the American options and only vary the optimization method used internally by RLSM to solve the $M$ recursive regression problems arising at the exercise dates. Depending on the chosen optimisation method only one, multiple or all datasets are taken into account to compute the parameters of the randomized neural network at the given exercise date.
We compare our regret-optimal (RO) optimisation method of Theorem <ref> to the following baselines.

* Local optimizer (LO-$i$). The local optimizers $\theta^\star_i$ solves the individual optimisation problem on the datasets $\mathcal{D}_i$ without considering the other datasets.

* Mean Local optimizer (MLO). The mean of the local optimizers $ \bar{\theta}^\star = \frac{1}{N} \sum_{i=1}^N \theta^\star_i$ is the bagging aggregation.

* Joint optimizer (JO). The joint optimizer ${\theta}^\star_{\mathcal{N}}$ with $\mathcal{N}=\{ 1, \dotsc, N\}$ solves the optimisation problem of the pooled regime $\mathcal{D} = \cup_{i=1}^N \mathcal{D}_i$.

* Joint Subset optimizer. This is the joint optimizer for any subset $\mathcal{I} \subset \mathcal{N}$ of the datasets $\mathcal{D}_{\mathcal{I}} = \cup_{i \in \mathcal{I}} \mathcal{D}_i$.

In this application we are only interested in the terminal parameters $\Theta(T)$. We know that our algorithm yields the parameters optimising (<ref>) for any $T\geq 1$ since there exists a closed-form solution. Hence, we choose $T=1$ which is most specialized on minimizing (<ref>), while the influence of the running costs (<ref>) can be controlled via the choice of $\lambda$ and $\beta$. Additionally, this minimizes the computation time for our algorithm.

As fKRR we use the same randomized neural networks as in Section <ref>.
The random weights of the randomized neural networks, as well as the randomly sampled paths both introduce some randomness into the pricing problem. Therefore we follow <cit.> by always running the experiments $n_{runs} = 100$ times with different seeds and taking means over the prices. For easier comparison, we then consider the relative performance ($RP$) on our main dataset, computed as the mean price of any given method divided by the mean price of the local optimizer on the main dataset.
Additionally, we compute asymptotically valid (by the central limit theorem) $95\%$-confidence-intervals as $RP \pm z_{0.975}\hat\sigma/\sqrt{n_{runs}}$, where $z_{0.975}$ is the $0.975$-quantile of the standard normal distribution and $\hat\sigma$ is the standard deviation  of $RP$ (computed as the sample standard deviation of the prices over the runs divided by the mean price of the local optimizer on the main dataset).
To make the results comparable between the optimization methods, we use the same seeds for all of them, such that run $i$ always uses the exact same paths and the exact same random weights in the randomized neural network, independently of the chosen method.
Moreover, we always use $\nu_2 \geq 50,000$ evaluation paths, such that the noise of the Monte Carlo approximation of the expectations is relatively small when evaluating the prices, even if we use much smaller training set sizes $\nu_1$.

In the first experiment we study the ability of our method to transfer knowledge of those datasets, out a large number of small equally-sized datasets, benefiting the main task.

In this experiment we generate $13$ different datasets, all of a Heston model with $\nu_1 = 100$ training samples but with different parameters.
The Heston model is defined as

\begin{equation}\label{equ:heston model}
\begin{split}
d X_t &= (r - \delta) X_t dt + \sqrt{v_t} X_t d W_t, \\
d v_t &= -k (v_t - v_{\infty}) dt + \sigma \sqrt{v_t}d B_t,
\end{split}
\end{equation}

where $W$ and $B$ are two Brownian motions with correlation $\rho \in [-1,1]$.
The main dataset $\mathcal{D}_1$ has the parameters $r=0.05$, $\delta=0.1$, $k=2$, $v_{\infty}=0.01$, $\sigma=0.2$, $\rho=-0.3$ and starting values $X_0=100$ and $v_0=v_\infty$.
The other datasets are define as all combinations of $r \in \{0.05, 0.5\}$, $\sigma \in \{ 0.15, 0.2, 0.25 \}$, $v_{\infty} \in \{ 0.005, 0.015 \}$ and otherwise the same parameters as the main dataset.
For all datasets we use a max call option with strike $K=100$ on $d=2$ i.i.d. stocks of the given model with $M=9$ equidistant exercise dates until the maturity $T_{mat} = 3$. The discounting factor is $e^{-r\,T_{mat}\,t/M}$ for $0\leq t \leq M$.
Moreover, we use a randomized neural network with one hidden layer with $300$ nodes and all methods use the ridge coefficient $\kappa=2$ and, where applicable, regret coefficients $\lambda=2, \beta=1$.

We report the results in Table <ref>.
For our regret-optimal method (RO) we use information sharing levels $\eta \in \{  10, 100, 500\}$ and see that $\eta=100$ leads to the best result outperforming the local optimizer of the main dataset (LO-1) by $9\%$. It also outperforms all other local optimizers (LO-$i$), the mean local optimizer (MLO) and the joint optimizer (JO).
For $\eta=10$ we have similar performance as LO-1 and for $\eta=500$ we also have significant outperformance, though smaller than for $\eta=100$.

The 6 datasets with $r=0.05$ (LO-2, $\dotsc$, LO-7) are relatively similar to the main dataset (with $RP \geq 0.96$), while the other 6 (LO-8, $\dotsc$, LO-13) are quite different (with $RP < 0.73$).
Therefore, we claim that our algorithm should mainly transfer knowledge from those 6 datasets with $r=0.05$ while avoiding to learn something wrong from the others.
In line with this we provide results of the joint optimizer on the 7 datasets with $r=0.05$, which performs better than the other baselines and outperforms LO-1 by $6.2\%$. In particular, this confirms that the datasets with $r=0.05$ are those to transfer knowledge from, since the joint optimizer for all datasets performs worse than LO-1. The joint subset optimizer is significantly outperformed by our regret-optimal algorithm RO ($\eta=100$) by about $3\%$, which shows that our method constitutes a better way for transfer learning than simply selecting the datasets most similar to $\mathcal{D}_1$.

We additionally compare our method to the local optimizer for the main dataset which uses either $700$ (the number of training samples in the 7 similar datasets with $r=0.05$) or $50,000$ training samples. Importantly, both are just references that would not be available in a real world setting where the datasets are limited.
The first one only leads to a slight increase in $RP$ of about $1\%$, which shows that our method extracts nearly as much knowledge out of the (same amount of) samples of the similar datasets with $r=0.05$.
The second one (which is considered a good proxy for the true American option price on the main dataset) has about $10\%$ better performance. This also shows the limitations of our method, in particular, it can only transfer as much knowledge as available in the other datasets.

Relative performance (RP) and $95\%$-confidence-intervals for different optimization methods in Experiment <ref>. We compare the local optimizers on the different datasets (LO-$n$), with the mean local optimizer (MLO), the joint optimizer (JO) and our regret-optimal method (RO). The “oracle” local optimizer on the main dataset with additional training samples (standard: $100$ samples per dataset) is included.

When the information sharing parameter $\eta$ is set to $100$, the regret-optimal (RO) algorithm outperforms all “non-oracle” baselines.

method      $RP$      $95\%$-CI

[head to column names, late after line=

[head to column names, late after line=

In the second experiment we test the ability of our method to transfer knowledge from a small dataset which is similar to the small main dataset, even if it is confronted with a very different and much larger (dominating in terms of samples) dataset.

In this experiment we generate 3 different datasets. The main dataset is generated from a rough Heston model

\begin{align*}
dX_t &=  (r - \delta) X_t dt +  \sqrt{v_t} X_t dW_t, \\
v_t &= v_0 + \int_0^t \frac{(t-s)^{H-1/2} }{\Gamma(H + 1/2)} \kappa (v_\infty - v_s)ds + \int_0^t \frac{(t-s)^{H-1/2} }{\Gamma(H + 1/2)} \sigma \sqrt{v_s}dB_s,
\end{align*}

which is similar to the Heston model except that the volatility is  a rough process with Hurst parameter $H\in(0,1/2]$, where the case $H=1/2$ coincides with the standard Heston model.
We use $H=0.1$ and otherwise the same parameters as for the main dataset in Experiment <ref>. The second and third dataset generated from standard Heston model (i.e.  $H=1/2$), with the same parameters, except for the drift $r=0.5$ for the second dataset.
In particular, dataset 3 is the closest non-rough dataset to the main dataset, while dataset 2 is quite different.
For datasets 1 and 3 we use $\nu_1=100$ training samples, while we use $\nu_1=50,000$ for dataset 2, such that it dominates the others in its size.
All the other choices are the same as in Experiment <ref>.

The results are reported in Table <ref>. The baselines MLO and  JO perform significantly worse than the LO on the main dataset. As expected, the JO suffers much more from the dominating dataset 2 than the MLO. Indeed, the MLO would suffer more from a large number of datasets that are quite different from the main dataset, no matter their sample size.
In contrast to this, our regret optimal method with $\eta \in\{50,100\}$ outperforms the LO-1 by $4\%$, showing that also a quite different dataset dominating in sample size is not a problem.

Relative performance (RP) and $95\%$-confidence-intervals for different optimization methods in Experiment <ref>. We compare the local optimizers on the different datasets (LO-$n$), the mean local optimizer (MLO), the joint optimizer (JO) and our regret-optimal method (RO). With $\eta=100$ RO outperforms all baselines.

method      $RP$      $95\%$-CI

[head to column names, late after line=

Both of these experiments show that our regret optimal method is well suited for automatic (up to the hyper-parameter $\eta$) dataset selection in transfer learning tasks that have potentially large numbers of datasets and do not allow a (manual) pre-selection.
In particular, neither a large number of “bad” datasets (of which the MLO suffers) nor some “bad” datasets that dominate in sample size (of which the JO suffers) constitute a problem for our regret-optimal method.
Moreover, the sample efficiency when transferring knowledge from similar datasets to the main task is very high (outperforming the JO on the respective subset of similar datasets and nearly being on par with the LO on the main dataset using a larger training sample), as we saw in Experiment <ref>.

§ CONCLUSION

In this work we presented a regret-optimal algorithm for federated transfer learning based on finite-rank kernel ridge-regressors (fKRR).
Besides the theoretical properties of this algorithm, we also provided experiments demonstrating the transfer learning capabilities of our method in the context of American option pricing.

One possible future research direction is to apply our transfer learning method to real world datasets, where the similarity between the datasets is a priori not easy to quantify.
From the theoretical side we would like to extend our method to regression models that do not necessarily permit a closed-form solution.

AK acknowledges financial support from an NSERC Discovery Grant No. RGPIN-2023-04482 and their McMaster Startup Funds. MG acknowledges support from an NSERC Discovery Grant No. RGPIN-2013-276671.

§ PROOFS AND TECHNICAL RESULTS

We keep our notation light and use $u^{\top}v$ to implement the Euclidean inner-product between any two vectors $u,v\in \mathbb{R}^k$ for any $k\in \mathbb{N}_+$.
We use $| \cdot |$ to denote the Euclidean norm of a vector or Fröbenius norm of a matrix, whichever is applicable.
For two square matrices $A$ and $B$, $A\leq B$ ($A<B$, resp.) means that $B-A$ is positive semi-definite (positive definite, resp.); see <cit.> for details.

Lemmas <ref> and <ref> give some elementary properties of matrix algebra, which will be used repeatedly in the subsequent proofs.

The Euclidean norm $|\cdot|$ of matrices has the following properties:

i). For any matrices $A\in \mathbb{R}^{n_1 \times n_2}$ and $B\in\mathbb{R}^{n_2\times n_3}$, we have that
$|AB| \leq |A|\cdot |B|$;

ii). For any matrices $A$, $B\in \mathbb{R}^{n_1 \times n_2}$, we have $|A+B| \leq |A|+|B|$.

Denote $A=(A_{ij})_{1\leq i\leq n_1 , 1\leq j \leq n_2}$ and $B=(B_{jk})_{1\leq j \leq n_2, 1\leq k \leq n_3}$.

By the Cauchy-Schwarz inequality, we have for each column vector $B_{\cdot k}$ of $B$ that

\begin{align}
|A B_{\cdot k}|^2 = \sum_{i=1}^{n_1} \big( \sum_{j=1}^{n_2} A_{ij} B_{jk} \big)^2
\leq \sum_{i=1}^{n_1} \big( \sum_{j=1}^{n_2} A_{ij}^2 \sum_{j=1}^{n_2} B_{jk}^2 \big)
= \big( \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} A_{ij}^2 \big) \sum_{j=1}^{n_2} B_{jk}^2
= |A|^2 \sum_{j=1}^{n_2} B_{jk}^2 .
\notag
\end{align}

It then follows that

\begin{align}
|AB|^2 = & |(A B_{\cdot 1}, A B_{\cdot 2}, ..., A B_{\cdot n_3} ) |^2
= |A B_{\cdot 1}|^2 + |A B_{\cdot 2}|^2 + \cdots + |A B_{\cdot n_3}|^2 \notag \\
\leq & |A|^2 \sum_{j=1}^{n_2} \sum_{k=1}^{n_3} B_{jk}^2 = |A|^2 |B|^2 , \notag
\end{align}

and assertion i) is proved.

For $A$, $B\in \mathbb{R}^{n_1\times n_1}$, we have

\begin{align}
|A+B|^2 = & \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} (A_{ij} + B_{ij})^2
=   \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} A_{ij}^2 + \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} B_{ij}^2
+ 2 \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} A_{ij} B_{ij} \notag \\
\leq & |A|^2 + |B|^2 + 2 |A|^2 |B|^2
= (|A|+|B|)^2 , \notag
\end{align}

which proves assertion ii).

Let $w_1$, $w_2$, ..., $w_n$ be positive numbers such that $\sum_{k=1}^n w_k = 1$.
For any $n$ matrices $A_1$, ..., $A_n\in \mathbb{R}^{n_1 \times n_2}$, it satisfies that
$\big| \sum_{k=1}^n w_k A_k \big|^2 \leq \sum_{k=1}^n w_k |A_k|^2$;
particularly, when $w_1=w_2=\cdots = w_n=1/n$, we have
$\big|\sum_{k=1}^n  A_k \big|^2 \leq n \sum_{k=1}^n |A_k|^2$ .

Denote $A_k = ( A_{k, ij})_{1\leq i\leq n_1,  1\leq j \leq n_2}$, $k=1$, ..., $n$.
By the Jensen's inequality, we have

\begin{align}
( w_1 A_{1, ij} + \cdots + w_n A_{n, ij})^2 \leq  w_1 A_{1, ij}^2 + \cdots + w_n A_{n, ij}^2  . \notag
\end{align}

It then follows that

\begin{align}
|w_1 A_1 + \cdots + w_n A_n|^2 = & \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} (w_1 A_{1, ij} + \cdots + w_n A_{n, ij})^2
\leq  \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} w_1 A_{1, ij}^2 + \cdots + w_n A_{n, ij}^2   \notag \\
= &   w_1 \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} A_{1, ij}^2 + \cdots + w_n \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} A_{n, ij}^2  \notag \\
= &   w_1 |A_1|^2 + \cdots + w_n |A_n|^2  .
\notag
\end{align}

We begin with the relatively short proof of Theorem <ref> establishing the optimality of the weights obtained in Algorithm <ref>.

For $i=1,\dots,N$, consider the probability measure

\[
\bar{\mathbb{P}}^{\eta}
\eqdef
\sum_{i=1}^N\,
\frac{
[s^1+\eta -s^i]_+
\sum_{k=1}^N\, [s^1 + \eta - s^k]_+
\,
\delta_{\theta^\star_i}
\]

on the discrete space $\{\theta^\star_i\}_{i=1}^N$ topologized by the discrete metric.

Any minimizer of (<ref>) is a minimizer of the following problem, and vice-versa

\begin{align}
\nonumber
\min_{
\mathbb{P}_w \in \mathcal{P}(\{\theta^\star_i\}_{i=1}^N)
\frac1{|\mathcal{D}_1|}\sum_{j=1}^{|\mathcal{D}_1|}\,
\big\|
\big\|_2^2
\eta\,
% \frac{(1-\eta)}{\eta}
\operatorname{D}_{\operatorname{KL}}\big(
\mathbb{P}_w
\|
\bar{\mathbb{P}}^{\eta}
\big)
% +\iota(\mathbb{P}_w \ll \bar{P}^{\eta})
\\
\nonumber
& =
\min_{
\mathbb{P}_w \in \mathcal{P}(\{\theta^\star_i\}_{i=1}^N)
\mathbb{E}_{\theta \sim \mathbb{P}_w}\,
\Biggl[
\frac1{|\mathcal{D}_1|}\sum_{j=1}^{|\mathcal{D}_1|}\,
\big\|
\big\|_2^2
\Biggr]
% \frac{(1-\eta)}{\eta}
\eta\,
\operatorname{D}_{\operatorname{KL}}\big(
\mathbb{P}_w
\|
\bar{\mathbb{P}}^{\eta}
\big)
\iota(\mathbb{P}_w \ll \bar{P}^{\eta})
\\
\label{eq:optimizingQ__regularized__Renormalized}
& =
\min_{
\mathbb{P}_w \in \mathcal{P}(\{\theta^\star_i\}_{i=1}^N)
\mathbb{E}_{\theta \sim \mathbb{P}_w}\,
\biggl[
\sum_{i=1}^N\,
I_{\theta = \theta^\star_i}
\biggl]
% \frac{(1-\eta)}{\eta}
\eta\,
\operatorname{D}_{\operatorname{KL}}\big(
\mathbb{P}_w
\|
\bar{\mathbb{P}}^{\eta}
\big)
\iota(\mathbb{P}_w \ll \bar{P}^{\eta})
\end{align}

where the map $\iota(\mathbb{P}_w \ll \bar{P}^{\eta})$ is an “indicator function in the sense of convex analysis”; that is, $\iota(\mathbb{P}_w \ll \bar{P}^{\eta})=0$ if $\mathbb{P}_w$ is absolutely continuous with respect to $\bar{P}^{\eta}$ and $\infty$ otherwise.
We may now apply <cit.>[This result is an extension of <cit.> which allows us to vary the influence of the regulative-entropy term in (<ref>).] to (<ref>) deducing that it admits a minimizer $\mathbb{P}_{w^{\star}}=\sum_{i=1}^N\,w_i^{\star}\,\delta_{\theta^\star_i}$ and its Radon-Nikodym derivative with respect to the measure $\bar{\mathbb{P}}^{\eta}$ is

\begin{align}
\label{eq:optimizingQ__regularized__Renormalized____density_____BEGIN}
\frac{
\bar{\mathbb{P}}^{\eta}
= &
\frac{
\sum_{l=1}^N \, I_{\theta^\star_l=\theta^\star_i}
\mathbb{E}_{\bar{\mathbb{P}^{\eta}}}\biggl[
\exp\big(
\sum_{i=1}^N\,
I_{\theta = \theta^\star_i}
\big)
\biggr]
\\
\nonumber
= &
\frac{
\sum_{k=1}^N\,
\frac{e^{s_k/\eta}[s^1+\eta-s^k]_+}{
\sum_{m=1}^N\,
\\
\label{eq:optimizingQ__regularized__Renormalized____density_____END}
= &
\Big(
\sum_{m=1}^N\,
\Big)
\,
\frac{
\sum_{k=1}^N\,
\end{align}

for $i=1,\dots,N$.  Since $\mathbb{P}_{w^{\star}}$ is a minimizer of (<ref>) then we may deduce the expression for the optimal weights $w^{\star}$ by first noting that

\begin{align}
\label{eq:optimizingQ__regularized___weights____P1}
\mathbb{P}_{w^{\star}}
\eqdef
\sum_{i=1}^N\, w^{\star}_i\delta_{\theta^\star_i}\,
\sum_{i=1}^N\,
\frac{d\mathbb{P}_{w^{\star}}}{d\bar{\mathbb{P}^{\eta}}}(\theta^\star_i)
\,
\left(
\frac{
\sum_{j=1}^N\,
\right)
\,
\delta_{\theta^\star_i}\,
\end{align}

Therefore, (<ref>)- (<ref>) and (<ref>) imply that, for each $i=1,\dots,N$, we have

\begin{align}
\nonumber
% \label{eq:optimizingQ__regularized___weights}
= &
\frac{d\mathbb{P}_{w^{\star}}}{d\bar{\mathbb{P}^{\eta}}}(\theta^\star_i)
\,
\left(
\frac{
\sum_{j=1}^N\,
\right)
\,
\\
\nonumber
= &
\Big(
\sum_{m=1}^N\,
\Big)
\,
\frac{
\sum_{k=1}^N\,
\,
\left(
\frac{
\sum_{j=1}^N\,
\right)
\\
\nonumber
= &
\frac{
\,
\sum_{k=1}^N\,
\end{align}

Next we now prove our main result, namely Theorem <ref>. The result is first rewritten, in the optimal control-theoretic notation introduced in Section <ref>.

Suppose that Assumptions <ref> and <ref> hold.  Then, the optimal control problem
(<ref>) and (<ref>) admits a unique solution for $t=0, 1, ..., T-1$ such that

\begin{align}
\boldsymbol\alpha(t)
= & -[(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
\big[ (\lambda I_{Np} + P(t+1)) \Theta(t)
- \lambda  \Theta^\star + S(t+1)  \big] ,   \label{optimala}
\end{align}

and the resulting cost is

\begin{align}
L^\star(\boldsymbol\alpha; \mathcal{D})
= \Theta^{\star \top} P(0)
\Theta^\star + 2 S^\top(0) \Theta^\star + r(0) .
\label{optimalL}
\end{align}

The $P(t)$, $S(t)$, and $r(t)$ in (<ref>) and (<ref>) are determined by the backward equation system on $t=1,\dots,T$:

\begin{align}
& \begin{cases}
P(t)  =  \beta I_{Np} - \beta^2 [(\lambda+\beta) I_{Np} + P(t+1) ]^{-1} , \label{eqnP(t)} \\
P(T) =  [w^\star_1 I_p ,\cdots, w^\star_N I_p]^{\top}
\big( \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}   u^i_j u^{i\top}_j \big)
[ w^\star_1 I_p, \cdots, w^\star_N I_p] ,
\end{cases}   \\
& \begin{cases}
S(t) =  \beta \big[ (\lambda+\beta)I_{Np} + P(t+1) \big]^{-1} (S(t+1) - \lambda \Theta^\star )
\label{eqnS(t)} ,     \\
S(T) =  - [ w^\star_1 I_p, \cdots, w^\star_N I_p]^{\top} \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  {u^i_j} y^i_j ,
\end{cases} \\
r(t) =  - ( S(t+1) - \lambda \Theta^\star )^\top
[(\lambda+\beta)I_{Np} + P(t+1)]^{-1} \cdot ( S(t+1) - \lambda \Theta^\star )   \\
\hspace{1.1cm}   + \lambda \| \Theta^\star \|_2^2 + r(t+1),
\label{eqnr(t)} \\ r(T) = \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} |y^i_j|^2 .
\end{cases}
\end{align}

Let $V(t, \Theta)$ be the value function associated with the optimal control problem
(<ref>) and (<ref>).
By the dynamic programming principle, $V(t, \Theta)$  satisfies the equation

\begin{align}
& V (t, \Theta) = \min_{\boldsymbol\alpha}
\Big[   \lambda \left| \Theta + \boldsymbol\alpha - \Theta^\star \right|^2
+ \beta \left| \boldsymbol\alpha \right|^2
+ V (t+1, \Theta + \boldsymbol\alpha ) \Big] , \quad t=0, 1, ..., T-1 ,   \label{dpVw} \\
& V(T, \Theta )
= l (\theta^{w^\star}; \mathcal{D} ) .   \label{V(T,Theta)}
\end{align}

The loss function can be written as a linear quadratic function of
$\Theta=[\theta_1^\top, \cdots, \theta_N^\top]^\top$ such that

\begin{align}
&  l (\theta^{w^\star}; \mathcal{D} ) \notag \\
= &
\sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  \Big[ u^{i\top}_j \sum_{k=1}^N w^\star_k \theta_k - y^i_j  \Big]^2 \notag \\
= &  \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}
\Big[ \sum_{k=1}^N w_k^{\star 2} \theta_k^\top u^i_j u^{i\top}_j \theta_k
-2 \sum_{k=1}^N w_k^\star y^i_j u^{i \top}_j \theta_k + |y^i_j|^2 \Big]
\notag \\
= &  \Big(\sum_{k=1}^N w^\star_k \theta_k^\top \Big) \Big( \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j  \Big) \Big(\sum_{k=1}^N w^\star_k  \theta_k\Big)
- 2 \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} y^i_j u^{i\top}_j \sum_{k=1}^N w^\star_k \theta_k  + \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} | y^i_j |^2 \notag \\
= & \Theta^\top [w^\star_1 I_p ,\cdots, w^\star_N I_p]^\top
\Big( \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j \Big)
[ w^\star_1 I_p, \cdots, w^\star_N I_p ] \Theta \notag \\
& - 2 \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} y^i_j u^{i \top}_j
[ w^\star_1 I_p, \cdots, w^\star_N I_p ] \Theta
+ \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} | y^i_j |^2 . \label{lwT}
\end{align}

We claim that $V(t, \Theta)$ takes the linear quadratic form

\begin{align}
V (t, \Theta)
=  \Theta^\top P(t) \Theta + 2 S^\top (t) \Theta + r(t) . \label{ansatzVw}
\end{align}

To see this, note from (<ref>) and (<ref>) that
$V(T, \Theta)$ takes the linear quadratic form of (<ref>). By solving the dynamic programming equation (<ref>) at $t=T-1$ for the optimal $\boldsymbol\alpha$ and substituting it back into the equation, we obtain $V(T-1, \Theta)$ also taking the linear quadratic form of (<ref>). By backward induction we have that $V(t,\Theta)$ takes the form of (<ref>) for all $t=0,...,T$.

Here $r(t)$ is a constant term independent of the “state variable” $\Theta$, which we note exists because the cost $L^\star$ defined by (<ref>) contains terms independent of $\Theta(t)$.

Combining (<ref>),
(<ref>), and (<ref>)  verifies that $P(T)$ and $S(T)$ satisfy the terminal condition as given by (<ref>) and (<ref>).
We substitute (<ref>) into (<ref>) and reorganize the terms to get

\begin{align}
& \Theta^\top P(t) \Theta + 2 S^\top(t) \Theta + r(t)  \notag \\
= &  \min_{\boldsymbol\alpha} \Big\{  \boldsymbol\alpha^\top [(\lambda + \beta ) I_{Np} + P(t+1)] \boldsymbol\alpha
+ 2 [ \Theta^\top P(t+1)
+ S^\top (t+1) +\lambda (\Theta^\top - \Theta^{\star \top})]  \boldsymbol\alpha \notag \\
& \hspace{1cm} + \Theta^\top ( \lambda I_{Np} + P(t+1) )  \Theta
+ 2 ( S^\top(t+1) - \lambda \Theta^\star ) \Theta + \lambda |\Theta^\star|^2 + r(t+1) \Big\} .
\label{DPVansatz}
\end{align}

By the first order condition[NB, we can argue this way, since the problem is convexified upon fixing $w^{\star}$.  In fact, this is the control-theoretic motivation for using Algorithm <ref> to decouple the optimization of $w$ and of $\theta_{\cdot}$ from one another, and not treating them as a single (non-convex if left coupled) control problem.] for the right-hand side of (<ref>) with respect to $\boldsymbol\alpha$, we have

\begin{align}
0 = &   [(\lambda+\beta) I_{Np} + P(t+1)] \boldsymbol\alpha
+  [ ( \lambda I_{Np} + P(t+1) ) \Theta
+ S(t+1) - \lambda  \Theta^\star ]  .
\label{1stOrdera}
\end{align}

By Assumption <ref>, we have
$P(T)+(\lambda+1)I_{Np} >0$, and therefore from (<ref>) we obtain the optimal $\boldsymbol{\alpha}$ taking the form of (<ref>).

We further substitute $\boldsymbol{\alpha}$ of the form (<ref>) into the right side of (<ref>) to obtain

\begin{align}
& \Theta^\top  P(t) \Theta  + 2 S^\top(t) \Theta  + r(t) \notag \\
= & - \big[ \Theta^\top (\lambda I_{Np} + P(t+1))
+ S^\top(t+1) - \lambda \Theta^{\star \top}
\big] \cdot \big[(\lambda + \beta )I_{Np} + P(t+1)\big]^{-1} \cdot \notag \\
& \big[  (\lambda I_{Np} + P(t+1)) \Theta  + S(t+1)
-  \lambda  \Theta^\star  \big]
+ \Theta^\top (\lambda I_{Np} + P(t+1))
\Theta \notag \\
&  + 2 [ S^\top(t+1) - \lambda \Theta^{\star \top}]\Theta + \lambda |\Theta^\star|^2 + r(t+1) .
\label{eqn:Vansatz}
\end{align}

Reorganizing the terms on right side of (<ref>), we have

\begin{align}
& \Theta^\top  P(t) \Theta  + 2 S^\top(t) \Theta  + r(t) \notag \\
= & - \Theta^\top   (\lambda I_{Np} + P(t+1) ) [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
(\lambda I_{Np} + P(t+1))   \Theta \notag \\
& + \Theta^\top ( \lambda I_{Np} + P(t+1) ) \Theta + 2 [ S^\top(t+1) - \lambda \Theta^{\star \top}]\Theta  \notag \\
& - 2 [S^\top(t+1) - \lambda \Theta^{\star \top}]
[(\lambda +\beta )I_{Np} + P(t+1) ]^{-1} (\lambda I_{Np} + P(t+1)) \Theta \notag \\
& - (S^\top(t+1) - \lambda \Theta^{\star\top}) [(\lambda + \beta )I_{Np} + P(t+1)]^{-1} (S(t+1) - \lambda \Theta^\star)
\notag \\
& + \lambda |\Theta^\star|^2 + r(t+1) .
\label{eqn:Vansatz-simplified}
\end{align}

Matching the coefficients of the quadratic terms of $\Theta$ on both sides of (<ref>), we obtain

\begin{align}
= &  -  (\lambda I_{Np} + P(t+1) ) [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
(\lambda I_{Np} + P(t+1))   \notag \\
& +   \lambda I_{Np} + P(t+1) .
\label{eqn:P-unsimplified}
\end{align}

The right side of
(<ref>) further simplifies to

\begin{align}
& -  (\lambda I_{Np} + P(t+1) ) [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
(\lambda I_{Np} + P(t+1))
+  \lambda I_{Np} + P(t+1)   \notag \\
= &   -  [ ( \lambda + \beta ) I_{Np} + P(t+1) - \beta I_{Np} ]  [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
(\lambda I_{Np} + P(t+1)) \notag \\
& +  \lambda I_{Np} + P(t+1)  \notag \\
= & -   ( \lambda  I_{Np} + P(t+1) )  + \beta [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
(\lambda I_{Np} + P(t+1))
+  \lambda I_{Np} + P(t+1) \notag \\
= &
\beta [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
(\lambda I_{Np} + P(t+1))  \notag \\
= &
\beta [(\lambda + \beta )I_{Np} + P(t+1)]^{-1}
[ ( \lambda + \beta ) I_{Np} + P(t+1) - \beta I_{Np} ] \notag  \\
= & \beta I_{Np} - \beta^2 [(\lambda + 1)I_{Np} + P(t+1)]^{-1} .
\label{eqn:P-simplified}
\end{align}

Combining (<ref>)
and (<ref>) gives (<ref>) for $t=T-1$.

Matching the coefficients of the linear terms of $\Theta$ on both sides of (<ref>) and taking the transpose, we obtain

\begin{align}
S(t) = & - (\lambda I_{Np} + P(t+1) ) [(\lambda + \beta )I_{Np} + P(t+1) ]^{-1} ( S(t+1) - \lambda \Theta^\star ) \notag \\
& +  S(t+1) - \lambda \Theta^\star .
\label{eqn:S-unsimplified}
\end{align}

The right side of (<ref>) further simplifies to

\begin{align}
& - (\lambda I_{Np} + P(t+1) ) [(\lambda +  \beta )I_{Np} + P(t+1) ]^{-1} ( S(t+1) - \lambda \Theta^\star )
+ S(t+1) - \lambda \Theta^\star
\notag \\
= &   - [ (\lambda + \beta ) I_{Np} + P(t+1) - \beta I_{Np} ] [(\lambda + \beta )I_{Np} + P(t+1) ]^{-1} ( S(t+1) - \lambda \Theta^\star ) \notag \\
& + S(t+1) - \lambda \Theta^\star
\notag \\
= & - ( S(t+1) - \lambda \Theta^\star )
+ \beta [(\lambda +\beta )I_{Np} + P(t+1) ]^{-1} ( S(t+1) - \lambda \Theta^\star )
+ S(t+1) - \lambda \Theta^\star
\notag \\
= &  \beta [(\lambda + \beta )I_{Np} + P(t+1) ]^{-1} ( S(t+1) - \lambda \Theta^\star ) . \label{eqn:S-simplified}
\end{align}

Combining (<ref>) and (<ref>) gives (<ref>) for $t=T-1$.

Matching the terms independent of $\Theta$ on both sides of (<ref>) gives

\begin{align}
r(t) = &  - ( S(t+1) - \lambda \Theta^\star )^\top
[(\lambda+ \beta )I_{Np} + P(t+1)]^{-1} \cdot ( S(t+1) - \lambda \Theta^\star ) \notag    \\
&  + \lambda |\Theta^\star|^2 + r(t+1) , \notag
\end{align}

which is (<ref>) for $t=T-1$.

$P(T) + (\lambda + \beta )I_{Np} > \beta I_{Np}$,
we have

\begin{align}
P(T-1) = \beta I_{Np}
- \beta^2 [(\lambda + \beta )I_{Np} + P(T)]^{-1} > 0 . \notag
\end{align}

By backward induction on $t$, we can show that
$\boldsymbol{\alpha}$ takes the form of (<ref>) for all $t=0, 1, ..., T-1$, and $P$ and $S$ satisfy the system (<ref>) and (<ref>).

Under Assumptions <ref> and <ref>, the solution $P$ and $S$ of the system (<ref>)-(<ref>) satisfies that for all $t=1$, $2$, ..., $T$,

\begin{align}
& | [(\lambda + \beta) I_{Np} + P(t)]^{-1} |\leq  (\lambda + \beta)^{-1},  \label{M(P)invbd}  \\
& |S(t)| \leq C \cdot N
\label{|S|bd} ,
\end{align}

where $C> 0$ is a constant depending on $\lambda$, $\beta$, $(w_1^\star, ..., w_N^\star)$, $\Theta^\star$, $K_x$, $K_y$ and $T$.

Let $\| \cdot \|$ be the spectral norm of a matrix such that $\|A\|\eqdef \sup_{x\neq 0} |Ax|/|x| $.
Since $P(t)$ is positive semi-definite, we have $(\lambda + \beta) I_{Np} + P(t) \geq (\lambda + \beta) I_{Np}$ and

\begin{align}
0\leq [(\lambda + \beta) I_{Np} + P(t) ]^{-1} \leq (\lambda + \beta)^{-1} I_{Np}   ,  \notag
\end{align}

which, by the Courant-Fisher Theorem (see <cit.>), implies that $\|  [(\lambda + \beta) I_{Np} + P(t) ]^{-1} \| \leq \| (\lambda+\beta)^{-1} I_{Np} \| = (\lambda + \beta)^{-1}$.
Since all norms of finite dimensional normed spaces are equivalent, and in particular, $|\cdot|\le \|\cdot\|$ on the space of $N_p\times N_p$-matrices, then we find that

\begin{align}
\big| [(\lambda + \beta) I_{Np} + P(t) ]^{-1} \big| \leq \|  [(\lambda + \beta) I_{Np} + P(t) ]^{-1} \| , \notag
\end{align}

and thus $\big| [(\lambda + \beta) I_{Np} + P(t) ]^{-1} \big| \leq (\lambda + \beta)^{-1}$,
which proves (<ref>).

Next, we prove (<ref>).
By Lemmas <ref> and <ref>, and (<ref>), we obtain

\begin{align}
| S(T) |^2
= & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\Big| \sum_{i=1}^N w^\star_i   \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j y^i_j  \Big|^2
\leq  \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w_i^\star  \Big| \sum_{j=1}^{|\mathcal{D}_i|}   u^i_j y^i_j \Big|^2  \notag  \\
\leq &  \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w_i^\star |\mathcal{D}_i|  \sum_{j=1}^{|\mathcal{D}_i|}   | u^i_j y^i_j |^2
\leq  \Big( \sum_{i=1}^N w^{\star 2}_i \Big) \sum_{i=1}^N w_i^\star |\mathcal{D}_i|^2  K_x^2 K_y^2  \notag \\
\leq &  \Big( \sum_{i=1}^N w^{\star 2}_i \Big) \sum_{i=1}^N w_i^\star N^2  K_x^2 K_y^2 ,
\label{eq:TBBounded}
\end{align}

and therefore the estimate (<ref>) holds for $t=T$.

Assume by induction that (<ref>) holds for $t=u$.
By (<ref>), Lemma <ref>, and the triangle inequality, we have

\begin{align}
|S(u-1)| \leq &  \beta \big| \big[ (\lambda+\beta)I_{Np} + P(u) \big]^{-1} \big| \cdot \big| S(u) - \lambda \Theta^\star  \big| \notag \\
\leq & \beta \big| \big[ (\lambda+\beta)I_{Np} + P(u) \big]^{-1} \big| \cdot ( |S(u)| + \lambda | \Theta^\star | ) .
\notag
\end{align}

The estimate (<ref>) for $t=u-1$ then follows from (<ref>) and the induction hypothesis.
Thus, by induction we have shown that (<ref>) holds for all $t=1$, $2$, ..., $T$.

We proceed to prove Theorem <ref> about the complexity of the regret-optimal algorithm.
To facilitate the proof, we introduce the following Lemma <ref> that summarizes complexities of matrix operations.

Suppose all elementary arithmetic operations (e.g. addition, subtraction, multiplication) of two real numbers are of constant (i.e. $\mathcal{O}(1)$) complexity. Then elementary matrix computations have the following complexities:

(i) Adding two $n_1 \times n_2$ matrices has a complexity of $O(n_1 n_2)$.

(ii) Multiplying an $n_1 \times n_2$ matrices by a scalar has a complexity of $O(n_1 n_2)$.

(iii) Multiplying an $n_1 \times n_2$ matrix with an $n_2 \times n_3$ matrix has a complexity of $O(n_1 n_2 n_3)$.
In particular, multiplying two $n\times n$ matrices has a complexity of $\mathcal{O}(n^{2.373})$.

Adding two $n_1\times n_2$ matrices involves adding each entry of one matrix with the corresponding entry of the other matrix, and each entry addition has a complexity of $\mathcal{O}(1)$. Hence the total complexity is $\mathcal{O}(n_1 n_2)$.

Multiplying an $n_1 \times n_2$ matrix by a scalar involves multiplying each entry by the scalar, and the scalar multiplication of each entry has a complexity of $\mathcal{O}(1)$. Hence the total complexity is $O(n_1 n_2)$.

Multiplication of an $n_1 \times n_2$ matrix and an $n_2 \times n_3$ matrix requires multiplying the $n_1$ row vectors of the first matrix by the $n_3$ column vectors of the second matrix, which needs a total of $n_1 n_3$ multiplications. Multiplication of two $n_2$-dimensional vectors has a complexity of $O(n_2)$. Therefore the total complexity is $O(n_1n_2n_3)$.
When multiplying two matrices of $n\times n$-dimension with an optimized CW-like algorithm the complexity is $\mathcal{O}(n^{2.373})$.

Complexity of Computing $P(T)$ is $\mathcal{O}\big(N^2p^3\big)$

Since $u^i_j \in \mathbb{R}^{p\times 1}$, computing each
$u^i_j u^{i\top}_j$ has $\mathcal{O}(p^2)$ complexity, by Lemma <ref>-(iii).
Multiplying the $p\times p$ matrix $u^i_j u^{i\top}_j$ by the scalar $w^\star_i$ has a complexity of $\mathcal{O}(p^2)$, by Lemma <ref>-(iii).

Then computing the $p\times p$ matrix $\sum_{i=1}^N w^\star_i \sum_{j=1}^{S_i}\, u^i_j u^{i\top}_j$ by adding the $\bar{N}$ matrices of dimension $p\times p$ has complexity $\mathcal{O}(\bar{N}\,p^2)$, by Lemma <ref>-(i).
Since $[w_1^\star I_p ,\cdots, w_N^\star I_p]^{\top}$ is $Np\times p$-dimensional, then computing the $np\times p$-dimensional product $[w_1^\star I_p ,\cdots, w_N^\star I_p]^{\top}\sum_{i=1}^N w^\star_i \sum_{j=1}^{S_i}\, u^i_j u^{i \top}_j$ has $\mathcal{O}(Np^3)$ complexity, by Lemma <ref>-(iii).

Likewise, the product

\[
\Big([w_1^\star I_p ,\cdots, w_N^\star I_p]^{\top}\sum_{i=1}^N  w^\star_i \sum_{j=1}^{S_i}\, u^i_j u^{i\top}_j \Big)[w_1^\star I_p ,\cdots, w_N^\star I_p]
\]

has a complexity of $\mathcal{O}(N^2p^3)$.  Therefore, computing $P(T)$ itself is of $\mathcal{O}\big(\bar{N}p^2+N^2p^3\big)$ complexity.
Assumption (<ref>) implies that $\bar{N}\leq N^2$; hence, the complexity of computing $P(T)$ is $\mathcal{O}\big(N^2p^3\big)$.

Complexity of Computing $S(T)$ is $\mathcal{O}\big(Np\max\{N,p\}\big)$

Since $u^i_j\in\mathbb{R}^{p\times 1}$ and $y^i_j\in \mathbb{R}$,
computing the product $u^i_j y_j^i$ has complexity $\mathcal{O}(p)$, by Lemma <ref>-(iii).
And multiplying the $p\times 1$ dimensional vector $u^i_j y^i_j$ by the scalar $w^\star_i$ has complexity $\mathcal{O}(p)$, by Lemma <ref>-(ii). Then the complexity of computing $\sum_{i=1}^N w^\star_i \sum_{j=1}^{S_i}\,u^i_j y_j^i$ by adding $\bar{N}$ vectors of dimension $p$ is $\mathcal{O}\big(\bar{N}p\big)$, by Lemma <ref>-(i).
Since $- [ w_1^\star I_p, \cdots, w_N^\star I_p]^{\top}$ is an $Np\times p$-dimensional matrix and $\sum_{i=1}^N w^\star_i \sum_{j=1}^{S_i}\, u^i_j y_j^i$ is a $p\times 1$-dimensional vector, then computing their $Np\times 1$ dimensional product has $\mathcal{O}(Np^2)$ complexity, according to Lemma <ref>-(iii).
Therefore, the computing $S(T)$ has a complexity of $\mathcal{O}(\bar{N}p + N P^2)$.
Since (<ref>) implies that $\bar{N}\le N^2$ then, computing $S(T)$ had a computational cost of $\mathcal{O}\big(Np\max\{N,p\}\big)$.

Complexity of Computing $P(t)$ and $S(t)$ for $t=T-1,\dots, 1$ is $\mathcal{O}\big(T(Np)^{2.373}\big)$

The complexity of computing $(\lambda +\beta )I_{Np} + P(t+1)$ for each $t$ is $\mathcal{O}(Np)$, since we only need to add scalars on the matrix's diagonal.  By <cit.>, computing the matrix inverse $[(\lambda+\beta)I_{Np}+P(t+1)]^{-1}$ with a CW-like algorithm has a complexity of $\mathcal{O}\big((Np)^{2.373}\big)$.
Since the addition $\beta I_{Np} - \beta^2 [(\lambda+ \beta )I_{Np} + P(t+1) ]^{-1}$ only requires us to multiply all entries of $ [(\lambda+\beta)I_{Np} + P(t+1) ]^{-1}$ by $-\beta^2$ and add along the diagonal, then it has a complexity of $\mathcal{O}\big((Np)^2 + Np\big)$, by Lemma <ref>-(i) and (ii).
Performing the matrix addition $S(t+1)-\lambda \Theta^\star$ has complexity $\mathcal{O}(Np)$ by Lemma <ref>-(i).
Computing the $Np\times 1$-dimensional product
$[(\lambda+1)I_{Np} + P(t+1)]^{-1} (S(t+1) - \lambda \Theta^\star)$
by multiplying an $Np\times Np$ matrix by an $Np\times 1$ matrix has complexity $\mathcal{O}(N^2 p^2)$, according to Lemma <ref>-(iii).
Therefore, computing each $P(t)$ and $S(t)$ for $t=T-1,\dots,1$, has complexity  $\mathcal{O}\big((Np)^{2.373}\big)$.  Since there are $T-1$ such matrices to compute, then the total cost of computing every $P(T-1),\dots,P(1)$ and $S(T-1)$, ..., $S(1)$ is $\mathcal{O}\big(T(Np)^{2.373}\big)$.

The Cost of Computing Each $\Theta(t)$ is $\mathcal{O}(T\,N^2p^2)$

To compute $\Theta(t)$ for all $t=1$, $\dots$, $T$, starting from $\Theta(0)=\Theta^\star$,
we need to compute $\mathbf{a}(t)$ given by (<ref>) for $t=0$, $1$, $\dots$, $T-1$.
Hence we need to quantify the complexity of computing $\mathbf{a}(t)$.

By (<ref>), we can obtain $-[(\lambda + \beta)I_{Np} + P(t+1)]^{-1}$ by computing
$\beta^{-2}(P(t)- \beta I_{Np})$, which involves adding a scalar to the diagonal of an $Np\times Np$ matrix, and then multiplying the matrix by a scalar.
The cost of computing the product
$(\lambda I_{Np} + P(t+1)) \Theta(t)$ is $\mathcal{O}(N^2 p^2)$ and the cost of computing all sums in $(\lambda I_{Np} + P(t+1)) \Theta(t)
- \lambda \Theta^\star + S(t+1)  $ are of lower order complexity.  Finally, the cost of computing the matrix-product between $-[(\lambda + \beta)I_{Np} + P(t+1)]^{-1} $ and $
\big[ (\lambda I_{Np} + P(t+1)) \Theta(t)
- \lambda \Theta^\star + S(t+1)  \big]$ is $\mathcal{O}(N^2p^2)$, by Lemma <ref>-(iii). Therefore, the cost of computing each $\mathbf{a}(t)$ is $\mathcal{O}(N^2p^2)$, and therefore the cost of computing $\Theta(t)$ for all $t=1$, ..., $t=T$ from $\Theta(0)=\Theta^\star$ is $\mathcal{O}(T N^2 p^2)$.

Tallying up Complexities: The Complexity of Computing $(\Theta(t))_{t=0}^T$ is
$\mathcal{O}( N^2 p^3 + T\, (N p)^{2.373})$.

Looking over all computations involved, we deduce that the cost of computing the sequence $(\Theta(t))_{t=1}^T$ is (in order) dominated by the computational complexity of computing the sequence $(P(t))_{t=1}^{T-1}$ and $P(T)$; whence it is
$\mathcal{O}( N^2 p^3 + T\, (N p)^{2.373})$.

The following Lemma <ref>,
Proposition <ref>, and Corollary <ref> are devoted to proving Theorem <ref>, the adversarial robustness of Algorithm <ref>.

We introduce the linear map $M(\cdot): \mathbb{R}^{Np\times Np} \to \mathbb{R}^{Np\times Np}$ such that $M(Q) = (\lambda + \beta) I_{Np} + Q$ for any $Q\in \mathbb{R}^{Np\times Np}$.

Given any two data sets $\mathcal{D}=\cup_{i=1}^N \mathcal{D}_i=\cup_{i=1}^N \cup_{j=1}^{|\mathcal{D}_i|}
\{(x^i_j, y^i_j) \}$ and
$\widetilde{\mathcal{D}}=\cup_{i=1}^N \widetilde{\mathcal{D}}_i = \cup_{i=1}^N \cup_{j=1}^{|\mathcal{D}_i|}  \{(\tilde{x}^i_j, \tilde{y}^i_j) \}$ satisfying Assumption <ref>, let $(P, S)$ and $(\widetilde{P}, \widetilde{S})$ be the solutions of (<ref>)-(<ref>) corresponding to $\mathcal{D}$ and $\widetilde{\mathcal{D}}$, respectively. Then, for all $t=1, ..., T$, we have

\begin{align}
\big| P(t) - \widetilde{P}(t) \big|
\le &
C \cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2 \Big]^{1/2} ,
\label{eqn:estP-tildeP} 	 \\
\big| [ M(P(t)) ]^{-1} - [ M(\widetilde{P}(t)) ]^{-1} \big|
\le&
C \cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2 \Big]^{1/2} ,
\label{eqn:estM-tildeM} \\
|S(t) - \widetilde{S}(t)|
\le &
C \cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i|  \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j |^2 \Big]^{1/2} ,
\label{eqn:estS-tildeS}
\end{align}

where $C> 0$ is a constant depending on $\lambda$, $\beta$, $N$, $\bar{N}$, $(w_1^\star, ..., w_N^\star)$, $\Theta^\star$, $K_x$, $K_y$ and $T$.

We let $C$ be a constant depending on $\lambda$, $\beta$, $N$, $\bar{N}$, $\Theta^\star$, $(w_1^\star, ..., w_N^\star)$, $K_x$, $K_y$, and $T$,
and allow $C$ to vary from place to place throughout the proof.
We employ backward induction starting from $t=T$ to prove (<ref>), (<ref>), and

From (<ref>), the terminal condition $| P(T) - \widetilde{P}(T) |^2$ can be written as

\begin{align}
\big| P(T) - \widetilde{P}(T) \big|^2
= &   \Big| [w_1^\star I_p ,\cdots, w_N^\star I_p ]^\top
\Big[ \sum_{i=1}^N w_i^\star  \sum_{j=1}^{|\mathcal{D}_i|}  ( u^i_j u^{i\top}_j  - \tilde{u}^i_j \tilde{u}^{i\top}_j  ) \Big]
[ w_1^\star I_p, \cdots,
w_N^\star I_p ]  \Big|^2 \notag \\
= & \sum_{k, l=1}^N
w_k^{\star 2}
w_l^{\star 2}
\Big| \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  ( u^i_j u^{i\top}_j  - \tilde{u}^i_j \tilde{u}^{i\top}_j  ) \Big|^2 . \notag
\end{align}

By Lemmas <ref> and <ref>, we have

\begin{align}
\big| P(T) - \widetilde{P}(T)  \big|^2
\leq &   \sum_{k, l=1}^N w_k^{\star 2} w_l^{\star 2}
\sum_{i=1}^N w_i^\star
\Big| \sum_{j=1}^{|\mathcal{D}_i|}
(  u^i_j u^{i\top}_j  - \tilde{u}^i_j \tilde{u}^{i\top}_j  ) \Big|^2
\notag \\
\leq &   \sum_{k, l=1}^N w_k^{\star 2} w_l^{\star 2}
\sum_{i=1}^N w_i^\star |\mathcal{D}_i|  \sum_{j=1}^{|\mathcal{D}_i|}
\big| u^i_j u^{i\top}_j  -
\tilde{u}^i_j \tilde{u}^{i\top}_j  \big|^2
\notag \\
= &   \sum_{k, l=1}^N w_k^{\star 2} w_l^{\star 2}
\sum_{i=1}^N
w_i^\star |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|}
\big| ( u^i_j - \tilde{u}^i_j  )u^{i\top}_j
+ \tilde{u}^i_j (u^{i\top}_j - \tilde{u}^{i\top}_j)   \big|^2  \notag \\
\leq &   \sum_{k, l=1}^N w_k^{\star 2} w_l^{\star 2}
\sum_{i=1}^N w_i^\star |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|}
2 \big( \big|  u^i_j - \tilde{u}^i_j   \big|^2 \cdot \big| u^i_j \big|^2
+ \big| \tilde{u}^i_j \big|^2 \cdot \big| u^i_j - \tilde{u}^i_j  \big|^2 \big)
\notag \\
\leq &   \sum_{k, l=1}^N w_k^{\star 2} w_l^{\star 2}
\sum_{i=1}^N w_i^\star |\mathcal{D}_i|  \sum_{j=1}^{|\mathcal{D}_i|} 4 K_x^2
\big| u^i_j - \tilde{u}^i_j   \big|^2
. \label{estP(T)-tildeP(T)}
\end{align}

Therefore (<ref>) holds for $t=T$.

Assume by induction that
(<ref>) holds for $t=u$. We show that (<ref>) also holds for $t=u-1$.
From (<ref>), we have

\begin{align}
P(u-1) - \widetilde{P}(u-1)
= & [ M( \widetilde{P}(u) ) ]^{-1}
- [ M(P(u)) ]^{-1}\notag  \\
= & [ M( P(u) ) ]^{-1} ( P(u) - \widetilde{P}(u) )
[ M( \widetilde{P}(u) ) ]^{-1} ,  \notag
\end{align}

and by Lemma <ref> we further have

\begin{align}
\big| P(u-1) - \widetilde{P}(u-1)  \big|
\leq & \big| [ M( P(u) ) ]^{-1}
\big| \cdot \big| P(u) - \widetilde{P}(u) \big| \cdot
\big| [ M( \widetilde{P}(u) ) ]^{-1} \big| . \label{|tildeP-P|recursion}
\end{align}

By Corollary <ref> and the induction hypothesis for $t=u$,
(<ref>) implies that (<ref>) holds for $t=u-1$.
We have shown by induction that (<ref>) holds for all $t=1$, ..., $T$.

\begin{align}
[M(P(t))]^{-1} - [M(\widetilde{P}(t))]^{-1}
=  [ M( P(t) ) ]^{-1} (  \widetilde{P}(t) - P(t) )
[ M( \widetilde{P}(t) ) ]^{-1} ,
\notag
\end{align}

by Lemma <ref> we have

\begin{align}
\big| [M(P(t))]^{-1} - [M(\widetilde{P}(t))]^{-1} \big|
\leq  \big| [ M( P(t) ) ]^{-1} \big| \cdot
\big|  \widetilde{P}(t) - P(t) \big| \cdot
\big| [ M( \widetilde{P}(t) ) ]^{-1} \big|  .
\label{eqn:M-tildeM<=P-tildeP}
\end{align}

The estimate (<ref>) then follows from (<ref>), (<ref>), and Corollary <ref>.

Now we establish the estimate (<ref>).
By (<ref>), the terminal condition $S(T)-\widetilde{S}(T)$ satisfies

\begin{align}
| S(T) - \widetilde{S}(T) |^2
= &
\Big| [w_1^\star, \cdots, w_N^\star]^\top \sum_{i=1}^N w_i^\star \sum_{j=1}^{|\mathcal{D}_i|}
( u^i_j y^i_j - \tilde{u}^i_j \tilde{y}^i_j ) \Big|^2 \notag \\
=& \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\Big| \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} ( u^i_j y^i_j - \tilde{u}^i_j \tilde{y}^i_j ) \Big|^2 .  \notag
\end{align}

By Lemmas <ref> and <ref>, we obtain

\begin{align}
| S(T) - \widetilde{S}(T) |^2
\leq & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w^\star_i  \Big| \sum_{j=1}^{|\mathcal{D}_i|} ( u^i_j y^i_j - \tilde{u}^i_j \tilde{y}^i_j ) \Big|^2  \notag  \\
\leq & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|}  \big|  u^i_j y^i_j - \tilde{u}^i_j \tilde{y}^i_j \big|^2  \notag   \\
= & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|}  \big|  (u^i_j - \tilde{u}^i_j) y^i_j + \tilde{u}^i_j (y^i_j - \tilde{y}^i_j) \big|^2  \notag \\
\leq & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w^\star_i |\mathcal{D}_i| \cdot 2  \sum_{j=1}^{|\mathcal{D}_i|} \big[ \big|  (u^i_j - \tilde{u}^i_j) y^i_j \big|^2
+ \big|  \tilde{u}^i_j (y^i_j - \tilde{y}^i_j) \big|^2 \big]   \notag \\
\leq & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w^\star_i |\mathcal{D}_i| \cdot 2  \sum_{j=1}^{|\mathcal{D}_i|} \big[ |  u^i_j - \tilde{u}^i_j |^2 \cdot | y^i_j |^2
+ |  \tilde{u}^i_j |^2 \cdot | y^i_j - \tilde{y}^i_j  |^2 \big]   \notag\\ \leq & \Big( \sum_{i=1}^N w^{\star 2}_i \Big)
\sum_{i=1}^N w^\star_i |\mathcal{D}_i| \cdot 2  \sum_{j=1}^{|\mathcal{D}_i|} \big[ K_y^2 \cdot |  u^i_j - \tilde{u}^i_j |^2
+ K_x^2 \cdot | y^i_j - \tilde{y}^i_j  |^2 \big]   \notag
\\
\label{estST-tildeST}
\leq & \Big( \sum_{i=1}^N w^{\star 2}_i \Big) 2 (K_x^2 + K_y^2)
\sum_{i=1}^N w^\star_i |\mathcal{D}_i|   \sum_{j=1}^{|\mathcal{D}_i|} \big[ |  u^i_j - \tilde{u}^i_j |^2
+  | y^i_j - \tilde{y}^i_j  |^2 \big] ,
\end{align}

and thus (<ref>) holds for $t=T$.
Assume by induction that (<ref>) holds for $t=u$, we show that (<ref>) also holds for $t=u-1$.
From (<ref>), the difference $S(t)-\widetilde{S}(t)$ can be written as

\begin{align}
S(u-1) - \widetilde{S}(u-1) = & [ M(\widetilde{P}(u)) ]^{-1}[S(u) - \widetilde{S}(u) ]  \notag \\
& + \big[ ( M(P(u)) )^{-1} -
( M(\widetilde{P}(u)))^{-1} \big]
( S(u) - \lambda \Theta^\star ) ,
\notag
\end{align}

and satisfies by Lemma <ref> and Corollary <ref> that

\begin{align}
| S(u-1) - \widetilde{S}(u-1) | \leq &
\big| [ M(\widetilde{P}(u)) ]^{-1} \big| \cdot |S(u) - \widetilde{S}(u) |  \notag \\
& + \big| ( M(P(u)) )^{-1} -
( M(\widetilde{P}(u)))^{-1} \big| \cdot
( |S(u)| + \lambda | \Theta^\star | ) .
\label{|tildeS-S|recursion}
\end{align}

By Corollary <ref>, (<ref>), and the induction hypothesis,
(<ref>) implies that
(<ref>) holds for $t=u-1$.
We have shown by induction that (<ref>) holds for all $t=1$, ..., $T$.

For two arbitrary data sets $\mathcal{D}=\cup_{i=1}^N \mathcal{D}_i=\cup_{i=1}^N \cup_{j=1}^{|\mathcal{D}_i|}
\{(x^i_j, y^i_j) \}$ and
$\widetilde{\mathcal{D}}=\cup_{i=1}^N \widetilde{\mathcal{D}}_i = \cup_{i=1}^N \cup_{j=1}^{|\mathcal{D}_i|}  \{(\tilde{x}^i_j, \tilde{y}^i_j) \}$ satisfying Assumption <ref>,
let $\boldsymbol\alpha =(\alpha_1^{\top}, \cdots, \alpha_N^{\top})^{\top}$ and
$\widetilde{\boldsymbol\alpha}=(\widetilde{\alpha}_1^{\top}, \cdots, \widetilde{\alpha}_N^{\top})^{\top}$ be the optimal controls for (<ref>) and (<ref>) corresponding to $\mathcal{D}$ and $\widetilde{\mathcal{D}}$, respectively. Then, we have

\begin{align}
\sum_{t=0}^{T-1} |\boldsymbol\alpha(t) - \widetilde{\boldsymbol\alpha}(t) |
\le C  \Big[ \sum_{i=1}^N w^{\star 2}_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j|^2  \Big]^{1/2} ,
\label{eqn:stabilityControl}
\end{align}

where $C > 0$ is a constant depending on $\lambda$, $\beta$, $N$, $\bar{N}$, $\Theta^\star$,
$(w^\star_1, \cdots, w^\star_N)$, $K_x$, $K_y$, and $T$.

Let $C$ be a constant depending on $\lambda$, $\beta$, $N$, $\bar{N}$, $\Theta^\star$,
$(w^\star_1, \cdots, w^\star_N)$, $K_x$, $K_y$, and $T$,
and $C$ is allowed to vary from place to place throughout the proof.
By (<ref>) and (<ref>), the optimal control $\boldsymbol\alpha(t)$ along the optimal trajectory can be written as

\begin{align}
\boldsymbol\alpha(t)
= & [ (M(P(t+1)))^{-1} - I_{Np} ]  \Big[ \Theta^\star
+ \sum_{u=0}^{t-1} \boldsymbol\alpha(u) \Big] \notag \\
& +  [ M(P(t+1)) ]^{-1} ( \lambda  \Theta^\star - S(t+1) ) ,   \notag
\end{align}

and the difference $\boldsymbol{\alpha}(t)-\widetilde{\boldsymbol\alpha}(t)$ can be written as

\begin{align}
\boldsymbol\alpha(t) - \widetilde{\boldsymbol\alpha}(t)
= &  \big[(M(P(t+1)))^{-1} - (M(\widetilde{P}(t+1)) )^{-1} \big]
\Big[ \Theta^\star + \sum_{u=0}^{t-1}  \boldsymbol\alpha(u)   + \lambda \Theta^\star - S(t+1) \Big]  \notag \\
& + \big[ ( M(\widetilde{P}(t+1)) )^{-1} - I_{Np} \big]
\sum_{u=0}^{t-1} (\boldsymbol\alpha(u) - \widetilde{\boldsymbol\alpha}(u)) \notag \\
& + (M(\widetilde{P}(t+1)))^{-1} (S(t+1) - \widetilde{S}(t+1)) .
\label{a-tildea}
\end{align}

For $t=0$, the difference $\boldsymbol\alpha(0)-\widetilde{\boldsymbol\alpha}(0)$ can be written as

\begin{align}
\boldsymbol\alpha(0) - \widetilde{\boldsymbol\alpha}(0)
= & \big[ (M(P(1)))^{-1} - (M(\widetilde{P}(1)))^{-1}  \big]
\big[ (\lambda + 1)   \Theta^\star - S(1)
\big] \notag \\
& - [ M(\widetilde{P}(1)) ]^{-1} (S(1)-\widetilde{S}(1)) , \notag
\end{align}

and satisfies by Lemma <ref> that

\begin{align}
\big| \boldsymbol\alpha(0) - \widetilde{\boldsymbol\alpha}(0) \big|
\leq &   \big| (M(P(1)))^{-1} - (M(\widetilde{P}(1)))^{-1}  \big| \cdot
\big[ (\lambda + 1)   | \Theta^\star| +  | S(1) |
\big] \notag \\
& +  \big| [ M(\widetilde{P}(1)) ]^{-1} \big| \cdot  \big| S(1)-\widetilde{S}(1) \big| .  \notag
\end{align}

It then follows from Corollary <ref>, (<ref>), and (<ref>) that

\begin{align}
\big| \boldsymbol\alpha(0) - \widetilde{\boldsymbol\alpha}(0) \big| \leq C
\cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j
- \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j|^2  \Big]^{1/2} . \notag
\end{align}

Assume by induction that

\begin{align}
\sum_{u=0}^{t-1} \big| \boldsymbol\alpha(u) - \widetilde{\boldsymbol\alpha}(u) \big| \leq C
\cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j |^2  \Big]^{1/2} .  \notag
\end{align}

By Lemma <ref>, (<ref>) implies that

\begin{align}
& | \boldsymbol\alpha(t) - \widetilde{\boldsymbol\alpha}(t) | \notag \\
\leq &  \big| (M(P(t+1)))^{-1} - (M(\widetilde{P}(t+1)) )^{-1} \big| \cdot
\Big[ |\Theta^\star| + \sum_{u=0}^{t-1}  | \boldsymbol\alpha(u) |
+ \lambda | \Theta^\star | +| S(t+1) | \Big]  \notag \\
& + \big| ( M(\widetilde{P}(t+1)) )^{-1} - I_{Np}  \big|
\sum_{u=0}^{t-1} | \boldsymbol\alpha(u) - \widetilde{\boldsymbol\alpha}(u) |
+ | (M(\widetilde{P}(t+1)))^{-1} | \cdot | S(t+1) - \widetilde{S}(t+1) | .  \notag
\end{align}

By (<ref>),  (<ref>), and the induction hypothesis,
the above inequality implies that that

\begin{align}
|\boldsymbol\alpha(t) - \widetilde{\boldsymbol\alpha}(t)|
\leq C  \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j|^2  \Big]^{1/2} , \notag
\end{align}

and therefore

\begin{align}
\sum_{u=0}^{t} \big| \boldsymbol\alpha(u) - \widetilde{\boldsymbol\alpha}(u) \big| \leq C
\cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j |^2  \Big]^{1/2} .  \notag
\end{align}

The estimate (<ref>) is then proved by induction.

Under the hypothesis of Proposition <ref>, the costs (<ref>) under the optimal controls $\mathbf{a}$
and $\widetilde{\mathbf{a}}$ satisfy the estimate

\begin{align}
\big| L^\star ( \boldsymbol\alpha ; \mathcal{D} )
-  L^\star (  \widetilde{\boldsymbol\alpha} ;  \widetilde{\mathcal{D}} )  \big|
\leq  C \cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{x}^i_j |^2  + |y^i_j - \tilde{y}^i_j |^2 \Big]^{1/2} .
\label{eqn:stabilityLast}
\end{align}

where $C>0$ is a constant depending on $\lambda$, $\beta$, $N$, $\bar{N}$, $\Theta^\star$,
$(w^\star_1, \cdots, w^\star_N)$, $K_x$, $K_y$, and $T$.

Throughout the proof $C$ is a constant depending on $\lambda$, $\beta$, $\bar{N}$, $\Theta^\star$,
$(w_1^\star, \cdots, w_N^\star)$, $K_x$, $K_y$ and $T$, and may vary from place to place.
By (<ref>) and (<ref>), the difference of the costs can be written as

\begin{align}
L^\star ( \boldsymbol\alpha ; \mathcal{D} )
-  L^\star (  \widetilde{\boldsymbol\alpha} ;  \widetilde{\mathcal{D}} )
= & \Theta^\top (0) ( P(0) - \widetilde{P}(0) ) \Theta(0)
+ 2 ( S(0) - \widetilde{S}(0) )^\top \Theta(0) + r(0)-\widetilde{r}(0) ,  \notag \\
= &
\Theta^{\star \top} ( P(0) - \widetilde{P}(0) ) \Theta^\star
+ 2 ( S(0) - \widetilde{S}(0) )^\top \Theta^\star + r(0)-\widetilde{r}(0) .
\notag
\end{align}

By the triangular inequality, the difference satisfies

\begin{align}
\big| L^\star (  \boldsymbol\alpha ; \mathcal{D} )
-  L^\star (  \widetilde{\boldsymbol\alpha} ;  \widetilde{\mathcal{D}} ) \big|
\leq  |  P(0) - \widetilde{P}(0) | \cdot |\Theta^\star|^2
+ 2 | S(0) - \widetilde{S}(0) | \cdot | \Theta^\star | + | r(0)-\widetilde{r}(0) | .
\label{estLast-tildeLast}
\end{align}

We claim that $r(0)-\widetilde{r}(0)$ satisfies the estimate

\begin{align}
| r(0)-\widetilde{r}(0) | \leq  C \cdot \Big[ \sum_{i=1}^N w^\star_i |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \tilde{u}^i_j |^2
+ |y^i_j - \tilde{y}^i_j |^2 \Big]^{1/2} .
\label{eqn:estr-tilder(0)}
\end{align}

Then the desired estimate (<ref>) is established by (<ref>), (<ref>),
(<ref>)  and (<ref>).

Now we prove the claim (<ref>).
From (<ref>), the difference $r-\widetilde{r}$ satisfies

\begin{align}
& r(t) - \widetilde{r}(t) =  r(t+1) - \widetilde{r}(t+1)
- ( S(t+1) - \widetilde{S}(t+1) ) \big\{ (M(P(t+1)))^{-1}
\big[ S(t+1) - 2 \lambda \Theta^\star \big]  \notag \\
& \hspace{2cm} + (M(\widetilde{P}(t+1)))^{-1} \widetilde{S}(t+1) \big\}  \notag \\
&\hspace{2cm} - \widetilde{S}^\top(t+1) \big[ (M(P(t+1)))^{-1}
- (M(\widetilde{P}(t+1)))^{-1} \big]  \big[ S(t+1) - 2 \lambda \Theta^\star \big]    \notag \\
&\hspace{2cm} - \Theta^{\star \top} \lambda  \big[ (M(P(t+1)))^{-1}
- (M(\widetilde{P}(t+1)))^{-1} \big]  \lambda  \Theta^\star , \notag  \\
& r(T) - \widetilde{r}(T) =   \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} ( |y^i_j|^2 - |\widetilde{y}^i_j|^2 ) . \label{eqnr(t)-tilder(t)}
\end{align}

Inductively, $r(0)-\widetilde{r}(0)$ can be written in terms of $r(T)-\widetilde{r}(T)$ as

\begin{align}
& r(0) - \widetilde{r}(0) \notag \\
= & r(T) - \widetilde{r}(T)
- \sum_{t=1}^T \Big\{ ( S(t) - \widetilde{S}(t) ) \big[ (M(P(t)))^{-1}
\big( S(t) - 2 \lambda \Theta^\star \big)
+ (M(\widetilde{P}(t)))^{-1} \widetilde{S}(t) \big]   \notag \\
& + \widetilde{S}^\top(t) \big[ (M(P(t)))^{-1}
- (M(\widetilde{P}(t)))^{-1} \big]  \big[ S(t) - 2 \lambda  \Theta^\star \big]
+ \Theta^{\star \top} \lambda  \big[ (M(P(t)))^{-1}
- (M(\widetilde{P}(t)))^{-1} \big]  \lambda \Theta^\star
\Big\} . \notag
\end{align}

By Lemma <ref>, we have

\begin{align}
&   | r(0) - \widetilde{r}(0) | \notag \\
& \leq  | r(T) - \widetilde{r}(T) |
+ \sum_{t=1}^T \Big\{ \big| S(t) - \widetilde{S}(t) \big|
\cdot \Big[ \big| (M(P(t)))^{-1} \big| \cdot
\big(  | S(t) | + 2 \lambda | \Theta^\star| \big)
+ \big| (M(\widetilde{P}(t)))^{-1} \big| \cdot \big| \widetilde{S}(t) \big|  \Big]   \notag \\
& \hspace{2cm} +\big| (M(P(t)))^{-1}
- (M(\widetilde{P}(t)))^{-1} \big|  \cdot
\big( \big| \widetilde{S}(t) \big| \cdot \big| S(t) - 2 \lambda  \Theta^\star \big|
+  \lambda | \Theta^\star |^2 \big)
\Big\} . \label{estr(0)-tilder(0)}
\end{align}

By (<ref>), Lemma <ref>, and the Cauchy-Schwarz inequality, the terminal condition satisfies

\begin{align}
|r(T) - \widetilde{r}(T)|
= &  \Big| \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}
( y^i_j - \tilde{y}^i_j ) (y^i_j + \tilde{y}^i_j ) \Big|
\leq \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  | (y^i_j - \tilde{y}^i_j) (y^i_j + \tilde{y}^i_j ) | \notag \\
\leq & \sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  |y^i_j + \widetilde{y}^i_j| \cdot
|y^i_ j - \widetilde{y}^i_j |
\leq  2 K_y
\sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|}  |y^i_j -  \widetilde{y}^i_j|
\notag \\
\leq & 2 K_y \Big( \sum_{i=1}^N w^\star_i  \sum_{j=1}^{|\mathcal{D}_i|}  |y^i_j - \widetilde{y}^i_j|^2 \Big)^{1/2} .
\label{estr(T)-tilder(T)}
\end{align}

The estimate (<ref>) then follows from
(<ref>), (<ref>), (<ref>), (<ref>), and (<ref>).

We may now derive the proof of Theorem <ref> as a direct consequence of Corollary <ref>.

Fix $\widetilde{\mathcal{D}} \in \mathbb{D}^{p,\varepsilon}$ and let $I\subseteq \{(i,j):i=1,\dots,N,\,j=1,\dots,|\mathcal{D}_i|\}$ consist of all indices $(i,j)$ for which $\{(i,j):\, u^i_j\neq \widetilde{u}^i_j \mbox{ and } y^i_j\neq \widetilde{y}^i_j\}$.  By Corollary <ref>, there is a constant $C>0$ which is independent of $\mathcal{D}$ such that

\begin{align}
\nonumber
\big| L^\star (  \boldsymbol\alpha ; \mathcal{D} )
-  L^\star (  \widetilde{\boldsymbol\alpha} ;  \widetilde{\mathcal{D}} )  \big|^2
\le &
\widetilde{C} \sum_{i=1}^N w_i^\star |\mathcal{D}_i| \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j - \widetilde{u}^i_j |^2  + |y^i_j - \widetilde{y}^i_j |^2
\\
\nonumber
\le &
\widetilde{C} \cdot \sum_{(i,j)\in I}\, w_i^\star |\mathcal{D}_i|  ( | u^i_j - \widetilde{u}^i_j |^2  + |y^i_j - \widetilde{y}^i_j |^2 )
\\
\nonumber
\le &
\widetilde{C} \cdot \sum_{(i,j)\in I}\, w_i^\star |\mathcal{D}_i| \cdot \max\{| u^i_j - \widetilde{u}^i_j |^2, |y^i_j - \widetilde{y}^i_j |^2\}
\\
\nonumber
= &
\widetilde{C} \cdot \sum_{(i,j)\in I}\, w_i^\star |\mathcal{D}_i| \cdot \max\{| u^i_j - \widetilde{u}^i_j |, |y^i_j - \widetilde{y}^i_j |\}^2
\\
\nonumber
\le &
\widetilde{C} \cdot \sum_{(i,j)\in I}\, w_i^\star |\mathcal{D}_i| \cdot 2\varepsilon^2
\\
\nonumber
\le &
\widetilde{C} \#I\, p \, 2\varepsilon^2
\end{align}

where we have set $\widetilde{C}\eqdef C^{1/2}$.

The computational complexity of the matrix-valued function $P(\cdot)$, as determined by (<ref>), does not primarily stem from the high dimensionality of the output matrices, but rather from their inherent lack of symmetries.
As we will see later, such symmetry will allow us to efficiently encode $P(\cdot)$ into low-dimensional structures, therefore completely mitigating the high-dimensionality effect.

We therefore introduce the following control, which is optimal for a “symmetrized version” of the control problem (<ref>) where $w^{\star}=(1/N,\dots,1/N)$.  This advantages this surrogate problem is that the associated $P$ function, which we denote by $\widehat{P}$, outputs highly symmetric[As formalized in Proposition <ref>.] Toeplitz matrices which are determined by exactly three factors. This allows us to encode the entire control problem into a control problem on $\mathbb{R}^{\tilde{O}(P)}$ whose dimension is independent of $N$.
Furthermore, the optimal controls of both problems can be related by stability estimates, depending only on the deviation of $(1/N,\dots,1/N)$ from $w^{\star}$, which allows us to infer the degree of sub-optimality of our $\widehat{P}$.  Therefore, we consider the control

\begin{align}
\widehat{\boldsymbol{\alpha}}(t)
= & -[(\lambda + \beta )I_{Np} + \widehat{P}(t+1)]^{-1}
\big[ (\lambda I_{Np} + \widehat{P}(t+1)) \Theta(t)
- \lambda  \Theta^\star_{(N)} + \widehat{S}(t+1)  \big]   \label{hata}
\end{align}

We will show that the algorithm under control (<ref>) enjoys lower computational complexity in comparison with the regret-optimal algorithm (<ref>), whose computational complexity is summarized through Theorem <ref>. Furthermore, we will then establish the near-regret optimality of (<ref>), stated in Theorem <ref>.

The associated matrix-valued functions $\widehat{P}$ and $\widehat{S}$ are uniquely determined by

\begin{align}
& \begin{cases}
\widehat{P}(t)  =  \beta I_{Np} - \beta^2 [(\lambda+\beta ) I_{Np} + \widehat{P}(t+1) ]^{-1} , \label{eqnhatP1/N(t)} \\
\widehat{P}(T) =  [ (1/N) I_p ,\cdots, (1/N) I_p]^{\top}
\Big( \sum_{i=1}^N (1/N) \sum_{j=1}^{|\mathcal{D}_i|} u^i_j u^{i\top}_j  \Big)
[ (1/N) I_p , \cdots, (1/N) I_p ] ,
\end{cases} \\
& \begin{cases}
\widehat{S}(t) =
\beta [(\lambda + \beta ) I_{Np} + \widehat{P}(t+1) ]^{-1}
\big( \widehat{S}(t+1) - \lambda \Theta^\star_{(N)} \big)
\label{eqnhatS1/N(t)} ,     \\
\widehat{S}(T)
- [(1/N)I_p, \cdots, (1/N)I_p ]^{\top} \sum_{i=1}^N(1/N)\sum_{j=1}^{|\mathcal{D}_i|}  u^i_j y^i_j
\end{cases}
\end{align}

where $[ I_p, \cdots, I_p ]$ denotes the $Np\times p$-dimensional matrix formed by concatenating $N$ copies of the $p\times p$-dimensional identity matrix $\frac{1}{N} I_p$.   The parameter $\Theta^\star_{(N)}$ is defined as

\begin{align}
\Theta^\star_{(N)} \eqdef [\theta_{(N)}^{\star \top}, \cdots, \theta_{(N)}^{\star \top}]^{\top} , \quad
\theta^\star_{(N)} \eqdef \frac{1}{N} \sum_{i=1}^N \theta^\star_i . \notag
\end{align}

According to (<ref>) in Assumption <ref>, $\widehat{P}(T)\geq 0$ and therefore (<ref>) admits a unique solution.

Similar to Corollary <ref>, we have $\big| [ (\lambda+\beta)I_{Np} + \widehat{P}(t) ]^{-1} \big| \leq (\lambda + \beta)^{-1}$. By replacing $(w_1^\star, ..., w_N^\star)$ in the proof of Corollary <ref> with $(1/N, ..., 1/N)$, we obtain
$|\widehat{S}(t)|\leq C \cdot N^{1/2}$ for all $t=1$, ..., $T$,
where $C>0$ is a constant depending on $\lambda$, $\beta$, $\Theta^\star$, $K_x$, $K_y$ and $T$.

Let $P$ and $\widehat{P}$ be the solutions of the systems (<ref>) and (<ref>), respectively. Then they satisfy that for all $t=1$, $2$, ..., $T$,

\begin{align}
\big| P(t) - \widehat{P}(t) \big|
\leq  C \cdot \bar{N} N \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big| ,
\label{estPw-P1/N} \\
\big| (M(P(t)))^{-1} - (M(\widehat{P}(t) ))^{-1} \big|
\leq  C \cdot  \bar{N} N \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big| ,
\label{estMPw-MP1/N}
\end{align}

where $C>0$ is a constant depending on $\lambda$, $\beta$, $K_x$.

By (<ref>), (<ref>), and Lemmas <ref> and <ref>, the difference
$\widehat{P}(T)-P(T)$ satisfies

\begin{align}
|\widehat{P}(T) - P(T)|^2 = & \sum_{k,l=1}^N
\Big| \frac{1}{N^2} \sum_{k=1}^N \frac{1}{N}
\sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j  - w^\star_k w^\star_l \sum_{i=1}^N w^\star_i  \sum_{j=1}^{|\mathcal{D}_i|} u^i_j u^{i\top}_j  \Big|^2 \notag \\
\leq & N^2 \cdot \max_{1\leq k, l, i \leq N} \Big| \frac{1}{N^3}
- w^\star_k w^\star_l
w^\star_i \Big|^2 \cdot \Big| \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|}  u^i_j u^{i\top}_j  \Big|^2  \notag \\
\leq &   N^2 \cdot \max_{1\leq k, l, i \leq N} \Big| \frac{1}{N^3}
-  w^\star_k w^\star_l
w^\star_i \Big|^2 \cdot  \bar{N} \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|}  | u^i_j |^4   \notag
\end{align}

\begin{align}
|\widehat{P}(T) - P(T)|
\leq & \bar{N} N K_x^2 \cdot \max_{1\leq k, l, i \leq N} \Big| \frac{1}{N^3}
-  w^\star_k w^\star_l
w^\star_i \Big| .
\label{|P(T)-hatP(T)|est}
\end{align}

Therefore (<ref>) holds for $t=T$.

Assume by induction that (<ref>) holds for $t=u$. We show that (<ref>) also holds for $t=u-1$.
From (<ref>) and (<ref>), we have

\begin{align}
\widehat{P}(u-1) - P(u-1)
= & [ M(P(u)) ]^{-1}
-  [ M( \widehat{P}(u) ) ]^{-1}  \notag  \\
= & [ M( P(u) ) ]^{-1} (\widehat{P}(u) - P(u))
[ M( \widehat{P}(u) ) ]^{-1} .
\label{hatP-Precursion}
\end{align}

By Lemmas <ref> and <ref>, we further have

\begin{align}
| \widehat{P}(u-1) - P(u-1) |
= & | [ M(P(u)) ]^{-1}
-  [ M( \widehat{P}(u) ) ]^{-1} | \notag  \\
\leq & | [ M( P(u) ) ]^{-1} | \cdot  | \widehat{P}(u) - P(u) | \cdot
| [ M( \widehat{P}(u) ) ]^{-1} | .
\label{|hatP-P|recursion}
\end{align}

By (<ref>), Corollary <ref>, Remark <ref>, and the induction hypothesis, we have
that (<ref>) holds for $t=u-1$. By induction we have shown that (<ref>) holds for all $1\leq t\leq T$.

From (<ref>), we have

\begin{align}
& \big| ( M(P(t)) )^{-1}
-  ( M( \widehat{P}(t) ) )^{-1} \big|
\leq  \big| ( M( P(t) ) )^{-1} \big| \cdot  \big| \widehat{P}(t) - P(t) \big| \cdot
\big| ( M( \widehat{P}(t) ) )^{-1}
\big| ,  \notag
\end{align}

which, together with (<ref>), Corollary <ref>, and Remark <ref>, establishes (<ref>).

Let $S$ and $\widehat{S}$ be the solutions of (<ref>) and (<ref>), respectively. Then they satisfy that
for all $t=1$, ..., $T$,

\begin{align}
|S(t) - \widehat{S}(t)| \leq &
C \cdot \Big\{ \bar{N} N^{3/2} \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big|
\notag \\
& \hspace{1cm} + \bar{N} N^{1/2} \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|
+ |\Theta^\star - \Theta^\star_{(N)}| \Big\} ,
\label{sup|S-hatS|est}
\end{align}

where $C>0$ is a constant depending on $\lambda$, $\beta$, $\Theta^\star$, $K_x$, $K_y$, and $T$.

By (<ref>), (<ref>), and and Lemmas <ref> and <ref>, the terminal condition $S(T)-\widehat{S}(T)$ satisfies

\begin{align} \Big|S(T)-\widehat{S}(T)\Big|^2
= & \sum_{k=1}^N \Big|
\frac{1}{N} \sum_{i=1}^N \frac{1}{N} \sum_{j=1}^{|\mathcal{D}_i|} u^i_j y^i_j - w^\star_k
\sum_{i=1}^N w^\star_i \sum_{j=1}^{|\mathcal{D}_i|} u^i_j y^i_j
\Big|^2 \notag \\
\leq & N \cdot \max_{1\leq k, l \leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|^2 \cdot \Big| \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|} u^i_j y^i_j \Big|^2  \notag \\
\leq & N \cdot \max_{1\leq k, l \leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|^2 \cdot \bar{N} \sum_{i=1}^N \sum_{j=1}^{|\mathcal{D}_i|} | u^i_j|^2 \cdot |y^i_j|^2 ,
\notag
\end{align}

which implies

\begin{align}
\big| S(T) - \widehat{S}(T) \big|
\leq & \bar{N} N^{1/2} \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big| \cdot   K_x K_y .
\label{|S(T)-hatS(T)|est}
\end{align}

Therefore (<ref>) holds for $t=T$.

Assume by induction that (<ref>) holds for $t=u$, we show that (<ref>) holds for $t=u-1$.
From (<ref>) and (<ref>), the difference $S(u-1)-\widehat{S}(u-1)$ can be written as

\begin{align}
S(u-1) - \widehat{S}(u-1) = & [ M(P(u)) ]^{-1}[S(u) - \widehat{S}(u) - \lambda (\Theta^\star - \Theta^\star_{(N)}) ] \notag \\
& - \big[ ( M(\widehat{P}(u)))^{-1} - ( M(P(u)) )^{-1} \big]
(\widehat{S}(u) - \lambda \Theta^\star_{(N)}) ,
\label{S(t)-hatS(t)}
\end{align}

and satisfies by
Lemma <ref> that

\begin{align}
|S(u-1) - \widehat{S}(u-1)|
\leq & | [M(P(u))]^{-1} | \cdot
\big[  |S(u) - \widehat{S}(u)| + \lambda |\Theta^\star - \Theta^\star_{(N)}|
\big] \notag \\
& + \big| ( M(\widehat{P}(u)))^{-1}
- ( M(P(u)) )^{-1} \big| \cdot
\big( |\widehat{S}(u)| + \lambda |\Theta^\star_{(N)}| \big) . \label{|hatS-S|recursion}
\end{align}

By (<ref>), (<ref>), Remark <ref>, and the induction hypothesis, we have that
(<ref>) holds for $t=u-1$. By induction we have shown that (<ref>) holds for all $t=1$, ..., $T$.

We now obtain Theorem <ref>; which is a direct consequence of the following Propositions <ref> and <ref>.

The controls
(<ref>) and (<ref>) satisfy $|\widehat{\boldsymbol{\alpha}}(t)| \leq C \cdot N^{1/2}$ for all $t=0$, $1$, ..., $T-1$, and

\begin{align}
\sum_{t=0}^{T-1} |\boldsymbol{\alpha}(t) - \widehat{\boldsymbol{\alpha}}(t)|
\leq & C\cdot \Big\{ \bar{N} N^{3/2} \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big| \notag \\
&\hspace{2cm} + \bar{N} N^{1/2} \cdot \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|
+ |\Theta^\star - \Theta^\star_{(N)}| \Big\} ,
\label{eqn:stabilityaw-a1/N}
\end{align}

where $C>0$ is a constant depending on $\lambda$, $\beta$, $\Theta^\star$, $K_x$, $K_y$, and $T$.

We let $C$ be a constant depending on $\lambda$, $\beta$, $\Theta^\star$, $K_x$, $K_y$, and $T$,
and $C$ is allowed to vary from place to place throughout the proof.
By (<ref>), (<ref>), and (<ref>), we can write $\boldsymbol{\alpha}$ and $\widehat{\boldsymbol\alpha}$ as

\begin{align}
\boldsymbol{\alpha}(t)
= & [ (M(P(t+1)))^{-1} - I_{Np} ]  \Big[ \Theta^\star
+ \sum_{u=0}^{t-1} \boldsymbol{\alpha}(u) \Big]
+  [ M(P(t+1)) ]^{-1} ( \lambda  \Theta^\star - S(t+1) ) ,   \notag \\
\widehat{\boldsymbol{\alpha}}(t)
= & [ (M(\widehat{P}(t+1)))^{-1} - I_{Np} ]  \Big[ \Theta^\star
+ \sum_{u=0}^{t-1} \widehat{\boldsymbol{\alpha}}(u) \Big]
+  [ M(\widehat{P}(t+1)) ]^{-1} ( \lambda  \Theta^\star_{(N)} - \widehat{S}(t+1) ) .   \notag
\end{align}

For $t=0$, we have

\begin{align}
|\widehat{\boldsymbol{\alpha}}(0)| \leq  & | (M(\widehat{P}(1)))^{-1} - I_{Np} | \cdot |  \Theta^\star |
+ \big| [ M(\widehat{P}(1)) ]^{-1} \big| ( \lambda  |\Theta^\star_{(N)}| + |\widehat{S}(1)| ) ,  \notag
\end{align}

which implies $|\widehat{\boldsymbol{\alpha}}(0)| \leq C N^{1/2}$.
Assume by induction that $|\widehat{\boldsymbol{\alpha}}(u)| \leq C N^{1/2}$ for $u=0$, ..., $t-1$, then

\begin{align}
\widehat{\boldsymbol{\alpha}}(t)
\leq & \big| (M(\widehat{P}(t+1)))^{-1} - I_{Np} \big|   \Big[ | \Theta^\star |
+ \sum_{u=0}^{t-1} | \widehat{\boldsymbol{\alpha}}(u) | \Big] \notag \\
& + \big| [ M(\widehat{P}(t+1)) ]^{-1} \big| ( \lambda | \Theta^\star_{(N)} | + | \widehat{S}(t+1) | ) \notag
\end{align}

gives that $|\widehat{\boldsymbol{\alpha}}(t)|\leq CN^{1/2}$. By induction we have shown that
$|\widehat{\boldsymbol{\alpha}}(t)|\leq CN^{1/2}$ for all $t=0$, $1$, ..., $T-1$.

We further write the difference $\boldsymbol{\alpha} - \widehat{\boldsymbol\alpha}$ as

\begin{align}
& \boldsymbol{\alpha}(t) - \widehat{\boldsymbol\alpha}(t) \notag \\
= &  \big[(M(P(t+1)))^{-1} - (M(\widehat{P}(t+1)) )^{-1} \big]
\Big[ (\lambda + 1)\Theta^\star + \sum_{u=0}^{t-1}  \boldsymbol{\alpha}(u)    - S(t+1) \Big]  \notag \\
& + \big[ ( M(\widehat{P}(t+1)) )^{-1} - I_{Np} \big]
\sum_{u=0}^{t-1} (\boldsymbol{\alpha}(u) - \widehat{\boldsymbol{\alpha}}(u))
+ [ M(\widehat{P}(t+1)) ]^{-1} (\widehat{S}(t+1) - S(t+1))  \notag \\
& + \lambda [ M(\widehat{P}(t+1)) ]^{-1} ( \Theta^\star
- \Theta^\star_{(N)} ) . \label{a-hata}
\end{align}

For $t=0$, the difference

\begin{align}
\boldsymbol{\alpha}(0) - \widehat{\boldsymbol{\alpha}}(0)
= & \big[ (M(P(1)))^{-1} - (M(\widehat{P}(1)))^{-1}  \big]
\big[ ( \lambda + 1 )   \Theta^\star_{(N)} - \widehat{S}(1)  \big] \notag \\
& + [ M(P(1)) ]^{-1} [ (\lambda+1) (\Theta^\star - \Theta^\star_{(N)} ) + \widehat{S}(1) - S(1) ]   \notag
\end{align}

\begin{align}
\big| \boldsymbol{\alpha}(0)
- \widehat{\boldsymbol{\alpha}}(0) \big|
\leq &
C  \cdot \Big\{ \bar{N} N^{3/2} \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big| \notag \\
& \hspace{1cm} + \bar{N} N^{1/2} \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|
+ |\Theta^\star - \Theta^\star_{(N)}| \Big\} , \notag
\end{align}

due to (<ref>), (<ref>), Lemma <ref>, Corollary <ref> and Remark <ref>.
Assume by induction for all $0 \leq u\leq t-1$,
$\boldsymbol{\alpha}(u-1) - \widehat{\boldsymbol{\alpha}}(u-1)$ satisfies

\begin{align}
\big| \boldsymbol{\alpha}(u-1)
- \widehat{\boldsymbol{\alpha}}(u-1) \big|
\leq & C
\cdot \Big\{ \bar{N} N^{3/2} \max_{1\leq i, k, l \leq N} \Big| \frac{1}{N^3} -  w^\star_k w^\star_l w^\star_i \Big| \notag \\
& \hspace{1cm} + \bar{N} N^{1/2} \max_{1\leq k, l\leq N} \Big| \frac{1}{N^2} - w^\star_k w^\star_l \Big|
+ |\Theta^\star - \Theta^\star_{(N)}| \Big\}  . \notag
\end{align}

By (<ref>) we have

\begin{align}
& | \boldsymbol{\alpha}(t) - \widehat{\boldsymbol{\alpha}}(t) | \notag \\
\leq &  \big| (M(P(t+1)))^{-1} - (M(\widehat{P}(t+1)) )^{-1} \big| \cdot
\Big[ \sum_{u=0}^{t-1}  | \boldsymbol{\alpha}(u) |
+ | ( \lambda +1 ) \Theta^\star |
+| S(t+1) | \Big]  \notag \\
& + \big| ( M(\widehat{P}(t+1)) )^{-1} - I_{Np}  \big|
\sum_{u=0}^{t-1} | \boldsymbol{\alpha}(u) - \widehat{\boldsymbol{\alpha}}(u) |