Adaptive Kernel Learning in Heterogeneous Networks

TL;DR

This paper introduces HALK, a decentralized algorithm for heterogeneous networks that adaptively learns kernel-based regression functions while satisfying nonlinear proximity constraints, achieving improved convergence rates.

Contribution

It proposes a novel decentralized stochastic primal-dual method with adaptive kernel projections for heterogeneous networks, enhancing convergence and constraint satisfaction.

Findings

Achieves (\u221a{T}) sub-optimality attenuation with constant step-sizes.

Ensures long-term satisfaction of nonlinear proximity constraints.

Improves upon existing rates for vector-valued decentralized learning.

Abstract

We consider learning in decentralized heterogeneous networks: agents seek to minimize a convex functional that aggregates data across the network, while only having access to their local data streams. We focus on the case where agents seek to estimate a regression \emph{function} that belongs to a reproducing kernel Hilbert space (RKHS). To incentivize coordination while respecting network heterogeneity, we impose nonlinear proximity constraints. To solve the constrained stochastic program, we propose applying a functional variant of stochastic primal-dual (Arrow-Hurwicz) method which yields a decentralized algorithm. To handle the fact that agents' functions have complexity proportional to time (owing to the RKHS parameterization), we project the primal iterates onto subspaces greedily constructed from kernel evaluations of agents' local observations. The resulting scheme, dubbed Heterogeneous Adaptive Learning with Kernels (HALK), when used with constant step-sizes, yields $\mathcal{O}(\sqrt{T})$ attenuation in sub-optimality and exactly satisfies the constraints in the long run, which improves upon the state of the art rates for vector-valued problems.

Tables

Table 1. TABLE I : Comparison of our work with the related works in the literature. ℝ d superscript ℝ 𝑑 {\mathbb{R}}^{d} denotes d 𝑑 d -dimensional Euclidean space, ℋ ℋ {\mathcal{H}} denotes an RKHS, and p 𝑝 p denotes the data dimension. Our work achieves null constraint violation on average with comparable average sub-optimality rate.

Reference	Average Sub-optimality	Average Constraint violation	Multi-agent	Function Class	Model Complexity
[7, 8]	$𝒪\left(T^{-1/2}\right)$	$𝒪\left(T^{-1/2}\right)$	✗	${ℝ}^d$	N.A.
[9, 10]	$𝒪\left(T^{-1/2}\right)$	$𝒪\left(T^{-1/4}\right)$	✓	${ℝ}^d$	N.A.
[9]	$𝒪\left(T^{-1/4}\right)$	zero	✗	${ℝ}^d$	N.A.
[11]	$𝒪\left(T^{-1/2}\right)$	$𝒪\left(T^{-1}\right)$	✗	${ℝ}^d$	N.A.
[12, 13]	$𝒪\left(T^{-1/2}\right)$	$𝒪\left(T^{-1/4}\right)$	✗	${ℝ}^d$	N.A.
[14, 15, 16, 17]	$𝒪\left(T^{-1/2}\right)$	$𝒪\left(T^{-1/4}\right)$	✗	$\mathscr{H}$	finite
This Work	$𝒪\left(T^{-1/2}+\alpha \right)$	zero	✓	$\mathscr{H}$	$𝒪\left({\alpha }^{-2p}\right)$

Equations

$$\underset{f\in \mathcal{H}}{\min }{\mathbb{E}}_{\mathbf{x},y}\left[\ell \left(f(\mathbf{x}),y\right)\right]+\frac{\lambda }{2}{\Vert f\Vert }_{\mathcal{H}}^2.$$

$$\displaystyle\min_{\{f_{i}\in{\mathcal{H}}\}}S({\mathbf{f}}):=\sum_{i\in{\mathcal{V}}}\Big{(}{\mathbb{E}}_{{\mathbf{x}}_{i},y_{i}}\Big{[}\ell_{i}(f_{i}\big{(}{\mathbf{x}}_{i}),y_{i}\big{)}\Big{]}+\frac{\lambda}{2}\|f_{i}\|^{2}_{{\mathcal{H}}}\Big{)}.$$

$${\mathbf{f}}^{\star }=\underset{\mathbf{f}\in {\mathcal{H}}^V}{\mathrm{argmin}}S(\mathbf{f})\text{s.t. }{H}_{ij}(f_i,f_j)\le {\gamma }_{ij},\text{for all}j\in n_i,$$

$$(i)⟨\tilde{f},\kappa ({\mathbf{x}}_i,\cdot )){⟩}_{\mathcal{H}}=\tilde{f}({\mathbf{x}}_i),(ii)\mathcal{H}=\overline{\text{span}\{\kappa ({\mathbf{x}}_i,\cdot )\}}$$

$$\tilde{f}({\mathbf{x}}_i)=\sum _{k=1}^N{w}_{i,k}\kappa ({\mathbf{x}}_{i,k},{\mathbf{x}}_i),$$

$$\underset{\mathbf{f}\in {\mathcal{H}}^V}{\min }\frac{1}{N}\sum _{k=1}^N{S}^k(\mathbf{f}).$$

$$\underset{\{{\mathbf{w}}_i\}\in {\mathbb{R}}^N}{\min }\frac{1}{N}\sum _{k=1}^N\sum _{i\in \mathcal{V}}{\ell }_i({\mathbf{w}}_i^T{\mathbf{\kappa }}_{{\mathbf{X}}_i}({\mathbf{x}}_{i,k}),y_{i,k})+\frac{\lambda }{2}{\mathbf{w}}_i^T{\mathbf{K}}_i{\mathbf{w}}_i,$$

$${\tilde{\mathbf{f}}}^{\star }$$

$$+\frac{\lambda }{2}\sum _{n,m\in \mathcal{I}}{w}_{i,n}{w}_{i,m}\kappa ({\mathbf{x}}_{i,m},{\mathbf{x}}_{i,n})].$$

$$({\mathbf{f}}^{\star },{\mathbf{\mu }}^{\star })=\text{arg}\underset{\mathbf{\mu }}{\max }\underset{\mathbf{f}}{\min }\tilde{\mathcal{L}}(\mathbf{f},\mathbf{\mu }).$$

$$({\mathbf{f}}^{\star },{\mathbf{\mu }}^{\star })=\text{arg}\underset{\mathbf{\mu }}{\max }\underset{\mathbf{f}}{\min }{\mathcal{L}}^e(\mathbf{f},\mathbf{\mu }),$$

$${\mathcal{L}}^e(\mathbf{f},\mathbf{\mu }):=\frac{1}{N}\sum _{k=1}^N\left[S^k(\mathbf{f})+⟨\mathbf{\mu },{\mathbf{H}}^k(\mathbf{f})-\mathbf{\gamma }⟩\right]$$

$${\mathbf{f}}_{\nu }^{\star }=\underset{\{f_i\in \mathcal{H}\}}{\mathrm{argmin}}S(\mathbf{f})\text{s.t.}{H}_{ij}(f_i,f_j)+\nu \le {\gamma }_{ij},\text{for all}j\in n_i,$$

$$\mathcal{L}(\mathbf{f},\mathbf{\mu })=$$

$${\hat{\mathcal{L}}}_t(\mathbf{f},\mathbf{\mu }):=S^t(\mathbf{f})+⟨\mathbf{\mu },{\mathbf{H}}^t(\mathbf{f})+\nu 1-\mathbf{\gamma }⟩-\frac{\delta \eta }{2}{\Vert \mathbf{\mu }\Vert }^2$$

$${\nabla }_{f_i}{\ell }_i(f_i({\mathbf{x}}_{i,t}),y_{i,t})=\frac{\partial {\ell }_i(f_i({\mathbf{x}}_{i,t}),y_{i,t})}{\partial f_i({\mathbf{x}}_{i,t})}\frac{\partial f_i({\mathbf{x}}_{i,t})}{\partial f_i}$$

$$\frac{\partial f_i({\mathbf{x}}_{i,t})}{\partial f_i}=\frac{\partial ⟨f_i,\kappa ({\mathbf{x}}_{i,t},\cdot ){⟩}_{\mathcal{H}}}{\partial f_i}=\kappa ({\mathbf{x}}_{i,t},\cdot ).$$

$${\nabla }_{f_i}{\hat{\mathcal{L}}}_t({\mathbf{f}}_t,{\mathbf{\mu }}_t)$$

$$+\sum _{j\in n_i}{\mu }_{ij}{h}_{ij}^{\prime }(f_i({\mathbf{x}}_{i,t}),f_j({\mathbf{x}}_{i,t}))\kappa ({\mathbf{x}}_{i,t},\cdot ).$$

$${\nabla }_{\mathbf{f}}{\hat{\mathcal{L}}}_t({\mathbf{f}}_t,{\mathbf{\mu }}_t)=\text{vec}[{\nabla }_{f_i}{\hat{\mathcal{L}}}_t({\mathbf{f}}_t,{\mathbf{\mu }}_t)].$$

$${\mathbf{f}}_{t+1}=$$

$$\displaystyle+\sum_{j\in n_{i}}\mu_{ij,t}h^{\prime}_{ij}(f_{i,t}({\mathbf{x}}_{i,t}),f_{j,t}({\mathbf{x}}_{i,t}))\Big{]}\kappa({\mathbf{x}}_{i,t},\cdot)\Big{)},$$

$${\mathbf{\mu }}_{t+1}=$$

$$(f_{i,t+1},{\mathbf{D}}_{i,t+1},{\mathbf{w}}_{i,t+1})=\text{KOMP}({\tilde{f}}_{i,t+1},{\tilde{\mathbf{D}}}_{i,t+1},{\tilde{\mathbf{w}}}_{i,t+1},ϵ)$$

$$f_{i,t+1}=$$

$$\displaystyle+\sum_{j\in n_{i}}\mu_{ij,t}h^{\prime}_{ij}(f_{i,t}({\mathbf{x}}_{i,t}),f_{j,t}({\mathbf{x}}_{i,t}))\Big{]}\kappa({\mathbf{x}}_{i,t},\cdot).$$

$${\nabla }_{{\mu }_{ij}}{\hat{\mathcal{L}}}_t({\mathbf{f}}_t,{\mathbf{\mu }}_t)=h_{ij}(f_{i,t}({\mathbf{x}}_{i,t}),f_{j,t}({\mathbf{x}}_{i,t}))-{\gamma }_{ij}+\nu -\delta \eta {\mu }_{ij}.$$

$$\displaystyle\mu_{ij,t+1}\!\!=\!\!\Big{[}\mu_{ij,t}(\!1\!-\!\delta\eta^{2})\!+\!\eta\!\Big{(}\!\!h_{ij}(f_{i,t}({\mathbf{x}}_{i,t}),f_{j,t}({\mathbf{x}}_{i,t}))\!-\!\gamma_{ij}\!\!+\!{\nu}\!\Big{)}\!\Big{]}_{+}.$$

$$f_{i,t}(\mathbf{x})=\sum _{n=1}^{t-1}{w}_{i,n}\kappa ({\mathbf{x}}_{i,n},\mathbf{x})={\mathbf{w}}_{i,t}^T{\mathbf{\kappa }}_{{\mathbf{X}}_{i,t}}(\mathbf{x}).$$

$${\mathbf{X}}_{i,t+1}=$$

$$[{\mathbf{w}}_{i,t+1}{]}_u=$$

$$f_{i,t+1}=$$

$$:=$$

$$\displaystyle+\sum_{j\in n_{i}}\mu_{ij,t}h^{\prime}_{ij}(f_{i,t}({\mathbf{x}}_{i,t}),f_{j,t}({\mathbf{x}}_{i,t}))\Big{]}\kappa({\mathbf{x}}_{i,t},\cdot)\Big{]}$$

Full text

I Introduction

In decentralized optimization, each agent $i\in{\mathcal{V}}$ in a network $\mathcal{G}=({\mathcal{V}},\mathcal{E})$ has a local objective but seek to cooperate with other agents to minimize the global network objective. The agents communicate only with their neighbors for solving the global objective. This global objective is the sum of local convex objectives available at different nodes of the network and depends upon the locally observed information. This framework has yielded, for instance, networked controllers [2], signal processing [3], robotics [4], and communications [5].

In this work, we focus on the case where the agents that comprise the interconnected network may be of different types, such as aerial and ground robots collaboratively gathering information [4], or wireless channel estimation when spatial covariates are present [5]. In such settings, local information may be distinct, but performance may still be boosted by information sharing among agents. This phenomenon may be mathematically encapsulated as convex non-linear proximity constraints. We focus on the case where each agent’s objective depends on a data stream, i.e., the online case, and the observations provided to the network are heterogeneous, when agents decisions are defined not by a standard parameter vector but instead a nonlinear regression function that belongs to a reproducing kernel Hilbert space (RKHS) [6].

Setting aside the constraints, the solution of stochastic programs, assuming no closed form exists, necessitates iterative tools. The simplest approach, gradient descent, requires evaluating an expectation which depends on infinitely many data realizations. This issue may be overcome through stochastic gradient descent (SGD) [18], which alleviates the dependence on the sample size by using stochastic gradients in lieu of true gradients, and hence is popular in large-scale learning [19]. However, its limiting properties are intrinsically tied to the statistical parameterization (decision variable) one chooses. For vector-valued problems, i.e., linear statistical models, the convergence of SGD is well-understood via convexity [20].

By contrast, optimization problems induced by those with richer descriptive capability (owing to universal approximation [21]), are more challenging. Dictionary learning [22] and deep networks [23] trade convexity for descriptive richness, which has led to a flurry of interest in non-convex stochastic optimization [24]. Generally, overcoming non-convexity requires adding noise that degrades parameter estimates [25], which may then prove inoperable for online systems. Instead, one may preserve convexity while obtaining nonlinearity (universality) through the “kernel trick” [26]. Owing to the Representer Theorem [27], one may transform the function variable to an inner product of weights and kernel evaluations at samples. Unfortunately, the representational complexity is proportional with the sample size $N$ [27], which for online settings $N\rightarrow\infty$. To address this issue, we employ hard-thresholding projections onto subspaces constructed greedily from the history of data observation via matching pursuit [28], which nearly preserves global convergence [29].

Now, we shift focus to multi-agent optimization. Typically, the global cost is additive across agents’ individual costs. Thus decentralized schemes may be derived by constraining agents’ decisions to be equal. One may solve such problems via primal-only schemes via penalty method [30], reformulating the consensus constraint in the dual domain [31], and primal-dual approaches [32] which alternate primal/dual descent/ascent steps on the Lagrangian. Approximate dual methods [33], i.e., ADMM, have also been used [34]. Beyond linear equality constraints, motivated by heterogeneous networks, only primal methods and exact primal-dual approaches are viable, since dual methods/ADMM require solving a nonlinear $\operatornamewithlimits{argmin}$ in the inner-loop which is prohibitively costly. Hence, we adopt a primal-dual approach to solving the proximity-constrained problem [10] over the more general RKHS setting [35], which is developed in detail in Sec. II. To do so, we generalize RKHS primal-dual method [15] to multi-agent optimization (Sec III), and obtain a new collaborative learning systems methodology which we call Heterogeneous Adaptive Learning with Kernels (HALK)(Sec III). Compared to [15], several technical contributions are unique to this work:

• •

We establish an explicit dependence of the number of dictionary elements (model complexity) that must be retained for the proposed non-parametric approach to yield solutions that are arbitrarily accurate. Thm. 1 characterizes that for a given choice of the parameter $\alpha$, the model order is upper bounded by ${\mathcal{O}}(\alpha^{-2p})$. Such a non-asymptotic characterization of model complexity has never been reported in the literature for any of the non-parametric algorithms.

• •

Relative to existing primal-dual algorithms for constrained stochastic optimization in RKHS [15, 35], we have introduced a regularization of the dual update in terms of problem constants that permits one to obtain average sub-optimality rate of ${\mathcal{O}}(T^{-1/2}+\alpha)$ which is similar to the tightest sub-optimality rates ${\mathcal{O}}(\sqrt{T})$ available in the literature, while ensuring feasibility on average, i.e., null constraint violation, in contrast to ${\mathcal{O}}(T^{3/4})$ rate for existing settings (Thm. 2). The dependence of the optimality gap on the model complexity ($\alpha$) along with the number of iterations ($T$) is established for the first time in the RKHS literature.

• •

In Sec. V, we demonstrate the algorithm’s effectiveness for spatio-temporal correlated Gaussian random field estimation (Sec. V-A). Moreover, we experimentally test it on a real ocean data set for monitoring ocean salinity and temperature on buoys at various depths (Sec. V-B). An extended experiment adapts bandwidths online [36]: each agent’s kernel adapts to its local data, and thus outperforms settings where local hyper-parameters are fixed.

Discussion of Rates Regarding the context of Thm. 2, we note that existing efforts to obtain strict feasibility only obtain optimality gap attenuation at ${\mathcal{O}}(T^{3/4})$ [9], or obtain ${\mathcal{O}}(\sqrt{T})$ sub-optimality with comparable constraint violation [7, 8], with the exception of a complicated barrier methods that do not generalize to the learning settings [11]. In this work, we achieve strict feasibility on average without degrading the sub-optimality gap. See Table I for details. Moreover, [12, 13] obtain a tunable tradeoff $\beta\in(0,1)$ between sub-optimality and constraint violation, which we do not consider for simplicity, and specifically considering $\beta=0.5$, we present their best result in Table I.

II Problem Formulation

In supervised learning, data takes the form of input-output examples, $({\mathbf{x}},y)$, which are independent identically distributed (i.i.d.) realizations from a distribution $\mathbb{P}({\mathbf{x}},y)$ with support ${\mathcal{X}}\times{\mathcal{Y}}$. Here, ${\mathcal{X}}\subset{\mathbb{R}}^{p}$ and ${\mathcal{Y}}\subset{\mathbb{R}}$. In this work, we focus on expected risk minimization where one seeks to compute the minimizer of a loss quantifying the merit of a nonlinear statistical model $f\in{\mathcal{H}}$ averaged over the data $\{({\mathbf{x}},y)\}$. Setting aside the choice of ${\mathcal{H}}$ for now, the merit of estimator $f$ is quantified by the convex loss function $\ell:{\mathcal{H}}\times{\mathcal{X}}\times{\mathcal{Y}}\rightarrow{\mathbb{R}}$ which is small when estimator $f({\mathbf{x}})$ evaluated at feature vector ${\mathbf{x}}$ is close to target variable $y$. Averaged over the unknown data distribution, we obtain the statistical loss ${\mathbb{E}}_{{\mathbf{x}},{\mathbf{y}}}[\ell(f({\mathbf{x}}),y)]$. It is commonplace to add a Tikhonov regularizer, yielding the regularized loss ${\mathbb{E}}_{{\mathbf{x}},{\mathbf{y}}}[\ell(f({\mathbf{x}}),y)]+(\lambda/2)\|f\|^{2}_{{\mathcal{H}}}$. Within the non-parametric setting, the regularizer permits application of the Representer Theorem [27], as discussed later in this paper. The optimal (centralized) function is defined as

[TABLE]

In this work, we focus on extensions of the formulation in (1) to the case where data is scattered across an interconnected network that represents, for instance, robotic teams, communication systems, or sensor networks. To do so, we define a symmetric, connected, and directed network ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ with $|{\mathcal{V}}|=V$ nodes and $|{\mathcal{E}}|=E$ edges. The neighborhood of agent $i$ is denoted by $n_{i}:=\{j:(i,j)\in{\mathcal{E}}\}$. Each agent $i\in{\mathcal{V}}$ observes a local data sequence as realizations $({\mathbf{x}}_{i,t},y_{i,t})$ from random pair $({\mathbf{x}}_{i},y_{i})\in{\mathcal{X}}\times{\mathcal{Y}}$ and seeks to learn a optimal regression function $f_{i}$. Local samples are assumed to be i.i.d. from a local distribution $({\mathbf{x}}_{i,t},y_{i,t})\sim\mathbb{P}_{i}({\mathbf{x}}_{i},y_{i})$ associated with node $i$. This setting may be encoded by associating to each node $i$ a convex loss functional $\ell_{i}:{\mathcal{H}}\times{\mathcal{X}}\times{\mathcal{Y}}\rightarrow{\mathbb{R}}$ that quantifies the merit of the estimator $f_{i}({\mathbf{x}}_{i})$ evaluated at feature vector ${\mathbf{x}}_{i}$, and defining the goal for each node as the minimization of the common global loss

[TABLE]

Subsequently, define the space ${\mathcal{H}}^{V}$ whose elements are stacked functions ${\mathbf{f}}(\cdot)=[f_{1}(\cdot);\cdots;f_{V}(\cdot)]$ that yield vectors of length $V$ when evaluated at local random vectors as ${\mathbf{f}}({\mathbf{x}})=[f_{1}({\mathbf{x}}_{1});\cdots;f_{V}({\mathbf{x}}_{V})]\in{\mathbb{R}}^{V}$. Moreover, define the stacked random vectors ${\mathbf{x}}=[{\mathbf{x}}_{1};\cdots;{\mathbf{x}}_{V}]\in{\mathcal{X}}^{V}\subset{\mathbb{R}}^{Vp}$ and ${\mathbf{y}}=[y_{1};\cdots y_{V}]\in{\mathbb{R}}^{V}$ that may represent $V$ labels or physical measurements. Note that under consensus constraints of the form $f_{i}=f_{j}$ for all $(i,j)\in{\mathcal{E}}$ and for a connected network, the solutions to (1) and (2) coincide.

However, as has been recently shown for the linear models, compelling all nodes to make common decisions may ignore local differences in their data streams, and in particular, yields a sub-optimal solution with respect to their distinct data [10]. In other words, it is indeed possible to boost the statistical accuracy of local estimates by relaxing the consensus constraint, while still maintaining information exchange between neighbors. In the present work, we achieve this enhancement by incentivizing the agents to coordinate without forcing their estimators to completely coincide.

To this end, we consider a convex (symmetric) proximity function $h_{ij}:{\mathbb{R}}\times{\mathbb{R}}\rightarrow{\mathbb{R}}_{+}$ for $(i,j)\in{\mathcal{E}}$. A node $i\in{\mathcal{V}}$ ensures proximity to its neighbors by requiring that $f_{j}$ for $j\in n_{i}$ be such that $H_{ij}(f_{i},f_{j}):={\mathbb{E}}_{{\mathbf{x}}_{i}}[h_{ij}(f_{i}({\mathbf{x}}_{i}),f_{j}({\mathbf{x}}_{i}))]\leq\gamma_{ij}$ 111Due to symmetry, we assume that the proximity function is of the form $h_{ij}=\frac{1}{2}\tilde{h}_{ij}$ for a convex $\tilde{h}_{ij}$ similarly defined. Subsequently, we suppress the tilde notation for notational simplicity, which also subsumes the scaling factor of $2$ due to constraints being doubly defined. for a given tolerance $\gamma_{ij}\geq 0$. The overall stochastic network optimization problem can be written as the

[TABLE]

Observe that if $h_{ij}(f_{i}({\mathbf{x}}_{i}),f_{j}({\mathbf{x}}_{i}))=|f_{i}({\mathbf{x}}_{i})-f_{j}({\mathbf{x}}_{i})|$ and $\gamma_{ij}=0$, the problem (3) specializes to online consensus optimization in RKHS, which was solved in an approximate manner using penalty methods in [35]. Here, we seek to obtain the exact optimal solutions to (3). In particular, we seek a solution to (3) whose optimality gap and constraint violation decay with the number of iterations $T$. In the subsequent section, we shift focus to doing so based upon Lagrange duality [20]. We specifically focus on distributed online settings where nodes do not know the distribution of the random pair $({\mathbf{x}}_{i},y_{i})$ but observe local independent samples $({\mathbf{x}}_{i,n},y_{i,n})$ sequentially.

II-A Function Estimation in RKHS

The optimization problem in (1), and hence (3), is intractable in general, since it defines a variational inference problem integrated over the unknown joint distribution $\mathbb{P}({\mathbf{x}},y)$. However, when ${\mathcal{H}}$ is equipped with a reproducing kernel $\kappa:{\mathcal{X}}\times{\mathcal{X}}\rightarrow{\mathbb{R}}$ (see [26]), a function estimation problem of the form (1) reduces to a parametric form via the Representer Theorem [37]. Thus, we restrict $\mathcal{H}$ to be a RKHS, i.e., for $\tilde{f}:{\mathcal{X}}\rightarrow{\mathbb{R}}$ in ${\mathcal{H}}$, it holds that

[TABLE]

for all ${\mathbf{x}}_{i}\in{\mathcal{X}}$. Here $\langle\cdot,\cdot\rangle_{{\mathcal{H}}}$ denotes the Hilbert inner product for ${\mathcal{H}}$. Further assume that the kernel is positive semidefinite, i.e. $\kappa({\mathbf{x}}_{i},{\mathbf{x}}_{i}^{\prime})\geq 0$ for all ${\mathbf{x}}_{i},{\mathbf{x}}_{i}^{\prime}\in{\mathcal{X}}$.

For kernelized empirical risk minimization (ERM), under suitable regularization, Representer Theorem [37] establishes that the optimal $\tilde{f}$ in hypothesized function class ${\mathcal{H}}$ admits an expansion in terms of kernel evaluations only over samples

[TABLE]

where ${\mathbf{w}}_{i}=[w_{i,1},\cdots,w_{i,N}]^{T}\in{\mathbb{R}}^{N}$ denotes a set of weights. Here $N$ in (5) is referred to as the model order. For ERM the model order and sample size are equal.

Remark 1

(Empirical Function Representation)Suppose, for the moment, that we have access to $N$ i.i.d. realizations of the random pairs $({\mathbf{x}}_{i},y_{i})$ for each agent $i$ such that the expectation in (3) is computable, and we further ignore the proximity constraint. Defining $S^{k}({\mathbf{f}}):=\sum_{i\in{\mathcal{V}}}\ell(f_{i}({\mathbf{x}}_{i,k}),y_{i,k})+\frac{\lambda}{2}\|f_{i}\|^{2}_{{\mathcal{H}}}$, the objective in (3) becomes:

[TABLE]

From the Representer Thm. [cf. (5)], (6) can be rewritten as

[TABLE]

where we have defined the kernel matrix ${\mathbf{K}}_{i}\in{\mathbb{R}}^{N\times N}$, with entries given by the kernel evaluations between ${\mathbf{x}}_{i,m}$ and ${\mathbf{x}}_{i,n}$ as $[{\mathbf{K}}_{i}]_{m,n}=\kappa({\mathbf{x}}_{i,m},{\mathbf{x}}_{i,n})$. We further define the vector of kernel evaluations $\boldsymbol{\kappa}_{{\mathbf{X}}_{i}}(\cdot)=[\kappa({\mathbf{x}}_{i,1},\cdot)\ldots\kappa({\mathbf{x}}_{i,N},\cdot)]^{T}$ related to the kernel matrix as ${\mathbf{K}}_{i}=[\boldsymbol{\kappa}_{{\mathbf{X}}_{i}}({\mathbf{x}}_{i,1})\ldots\boldsymbol{\kappa}_{{\mathbf{X}}_{i}}({\mathbf{x}}_{i,N})]$, whose dictionary of associated training points is defined as ${\mathbf{X}}_{i}=[{\mathbf{x}}_{i,1},\ \ldots\ ,{\mathbf{x}}_{i,N}]$. The Representer Theorem allows us to transform a nonparametric infinite dimensional optimization problem in ${\mathcal{H}}^{V}$ (6) into a finite $NV$-dimensional parametric problem (7). Thus, for ERM, the RKHS permits solving nonparametric regression problems as a search over ${\mathbb{R}}^{VN}$ for a set of coefficients.

However, to solve problems of the form (6) when training examples $({\mathbf{x}}_{i,k},y_{i,k})$ become sequentially available or their total number $N$ is not finite, the objective in (6) becomes an expectation over random pairs $({\mathbf{x}}_{i},y_{i})$ as [26]

[TABLE]

The Representer Theorem, generalized to expected value minimization, i.e., in [38], implies that the basis expansion is over a countably infinite index set ${\mathcal{I}}$. That is, as the data sample size $N\rightarrow\infty$, the representation of $f_{i}$ becomes infinite as well. Our goal is to solve (1) in an approximate manner such that each $\tilde{f}_{i}$ admits a finite representation near $\tilde{f}_{i}^{\star}$, while satisfying the proximity constraints as mentioned in (3), omitted for the sake of discussion between (6) - (1).

One wrinkle in the story is that the Representer Theorem in its vanilla form [37] does not apply to constrained problems (3). However, recently, it has been generalized to the Lagrangian of constrained problems in RKHS [15][Thm. 1]. To this end, some preliminaries are required. Let us collect the thresholds $\gamma_{ij}$ into a vector ${\boldsymbol{\gamma}}\in{\mathbb{R}}^{\sum_{i}|n_{i}|}$ where $|n_{i}|$ is the number of neighbors of node $i$. Likewise, let us compact the the functions $\{H_{ij}\}$ into a vector-valued function ${\mathbf{H}}:{\mathcal{H}}^{V}\rightarrow{\mathbb{R}}^{\sum_{i}|n_{i}|}$. Associate a non-negative Lagrange multiplier $\mu_{ij}$ for each constraint in (3), and collect these Lagrange multipliers into a vector $\boldsymbol{\mu}\in{\mathbb{R}}^{\sum_{i}|n_{i}|}$. The Lagrangian of (3) is therefore given by $\tilde{\mathcal{L}}({\mathbf{f}},\boldsymbol{\mu})=S({\mathbf{f}})+\langle\boldsymbol{\mu},{\mathbf{H}}({\mathbf{f}})-{\boldsymbol{\gamma}}\rangle$. Throughout, we assume Slater’s condition, implying strong duality [39]. Hence the optimal of (3) is identical to that of the primal-dual optimal pair $({\mathbf{f}}^{\star},\boldsymbol{\mu}^{\star})$ of the saddle-point problem

[TABLE]

Function Representation Now, we establish the Representer Theorem for applying to a variant of the Lagrangian (9). Consider the empirical approximation of $\tilde{\mathcal{L}}({\mathbf{f}},\boldsymbol{\mu})$ with training set of node $i$ defined as ${\mathcal{S}}_{i}=\{({\mathbf{x}}_{i,1},{\mathbf{y}}_{i,1}),\dots,({\mathbf{x}}_{i,N},{\mathbf{y}}_{i,N})\}$ over $N$ samples. The empirical version of (9) over samples ${\mathcal{S}}:=\{{\mathcal{S}}_{1},\dots,{\mathcal{S}}_{V}\}$ is written as:

[TABLE]

where ${\mathcal{L}}^{e}({\mathbf{f}},\boldsymbol{\mu})$ the empirical form of the Lagrangian:

[TABLE]

where ${\mathbf{H}}^{k}({\mathbf{f}})$ collects the functions $H_{ij}^{k}(f_{i},f_{j}):=h_{ij}(f_{i}({\mathbf{x}}_{i,k}),f_{j}({\mathbf{x}}_{i,k}))$ over all $i\in{\mathcal{V}}$ and $j\in n_{i}$. Now, with this empirical formulation, we generalize the Representer Theorem for constrained settings to the multi-agent problem (11) as a corollary of [15][Thm. 1] with the proof provided in Appendix C of the supplementary material.

Corollary 1

Let ${\mathcal{H}}$ be a RKHS equipped with kernel $\kappa$ and ${\mathcal{S}}$ be the training data of the network. Each function $i$ that is a primal minimizer of (11) takes the form $f_{i}^{\star}=\sum_{k=1}^{N}w_{i,k}\kappa({\mathbf{x}}_{i,k},.)$ where $w_{i,k}\in{\mathbb{R}}$ are coefficients.

Next, we shift to solving (3) in distributed online settings where nodes do not know the distribution of the random pair $({\mathbf{x}}_{i},y_{i})$ but observe local samples $({\mathbf{x}}_{i,k},y_{i,k})$ sequentially, through use of the Representer Theorem as stated in Corollary 1 that makes the function parameterization computable.

With the Representer theorem in place as explained above, we are ready to present the conservative version of (3) to vanish the long term constraint violation. To do so, we add $\nu$ to the constraint in (3) to reformulate the problem as follows

[TABLE]

By modifying (12), we consider a stricter constraint than (3). Doing so permits us to establish that approximate algorithmic solutions to (12) ensure the constraints in (3) may be exactly satisfied. Moreover, we are able to do so while tightening existing bounds on the the sub-optimality in [9] in Sec. IV.

III Algorithm Development

In this section, we develop an online and decentralized algorithm for (12) when $\{f_{i}\}_{i\in{\mathcal{V}}}$ belong to a RKHS. Begin by defining the augmented Lagrangian relaxation of (12):

[TABLE]

where ${\mathbf{1}}$ denotes the all-one vector of appropriate size and the parameters $\delta,\nu>0$. The dual regularization term in (13) is included in the design to control the violation of non-negative constraints on the dual variable over time $t$. For future reference, we denote ${\mathcal{L}}^{s}({\mathbf{f}},\boldsymbol{\mu})$ as the standard Lagrangian, which is (13) with $\delta=0$. We consider the stochastic approximation of (13) evaluated at sample $({\mathbf{x}}_{i,t},y_{i,t})$,

[TABLE]

With this definition, we propose applying primal-dual method to (14) – see [33]. To do so, we first require the functional stochastic gradient of (14) evaluated at a sample point $({\mathbf{x}}_{t},{\mathbf{y}}_{t})$. Begin by considering the local loss term in (14) :

[TABLE]

where we have applied the chain rule. Now, define the short-hand notation $\ell_{i}^{\prime}(f_{i}({\mathbf{x}}_{i,t}),y_{i,t}):={\partial\ell_{i}(f_{i}({\mathbf{x}}_{i,t}),y_{i,t})}/{\partial f_{i}({\mathbf{x}}_{i,t})}$ for the derivative of $\ell_{i}(f({\mathbf{x}}_{i,t}),y_{i,t})$ with respect to its first scalar argument $f_{i}({\mathbf{x}}_{i,t})$ evaluated at ${\mathbf{x}}_{i,t}$.

To evaluate the second term on the right-hand side of (15), differentiate both sides of the expression and use the reproducing property of the kernel with respect to $f_{i}$ to obtain

[TABLE]

Substitute the kernel at ${\mathbf{x}}_{i,t}$ on the right of (16) into (15):

[TABLE]

where we apply analogous logic as that which yields (15) to (17). To simplify, define the $V$-fold stacking of (17) as

[TABLE]

With these definitions, saddle point method on the augmented Lagrangian (13), which operates by alternating primal/dual stochastic gradient descent/ascent steps, is given as [33]:

[TABLE]

We note that local losses $\ell$ are assumed differentiable; however, the technical development easily extends to non-differentiable cases and the theoretical analysis ends up being very similar [40, Eq. (9)]. If a unique gradient is not available is practice, $\partial$ denotes a subgradient and we choose it arbitrarily from the set of available subgradients to perform the update in (19). Moreover, we require the step-size $\eta<1/\lambda$ for regularizer $\lambda>0$ in (1). Observe that (19) decouples by agent $i\in{\mathcal{V}}$:

[TABLE]

Note that the dual update in (20) is vector-valued, and defined for each edge $(i,j)\in{\mathcal{E}}$. Since the constraints involve only pairwise interactions between nodes $i$ and neighbors $j\in n_{i}$, the dual update separates along each edge $(i,j)$:

[TABLE]

The update of $\boldsymbol{\mu}_{t}$ is carried out by substituting (22) in (20) and using the fact that vector-wise projection is applied entry-wise and thus the individual local updates $\mu_{ij}$ can be written as

[TABLE]

The sequence of $({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ is initialized by ${\mathbf{f}}_{0}=0\in{\mathcal{H}}^{V}$ and $\boldsymbol{\mu}=0\in{\mathbb{R}}^{E}_{+}$. Using the Representer theorem, $f_{i,t}$ can be written in terms of kernels evaluated at past observations as

[TABLE]

We define ${\mathbf{X}}_{i,t}=[{\mathbf{x}}_{i,1},\ldots,{\mathbf{x}}_{i,t-1}]\in{\mathbb{R}}^{p\times(t-1)}$, $\boldsymbol{\kappa}_{{\mathbf{X}}_{i,t}}(\cdot)=[\kappa({\mathbf{x}}_{i,1},\cdot),\ \ldots\ ,\kappa({\mathbf{x}}_{i,t-1},\cdot)]^{T}$, and ${\mathbf{w}}_{i,t}=[w_{i,1},\ \ldots\,w_{i,t-1}]^{T}\in{\mathbb{R}}^{t-1}$ on the right-hand side of (24). Combining the update in (21) along with the kernel expansion in (24), implies that the primal functional stochastic descent step in ${\mathcal{H}}^{V}$ results in the following $V$ parallel parametric updates on both kernel dictionaries ${\mathbf{X}}_{i}$ and ${\mathbf{w}}_{i}$:

[TABLE]

From (25) we note that each time one more column gets added to the columns in ${\mathbf{X}}_{i,t}$, an instance of the curse of kernelization [41]. We define the number of data points, i.e., the number of columns of ${\mathbf{X}}_{i,t}$ at time $t$ as the model order. We note that for the update in (21), the model order is $t-1$ and it grows unbounded with iteration index $t$. This challenge often appears in connecting nonparametric statistics and optimization methods [26]. Next, motivated by [29], we now define compressive subspace projections of the function sequence defined by (19) to trade off memory and optimality.

Complexity Control via Subspace Projections To alleviate the aforementioned memory bottleneck, we project the function sequence (21) onto a lower dimensional subspace such that ${\mathcal{H}}_{\mathbf{D}}\subseteq{\mathcal{H}}$, where ${\mathcal{H}}_{\mathbf{D}}$ is represented by a dictionary ${\mathbf{D}}=[{\mathbf{d}}_{1},\ \ldots,\ {\mathbf{d}}_{M}]\in{\mathbb{R}}^{p\times M}$. Being specific, ${\mathcal{H}}_{\mathbf{D}}$ has the form ${\mathcal{H}}_{\mathbf{D}}=\{f\ :\ f(\cdot)=\sum_{n=1}^{M}w_{n}\kappa({\mathbf{d}}_{n},\cdot)={\mathbf{w}}^{T}\boldsymbol{\kappa}_{{\mathbf{D}}}(\cdot)\}=\text{span}\{\kappa({\mathbf{d}}_{n},\cdot)\}_{n=1}^{M}$, and $\{{\mathbf{d}}_{n}\}\subset\{{\mathbf{x}}_{u}\}_{u\leq t}$. For convenience we define $\boldsymbol{\kappa}_{{\mathbf{D}}}(\cdot)=[\kappa({\mathbf{d}}_{1},\cdot),\ldots,\kappa({\mathbf{d}}_{M},\cdot)]^{T}$, and ${\mathbf{K}}_{{\mathbf{D}},{\mathbf{D}}}$ as the resulting kernel matrix from this dictionary. In a similar manner, we define dictinaries ${\mathbf{D}}_{i,t}$ and subspace ${\mathcal{H}}_{{\mathbf{D}}_{i,t}}$ for each agent at time $t$. Similarly, the model order (i.e., the number of columns) of the dictionary ${\mathbf{D}}_{i,t}$ is denoted by $M_{i,t}$. We enforce function parsimony by selecting dictionaries ${\mathbf{D}}_{i}$ with $M_{i,t}<<{\mathcal{O}}(t)$ for each $i$ [29].

Now, we propose projecting the update in (21) to a lower dimensional subspace ${\mathcal{H}}_{{\mathbf{D}}_{i,t+1}}=\text{span}\{\kappa({\mathbf{d}}_{i,n},\cdot)\}_{n=1}^{M_{t+1}}$ as

[TABLE]

where we define the projection operator ${\mathcal{P}}_{{\mathbf{D}}_{i,t+1}}$ onto subspace ${\mathcal{H}}_{{\mathbf{D}}_{i,t+1}}\subset{\mathcal{H}}$ by the update (26).

Coefficient update. The update (26) may be rewritten in terms of updates for the dictionary ${\tilde{\mathbf{D}}}_{i,t+1}$ and coefficients ${\tilde{\mathbf{w}}}_{i,t+1}$:

[TABLE]

from the un-projected functional update step in (21). We denote the un-projected functional update as $\tilde{f}_{i,t+1}={f}_{i,t}-\eta\nabla_{f_{i}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ whose stacked functional version, using the stacked functional stochastic gradient (18), can be written as

[TABLE]

The number of columns of dictionary ${\tilde{\mathbf{D}}}_{i,t+1}$ is $M_{i,t}+1$, which is also the length of ${\tilde{\mathbf{w}}}_{i,t+1}$. For now, to simplify notation, we denote $\tilde{M}_{i,t+1}:=M_{i,t}+1$. For a given dictionary ${\mathbf{D}}_{i,t+1}$, projecting $\tilde{f}_{i,t+1}$ onto the subspace ${\mathcal{H}}_{{\mathbf{D}}_{i,t+1}}$ is equivalent to finding coefficients ${\mathbf{w}}_{i,t+1}$ w.r.t. dictionary ${\mathbf{D}}_{i,t+1}$ given as

[TABLE]

where we define the cross-kernel matrix ${\mathbf{K}}_{{\mathbf{D}}_{i,t+1},{\tilde{\mathbf{D}}}_{i,t+1}}$ whose $(n,m)^{\text{th}}$ entry is given by $\kappa({\mathbf{d}}_{i,n},{\tilde{\mathbf{d}}}_{i,m})$. The other kernel matrices ${\mathbf{K}}_{{\tilde{\mathbf{D}}}_{i,t+1},{\tilde{\mathbf{D}}}_{i,t+1}}$ and ${\mathbf{K}}_{{\mathbf{D}}_{i,t+1},{\mathbf{D}}_{i,t+1}}$ are defined similarly. The number of columns in ${\mathbf{D}}_{i,t+1}$ is ${M}_{i,t+1}$, while the number of columns in ${\tilde{\mathbf{D}}}_{t+1}$ [cf. (27)] is $\tilde{M}_{i,t+1}=M_{i,t}+1$. Next we see, how the dictionary ${\mathbf{D}}_{i,t+1}$ is obtained from ${\tilde{\mathbf{D}}}_{i,t+1}$.

Dictionary Update. The dictionary ${\mathbf{D}}_{i,t+1}$ is selected based upon greedy compression [28], i.e., ${\mathbf{D}}_{i,t+1}$ is formed from ${\tilde{\mathbf{D}}}_{i,t+1}$ by selecting a subset of $M_{i,t+1}$ columns from $\tilde{M}_{i,t+1}$ number of columns of ${\tilde{\mathbf{D}}}_{i,t+1}$ that best approximate $\tilde{f}_{i,t+1}$ in terms of Hilbert norm error, i.e., $\|f_{i,t+1}-\tilde{f}_{i,t+1}\|_{{\mathcal{H}}}\leq\epsilon$, where $\epsilon$ is the error tolerance, which may be done by kernel orthogonal matching pursuit (KOMP) [42] written as $(f_{i,t+1},{\mathbf{D}}_{i,t+1},{\mathbf{w}}_{i,t+1})=\text{KOMP}(\tilde{f}_{i,t+1},\tilde{{\mathbf{D}}}_{i,t+1},\tilde{{\mathbf{w}}}_{i,t+1},\epsilon)$.

We use a destructive variant of KOMP with pre-fitting as done in [29] and present it in Algorithm 2. This algorithm starts with the full dictionary ${\tilde{\mathbf{D}}}_{i,t+1}$ and sequentially removes the dictionary elements till the condition $\|f_{i,t+1}-\tilde{f}_{i,t+1}\|_{{\mathcal{H}}}\leq\epsilon$ is violated. In order to ensure the boundedness of the primal iterates in the subsequent section, we consider a variant of KOMP that explicitly enforces the projection to be contained within a finite Hilbert norm ball, which has the practical effect of thresholding the coefficient vector if it climbs above a certain large finite constant. Next we delve into the convergence of updates (26) and (23), summarized as Algorithm 1.

IV Convergence Analysis

In this section, we establish the convergence of Algorithm 1 by characterizing both objective sub-optimality and constraint violation in expectation. Before doing so, we define terms to clarify the analysis. Specifically, the projected functional stochastic gradient associated with (26) is defined as

[TABLE]

Using (30), the update (26) can be rewritten as $f_{i,t+1}=f_{i,t}-\eta\tilde{\nabla}_{f_{i}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$. We stack the projected stochastic functional gradient $\tilde{\nabla}_{f_{i}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ and define the stacked version

$\tilde{\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})=[\tilde{\nabla}_{f_{1}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t}),\dots,\tilde{\nabla}_{f_{V}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})]$. Using this stacked gradient, the update (26) then takes the form

[TABLE]

Next, we state the assumptions required for the convergence.

Assumption 1

The feature space ${\mathcal{X}}\subset{\mathbb{R}}^{p}$ and target domain ${\mathcal{Y}}\subset{\mathbb{R}}$ are compact, and the kernel map may be bounded as

[TABLE]

Assumption 2

The local losses $\ell_{i}(f_{i}({\mathbf{x}}),y)$ are convex and non-differentiable with respect to the first (scalar) argument $f_{i}({\mathbf{x}})$ on ${\mathbb{R}}$ for all ${\mathbf{x}}\in{\mathcal{X}}$ and $y\in{\mathcal{Y}}$. Moreover, the instantaneous losses $\ell_{i}:{\mathcal{H}}\times{\mathcal{X}}\times{\mathcal{Y}}\rightarrow{\mathbb{R}}$ are $C_{i}$-Lipschitz continuous

[TABLE]

Further denote $C:=\max_{i}C_{i}$ (largest modulus of continuity).

Assumption 3

The constraint functions $h_{ij}$ for all $(i,j)\in{\mathcal{E}}$ are all uniformly $L_{h}$-Lipschitz continuous in its first (scalar) argument; i.e., for any $z$, $z^{\prime}\in{\mathbb{R}}$, there exist constant $L_{h}$, s.t.

[TABLE]

and is also convex w.r.t the first argument $z$.

Assumption 4

There exists ${\mathbf{f}}^{\dagger}$ such that for all $(i,j)\in{\mathcal{E}}$, we have $h_{ij}(f^{\dagger}_{i},f^{\dagger}_{j})+\xi\leq\gamma_{ij}$, for some $\xi>0$, which implies that the constraint is strictly satisfied.

Assumption 5

The functions $f_{i,t+1}$ output from KOMP have Hilbert norm bounded by $R_{{\mathcal{B}}}<\infty$. Also, the optimal $f_{i}^{\star}$ lies in the ball ${\mathcal{B}}$ with radius $R_{{\mathcal{B}}}$.

Often, Assumption 1 holds by data domain itself. Assumptions 2 and 3 ensure that the constrained stochastic optimization problem is convex and the objective as well as constraint function satisfies the Lipschitz property, which are typical of first-order methods [40]. Moreover, the analysis in this work easily extends to the case where $\ell$ is non-differentiable. Assumption 4, i.e., Slater’s condition, ensures the feasible set of (3) is non-empty, and is standard in primal-dual methods [39]. Assumption 5 ensures that the algorithm iterates and the optimizer are finite, and their domains overlap. It may be explicitly enforced by dividing the norm of the coefficient vector output from KOMP by a large constant.

We begin with establishing an upper bound on the memory order of function $f_{i,t}$ obtained from Algorithm 1. Hence, using Assumption 2 and 3 we present the model order theorem.

Theorem 1

Let $f_{i,t}$ denote the function sequence of agent $i$ at $t$th instant generated from Algorithm 1 with dictionary ${\mathbf{D}}_{i,t}$. Denote $M_{i,t}$ as the model order representing the number of dictionary elements in ${\mathbf{D}}_{i,t}$. Then with constant step size $\eta$ and compression budget $\epsilon$, for a Lipschitz Mercer kernel $\kappa$ on a compact set ${\mathcal{X}}\subset{\mathbb{R}}^{p}$, there exists a constant $\beta$ such that for any training set $\{{\mathbf{x}}_{i,t}\}_{t=1}^{\infty}$, $M_{i,t}$ satisfies

[TABLE]

where $\alpha=\epsilon/\eta$, ${R_{M_{i,t}}}=C+L_{h}ER_{i,t}$ and $R_{i,t}=\max_{j\in n_{i}}|\mu_{ij,t}|$. The model order of the multi-agent system is then $M_{t}=\sum_{i=1}^{N}M_{i,t}$.

Thm. 1 establishes an explicit non-asymptotic bound on the model complexity as a function of problem parameters as well as algorithm parameters $\epsilon$ and $\eta$. The proof of Thm. 1 is provided in Appendix A. Observe from (35) that the model order depends on the algorithm parameters only through the ratio $\alpha=\eta/\epsilon$. In other words, by keeping $\alpha$ large, it is possible to ensure that the model order remains small. On the other hand, if we keep $\alpha$ small, the bound in (35) would be large. It is remarked that the non-asymptotic characterization of the model complexity in the distributed non-parametric setting in Thm. 1 is the first of its kind. Existing results, such as in [29], only provide asymptotic guarantees. As an implication of Thm. 1, observe that to ensure that the model order remains bounded by a given $M$, we must choose $\alpha={\mathcal{O}}(M^{-1/2p})$. Intuitively, if the model order is fixed a priori, the optimality gap cannot decay to zero in general, regardless of the number of iterations the algorithm is run. The resulting error floor introduced due to such a budget restriction will be characterized subsequently.

Via Lemma 1–4 (see the supplementary material), we establish our central result, which is the mean convergence of Algorithm 1 in terms of sub-optimality and feasibility.

Theorem 2

Suppose Assumptions 1-5 hold and $({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ be the primal-dual sequence of Algorithm 1 under constant step-size $\eta=T^{-1/2}$, and $\nu=\zeta T^{-1/2}+\Lambda\alpha$, where \zeta\geq\frac{1}{2}\bigg{[}R_{{\mathcal{B}}}^{2}+(1+\delta)\bigg{(}2+2\Big{(}\frac{{4}VR_{{\mathcal{B}}}(CX+\lambda R_{{\mathcal{B}}})}{\xi}\Big{)}^{2}\bigg{)}+K\bigg{]} and $\Lambda\geq 4VR_{{\mathcal{B}}}$ and ${K=8VX^{2}C^{2}+4V\lambda^{2}\cdot R_{{\mathcal{B}}}^{2}+2EK_{1}+2EL_{h}^{2}X^{2}\cdot R_{{\mathcal{B}}}^{2}}$.

1. (i)

The average expected sub-optimality decays as

[TABLE] 2. (ii)

Moreover, the average of aggregate constraint is met, i.e.,

[TABLE]

The full proof of Thm. 2 is provided in Appendix B. Thm. 2 can be viewed an explicit non-asymptotic characterization of the optimality gap as a function of the iterations $T$ and the model complexity parameter $\alpha$. Recall from Thm. 1 that if $M$ is fixed a priori, it amounts to fixing $\alpha={\mathcal{O}}(M^{-1/2p})$. In other words, a fixed model order $M$ translates to a floor of ${\mathcal{O}}(M^{-1/2p})$ on the optimality gap.

In the extreme case when the model order is allowed to grow indefinitely, the optimality gap does decay as ${\mathcal{O}}(T^{-1/2})$ while the constraint violation is zero. The rate provided here therefore improves upon the existing results in [9, 10, 8], where the optimality gap is shown to decay as ${\mathcal{O}}(T^{-1/4})$. Further, the rate provided here matches the best possible rate achievable for unconstrained general convex optimization problems. The specific approach responsible for allowing us to achieve zero constraint violation can be intuitively explained as follows: instead of the original problem, we actually solve a $\nu$-tightened problem with a smaller constraint set. As long as the original problem is strongly feasible and we set $\nu$ appropriately, such a tightening only leads to ${\mathcal{O}}(T^{-1/2})$ suboptimality, so that the overall optimality gap only changes by a constant factor.

It is instructive to express the results in Thm. 2 in terms of the iteration and model complexity parameters. Suppose that the goal is to achieve an desired optimality gap of $\varepsilon$. The following corollary, whose proof is provided in Appendix I of the supplementary material, characterizes the iteration and storage complexities of the proposed algorithm.

Corollary 2

For step-size $\eta=T^{-1/2}$ and dual regularizer $\nu=\zeta T^{-1/2}+\Lambda\alpha$, under Assumptions 1-5 hold, the optimality gap to be less than $\varepsilon$, i.e., $\frac{1}{T}\sum_{t=1}^{T}{\mathbb{E}}[S({\mathbf{f}}_{t})-S({\mathbf{f}}^{\star})]\leq\varepsilon$ requires ${\mathcal{O}}(1/\varepsilon^{2})$ number of iterations and the model complexity bound of ${\mathcal{O}}(1/\varepsilon^{2p})$ with the parameter $\alpha\leq\frac{\varepsilon}{2(4VR_{{\mathcal{B}}}+Q\Lambda)}$, where $Q=4VR_{{\mathcal{B}}}(CX+\lambda R_{{\mathcal{B}}})/\xi$.

As expected, the iteration complexity matches the state-of-the-art first order oracle complexity of SGD for the unconstrained convex setting. In addition, the model order complexity bound characterizes the corresponding storage requirements.

V Numerical Results

In this section, we apply the proposed algorithm to solve a spatial temporal random field estimation problem and an another problem of inferring from oceanographic data.

V-A Spatio-temporal Random Field Estimation

The estimation of a spatio-temporal field using a set of sensors spread across a region with required level of accuracy is an central challenge in wireless sensor networks (WSNs) [5]. We model this problem by considering the problem of estimating a temporally varying spatial planar correlated Gaussian random field in a given region ${\mathcal{G}}\subset{\mathbb{R}}^{2}$ space. A spatial temporal random field is a random function of the spatial components $u$ (for x-axis) and $z$ (for y-axis) across a region ${\mathcal{G}}$ and time. Moreover, the random field is parameterized by its correlation matrix ${\mathbf{R}}_{s}$, which depends on the location of the sensors. Each element of $[{\mathbf{R}}_{s}]_{ij}$ is assumed to have a structure of the form $\Omega(l_{i},l_{j})=e^{-\|l_{i}-l_{j}\|}$, where $l_{i}$ and $l_{j}$ are the respective locations of sensor $i$ and $j$ in region ${\mathcal{G}}$ [5]. From this correlation, note that the nodes close to each other have high correlation whereas nodes located far away are less correlated, meaning that observations collected from the nearby nodes are more relevant than observations from distant nodes.

We experimented with a sensor network with $V=40$ nodes spatially distributed in a $100\times 100$ meter square area. Each node $i$ collects the observation $y_{i,t}$ at time instant $t$. In the collected data, $y_{i,t}$ denotes the noisy version of the original field $s_{i,t}$ at node $i$ for time instant $t$. The observation model is given by $y_{i,t}=s_{i,t}+n_{i,t}$, where $n_{i,t}\sim{\mathcal{N}}(0,0.5)$ is i.i.d. Each node seeks to sequentially minimize its local loss, i.e.,

[TABLE]

where $\hat{s}_{i,t}$ is the estimated value of actual field $s_{i,t}$ and $\phi$ is the tuning parameter. In the experiments, we choose $\delta=10^{4}$ to ensure that the local loss function is Lipschitz and the gradient is bounded by $2\delta$. The instantaneous observation ${\mathbf{s}}_{t}$ across the network is given by {\mathbf{s}}_{t}=\boldsymbol{\pi}+{\mathbf{C}}^{T}\big{(}\mathbf{1}\sin(\omega t)+{\mathbf{v}}_{t}\big{)}, where $\mathbf{1}$ is a vector of ones of length $V$, $\sin(\omega t)$ is a sinusoidal with angular frequency $\omega=2$, $\boldsymbol{\pi}=\{1/V,2/V,\dots,1\}$ is a fixed mean vector of length $V$, ${\mathbf{C}}$ is the Cholesky factorization of the correlation ${\mathbf{R}}_{s}$, and ${\mathbf{v}}_{t}\sim{\mathcal{N}}(\mathbf{0},0.1\mathbf{I})$, where $I$ denotes $40\times 40$ identity matrix. We select tolerance parameter to be $\gamma_{ij}=\Omega(l_{i},l_{j})$. We solve the problem (3) of minimizing the regularized loss function over $f_{i}$ where we learn the function $f_{i}(t)$, which is the function approximation to the actual field, $s_{i,t}$, i.e., we solve an online decentralized regression (curve fitting) problem.

We run Algorithm 1 to generate local functions $f_{i}(t)$ to track ground-truth $y_{i,t}$. We select parameters: parsimony constant $P=8$ [Thm. 2], Gaussian kernel with bandwidth $\sigma=0.05$ such that we capture the variation of sinusoidal function with angular frequency of $2$, and primal/dual regularizers $\lambda=10^{-5}$ and $\delta=10^{-5}$. We run the algorithm for $1500$ iterations with step-size $\eta=0.01$. For comparison purposes, we consider two other comparable techniques: consensus with kernels via penalty method [35] and the simplification of (3) to linear statistical models, i.e., $f_{i}({\mathbf{x}}_{i})={\mathbf{w}}_{i}^{T}\phi({\mathbf{x}})$ for some $d$ fixed-dimensional parameter vectors $\{{\mathbf{w}}_{i}\}_{i=1}^{N}$, which we call linear method [10]. For comparison with the penalty method, the penalty coefficient is $0.08$ which is tuned for the best performance, with all other parameters held fixed to Algorithm 1. For linear method, we consider three parametric models (which imply a different structural form for the feature map $\phi({\mathbf{x}})$): (a) Quadratic polynomial; (b) Cubic polynomial and (c) Sine polynomial (i.e., of the form $at+b\text{sin}(\omega t)$ where $a$ and $b$ are the model parameters and $\omega=1$ is the angular frequency).

Fig. 1 displays the results of this comparison. In particular, Fig. 1(a) compares the global network wide averaged loss for three different methods, which shows that the linear method is unable to effectively track the target variable the model behavior. The linear sine polynomial model (denoted as “Sine”) is closer to the target than the quadratic (denoted as “Degree 2”) or cubic polynomial (denoted as “Degree 3”). Algorithm 1 attains superior performance to the RKHS-based penalty method [35]. Moreofer, Fig. 1(b) demonstrates that Algorithm 1 achieves tighter constraint satisfaction relative to penalty method, and is comparable to primal-dual schemes for linear models. Doing so allows nodes estimates tune their closeness to neighbors through proximity tolerances $\gamma_{ij}$.

Fig. 1(c) plots the model order for primal-dual method and penalty method, which omits linear method plots because its a parametric method with fixed complexity equal to the parameter dimension $d$. Early on, primal-dual method (being an exact method) has higher complexity than penalty method. In steady state, Algorithm 1 and penalty method have comparable complexity. With a similar model complexity, we attain near-exact constraint satisfaction via primal-dual method as compared to penalty method. Overall, the model complexity of $30$ is orders of magnitude smaller than sample size $1500$.

V-B Inferring Oceanographic Data

Wireless sensor networks may also be used to monitor various environmental parameters, especially in oceanic settings. To this end, we associate each node in the network to an oceanic buoy tasked with estimating salinity and temperature when deployed at standard depths. Decentralization is advantageous here due to the fact that server stations are impractical at sea, and centralization may exceed the cost of computation per node [5]. Thus, we run Algorithm 1 on the World Oceanic Database [43], obtained from multiple underwater sensors in the Gulf of Mexico. In this Regional Climatology data set, temperature, salinity, oxygen, phosphate, silicate, and nitrate are recorded at various depth levels.

We restrict focus to temperature and salinity parameters at different locations with varying depths, during the winter time-period. The readings of the climatological fields are obtained for a particular latitude and longitude at standard depths starting from [math] meters to $5000$ meters. The latitude and longitude specifies the node (sensor) location. Similar readings are obtained for various locations spanning the water body. The experiment is carried out considering $50$ nodes, where edges are determined by measuring the distance between two nodes, and drawing an edge to a particular node if its distance is less than $1000$ kilometers away. The proximity parameter $\gamma_{ij}$ is obtained by evaluating $\exp(-\text{dist}(i,j)/1000)$, where $\text{dist}(i,j)$ denotes distance nodes in kilometers.

We use Algorithm 1 to estimate the value of climatological field $y_{i}$ at a depth $d_{i}$ such as salinity or temperature at each node $i$. We solve problem (3) by minimizing the regularized quadratic loss between estimated climatological field and observed climatological field $y_{i}$ over function $f_{i}$. A key benefit of doing so is the ability to interpolate missing values, which arise due to, e.g., limited battery or bandwidth.

Bandwidth Adaptation We consider an extended implementation of Algorithm 1 for the purpose of experimentation, where kernel hyper-parameters may be selected adaptively via online maximum likelihood, rather than fixed a priori [36]. The motivation is that in this oceanic setting, spatial fields exhibit scale heterogeneity that can be exploited for boosting accuracy. In the decentralized online setting, we run a bandwidth at step $(t+1)$ for the Gaussian kernel of agent $i$ as:

[TABLE]

where ${\mathbf{d}}_{l}$ and ${\mathbf{d}}_{k}$ are the dictionary elements of ${\tilde{\mathbf{D}}}_{i,t+1}$ [cf (27)]. The utility of this adaptive bandwidth selection is investigated next. For better interpretability, we denote the centralized method with fixed kernel bandwidth as “N=1, fix” and distributed method with fixed and adaptive bandwidth as “N=50, fix” and “N=50, adap”.

V-B1 Temperature

Here we use Algorithm 1 for predicting the statistical mean of the temperature field of different nodes at varying depths. The real data obtained from the World Oceanic database has statistical mean of the temperature field. We run Algorithm 1 for $T=2430$ iterations with constant step-size $\eta=0.01$ and regularizers $\lambda=10^{-5}$, $\delta=10^{-5}$ with a Gaussian kernel with fixed bandwidth parameter $\sigma=50$. The adaptive scheme employs (39) with the same bandwidth initialization. The parsimony constant is fixed at two values, $P\in\{0.4,40\}$. The adaptive bandwidth case is studied only for $P=40$. The parsimony constant is set to $P=0.001$ for [29] to ensure comparably sized models across cases.

Fig. 2 shows the numerical experiment results for the ocean data. Fig. 2(a) demonstrates that the prediction error for test cases reduces with increasing samples, and illustrates that centralization is inappropriate here: local spatial variability of the field causes [29] to be outperformed by Algorithm 1 under both fixed and variable bandwidths. The variable bandwidth along with appropriately chosen parsimony constant and step size performs best in terms of model fitness. A similar trend may be gleaned from the plot of constraint violation over time in Fig. 2(b). Interestingly, the adaptive bandwidth scheme obtains more accurate model with comparable constraint violation, both of which are superior to centralized approaches, thus substantiating the experimental merits of (39). The evolution of a few random agents’ bandwidths is visualized in Fig. 2(c) – since (39) is a stochastic fixed point iteration, each agent’s bandwidth converges to a neighborhood.

V-B2 Salinity

Next we consider Algorithm 1 for predicting the mean salinity from various oceanic locations and depths. We set the primal and dual regularizer $\lambda=\delta=10^{-5}$, and run it for $T=2500$ iterations with constant step-size of $0.01$. We select two values of the parsimony constant $P\in\{0.4,40\}$. For the centralized case, we fix $P=0.001$ so that its complexity is comparable to the distributed approach to ensure a fair comparison. The bandwidth of the Gaussian kernel is fixed at $50$ and also considered adaptive for the simulation with initial bandwidth value set at $50$.

We display these results in Fig. 3. Fig. 3(a) demonstrates that learning a single function to fit all data is unable to filter out correlation effects and hence gives poor model fitness, as compared to fitting multiple $f_{i}$’s to different nodes considering proximity constraints. Moreover, the average model order for a single node for the distributed case is $39$ for $P=0.4$, thus giving an aggregate complexity of $1950$ for $50$ nodes. This is less than the centralized model order of $2372$ for a single function. Thus the distributed approach yields improved accuracy with reduced complexity. We note that increasing the parsimony constant results in worse model fit but saves complexity, yielding a tunable tradeoff between fitness and complexity. Similar to the temperature data, in Fig. 3(a) also we observe improvement in performance for adaptive bandwidth case with smaller average model order. In Fig. 3(b) we may observe that Algorithm 1 incurs attenuating constraint violation for both fixed and adaptive bandwidth case. Overall, complexity settles to around $39$, which is orders of magnitude smaller than the $2500$ sample size.

V-B3 Radial basis function (RBF) network

To evaluate the merit of the proposed decentralized projection scheme, we considered a variant of our approach which employs a fixed feature representation, i.e., linear regression over a nonlinear RBF basis whose points are selected according to a uniform grid across the feature space. Doing so may be interpreted as a compatible comparison to [10]. We call this scheme RBF, as it is a multi-agent proximity-constrained problem for learning the coefficients of a one-layer RBF network. We implement this scheme on the real ocean data set considered in the paper with the same parameters as considered for HALK.

The fixed model order for a single agent for the RBF case is fixed to the final settled model order obtained in case of our paper (for parsimony constant $P=40$), i.e., 26 and 31 for temperature and salinity field estimation. The fixed dictionary points for the RBF case are randomly selected from a uniform distribution. It can be observed from both Fig. 2(a) and Fig. 3(a) that for salinity and temperature field estimation, the RBF based algorithm performs inferior to our algorithm in the long run. IEarly on, RBF performs better since it has fewer parameters to learn due to its fixed basis, but over time the quality of its fixed basis is exceeded by the dynamically chosen one of HALK. The constraint violation for both the field estimation parameters are presented in Fig. 2(b) and Fig. 3(b). The important observation here is that our proposed algorithm discerns the number of dictionary points (model order) dynamically and it evolves with the system. However, in case of RBF based algorithm, we require to apriori fix the size of the RBF network for the algorithm to initialize and iterate, which in many cases is difficult to fix in advance.

VI Conclusion

We proposed learning in heterogeneous networks via constrained functional stochastic programming. We modeled heterogeneity using proximity constraints, which allowed agents to make decisions which are close but not necessarily equal. Moreover, motivated by their universal function approximation properties, we restricted focus to the case where agents decisions are defined by functions in RKHS. We formulated the augmented Lagrangian, and proposed a decentralized stochastic saddle point method to solve it. Since decision variables were functions belonging to RKHS, not vectors, we required generalizing the Representer Theorem to this setting, and further projecting the primal iterates onto subspaces greedily constructed from subsets of past observations.

The algorithm, Heterogeneous Adaptive Learning with Kernels (HALK), converges in terms of primal sub-optimality and satisfies the constraints on average. We further established a controllable trade-off between convergence and complexity. We validated HALK for estimating a spatial temporal Gaussian random field in a heterogeneous sensor network, and employed it to predict oceanic temperature and salinity from depth. We also considered a generalization where the kernel bandwidth of each agent’s function is allowed to vary, and observed better experimental performance as compared to the fixed bandwidth. In future work, we hope to relax communications requirements and allow asynchronous updates [44, 45], permit the learning rates to be distinct among agents, and obtain tighter dependence of the convergence results on the network data.

Appendix A Proof of Thm. 1

The proof is motivated from the derivation presented in [29, Thm. 4] and is presented here for the proposed algorithm. Consider the function iterates $f_{i,t}$ and $f_{i,t+1}$ of agent $i$ generated from Algorithm 1 at $t$-th and $(t+1)$-th instant. The function iterates $f_{i,t}$ and $f_{i,t+1}$ are parametrized by dictionary ${\mathbf{D}}_{i,t}$ and ${\mathbf{D}}_{i,t+1}$ and weights ${\mathbf{w}}_{i,t}$ and ${\mathbf{w}}_{i,t+1}$, respectively. The dictionary size corresponding to $f_{i,t}$ and $f_{i,t+1}$ in dictionary ${\mathbf{D}}_{i,t}$ and ${\mathbf{D}}_{i,t+1}$ are denoted by $M_{i,t}$ and $M_{i,t+1}$, respectively. The kernel dictionary ${\mathbf{D}}_{i,t+1}$ is formed from ${\tilde{\mathbf{D}}}_{i,t+1}=[{\mathbf{D}}_{i,t},\;\;{\mathbf{x}}_{i,t}]$ by selecting a subset of $M_{i,t+1}$ columns from $\tilde{M}_{i,t+1}=M_{i,t}+1$ number of columns of ${\tilde{\mathbf{D}}}_{i,t+1}$ that best approximate $\tilde{f}_{i,t+1}$ in terms of Hilbert norm error, i.e., $\|f_{i,t+1}-\tilde{f}_{i,t+1}\|_{{\mathcal{H}}}\leq\epsilon$, where $\epsilon$ is the error tolerance. Suppose the model order of function $f_{i,t+1}$ is less than equal to that of $f_{i,t}$, i.e., $M_{i,t+1}\leq M_{i,t}$, which holds when the stopping criteria of KOMP is violated for dictionary ${\tilde{\mathbf{D}}}_{i,t+1}$:

[TABLE]

where $\gamma_{j}$ is the minimal approximation error with dictionary element ${\mathbf{d}}_{i,j}$ removed from dictionary ${\tilde{\mathbf{D}}}_{i,t+1}$ defined as

[TABLE]

where ${\mathcal{I}}=\{1,\dots,M_{i,t}+1\}$.

Observe that (40) lower bounds the approximation error $\gamma_{M_{i,t}+1}$ of removing the most recently added feature vector ${\mathbf{x}}_{i,t}$. Thus if $\gamma_{M_{i,t}+1}\leq\epsilon$, then (40) is satisfied and the relation $M_{i,t+1}\leq M_{i,t}$ holds, implying the model order does not grow. Hence it is adequate to consider $\gamma_{M_{i,t}+1}$.

Using the definition of $\tilde{f}_{i,t+1}$ and denoting ${\mathcal{I}}^{\prime}\coloneqq{\mathcal{I}}\setminus\{M_{i,t}+1\}$, we write $\gamma_{M_{i,t}+1}$ as

[TABLE]

The minimizer of (42) is obtained for ${\mathbf{u}}^{*}$ is obtained via a least-squares computation, and takes the form:

[TABLE]

Further, we define $\ell_{i}^{\prime}(f_{i,t})$ as

[TABLE]

Substituting ${\mathbf{u}}^{*}$ from (43) into (42) and applying Cauchy-Schwartz with (A), we obtain $\gamma_{M_{i,t}+1}$ as

[TABLE]

It can be observed from the right hand side of (45) that the Hilbert-norm term can be replaced by using the definition of subspace distance from (117). Thus, we get

[TABLE]

Now for $\gamma_{M_{i,t}+1}\leq\epsilon$, the right hand side of (46) should also be upper bounded by $\epsilon$ and thus can be written as

[TABLE]

where we have divided both the sides by $|{\ell_{i}^{\prime}(f_{i,t})}|$. Note that if (47) holds, then $\gamma_{M_{i,t}+1}\leq\epsilon$ and since $\gamma_{M_{i,t}+1}\geq\min_{j}\gamma_{j}$, we may conclude that (40) is satisfied. Implying the model order at the subsequent steps does not grow, i.e., $M_{i,t+1}\leq M_{i,t}$.

Now, let’s take the contrapositive of the expression in (47) to observe that growth in the model order ($M_{i,t+1}=M_{i,t}+1$) implies that the condition

[TABLE]

holds. Therefore, every time a new point is added to the model, the corresponding kernel function, i.e., $\kappa({\mathbf{x}}_{i,t},\cdot)$ is at least $\frac{\epsilon}{\eta|{\ell_{i}^{\prime}(f_{i,t})}|}$ distance far away from every other kernel function in the current model defined by dictionary ${\mathbf{D}}_{i,t}$.

Now, to have a bound on the right-hand side term of (48), we bound the denominator of the right-hand side of (48). Thus, we upper bound $|{\ell_{i}^{\prime}(f_{i,t})}|$ as

[TABLE]

where in the first equality we have used the definition of $\ell_{i}^{\prime}(f_{i,t})$ from (A) and the second inequality is obtained by using triangle inequality. To obtain the last inequality in (A), we use Assumption 2 and 3, with $R_{i,t}=\max_{j\in n_{i}}|\mu_{ij,t}|$ and upper-bounded $|n_{i}|$ by the total number of edges $E$. Subsequently, we denote the right hand side of (A) as ${R_{M_{i,t}}}\coloneqq C+L_{h}ER_{i,t}$. Substituting into (48) yields

[TABLE]

Hence, the stopping criterion for the newest point is violated whenever it satisfies the condition, $\|\phi({\mathbf{x}}_{i,t})-\phi({\mathbf{d}}_{i,k})\|_{2}\leq\frac{\epsilon}{\eta{R_{M_{i,t}}}}$ for $k\in\{1,\dots,M_{i,t}\}$ meaning $\phi({\mathbf{x}}_{i,t})$ can be well approximated by $\phi({\mathbf{d}}_{i,k})$ with already existing point ${\mathbf{d}}_{i,k}$ in the dictionary. Now for the finite model order proof, we proceed in a manner similar to the proof of [46, Thm. 3.1]. Since feature space ${\mathcal{X}}$ is compact and $\phi({\mathbf{x}})=\kappa({\mathbf{x}},\cdot)$ is continuous (as $\kappa$ is continuous), we can deduce that $\phi({\mathcal{X}})$ is compact. Therefore, the covering number (number of balls with radius $\varrho=\frac{\epsilon}{\eta{R_{M_{i,t}}}}$ to cover the set $\phi({\mathcal{X}})$) of set $\phi({\mathcal{X}})$ is finite [47]. The covering number of a set is finite if and only if its packing number (maximum number of points in $\phi({\mathcal{X}})$ separated by distance larger than $\varrho$) is finite. This means that the number of points in $\phi({\mathcal{X}})$ separated by distance $\varrho$ is finite. Note that KOMP retains points which satisfy $\|\phi({\mathbf{d}}_{i,j})-\phi({\mathbf{d}}_{i,k})\|_{2}>\varrho$, i.e., the dictionary points are $\varrho$ separated. Thus, when the packing number of $\phi({\mathcal{X}})$ with scale $\varrho$ is finite, the number of dictionary points is also finite.

From [46, Proposition 2.2], we know that for a Lipschitz continuous Mercer kernel $\kappa$ on a compact set ${\mathcal{X}}\subset{\mathbb{R}}^{p}$, there exists a constant $\beta$ depending upon the ${\mathcal{X}}$ and the kernel function such that for any training set $\{{\mathbf{x}}_{i,t}\}_{t=1}^{\infty}$ and any $\vartheta>0$, the number of elements in the dictionary satisfies

[TABLE]

From (50), observe $\vartheta=\frac{\epsilon}{\eta{R_{M_{i,t}}}}$. With (51) we may write $M_{i,t}\leq\beta({\eta{R_{M_{i,t}}}}/{\epsilon})^{2p}$, which is the required result stated in (35) of Thm. 1. To obtain the model order $M_{t}$ of the multi-agent system, we sum the model order of individual nodes across network, i.e., $M_{t}=\sum_{i=1}^{N}M_{i,t}$ as stated in Thm. 1. $\blacksquare$

Appendix B Proof of Thm. 2

Before discussing the proof, we introduce the following compact notations to make the analysis clear and compact. We further use the following short-hand notations to denote the expressions involving $h_{ij}(\cdot,\cdot)$ as

[TABLE]

Moreover, the $E$-fold stacking of the constraints across all the edges is denoted as ${\mathbf{G}}^{t}({\mathbf{f}})\coloneqq\text{vec}[g_{ij}(f({\mathbf{x}}_{t}))]$ and ${\mathbf{G}}({\mathbf{f}})\coloneqq{\mathbb{E}}_{{\mathbf{x}}_{t}}[{\mathbf{G}}^{t}({\mathbf{f}})]$. The intermediate results required for the proof are stated in Lemma 1-4 (detailed in the supplementary material). Consider the statement of Lemma 4 (c.f. (4)), expand the left hand side of (4) using the definition of (14), further utilizing the notation of $S^{t}({\mathbf{f}}_{t})$ stated in (3) and ${\mathbf{G}}^{t}(\cdot)$ in (52), we can write

[TABLE]

[TABLE]

where $\Delta_{t}=(\|{\mathbf{f}}_{t}-{\mathbf{f}}\|_{{\mathcal{H}}}^{2}-\|{\mathbf{f}}_{t+1}-{\mathbf{f}}\|_{{\mathcal{H}}}^{2}+\|\boldsymbol{\mu}_{t}-\boldsymbol{\mu}\|^{2}-\|\boldsymbol{\mu}_{t+1}-\boldsymbol{\mu}\|^{2})$. Next, let us take the total expectation on both sides of (B). From Lemma 2, substituting the upper bounds of $\|{\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})\|_{{\mathcal{H}}}^{2}$ and $\|\nabla_{\boldsymbol{\mu}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\mathbf{\boldsymbol{\mu}}_{t})\|^{2}$ to (B), the right hand side of (B) can be written as

[TABLE]

Since each individual $f_{i,t}$ and $f_{i}$ for $\in\{i,\dots,V\}$ in the ball ${\mathcal{B}}$ have finite Hilbert norm and is bounded by $R_{{\mathcal{B}}}$, the term $\|{\mathbf{f}}_{t}-{\mathbf{f}}\|_{{\mathcal{H}}}$ can be upper bounded by $2\sqrt{V}R_{{\mathcal{B}}}$. Next, we define $K\coloneqq 8VX^{2}C^{2}+4V\lambda^{2}\cdot R_{{\mathcal{B}}}^{2}+2EK_{1}+2EL_{h}^{2}X^{2}\cdot R_{{\mathcal{B}}}^{2}$. Now using the the bound of $\|{\mathbf{f}}_{t}-{\mathbf{f}}\|_{{\mathcal{H}}}$ and the definition of $K$, we can upper bound the expression in (B), and then collectively writing the left and right hand side terms together, we get

[TABLE]

where $C(\delta):=8VX^{2}L_{h}^{2}E+2E\delta^{2}\eta^{2}-\delta$. Next, we select the constant parameter $\delta$ such that $C(\delta)\leq 0$, which then allows us to drop the term involving $\|\boldsymbol{\mu}_{t}\|^{2}$ from the second expected term of right-hand side of (55). Further, take the sum of the expression in (55) over times $t=1,\dots,T$, assume the initialization ${\mathbf{f}}_{1}=0\in{\mathcal{H}}^{V}$ and $\boldsymbol{\mu}_{1}=0\in{\mathbb{R}}_{+}^{E}$, we obtain

[TABLE]

where we drop the negative terms remaining after the telescopic sum since $\|{\mathbf{f}}_{T+1}-{\mathbf{f}}\|_{{\mathcal{H}}}^{2}$ and $\|\boldsymbol{\mu}_{T+1}-\boldsymbol{\mu}\|^{2}$ are always positive. It can be observed from (56) that the right-hand side of this inequality is deterministic. We now take ${\mathbf{f}}$ to be the solution ${\mathbf{f}}_{\nu}^{\star}$ of (12), which in turn implies ${\mathbf{f}}_{\nu}^{\star}$ must satisfy the inequality constraint of (12). This means that ${\mathbf{f}}_{\nu}^{\star}$ is a feasible point, such that $\sum_{t=1}^{T}\sum_{(i,j)\in{\mathcal{E}}}\mu_{ij,t}(g_{ij}(f^{\star}_{i}({\mathbf{x}}_{i,t}),f^{\star}_{j}({\mathbf{x}}_{i,t}))+\nu)\leq 0$ holds. Thus we can simply drop this term in (56) and collecting the terms containing $\|\boldsymbol{\mu}\|^{2}$ together, we obtain

[TABLE]

where $z(\eta,T):=\frac{\delta\eta T}{2}+\frac{1}{2\eta}$. Next, we maximize the left-hand side of (B) over $\boldsymbol{\mu}$ to obtain the optimal Lagrange multiplier which controls the growth of the long-term constraint violation, whose closed-form expression is given by

[TABLE]

Therefore, selecting $\boldsymbol{\mu}=\bar{\boldsymbol{\mu}}$ in (B), we obtain

[TABLE]

Firstly, consider the objective error sequence {\mathbb{E}}\big{[}{S({\mathbf{f}}_{t})}-{S({\mathbf{f}}_{\nu}^{\star})}\big{]}, we observe from (B) that the second term present on the left-side of the inequality can be dropped without affecting the inequality owing to the fact that it is positive. So we obtain

[TABLE]

Using Lemma 1 and summing over $t=1,\dots,T$, we have

[TABLE]

Adding these inequalities, and then taking the average over $T$, we obtain

[TABLE]

Since $\epsilon\leq{2}R_{{\mathcal{B}}}$, we upper bound $\epsilon$ present in the second term of the left hand side of (B) by ${2}R_{{\mathcal{B}}}$. Next, using the definition $\alpha=\epsilon/\eta$ from Thm. 1, and setting step size $\eta=1/\sqrt{T}$ and $\nu=\zeta T^{-1/2}+\Lambda\alpha$ for some $(\zeta,\Lambda)>0$, we obtain

[TABLE]

which is the required result in (36).

Next, we establish the bound on the growth of the constraint violation. For this, we first denote ${\mathcal{L}}^{s}$ as the standard Lagrangian for (12) and write it for ${{\mathbf{f}}}$ and $\boldsymbol{\mu}$ as, ${\mathcal{L}}^{s}({{\mathbf{f}}},\boldsymbol{\mu})={S({\mathbf{f}})}+\langle\boldsymbol{\mu},{{\mathbf{G}}({\mathbf{f}})}+\nu{\mathbf{1}}\rangle$. The standard Lagrangian for ${({\mathbf{f}}_{t},\boldsymbol{\mu})}$ and $({\mathbf{f}},\boldsymbol{\mu}_{t})$ are defined similarly. Now using the expressions of ${\mathcal{L}}^{s}({\mathbf{f}}_{t},\boldsymbol{\mu})$ and ${\mathcal{L}}^{s}({\mathbf{f}},\boldsymbol{\mu}_{t})$, we rewrite (56) as

[TABLE]

Since $({\mathbf{f}}_{\nu}^{*},\boldsymbol{\mu}_{\nu}^{*})$ is the optimal pair for standard Lagrangian ${\mathcal{L}}^{s}({\mathbf{f}},\boldsymbol{\mu})$ of (12) and assuming $\mathbf{1}_{i}$ to be a vector of all zeros except the $i$-th entry which is unity, write for $\boldsymbol{\mu}=\mathbf{1}_{i}+\boldsymbol{\mu}_{\nu}^{*}$:

[TABLE]

The first equality in (B) is written by using the definition of standard Lagrangian ${\mathcal{L}}^{s}$ and ${\mathbf{G}}$ denotes the stacking of the constraints of all edges as defined in the paragraph following (52). In the second equality, we split $\langle\mathbf{1}_{i}+\boldsymbol{\mu}_{\nu}^{*},{\mathbf{G}}({\mathbf{f}}_{t})+\nu\mathbf{1}\rangle$ into $\langle\boldsymbol{\mu}_{\nu}^{*},{\mathbf{G}}({\mathbf{f}}_{t})+\nu\mathbf{1}\rangle$ and $\langle\mathbf{1}_{i},{\mathbf{G}}({\mathbf{f}}_{t})+\nu\mathbf{1}\rangle$ using additivity of the inner product. Via the definition of ${\mathcal{L}}^{s}$, we rewrite ${S}({\mathbf{f}}_{t})+\langle\boldsymbol{\mu}_{\nu}^{*},{\mathbf{G}}({\mathbf{f}}_{t})+\nu\mathbf{1}\rangle$ in the second equality as ${\mathcal{L}}^{s}({\mathbf{f}}_{t},\boldsymbol{\mu}_{\nu}^{*})$ in the third equality. The last term in the third equality is written using the definition of $\mathbf{1}_{i}$ and $G_{i}$ denotes the $i$-th constraint of the stacked constraint vector ${\mathbf{G}}$. Since $({\mathbf{f}}_{\nu}^{*},\boldsymbol{\mu}_{\nu}^{*})$ is a saddle point of ${\mathcal{L}}^{s}$, it holds that

[TABLE]

From the relation of (65), we know the first term in the third equality of (B) is positive, and thus drop it to write:

[TABLE]

where in the second equality we have added and subtracted {\mathbb{E}}\big{[}{\mathcal{L}}^{s}({\mathbf{f}}_{\nu}^{*},\boldsymbol{\mu}_{t})\big{]} and the last inequality comes from the fact that {\mathbb{E}}\big{[}{\mathcal{L}}^{s}({\mathbf{f}}_{\nu}^{*},\boldsymbol{\mu}_{t})\big{]}-{\mathbb{E}}\big{[}{\mathcal{L}}^{s}({\mathbf{f}}_{\nu}^{*},\boldsymbol{\mu}_{\nu}^{*})\big{]}\leq 0 from the relation (65). Now, summing (B) over $t=1,\dots,T$ and using (B) we get:

[TABLE]

where $\alpha=\epsilon/\eta$ defined in Thm. 1. Next, using Assumption 5, the first term in (B) by $R_{{\mathcal{B}}}^{2}$, and the second term on the right hand side of (B) is bounded as

[TABLE]

In the second inequality of (B), we have used the fact that each $\mu_{i,\nu}^{*}$ is positive thus allowing us to use the inequality $(a^{2}+b^{2}+c^{2})\leq(a+b+c)^{2}$ where $a,b$ and $c$ are positive. In the last inequality of (B) we have used the upper bound of (85) (in the supplementary) to bound the \big{(}\sum_{i=1}^{E}\mu_{i,\nu}^{*}\big{)}^{2} term. Now setting step size $\eta=1/\sqrt{T}$ and using Assumption 5 and the upper bound of (B), we get

[TABLE]

where \Gamma\coloneqq\frac{1}{2}\bigg{[}R_{{\mathcal{B}}}^{2}+(1+\delta)\bigg{(}2+2\Big{(}\frac{{4}VR_{{\mathcal{B}}}(CX+\lambda R_{{\mathcal{B}}})}{\xi}\Big{)}^{2}\bigg{)}+K\bigg{]}. If $\epsilon\leq 2R_{{\mathcal{B}}}$, then it follows from (69) that

[TABLE]

Setting $\nu=\zeta T^{-1/2}+\Lambda\alpha$, where $\zeta\geq\Gamma$ and $\Lambda\geq 4VR_{{\mathcal{B}}}$ ensures the aggregation of constraints gets satisfied on long run, i.e., \sum_{t=1}^{T}{\mathbb{E}}\Big{[}G_{i}({\mathbf{f}}_{t})\Big{]}\leq 0, as stated in (37). Analogous logic applies for all edges $i\in{\mathcal{E}}$. $\blacksquare$

Appendix C Proof of Corollary 1

The proof generalizes that of the classical Representer Theorem. The inner minimization in (10) with respect to ${\mathbf{f}}$ can be written as

[TABLE]

Let the subspace of functions spanned by the kernel function $\kappa({\mathbf{x}}_{i,t},.)$ for ${\mathbf{x}}_{i,k}\in{\mathcal{S}}_{i}$ be denoted as ${\mathcal{F}}_{\kappa,{\mathcal{S}}_{i}}$, i.e.,

[TABLE]

We denote the projection of $f_{i}$ on the subspace ${\mathcal{F}}_{\kappa,{\mathcal{S}}_{i}}$ as $f_{ip}$ and the component perpendicular to the subspace as $f_{i\perp}$, which can be written as $f_{i\perp}=f_{i}-f_{ip}$. Now we can write

[TABLE]

Thus the evaluation of $f_{i}$ at any arbitrary training point ${\mathbf{x}}_{ik}$ is independent of $f_{i\perp}$. Using this fact, we can now write (71) as,

[TABLE]

Thus from (74), we can say that ${\mathcal{E}}({\mathbf{f}};{\mathcal{S}},\boldsymbol{\mu})$ is independent of $f_{i\perp}$. As we are minimizing (11) with respect to $f_{i}$, the evaluation of $f_{j}$ at the training point of node $i$ can be treated as a constant in ${\mathcal{E}}({\mathbf{f}};{\mathcal{S}},\boldsymbol{\mu})$ which is the first part in (11). Additionally, note that $\lambda\cdot\|f_{i}\|_{{\mathcal{H}}}^{2}\cdot 2^{-1}\geq\lambda\cdot\|f_{ip}\|_{{\mathcal{H}}}^{2}\cdot 2^{-1}$. Therefore, given any $\boldsymbol{\mu}$, the quantity ${\mathcal{E}}({\mathbf{f}};{\mathcal{S}},\boldsymbol{\mu})+\sum_{i=1}^{V}\lambda\cdot\|f_{i}\|_{{\mathcal{H}}}^{2}\cdot 2^{-1}$ is minimized at some $f_{i}^{*}(\boldsymbol{\mu}_{i})$ such that $f_{i}^{*}(\boldsymbol{\mu}_{i})$ lies in ${\mathcal{F}}_{k,{\mathcal{S}}_{i}}$. This holds specifically for $\boldsymbol{\mu}_{i}^{*}$ where $f_{i}^{*}=f_{i}^{*}(\boldsymbol{\mu}_{i}^{*})$, there by completing the proof. $\blacksquare$

Appendix D Statement and Proof of Lemma 1

Using Assumption 4, we bound the gap between optimal of problem (3) and (12) and is presented as Lemma 1.

Lemma 1

Under Assumption 2, 4 and 5, for $0\leq\nu\leq\xi/2$, it holds that:

[TABLE]

Proof.

Let $({\mathbf{f}}^{*},\boldsymbol{\mu}^{*})$ be the solution to (3) and $({\mathbf{f}}_{\nu}^{*},\boldsymbol{\mu}_{\nu}^{*})$ be the solution to (12). As $\nu\leq\frac{\xi}{2}\leq\xi$, there exists a strictly feasible primal solution ${\mathbf{f}}^{\dagger}$ such that $G({\mathbf{f}}^{\dagger})+\mathbf{1}\nu\leq G({\mathbf{f}}^{\dagger})+\mathbf{1}\xi$, where $\mathbf{1}$ denotes the vector of all ones and $G$ denotes the stacked vector of constraints as defined in the proof of Theorem 2. Hence strong duality holds for (12). Therefore using the definition of $S({\mathbf{f}})$ from (2), we have

[TABLE]

where the inequality in (76) comes from from the optimality of ${\mathbf{f}}_{\nu}^{*}$ and (77) comes from the fact that $G({\mathbf{f}}^{*})\leq 0$. Next using Assumption 4, we have strict feasibility of ${\mathbf{f}}^{\dagger}$, so using (76) we can write:

[TABLE]

Thus from (78), we can equivalently write,

[TABLE]

Now we upper bound the difference of $S({\mathbf{f}}^{\dagger})-S({\mathbf{f}}_{\nu}^{*})$. Using the definition of $S({\mathbf{f}})$, we write the difference of $S({\mathbf{f}}^{\dagger})-S({\mathbf{f}}_{\nu}^{*})$ as

[TABLE]

Next, we bound the sequence in (80) as

[TABLE]

where using triangle inequality we write the first inequality and then using Assumption (2) of Lipschitz-continuity condition we write the second inequality. Further, using reproducing property of $\kappa$ and Cauchy-Schwartz inequality, we simplify |f^{\dagger}_{i}\big{(}{\mathbf{x}}_{i,t})-f_{i,\nu}^{\star}\big{(}{\mathbf{x}}_{i,t})| in (D) as

[TABLE]

where the last inequality comes from Assumption 1 and 5. Now, we consider the $|\|f^{\dagger}_{i}\|^{2}_{{\mathcal{H}}}-\|f_{i,\nu}^{\star}\|^{2}_{{\mathcal{H}}}|$ present in the right-hand side of (D),

[TABLE]

Substituting (82) and (83) in (D), we obtain

[TABLE]

Now using (84), we rewrite (79) as

[TABLE]

Finally, we use (85) in (77) and get the required result:

[TABLE]

∎

The importance of Lemma 1 is that it establishes the fact that the gap between the solutions of the problem (3) and (12) is ${\mathcal{O}}(\nu)$.

Appendix E Statement and Proof of Lemma 2

We bound the primal and dual stochastic gradients used for (26) and (23), respectively in the following lemma.

Lemma 2

Using Assumptions 1-5, the mean-square-magnitude of the primal and dual gradients of the stochastic augmented Lagrangian $\hat{{\mathcal{L}}}_{t}({\mathbf{f}},\boldsymbol{\mu})$ defined in (14) are upper-bounded as

[TABLE]

for some $0<K_{1}<\infty$.

Proof.

In this proof for any $({\mathbf{f}},\boldsymbol{\mu})\in{\mathcal{H}}^{V}\times{\mathbb{R}}^{E}_{+}$ we upper bound the mean-square-magnitude of primal gradient as

[TABLE]

where for the first equality we have used the fact that the functional gradient is a concatenation of functional gradients associated with each agent. The second inequality is obtained by considering the worst case estimate across the network. In the right-hand side of (89) we substitute the value of $\nabla_{fi}\hat{{\mathcal{L}}}_{t}({\mathbf{f}},\boldsymbol{\mu})$ from (17) to obtain,

[TABLE]

In (90), we have used the fact that $\|a+b\|_{{\mathcal{H}}}^{2}\leq 2\cdot(\|a\|_{{\mathcal{H}}}^{2}+\|b\|_{{\mathcal{H}}}^{2})$ for any $a,b\in{\mathcal{H}}$, i.e., the sum of squares inequality. Next we again use the sum of squares inequality for the first bracketed term in the right hand side of (90) and also used Assumption 5 to upper bound $\|f_{i}\|_{{\mathcal{H}}}^{2}$ by $R_{\mathcal{B}}^{2}$ and get,

[TABLE]

where $c(\lambda):=\!\!2V\!\lambda^{2}\!\cdot\!R_{{\mathcal{B}}}^{2}$. Using Cauchy-Schwartz inequality, the first term on the right-hand side of (91) can be written as

[TABLE]

Then using Assumptions 1 and 2, we bound $\|\ell_{i}^{\prime}(f_{i}({\mathbf{x}}_{i,t}),y_{i,t})\|^{2}$ by $C^{2}$ and $\|\kappa({\mathbf{x}}_{i,t},\cdot)\|^{2}$ by $X^{2}$. Similarly we use Cauchy-Schwartz inequality for the second term in (91) and bound $\|\kappa({\mathbf{x}}_{i,t},\cdot)\|^{2}$ by $X^{2}$. Now using these, (91) can be written as,

[TABLE]

where $V^{\prime}:=VX^{2}$. Using Assumption 3, we bound $h^{\prime}_{ij}(f_{i}({\mathbf{x}}_{i,t}),f_{j}({\mathbf{x}}_{i,t}))$ present in the second term on the right-hand side of (92) by $L_{h}$ and then taking the constant $L_{h}$ out of the summation, we get

[TABLE]

Here, $|n_{i}|$ denotes the number of neighborhood nodes of agent $i$. Then we have used the fact $\|\sum_{j=1}^{|n_{i}|}\mu_{ij}\|^{2}\leq{|n_{i}|}\sum_{j=1}^{|n_{i}|}|\mu_{ij}|^{2}$ and got

[TABLE]

Next we upper bound ${|n_{i}|}$ and $\sum_{j=1}^{|n_{i}|}|\mu_{ij}|^{2}$ by ${E}$ and $\|\boldsymbol{\mu}\|^{2}$ and write (94) as

[TABLE]

Thus (95) which establishes an the upper bound on $\mathbb{E}[\|\nabla_{\mathbf{f}}\hat{L}_{t}({\mathbf{f}},\boldsymbol{\mu})\|^{2}_{{\mathcal{H}}}]$ is valid.

With this in hand, we now shift focus to deriving a similar upper-bound on the magnitude of the dual stochastic gradient of the Lagrangian \mathbb{E}\Big{[}\|\nabla_{\boldsymbol{\mu}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}},\boldsymbol{\mu})\|^{2}_{{\mathcal{H}}}\Big{]} as

[TABLE]

In the first equality we write the concatenated version of the dual stochastic gradient associated with each agent, whereas the second inequality is obtained by considering the worst case bound. In the third inequality, we use the fact $|a-b-c|^{2}\leq|a-c|^{2}$ owing to the fact that the right hand side of the inequality is a scalar. Next, applying $\|a+b\|_{{\mathcal{H}}}^{2}\leq 2\cdot(\|a\|_{{\mathcal{H}}}^{2}+\|b\|_{{\mathcal{H}}}^{2})$ for any $a,b\in{\mathcal{H}}$, we get

[TABLE]

Here we have ignored the $\nu^{2}$ term as $\nu<1$ and can be subsumed within the first term. Then we bound the first term in (97) using Assumption 3 and the second term is upper bounded by $\|\boldsymbol{\mu}\|^{2}$

[TABLE]

Next, we use $|f_{i}({\mathbf{x}}_{i,t})|^{2}=|\langle f_{i},\kappa({\mathbf{x}}_{i,t},\cdot)\rangle_{\mathcal{H}}|^{2}\leq\|f_{i}\|_{\mathcal{H}}^{2}\cdot\|\kappa({\mathbf{x}}_{i,t},\cdot)\|_{\mathcal{H}}^{2}$ and then we have upper bounded $\|f_{i}\|_{\mathcal{H}}^{2}$ and $\|\kappa({\mathbf{x}}_{i,t},\cdot)\|_{\mathcal{H}}^{2}$ by $R_{{\mathcal{B}}}^{2}$ and $X^{2}$, and we obtain

[TABLE]

∎

Appendix F Statement and Proof of Lemma 3

The following lemma bounds the difference of projected stochastic functional gradient and un-projected stochastic functional gradient.

Lemma 3

The difference between the stochastic functional gradient defined by ${\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ and projected stochastic functional gradient $\tilde{\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$, is bounded as

[TABLE]

for all $t>0$. Here, $\eta>0$ is the algorithm step-size and $\epsilon>0$ is the error tolerance parameter of the KOMP.

Proof.

Considering the squared-Hilbert- norm difference of the left hand side of (100)

[TABLE]

In (101c) we used the stacked version of $\tilde{f}_{i,t+1}$ to substitute $\tilde{{\mathbf{f}}}_{t+1}$ in place of ${\mathbf{f}}_{t}-\eta\tilde{\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$. In (101d) we used the error tolerance parameter of the KOMP update. Then taking the square root of (101d) gives the inequality stated in (100) and concludes the proof . ∎

Appendix G Definition and Proof of Lemma 4

Next, Lemma 4 characterizes the instantaneous Lagrangian difference $\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu})-\hat{{\mathcal{L}}}_{t}({\mathbf{f}},\boldsymbol{\mu}_{t})$.

Lemma 4

Under Assumptions 1-5 and the primal and dual updates generated from Algorithm 1, the instantaneous Lagrangian difference satisfies the following decrement property

[TABLE]

Proof.

Considering the squared hilbert norm of the difference between the iterate ${\mathbf{f}}_{t+1}$ and any feasible point ${\mathbf{f}}$ with each individual $f_{i}$ in the ball ${\mathcal{B}}$ and exoanding it using the (31), we get

[TABLE]

where we have added and subtracted $2\eta\langle{\mathbf{f}}_{t}-{\mathbf{f}},{\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})\rangle$ and gathered like terms on the right-hand side. Now to handle the second term on the right hand side of (G), we use Cauchy Schwartz inequality along with the Lemma 3 to replace the directional error associated with sparse projections with the functional difference defined by the KOMP stopping criterion:

[TABLE]

Now to bound the norm of $\tilde{\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$, the last term in the right hand side of (G), we add and subtract ${\nabla}_{{\mathbf{f}}}\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ and then use the identity $\|a+b\|_{{\mathcal{H}}}^{2}\leq 2.(\|a\|_{{\mathcal{H}}}^{2}+\|b\|_{{\mathcal{H}}}^{2})$ and further use Lemma 3 and finally get,

[TABLE]

Now we substitute the expressions in (104) and (105) in for the second and fourth terms in (G) which allows us to write

[TABLE]

By re-ordering the terms of the above equation, we get

[TABLE]

Using the first order convexity condition for instantaneous Lagrangian $\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$, since it is convex with respect to ${\mathbf{f}}_{t}$ and write

[TABLE]

Next we use (108) in (107) and get

[TABLE]

Similarly, we consider the squared difference of dual variable update $\boldsymbol{\mu}_{t+1}$ in (20) and an arbitrary dual variable $\boldsymbol{\mu}$,

[TABLE]

The above inequality in (110) comes from the non-expansiveness of the projection operator $[.]_{+}$. Next we expand the square of the right-hand side of (110) and get,

[TABLE]

We re-arrange the terms in the above expression and get,

[TABLE]

Since the instantaneous Lagrangian $\hat{{\mathcal{L}}}_{t}({\mathbf{f}}_{t},\boldsymbol{\mu}_{t})$ is concave with respect to the dual variable $\boldsymbol{\mu}_{t}$, i.e.,

[TABLE]

Next we use the left-hand side of the inequality (113) in (112) and get the expression,

[TABLE]

We subtract (114) from (109) to obtain the final expression,

[TABLE]

∎

Appendix H Definition and proof of Lemma 5

In this section, we present the proof of Theorem 1, where we upper bound the growth of the dictionary. But before going into the proof of Theorem 1, we present Lemma 5 which defines the notion of measuring distance of a point from subspace which will be subsequently used in the proof of Theorem 1. Using Lemma 5, we establish the relation between the stopping criteria of the compression procedure to a Hilbert subspace distance.

Lemma 5

Define the distance of an arbitrary feature vector ${\mathbf{x}}$ obtained by the feature transformation $\phi({\mathbf{x}})=\kappa({\mathbf{x}},\cdot)$ to the subspace of the Hilbert space spanned by a dictionary ${\mathbf{D}}$ of size $M$, i.e., ${\mathcal{H}}_{{\mathbf{D}}}$ as

[TABLE]

This set distance gets simplified to the following least-squares projection when dictionary, ${\mathbf{D}}\in{\mathbb{R}}^{p\times M}$ is fixed

[TABLE]

Proof.

The distance to the subspace ${\mathcal{H}}_{{\mathbf{D}}}$ is defined as

[TABLE]

where the second equality comes from the fact that as ${\mathbf{D}}$ is fixed so minimizing over ${\mathbf{f}}$ translates down to minimizing over ${\mathbf{v}}$ since it is the only free parameter now. Now we solve (118) and obtain ${\mathbf{v}}^{*}={\mathbf{K}}_{{\mathbf{D}},{\mathbf{D}}}^{-1}\boldsymbol{\kappa}_{{\mathbf{D}}}({\mathbf{x}})$ minimizing (118) in a manner similar to logic which yields (29). Now using ${\mathbf{v}}^{*}$ we obtain the required result given in (117), thereby concluding the proof. ∎

Appendix I Proof of Corollary 2

Considering the optimality gap bound in (B) and using the fact that $\epsilon\leq{2}R_{{\mathcal{B}}}$, we can write (B) with step size $\eta=1/\sqrt{T}$ as

[TABLE]

Now denoting $Q\coloneqq 4VR_{{\mathcal{B}}}(CX+\lambda R_{{\mathcal{B}}})/\xi$ and using $\nu=\zeta T^{-1/2}+\Lambda\alpha$, we write (119) as

[TABLE]

Now for the optimality gap in (I) to be less than $\varepsilon$, it is sufficient if we consider that both the terms in the right hand side of the inequality (I) are bounded by $\varepsilon/2$, i.e.,

[TABLE]

Thus, from (121) and (122) we can deduce that for an optimality gap to be less than $\varepsilon$ we require ${\mathcal{O}}(1/\varepsilon^{2})$ iterations with $\alpha$ satisfying (121). Now using the bound of $\alpha$ from (121) in (35) of Theorem 1, we get model order complexity bound of ${\mathcal{O}}(1/\varepsilon^{2p})$ required to achieve an optimality gap of $\varepsilon$.