Nonparametric Compositional Stochastic Optimization for Risk-Sensitive Kernel Learning

TL;DR

This paper introduces COLK, a memory-efficient stochastic algorithm for risk-sensitive kernel learning in compositional optimization problems within RKHS, demonstrating convergence and robustness in supervised learning tasks.

Contribution

It presents the first memory-efficient stochastic method for compositional problems in RKHS, with proven convergence and practical robustness.

Findings

COLK reliably converges in experiments.

It attains consistent performance across training runs.

It overcomes overfitting in supervised learning.

Abstract

In this work, we address optimization problems where the objective function is a nonlinear function of an expected value, i.e., compositional stochastic {strongly convex programs}. We consider the case where the decision variable is not vector-valued but instead belongs to a reproducing Kernel Hilbert Space (RKHS), motivated by risk-aware formulations of supervised learning and Markov Decision Processes defined over continuous spaces. We develop the first memory-efficient stochastic algorithm for this setting, which we call Compositional Online Learning with Kernels (COLK). COLK, at its core a two-time-scale stochastic approximation method, addresses the fact that (i) compositions of expected value problems cannot be addressed by classical stochastic gradient due to the presence of the inner expectation; and (ii) the RKHS-induced parameterization has complexity which is proportional to the iteration index which is mitigated through greedily constructed subspace projections. We establish almost sure convergence of COLK with attenuating step-sizes, and linear convergence in mean to a neighborhood with constant step-sizes, as well as the fact that its complexity is at-worst finite. The experiments with robust formulations of supervised learning demonstrate that COLK reliably converges, attains consistent performance across training runs, and thus overcomes overfitting.

Figures

Tables

Table 1. TABLE I : Summary of recent related results of model complexity lower bounds and convergence rate upper bounds where we differentiate as follows: objective structure is compositional (comp.) or non-compositional., objective function is strongly convex (SC) or non-convex (NC). In [ 35 ] , E 𝐸 E denotes the average error in the gradient, and in [ 51 ] , d λ subscript 𝑑 𝜆 d_{\lambda} denotes the degree of freedom for a given regularization parameter λ 𝜆 \lambda . In the table, δ 𝛿 \delta denotes the optimal accuracy parameter.

Refs.	Structure	Obj.	Model compl.	Convergence rate
[35]	Non-comp.	SC	$\Omega \left(1\right)$	$𝒪\left(1/{\left(\delta -E\right)}^2\right)$
[51]	Non-comp.	SC	$\Omega \left(T\right)$	$𝒪\left(1/\delta \right)$
[17]	Non-comp.	SC	$\Omega \left(T^p\right)$	$𝒪\left(\left(1/{\delta }^2\right)\right)$
[52]	Non-comp.	SC	$\Omega \left(\sqrt{T}\log \left(d_{\lambda }\right)\right)$	$𝒪\left(1/{\delta }^2\right)$
Our	Comp.	SC	$\Omega \left({\left(1/\delta \right)}^p\right)$	$𝒪\left(\left(1/\delta \right)\log \left(1/\delta \right)\right)$
Our	Comp.	NC	$\Omega \left({\left(1/\delta \right)}^p\right)$	$𝒪\left(1/{\delta }^2\right)$

Equations

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

$$f^{\star }=\underset{f\in \mathcal{H}}{\mathrm{argmin}}R(f):=J(f)+\frac{\lambda }{2}\Vert f{\Vert }_{\mathcal{H}}^2,$$

$$(i)⟨f,\kappa (\mathbf{u},\cdot ){⟩}_{\mathcal{H}}=f(\mathbf{u})\text{for all }\mathbf{u}\in \mathcal{U},$$

$$(ii)\mathcal{H}=\overline{\text{span}\{\kappa (\mathbf{u},\cdot )\}}\text{for all }\mathbf{u}\in \mathcal{U}.$$

$$f(\mathbf{u})=\sum _{n=1}^N{w}_n\kappa ({\mathbf{\xi }}_n,\mathbf{u}).$$

$$f^{\star }=\underset{f\in \mathcal{H}}{\mathrm{argmin}}\mathbb{E}[l(f({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})],$$

$${\hat{\mathbf{y}}}^{\star }=\underset{\hat{\mathbf{y}}\in {\mathcal{Y}}^{\mathcal{X}}}{\mathrm{argmin}}\mathbb{E}[l(\hat{\mathbf{y}}({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})],$$

$$\underset{\text{Estimation Error}}{\underbrace{\mathbb{E}[l(\hat{f}({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})]-\underset{f\in \mathcal{H}}{\min }\mathbb{E}[l(f({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})]}}$$

$$+\underset{\text{Approximation Error}}{\underbrace{\underset{f\in \mathcal{H}}{\min }\mathbb{E}[l(f({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})]-\underset{\hat{\mathbf{y}}\in {\mathcal{Y}}^{\mathcal{X}}}{\min }\mathbb{E}[l(\hat{\mathbf{y}}({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})]}},$$

$$\!\!\widetilde{\text{Var}}[l(f({\mathbf{x}}_{\boldsymbol{\theta}}),{\mathbf{y}}_{\boldsymbol{\theta}})]=\mathbb{E}\Big{\{}\!\big{(}l(f({\mathbf{x}}_{\boldsymbol{\theta}}),{\mathbf{y}}_{\boldsymbol{\theta}})-{\mathbb{E}[l(f({\mathbf{x}}_{\boldsymbol{\xi}}),{\mathbf{y}}_{{\boldsymbol{\xi}}})]}\big{)}_{+}^{2}\Big{\}}\;,$$

$${\ell }_{\mathbf{\theta }}(\mathbf{u})$$

$$[h_{\mathbf{\xi }}(f(\mathbf{\theta })){]}_1$$

$$[h_{\mathbf{\xi }}(f(\mathbf{\xi })){]}_2$$

$$f^{\star }=\underset{f\in \mathcal{H}}{\mathrm{argmin}}\mathbb{E}[l(f({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})]+\eta \widetilde{\text{Var}}[l(f({\mathbf{x}}_{\mathbf{\theta }}),{\mathbf{y}}_{\mathbf{\theta }})]+\frac{\lambda }{2}\Vert f{\Vert }_{\mathcal{H}}^2.$$

$$\!\!\!\!V^{\pi}({\mathbf{x}})=\mathbb{E}_{{\mathbf{y}}}\Big{[}\sum_{t=0}^{\infty}\gamma^{t}r({\mathbf{x}}_{t},{\mathbf{a}}_{t},{\mathbf{y}}_{t}){\,\big{|}\,}{\mathbf{x}}_{0}={\mathbf{x}},\{{\mathbf{a}}_{t}=\pi({\mathbf{x}}_{t})\}_{t=0}^{\infty}\Big{]}\;.$$

$$V^{\pi}({\mathbf{x}})=\int_{{\mathcal{X}}}[r({\mathbf{x}},\pi({\mathbf{x}}),{\mathbf{y}})+\gamma V({\mathbf{y}})]\mathbb{P}(d{\mathbf{y}}{\,\big{|}\,}{\mathbf{x}},\pi({\mathbf{x}})),\;$$

$$\displaystyle J\!(V)\!\!:=\!\mathbb{E}_{{\mathbf{x}},\pi({\mathbf{x}})}\!\Big{\{}\!\frac{1}{2}\big{(}\mathbb{E}_{{\mathbf{y}}}\big{[}r({\mathbf{x}},\pi({\mathbf{x}}),{\mathbf{y}})\!+\!\gamma V({\mathbf{y}})\!-\!V({\mathbf{x}}){\,\big{|}\,}{\mathbf{x}},\pi({\mathbf{x}})\big{]}\big{)}^{2}\!\Big{\}}.$$

$$⟨\mathbb{E}\left[{\nabla }_f{h}_{\mathbf{\xi }}(f(\mathbf{\xi }))\right],\mathbb{E}\left[\nabla {\ell }_{\mathbf{\theta }}(\mathbb{E}[h_{\mathbf{\xi }}(f(\mathbf{\xi }))])\right]⟩$$

$$⟨{\nabla }_f{h}_{{\mathbf{\xi }}_t}(f({\mathbf{\xi }}_t)),\nabla {\ell }_{{\mathbf{\theta }}_t}(\mathbb{E}[h_{\mathbf{\xi }}(f(\mathbf{\xi }))])⟩.$$

$$f_{t+1}=(1-\lambda \alpha )f_t-\alpha ⟨{\nabla }_f{h}_{{\mathbf{\xi }}_t}(f_t({\mathbf{\xi }}_t)),\nabla {\ell }_{{\mathbf{\theta }}_t}(\mathbb{E}[h_{\mathbf{\xi }}(f_t({\mathbf{\xi }}_t))])⟩.$$

$${\mathbf{g}}_{t+1}=(1-\beta ){\mathbf{g}}_t+\beta h_{{\mathbf{\xi }}_t}(f_t({\mathbf{\xi }}_t))$$

$${\mathbf{g}}_{t+1}$$

$$f_{t+1}=(1-\lambda \alpha )f_t-\alpha ⟨{\nabla }_f{h}_{{\mathbf{\xi }}_t}(f_t({\mathbf{\xi }}_t)),\nabla {\ell }_{{\mathbf{\theta }}_t}({\mathbf{g}}_{t+1})⟩,$$

$$⟨{\nabla }_f$$

$$=\sum _{i=1}^m\frac{\partial h_{{\mathbf{\xi }}_t}^i(\omega )}{\partial \omega }{|}_{\omega =f({\mathbf{\xi }}_t)}\times \frac{\partial {\ell }_{{\mathbf{\theta }}_t}(\mathbf{u})}{\partial u_i}{|}_{\mathbf{u}={\mathbf{g}}_{t+1}}\kappa ({\mathbf{\xi }}_t,\cdot )$$

$$=⟨h_{{\mathbf{\xi }}_t}^{\prime }(f_t({\mathbf{\xi }}_t)),{\ell }_{{\mathbf{\theta }}_t}^{\prime }({\mathbf{g}}_{t+1})⟩\kappa ({\mathbf{\xi }}_t,\cdot ).$$

$$f_{t+1}=(1-\lambda \alpha )f_t-\alpha ⟨h_{{\mathbf{\xi }}_t}^{\prime }(f_t({\mathbf{\xi }}_t)),{\ell }_{{\mathbf{\theta }}_t}^{\prime }({\mathbf{g}}_{t+1})⟩\kappa ({\mathbf{\xi }}_t,\cdot ).$$

$$f_t(\mathbf{u})=\sum _{n=1}^{t-1}{w}_n\kappa ({\mathbf{\xi }}_n,\mathbf{u})={\mathbf{w}}_t^T{\kappa }_{{\mathbf{U}}_t}(\mathbf{u}),$$

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

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

$$f_{t+1}$$

$$:={\mathcal{P}}_{{\mathcal{H}}_{{\mathbf{D}}_{t+1}}}\left[(1-\lambda \alpha )f_t-\alpha ⟨h_{{\mathbf{\xi }}_t}^{\prime }(f_t({\mathbf{\xi }}_t)),{\ell }_{{\mathbf{\theta }}_t}^{\prime }({\mathbf{g}}_{t+1})⟩\kappa ({\mathbf{\xi }}_t,\cdot )\right].$$

$${\mathbf{g}}_{t+1}$$

$${\tilde{f}}_{t+1}=(1-\lambda \alpha )f_t-\alpha ⟨h_{{\mathbf{\xi }}_t}^{\prime }(f_t({\mathbf{\xi }}_t)),{\ell }_{{\mathbf{\theta }}_t}^{\prime }({\mathbf{g}}_{t+1})⟩\kappa ({\mathbf{\xi }}_t,\cdot )$$

Full text

I Introduction

In this work, we focus on compositional stochastic programming, a setting where the objective function is an expectation over a set of random convex functions, each of which depends on the expected value of a different random convex function. This problem setting has received recent attention in operations research [3, 4] and machine learning [5] due to its ability to gracefully address technicalities that arise in the theory of Markov Decision Problems (MDPs) [6] and bias-variance issues in supervised learning [7]. Our goal is to solve this class of problems when the decision variable is not vector-valued, as in [3], but is instead itself a function. This setting arises intrinsically in MDPs defined over the continuous state and action spaces [8] or when accounting for risk [9] in supervised learning with nonlinear interpolators [10].

The theory of optimization in function space began with variational calculus [11] and Hamilton’s Principle [12]. However, in modern applications, we require solutions to such problems in situations where classical methods no longer apply. Two different issues arise: (1) how to evaluate the expectation (integral) and (2) how to parameterize the function so that tractable updates may be obtained. Setting aside (1) for now, to address (2), i.e., to handle the intractability of general functional optimization, one must restrict the function we seek to not only yield a computationally tractable formulation, but also one be rich enough to address common experimental settings. In learning theory, for instance, we typically restrict the function to be a neural network [13] or a nonparametric basis expansion in terms of data [14], whereas in control systems, polynomial interpolation [15] and kriging [16] are popular. In this work, we address the case where the function class is a Reproducing Kernel Hilbert Space (RKHS), motivated by a recently developed memory-efficient parameterization of a function that is infinite dimensional [17]. This approach subsumes polynomial interpolation [15], avoids the memory explosion associated with large sample-size kriging [16], and preserves convexity, thus avoiding convergence to poor stationary points rampant in neural network training [18].

With the function class specified, we turn to discuss how to solve the associated functional stochastic program: doing so requires iterative stochastic methods [19, 20], since deterministic approaches [21] require computing gradients that depend on infinitely many realizations of a random variable, thus exhibiting prohibitive complexity. Unfortunately, standard stochastic gradient descent (SGD) is inapplicable to the compositional setting, because, for a single stochastic descent direction, one requires the evaluation of an additional inner expectation, an observation that was popularized in reinforcement learning as “the double sampling problem” [6]. To ameliorate this issue, we develop a functional nonparametric extension of stochastic quasi-gradient (SQG) method [22, 3], which uses two time-scale stochastic approximation: one uses a quasi-stationary estimate of the inner expectation, whereas the other executes stochastic descent [23].

Different from the estimator for the inner expectation in [3], in this work we consider a momentum (gradient tracking) scheme that uses a difference of instantaneous costs, which permits us to tighten existing rate analyses for compositional problems with minimal additional computational overhead. Gradient tracking has been investigated extensively in the distributed optimization literature to have agents’ track the gradient of the global objective – see for instance, [24, 25, 26] and more recently [27, 28, 29, 30]. Here we adopt the spirit of using the previous iterate, but rather than track the gradient itself, we use history information to track the evaluation of the inner objective that appears in compositional problems, which to the best of our knowledge, is the first time it has been employed in this context, and permits us to refine existing rates for compositional stochastic programming. This contrasts with its use for improving the vanilla gradient estimation error that arises in non-convex optimization as studied by [31].

However, with the choice of ${\mathcal{H}}$ as an RKHS, and sequential application of the Representer Theorem [32], means that the function parameterization grows with the iteration index [10], and thus becomes untenable for expected value problems. While memory-reduction methods for RKHS optimization exist using, for instance, random feature approximations [33], forgetting factors [10], random dropping [34], or projections onto subspaces of fixed size [35], none of the aforementioned approaches can ensure that the stochastic descent properties of the algorithm are nearly preserved with probability $1$, and hence fail to (approximately) preserve almost sure convergence of the optimization method to which they are applied. Thus, it is an open challenge how to generalize the results of [3] to optimization over RKHS. Thus, our main contributions are to:

• •

extend SQG to RKHS, which we compress with matching pursuit [36] (Sec. III). We tailor the compression to ensure valid descent [17, 37]. We call this method Compositional Online Learning with Kernels (COLK).

• •

characterize the number of iterations and the memory size requirement mathematically to achieve $\delta$-suboptimal solutions for strongly convex (Theorem 2) and non-convex (Theorem 3) compositional objectives in terms of mean-square error to the optimizer and the expected gradient norm, respectively.

• •

establish that the worst-case complexity of the function sequence is finite (Theorem 1), with an explicit dependence on the parameter dimension, algorithm step-size, and compression budget. This result establishes a trades-off between convergence accuracy and model complexity.

• •

experimentally (Sec. V) validate this method on a problem instantiation defined by robust supervised learning. Doing so yields nonlinear statistical models whose bias and variance are small, first on a synthetic data regression outliers which has a heavier tailed distribution, i.e., more outliers are present, and then on benchmark data: lidar [38]. We observe that COLK yields consistently accurate performance across training realizations, meaning that it does not overfit, in contrast to other methods that cannot minimize risk functionals [7].

Overall, then, this work differs from [3] in the following ways: (i) [3] considers vector-valued decision variables, whereas in this case, we focus on the RKHS setting. This more general setting implies that the function class to which statistical models belong is universal, but this universality comes at the cost of complexity. In the vector setting, stochastic compositional gradients may be used for updates without any memory bottleneck, whereas for RKHS one must use approximate gradients to ensure tractable implementation, but this approximation incurs errors in the gradients, which we analyze. (ii) We consider a modification of stochastic quasi-gradient methods that employ a momentum scheme in the faster time-scale in order to improve the overall convergence rates. We further note that all analysis in [3] is for attenuating step-sizes, whereas our focus is on constant step-sizes. (iii) We further characterize the dependency of the algorithm’s non-asymptotic parameterization complexity on the underlying dimension of the feature space under constant step-sizes (Theorem 1), an attribute that is unique to nonparametric statistics which does not appear in [3]. Moreover, we illuminate the differences between these works experimentally in Section V.

Notation: All the scalars are denoted by letters in regular font $s$, vectors are bold ${\mathbf{s}}$, and matrices are capitalized ${\mathbf{S}}$. The notation $\left\|\cdot\right\|_{\mathcal{H}}$ represents the norm in the RKHS symbolized by $\mathcal{H}$. Subspaces of the RKHS are denoted by ${\mathcal{H}}_{{\mathbf{S}}}$. The inner product operator between two functions $f\in\mathcal{H}$ and $f^{\prime}\in\mathcal{H}$ is denoted by $\boldsymbol{\langle}f,f^{\prime}\boldsymbol{\rangle}_{\mathcal{H}}$. The expectation operator is symbolized by $\mathbb{E}$. The notation $\circ$ describes the compositional operator, for instance, $f\circ g(x)$ denotes the composition as $f(g(x))$. A kernel function is denote by $\kappa(\cdot,\cdot)$ which takes two vectors as arguments. Stacking of kernel evaluations is denoted as the empirical kernel map $\boldsymbol{\kappa}_{{\mathbf{U}}}(\cdot)$ and the kernel matrix ${\mathbf{K}}_{{\mathbf{U}},{\mathbf{V}}}$ is the matrix whose $(i,j)$ entry is given by $\kappa({\mathbf{u}}_{i},{\mathbf{v}}_{j})$. ${\mathcal{P}}_{\mathcal{C}}$ denotes the orthogonal projection operator onto a convex set ${\mathcal{C}}$ . Gradients are denoted in the usual way, and partial derivative, for instance, of scalar valued function $\ell({\mathbf{u}})$ with respect to $u_{i}$ at ${\mathbf{v}}$ , is denoted by, $\frac{\partial\ell({\mathbf{u}})}{\partial u_{i}}|_{{\mathbf{u}}={\mathbf{v}}}$.

II Compositional Stochastic Programming in RKHS

In this work, we focus on solving functional optimization problems whose objective is a nonlinear function of an expected value. That is, the objective function is a composition of two functions, each of which is an expected value over a set of functions parameterized by a pair of random variables. More specifically, there are two sets of random variables $\{\boldsymbol{\xi}_{t}\}\subset{\mathbb{R}}^{p}$ and $\{\boldsymbol{\theta}_{t}\}\subset\mathbb{R}^{p}$. In general both the random variables are allowed to be dependent, but for the ease of analysis and understanding, we assume that $\boldsymbol{\xi}\in\boldsymbol{\Xi}\subset{\mathbb{R}}^{p}$, $\boldsymbol{\theta}\in\boldsymbol{\Theta}\subset{\mathbb{R}}^{p}$ and $\boldsymbol{\xi}$, $\boldsymbol{\theta}$ are independent of each other. Considering these random pairs, the cost takes the form $J(f):=(L\circ{\mathbf{H}})(f)$, where ${\mathbf{H}}(f)=\mathbb{E}_{\boldsymbol{\xi}}\left[\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}}(f(\boldsymbol{\xi}))\right]$ is a map ${\mathbf{H}}:{\mathcal{H}}\rightarrow{\mathbb{R}}^{m}$ that is an expectation over a set of random functions $\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}}(f(\boldsymbol{\xi}))$. Specific instances of $\mathscr{h}$ will be discussed in Examples 1-2 shortly. Often in applications, $\mathscr{h}$ quantifies statistics of some loss function, possibly $\ell$. Similarly, $L({\mathbf{u}})=\mathbb{E}_{\boldsymbol{\theta}}\left[\ell_{\boldsymbol{\theta}}({\mathbf{u}})\right]$ is a map $L:\mathbb{R}^{m}\rightarrow{\mathbb{R}}$ that is an expected value over a random collection variable. Further, ${\mathcal{H}}$ is a function space to be subsequently specified. In this work, we focus on the functional compositional stochastic program:

[TABLE]

We consider both the cases that $J(f)$ is convex with respect to function $f$ with Tikhonov regularizer $\frac{\lambda}{2}\|f\|_{{\mathcal{H}}}$ to ensure strong convexity, defining the regularized loss [39]

[TABLE]

and the case where $R(f)$ is non-convex111In which case the regularizer is inconsequential, and may be set as $\lambda=0$.. The feasible set ${\mathcal{H}}$ of (1), the domain of $H$, and hence $J$, is not Euclidean space ${\mathbb{R}}^{p}$, as in [3], but instead is a Hilbert space equipped with a unique kernel function, $\kappa:{\mathcal{U}}\times{\mathcal{U}}\rightarrow{\mathbb{R}}$, such that:

[TABLE]

where $\langle\cdot,\cdot\rangle_{{\mathcal{H}}}$ denotes the Hilbert inner product for ${\mathcal{H}}$ and ${\mathcal{U}}:=\boldsymbol{\Xi}\cup\boldsymbol{\Theta}$ denotes the union of data domains $\boldsymbol{\Xi}$ and $\boldsymbol{\Theta}$, whose elements ${\mathbf{u}}$ are random variables $\boldsymbol{\xi}$ or $\boldsymbol{\theta}$. We assume that the kernel is nonnegative, i.e. $\kappa({\mathbf{u}},{\mathbf{u}}^{\prime})\geq 0$ for all ${\mathbf{u}},{\mathbf{u}}^{\prime}\in{\mathcal{U}}$ so that it is a Mercer kernel. Function spaces of this type are called reproducing kernel Hilbert spaces (RKHS) [40].

In (II), property (i) is called the reproducing property of the kernel and comes from the Riesz Representation Theorem [41]. Replacing $f$ by $\kappa({\mathbf{u}}^{\prime},\cdot)$ in (II) (i) yields the expression $\langle\kappa({\mathbf{u}}^{\prime},\cdot),\kappa({\mathbf{u}},\cdot)\rangle_{{\mathcal{H}}}=\kappa({\mathbf{u}},{\mathbf{u}}^{\prime})$, which is why $\kappa$ is called “reproducing.” This property provides a practical means by which to access a nonlinear transformation of the input space ${\mathcal{U}}$. Specifically, denote by $\phi(\cdot)$ a nonlinear map of the feature space that assigns to each ${\mathbf{u}}$ the kernel function $\kappa(\cdot,{\mathbf{u}})$. Then the reproducing property of the kernel allows us to write the inner product of the image of distinct feature vectors ${\mathbf{u}}$ and ${\mathbf{u}}^{\prime}$ under the map $\phi$ in terms of kernel evaluations only: $\langle\phi({\mathbf{u}}),\phi({\mathbf{u}}^{\prime})\rangle_{{\mathcal{H}}}=\kappa({\mathbf{u}},{\mathbf{u}}^{\prime})$. This is the kernel trick, and it provides a principled method for function estimation.

Moreover, property (II) (ii) states that any function $f\in{\mathcal{H}}$ may be written as a linear combination of kernel evaluations. For kernelized and regularized empirical risk minimization (i.e., the sample average approximation of (1) for some fixed $N$ realizations of $\boldsymbol{\xi}$ and $\boldsymbol{\theta}$), the Representer Theorem [40, 32] establishes that the optimal $f$ in function class ${\mathcal{H}}$ may be written as an expansion of kernel evaluations only at elements of the training set as222In the non-convex setting, this holds when $R$ is either quasi-convex [42] or strongly convex in a neighborhood of a stationary point.

[TABLE]

where ${\mathbf{w}}=[w_{1},\cdots,w_{N}]^{T}\in{\mathbb{R}}^{N}$ denotes the weight vector. The upper summand index $N$ in (4) is henceforth referred to as the model order. Common choices $\kappa$ include the polynomial kernel and the radial basis kernel, i.e., $\kappa({\mathbf{u}},{\mathbf{u}}^{\prime})=\left({\mathbf{u}}^{T}{\mathbf{u}}^{\prime}+b\right)^{c}$ and $\kappa({\mathbf{u}},{\mathbf{u}}^{\prime})=\exp\left\{-\frac{\lVert{\mathbf{u}}-{\mathbf{u}}^{\prime}\rVert_{2}^{2}}{2c^{2}}\right\}$, respectively, with ${\mathbf{u}},{\mathbf{u}}^{\prime}\in{\mathcal{U}}$.

Then, one may use the Representer Theorem (4) to transform the sample average approximation of (1) over all of ${\mathcal{H}}$ into the parametric problem of two $N$-dimensional weight vector ${\mathbf{w}}$. However, as $N\rightarrow\infty$, the function representation becomes infinite as well. In this work, we seek to find functions that are close-to-optimal solutions to (1) but also admit a finite-memory representation. Before turning to develop an algorithmic tool which does so, we note that the problem setting (1) arises in diverse applications. Here we mention two, the first of which is the focus of this work. In both instances, we focus on the case where samples from the data distribution are revealed incrementally, the input-ouput relationship between features and targets is not necessarily linearly separable, and hence feature selection must done on the fly during training. This setting makes RKHS parameterizations advantageous in terms of balancing descriptive richness with preserving convexity, and hence defining training schemes that may find the optimal representation within the defined function class.

Example 1** **(Robust Supervised Learning)

Consider a random pair $({\mathbf{x}},{\mathbf{y}})\in{\mathcal{X}}\times{\mathcal{Y}}$, realizations of which are training examples $({\mathbf{x}}_{n},{\mathbf{y}}_{n})$, and ${\mathcal{X}}\subset{\mathbb{R}}^{p}$, the $p$-dimensional Euclidean space. In comparison to formulation in (1), we have $\boldsymbol{\theta}=\boldsymbol{\xi}={\mathbf{x}}$ and ${\mathbf{y}}^{\boldsymbol{\theta}}={\mathbf{y}}^{\boldsymbol{\xi}}={\mathbf{y}}$ for this example which represents the corresponding target values. In the case of classification with $C$ classes, ${\mathcal{Y}}=\{1,\dots,C\}$, whereas in the case of regression ${\mathcal{Y}}\subset{\mathbb{R}}^{q}$. In supervised learning, we learn an estimator $f({\mathbf{x}})$ according to its ability to minimize a loss function $l:{\mathcal{H}}\times{\mathcal{X}}\times{\mathcal{Y}}\rightarrow{\mathbb{R}}$ averaged over data:

[TABLE]

where we define $L(f)=\mathbb{E}[l(f({\mathbf{x}}_{\boldsymbol{\theta}}),{\mathbf{y}}_{\boldsymbol{\theta}})]$, and ignore the regularizer for the moment. The loss $l$ quantifies the merit of the estimator $f({\mathbf{x}})$ with respect to its target ${\mathbf{y}}$. However, as it is well known in statistics [39], solving (5) is really only an approximation of the Bayes optimal estimator

[TABLE]

where ${\mathcal{Y}}^{{\mathcal{X}}}$ denotes the space of all functions $\hat{{\mathbf{y}}}:{\mathcal{X}}\rightarrow{\mathcal{Y}}$ that map data ${\mathbf{x}}$ to target variables ${\mathbf{y}}$. Suppose we obtain some estimate $\hat{f}$ by approximately minimizing (5). Then the performance difference associated with $\hat{f}$ and the Bayes optimal $\hat{{\mathbf{y}}}^{\star}$ (6) is given as

[TABLE]

where we add and subtract the optimal supervised cost $\min_{f\in{\mathcal{H}}}\mathbb{E}[l(f({\mathbf{x}}_{\boldsymbol{\theta}}),{\mathbf{y}}_{\boldsymbol{\theta}})]$ to obtain that this discrepancy decomposes into two terms: the estimation error, and approximation error. In most cases, the estimation error may be identified with bias, and the approximation error may be identified with variance plus noise – we subsequently identify the decomposition (1) with bias and variance, as bias and estimation error are both defined as subsampling errors, whereas variance and approximation error are both defined by quality of model fit across training runs. See [43][Sections 3.1 - 3.2].

In supervised learning, typically we try to make the model bias small as the number of data points goes to infinity, resigning ourselves to the fact that the variance is an intrinsic penalty we suffer for selecting a particular modeling hypothesis: the function class to which $f$ belongs. However, authors in operations research [9] and applied probability [44] have proposed to optimize both the expected loss over all data plus a measure of the dispersion of the estimate with respect to its target variable as a way of accounting for the unknown approximation error of the modeling hypothesis in (1). Many measures of dispersion are possible, but one which yields a convex formulation is the semivariance.

[TABLE]

where $a_{+}=\max(a,0)$ denotes the positive projection. Note that when we omit the positive projection, (8) reduces to the variance of the instantaneous loss $l(f({\mathbf{x}}_{\boldsymbol{\theta}}),{\mathbf{y}}_{\boldsymbol{\theta}})$. To see that variance is a composite function observe that when $\boldsymbol{\theta}$ and $\boldsymbol{\xi}$ are independent, we have that

[TABLE]

where note that $u_{1}$ and $u_{2}$ are also random quantities so the first expectation cannot be dropped.

However, the subtraction of the second moment of $l(f({\mathbf{x}}_{\boldsymbol{\theta}}),{\mathbf{y}}_{\boldsymbol{\theta}})$ without positive projection makes the problem non-convex. Thus, using positive projections yields semivariance rather than true variance [9]. With the measure of dispersion in (8), one robust formulation of supervised learning over an RKHS ${\mathcal{H}}$ which accounts for approximation error is

[TABLE]

where $\eta$ is a scaling parameter that tunes the importance of estimation or approximation error. Note that we have added the regularizer back into (12). Solutions of (12), as compared with (5), are better attuned to data points associated with high variance objective evaluation, which practically may be interpreted as outliers or in the case of classification, situations where training examples possess characteristics corresponding to multiple classes. Due to the fact that our analysis requires Lipschitz gradients (Section IV), numerically we approximate the positive projection in (8) with the softmax in Section V – see [21]. An alternative risk measure popular in finance is the conditional value-at-risk (CVaR) [45], which quantifies the loss function at different quantiles of its distribution.

Example 2** **(Policy Evaluation in MDPs)

Another instantiation of (1) is the task of policy evaluation in a continuous Markov Decision Problem (MDP) [46], a quintuple $({\mathcal{X}},{\mathcal{A}},\mathbb{P},r,\gamma)$, where $\mathbb{P}$ is the action-dependent transition probability of the process: when the agent starts in state ${\mathbf{x}}_{t}\in{\mathcal{X}}\subset{\mathbb{R}}^{p}$ at time $t$ and takes an action ${\mathbf{a}}_{t}\in{\mathcal{A}}$, a transition to next state ${\mathbf{y}}_{t}\in{\mathcal{X}}$ is distributed according to {\mathbf{y}}_{t}\sim\mathbb{P}(\cdot{\,\big{|}\,}{\mathbf{x}}_{t},{\mathbf{a}}_{t}). After the agent transitions to a particular ${\mathbf{y}}_{t}$, the MDP provides to it an instantaneous reward $r({\mathbf{x}}_{t},{\mathbf{a}}_{t},{\mathbf{y}}_{t})$, where the reward function is a map $r:{\mathcal{X}}\times{\mathcal{A}}\times{\mathcal{X}}\rightarrow{\mathbb{R}}$.

In policy evaluation, control decisions ${\mathbf{a}}_{t}$ are chosen according to a stochastic stationary policy $\pi:{\mathcal{X}}\rightarrow\rho({\mathcal{A}})$, where $\rho({\mathcal{A}})$ denotes the set of probability distributions over ${\mathcal{A}}$, and one seeks to compute the value of a policy when starting in state ${\mathbf{x}}$, quantified by the discounted expected sum of rewards, or value function $V^{\pi}({\mathbf{x}})$:

[TABLE]

For a single trajectory through the state space ${\mathcal{X}}$, ${\mathbf{y}}_{t}={\mathbf{x}}_{t+1}$. The discount factor $\gamma\in(0,1)$ determines the agent’s farsightedness. From the definition of the value function in (13), one may derive the Bellman evaluation equation [46]:

[TABLE]

for all ${\mathbf{x}}\in{\mathcal{X}}$. The functional fixed point problem (14) defined by Bellman’s equation may be reformulated as a nested stochastic program. To do so, rewrite the integral as an expectation, subtract the value function $V^{\pi}({\mathbf{x}})$ that satisfies the fixed point relation from both sides, and then pull it inside the expectation. Then, to solve (14) in an initialization-independent manner, integrate out ${\mathbf{x}}$, the starting point of the trajectory defining the value function (13), as well as policy $\pi({\mathbf{x}})$. Doing so, i.e., considering the Galerkin relaxation of (14) followed by squaring the resulting expression, yields the compositional objective

[TABLE]

However, since it is intractable to optimize over all bounded functions ${\mathcal{B}}({\mathcal{X}})$, one may restrict the minimization of $J(V)$ to an RKHS ${\mathcal{H}}$. This hypothesis, however, requires the introduction of regularization. Assuming that the Bellman fixed point $V^{\pi}$ is a continuous function, the RKHS approximation may be made close to the true $V^{\pi}$ when used with a universal kernel [47] – see [37]. We note that the problem formulation in (15) is slight generalization of the problem formulation in (12). In particular, observe that the inner random variable ${\mathbf{y}}$ in (15) depends on the outer random variables ${\mathbf{x}}$ and $\pi({\mathbf{x}})$, necessitating the use of the conditional expectation operator. However, this generalization can be readily handled as detailed in [3] and the subsequent analysis considers (12) for brevity.

With the problem setting clarified, we next shift focus to developing an iterative numerical method to solve (1).

III Algorithm Development

Now we turn to solving the stochastic compositional optimization (1) over the RKHS ${\mathcal{H}}$. We focus on the development of stochastic approximation methods such that we may minimize $R$ over ${\mathcal{H}}$ with only sequentially revealed independent and identically distributed realizations of $\boldsymbol{\theta}$ and $\boldsymbol{\xi}$ without knowledge of their underlying distributions. Related ideas are developed for the specialized objective of Example 2 in [37] and in the vector-valued case in [3]. A key distinguishing feature of our method is the use of a gradient tracking scheme to reduce the estimation error of the inner expectation – this is the first time momentum-style updates have been employed for compositional problems, and their use permits us to improve upon existing rate analyses. The fundamental building blocks of our proposed algorithm are a functional generalization of the stochastic quasi-gradient method (Section III-A) operating in tandem with low-dimensional subspace projections that are greedily constructed using matching pursuit (Section III-B), which we detail next.

III-A Functional Stochastic Quasi-gradient Descent

Note that the functional gradient of the objective function $J(f)$ is given by

[TABLE]

where $\nabla_{f}\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}}(f(\boldsymbol{\xi}))\in\mathcal{H}\times\mathbb{R}^{m}$ and $\nabla\ell_{\boldsymbol{\theta}}\in\mathbb{R}^{m}$ which is defined as $\left[\nabla\ell_{\boldsymbol{\theta}}({\mathbf{u}})\right]_{i}=\frac{\partial\ell_{\boldsymbol{\theta}}({\mathbf{u}})}{\partial u_{i}}$. Observe that the stochastic version of the functional gradient is obtained by dropping the outer expectation as follows

[TABLE]

In vanilla stochastic gradient descent, one descends along the stochastic gradient (17), i.e., performs the update given by

[TABLE]

However, in (18), the stochastic gradient at a specific random variable $\boldsymbol{\xi}_{t}$, $\boldsymbol{\theta}_{t}$ is not available due to the expectation involved in the argument of $\nabla\ell_{\boldsymbol{\theta}_{t}}(\mathbb{E}[\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}}(f(\boldsymbol{\xi}))])$. This issue precludes use of vanilla stochastic gradient method for solving (1).

Thus, we propose using a two time-scale stochastic approximation strategy called stochastic quasi-gradient method [48, 22]. To estimate the inner-expectation, one may recursively average an instantaneous approximation of the inner expectation $\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f(\boldsymbol{\xi}_{t}))$ evaluated at $\boldsymbol{\xi}_{t}$ as [3]

[TABLE]

with the intent of estimating the expectation $\mathbb{E}_{\boldsymbol{\xi}}\left[\boldsymbol{\mathscr{h}}(f(\boldsymbol{\xi}))\right]$. By contrast, in this work, we consider a momentum analogue of this update instead

[TABLE]

Observe here that only a single observation $\xi_{t}$ is required per update, but two function evaluations of $\mathscr{h}$. In (19), $\beta$ is a scalar learning rate chosen from the unit interval $(0,1)$. This update may be interpreted as belonging to the family of gradient tracking schemes – see [29][p. 500-501] for illuminating discussion. Then, we define a function sequence $f_{t}\in{\mathcal{H}}$ initialized as null $f_{0}=0$ is sequentially updated using stochastic quasi-gradient descent:

[TABLE]

where $\alpha$ is a step-size parameter we fix as constant throughout. Further note that the term $\langle\nabla_{f}\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}({f_{t}(\boldsymbol{\xi}_{t})})$ is a function in ${\mathcal{H}}$, and thus infinite dimensional. However, by applying the chain rule and the reproducing property of the kernel [cf. (II)(i)], we obtain

[TABLE]

In (23), we have used the vector inner product notation to denote the summation in (22). Note that the kernel function $\kappa(\boldsymbol{\xi}_{t},\cdot)$ in (22) is common and therefore outside the inner product in (23). Utilizing this notation, the function update equation of (21) may be written as

[TABLE]

Observe that in (24), the vector $\boldsymbol{\mathscr{h}}^{\prime}({f_{t}(\boldsymbol{\xi}_{t})})$ associates with the sample point $\boldsymbol{\xi}_{t}$ evaluated by the kernel, whereas $\ell^{\prime}({\mathbf{g}}_{t+1})$ associates with the tracking parameter ${\mathbf{g}}_{t+1}$. Moreover, the function sequence in (24) belongs to a RKHS defined over the span of kernels $\{\kappa(\boldsymbol{\xi},\cdot)\}_{\boldsymbol{\xi}\in\Xi}$. Specifically, using the fact that $f_{0}=0\in{\mathcal{H}}$, one may obtain through induction that the function $f_{t}$ at time $t$ admits an expansion via kernel evaluations of past data realizations $\boldsymbol{\xi}_{n}$ and scalar weights $w_{n}$ for $n<t$:

[TABLE]

where we define weight vector ${\mathbf{w}}_{t}:=[w_{1},\cdots,w_{t-1}]^{T}$, kernel dictionary ${\mathbf{U}}_{t}=[\boldsymbol{\xi}_{1};\cdots;\boldsymbol{\xi}_{t-1}]$, and the empirical kernel map $\mathbf{\kappa}_{{\mathbf{U}}_{t}}({\mathbf{u}}):=[\kappa(\boldsymbol{\xi}_{1},{\mathbf{u}}),\cdots,\kappa(\boldsymbol{\xi}_{t-1},{\mathbf{u}})]$. Thus, performing the stochastic quasi-gradient iteration in the RKHS amounts to the following parametric updates on the coefficient vector ${\mathbf{w}}$ and kernel dictionary ${\mathbf{U}}$, given by

[TABLE]

In (III-A), observe the kernel dictionary parameterizing function $f_{t}$ is a matrix ${\mathbf{U}}_{t}\in{\mathbb{R}}^{p\times(t-1)}$ which stacks past realizations of random variable $\boldsymbol{\xi}$, and the coefficient vector ${\mathbf{w}}_{t}\in{\mathbb{R}}^{t-1}$ as the associated scalars in the kernel expansion (25) which are updated according to (III-A). Observe that the function update of (24) implies that the complexity of $f_{t}$ is ${\mathcal{O}}(t)$, due to the fact that the number of columns in ${\mathbf{U}}_{t}$, or model order $M_{t}$, is $(t-1)$, and thus is unsuitable for settings where the total number of data samples is not finite, or are arriving sequentially and repeatedly. This is an inherent challenge of extending [3] to optimizing over nonlinear functions that belong to RKHS. To address this, we consider projections of (24) onto low-dimensional subspaces, inspired by [17], which we detail in the following subsection.

III-B Subspace Projections for Complexity Control

In this subsection, we turn to address the untenable growth of the function representational complexity discussed in the previous section, namely, that the model order is $M_{t}=(t-1)$, and grows without bound with the iteration index $t$. To do so, we adopt the idea of bias-inducing proximal projections onto low-dimensional subspaces developed in [17].

Specifically, we construct an approximate sequence of functions by orthogonally projecting functional stochastic gradient updates onto subspaces ${\mathcal{H}}_{\mathbf{D}}\subseteq{\mathcal{H}}$ that consist only of functions that can be represented using some dictionary ${\mathbf{D}}=[{\mathbf{d}}_{1},\ \ldots,\ {\mathbf{d}}_{M}]\in{{\mathbb{R}}^{p\times M}}$, i.e., ${\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}$. Here we define ${{\mathbf{d}}_{n}\in{\mathbb{R}}^{p}}$ as a model point which stacks exemplar realizations of $\boldsymbol{\xi}$, i.e., ${\mathbf{d}}_{n}=\boldsymbol{\xi}_{n}$. Further define $\boldsymbol{{\kappa}}_{{\mathbf{D}}}(\cdot)=[{\kappa}({\mathbf{d}}_{1},\cdot)\ldots{\kappa}({\mathbf{d}}_{M},\cdot)]$, and ${\mathbf{K}}_{{\mathbf{D}},{\mathbf{D}}}$ as the resulting kernel matrix from dictionary ${\mathbf{D}}_{t}$. The dictionary ${\mathbf{D}}_{t}$ is updated with the new sample $\boldsymbol{\xi}_{t}$. Then to enforce projection, we replace the update in (21) with the following one which performs projection on to the subspace $\mathcal{H}_{{\mathbf{D}}_{t+1}}$:

[TABLE]

where $\mathcal{P}_{{\mathcal{H}}_{{\mathbf{D}}_{t+1}}}$ denotes the projection on to the subspace ${\mathcal{H}}_{{\mathbf{D}}_{t+1}}$.

In order to perform the projection in (III-B), let us denote the unprojected version of the function as $\tilde{f}_{t+1}$, dictionary $\tilde{{\mathbf{D}}}_{t+1}$ and the weights $\tilde{{\mathbf{w}}}_{t+1}$ given by

[TABLE]

This update may be represented parametrically as

[TABLE]

Given that we need to project $\tilde{f}_{t+1}$ onto the stochastic subspace ${\mathcal{H}}_{{\mathbf{D}}_{t+1}}$, for a fixed dictionary ${\mathbf{D}}_{t+1}$, the stochastic projection in (III-B) amounts to update the coefficient vector ${\mathbf{w}}_{t+1}$ as

[TABLE]

where ${\mathbf{K}}_{{\mathbf{D}}_{t+1},{\mathbf{D}}_{t+1}}$ and ${\mathbf{K}}_{{\mathbf{D}}_{t+1},{\tilde{\mathbf{D}}}_{t+1}}$ are the cross kernel matrices between $\{{\mathbf{D}}_{t+1},{\mathbf{D}}_{t+1}\}$ and $\{{\mathbf{D}}_{t+1},{\tilde{\mathbf{D}}}_{t+1}\}$, respectively. So the basic idea is to obtain a compressed version of $\tilde{f}_{t+1}$ and corresponding dictionary ${\mathbf{D}}_{t+1}$ and the coefficient vector ${\mathbf{w}}_{t+1}$. As previously noted, numerous approaches are possible for seeking this sparse representation. We make use of kernel orthogonal matching pursuit (KOMP) [49, Sec. 2.3] with allowed error tolerance $\epsilon$ to find a kernel dictionary matrix ${\mathbf{D}}_{t+1}$ based on the one which adds the latest sample point ${\tilde{\mathbf{D}}}_{t+1}$. This choice is due to the fact that we can tune its stopping criterion to guarantee a decrement property in expectation [17], as well as ensure the function complexity remains finite – see Sec. IV. We now describe the variant of KOMP we propose using, called Destructive KOMP with Pre-Fitting (see [49], Section 2.3), (Algorithm 2), which takes a candidate function $\tilde{f}$ of model order $\tilde{M}$ parameterized by its kernel dictionary ${\tilde{\mathbf{D}}}\in{\mathbb{R}}^{p\times\tilde{M}}$ and coefficient vector ${\tilde{\mathbf{w}}}\in{\mathbb{R}}^{\tilde{M}}$, and approximates $\tilde{f}$ by a function $f\in{\mathcal{H}}$ with a lower model order. Initially, this sparse approximation is the original function $f=\tilde{f}$ so that its dictionary is initialized with that of the original function ${\mathbf{D}}={\tilde{\mathbf{D}}}$, with corresponding coefficients ${\mathbf{w}}={\tilde{\mathbf{w}}}$. Then, the algorithm sequentially removes dictionary elements from dictionary ${\tilde{\mathbf{D}}}$, yielding a sparse approximation $f$ of $\tilde{f}$, until the error threshold $\|f-\tilde{f}\|_{{\mathcal{H}}}\leq\epsilon$ is violated, in which case it terminates.

Now we explain the execution of the KOMP algorithm in detail. At each stage of KOMP, a pair of dictionary element associated with index $j$ of ${\mathbf{D}}$ is selected to be removed which contributes the least to the Hilbert-norm approximation error $\min_{f\in{\mathcal{H}}_{{\mathbf{D}}\setminus\{j\}}}\|\tilde{f}-f\|_{{\mathcal{H}}}$ of the original function $\tilde{f}$, when dictionary ${\mathbf{D}}$ is used. Since at each stage, the kernel dictionary is fixed, this amounts to a computation involving weights ${\mathbf{w}}\in{\mathbb{R}}^{(M-1)}$ only; that is, the error of removing dictionary point ${\mathbf{d}}_{j}$ is computed for each $j$ as $\gamma_{j}={\min_{{\mathbf{w}}_{{\mathcal{I}}\setminus\{j\}}\in{\mathbb{R}}^{({M}-1)}}}\|\tilde{f}(\cdot)-\sum_{k\in{\mathcal{I}}\setminus\{j\}}w_{k}{\kappa}({\mathbf{d}}_{k},\cdot)\|.$ We use the notation ${\mathbf{w}}_{{\mathcal{I}}\setminus\{j\}}$ to denote the entries of ${\mathbf{w}}\in{\mathbb{R}}^{M}$ restricted to the sub-vector associated with indices ${\mathcal{I}}\setminus\{j\}$. Then, we define the dictionary element which contributes the least to the approximation error as $j^{\star}=\operatornamewithlimits{argmin}_{j}\gamma_{j}$. If the error associated with removing this kernel dictionary element exceeds the given approximation budget $\gamma_{j^{\star}}>\epsilon$, the algorithm terminates. Otherwise, this dictionary elements associated with ${\mathbf{d}}_{j^{\star}}$ are removed, the weights ${\mathbf{w}}$ are revised based on the pruned dictionary as ${\mathbf{w}}={\operatornamewithlimits{argmin}_{{\mathbf{w}}\in{\mathbb{R}}^{{M}}}}\lVert\tilde{f}(\cdot)-{\mathbf{w}}^{T}\boldsymbol{{\kappa}}_{{\mathbf{D}}}(\cdot)\rVert_{{\mathcal{H}}}$, and the process repeats as long as the current function approximation is defined by a nonempty dictionary. 333We assume that the output of Algorithm 2 has bounded Hilbert norm, which may be explicitly enforced, for instance, by thresholding the norm of the weight vector if it climbs above some large constant.

With the subspace projection procedure in place, we may summarize the key steps of the proposed method in Algorithm 1 for solving (1) while maintaining a finite model order, thus breaking the “curse of kernelization” for compositional stochastic programming over an RKHS. The method, Compositional Online Learning with Kernels (COLK), executes the stochastic subspace projection of a functional stochastic quasi-gradient step onto sparse subspaces ${\mathcal{H}}_{{\mathbf{D}}_{t+1}}$ stated in (III-B). The initial function is set to null $f_{0}=0$, meaning that it has empty kernel dictionary ${\mathbf{D}}_{0}=[]$ and coefficient vector ${\mathbf{w}}_{0}=[]$. The notation $[]$ is used to denote the empty matrix or vector respective size $p\times 0$ or [math]. Then, at each step, given an independent training example $\boldsymbol{\xi}_{t}$ and step-size $\eta_{t}$, we update the auxiliary variable ${\mathbf{g}}_{t+1}=(1-\beta)({\mathbf{g}}_{t}-\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f_{t-1}(\boldsymbol{\xi}_{t})))+\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f_{t}(\boldsymbol{\xi}_{t}))$. Then, this updated scalar is used to compute the unconstrained functional stochastic gradient iterate $\tilde{f}_{t+1}={(1-\lambda\alpha)f_{t}-{\alpha\langle{\boldsymbol{\mathscr{h}}^{\prime}_{\boldsymbol{\xi}_{t}}(f_{t}(\boldsymbol{\xi}_{t}))},\ell^{\prime}_{\boldsymbol{\theta}_{t}}({\mathbf{g}}_{t+1})\rangle\kappa(\boldsymbol{\xi}_{t},\cdot)}}$, which admits the parametric representation ${\tilde{\mathbf{D}}}_{t+1}$ and ${\tilde{\mathbf{w}}}_{t+1}$ as stated in (31). These parameters are then fed into KOMP with approximation budget $\epsilon$, such that

[TABLE]

Next, we establish conditions and parameter selections under which Algorithm (1) solves the compositional problem in (1) to near-optimality while ensuring that the memory of the learned function remains finite.

IV Convergence Analysis

We establish the convergence of Algorithm 1, which is characterized in terms of the behavior of the function iterates $f_{t}$ with respect to the minimizer of (1), as well as its finite-memory properties. Our analysis is inspired by the notion of coupled supermartingales illuminated in [3]. There are two key points of departure in our analysis: (i) the projection step introduces a finite error to the gradient estimation at each time step in pursuit of low-dimensional function representations which then determines the $\epsilon$-factor in the radius/rate of convergence; (ii) the use of a difference-of-iterates scheme in the faster time-scale (20) yields tighter dependences on problem constants in the convex setting, and improved rates to stationarity in non-convex settings. The specific effect of (ii) is summarized in Lemma 1, which is key to deriving these results.

To proceed with the analysis, a few quantities are first defined to simplify the exposition. Firstly, define the functional stochastic quasi-gradient of the objective function (cf. (2))

[TABLE]

and corresponding compressed projected version is written as

[TABLE]

Utilizing this notation, the main function update step in (30) can be written as

[TABLE]

which is an update for the function $f_{t+1}$ re-written as a step involving projected gradients. These definitions are needed to clarify the technical setting, which we do next. Begin by defining as ${\mathcal{F}}_{t}$ the filtration, i.e., the time-dependent sigma algebra containing the algorithm history ${\mathcal{F}}_{t}\supset(\{f_{u},g_{u}\}_{u=0}^{t}\cup\{\boldsymbol{\theta}_{s},\boldsymbol{\xi}_{s}\}_{s=0}^{t-1})$. With this definition of algorithm history contained by ${\mathcal{F}}_{t}$, we state the assumptions as follows.

Assumption 1

The reproducing kernel map is bounded as

[TABLE]

Assumption 2

At each time instant $t$, the second moment of the derivative of inner function ${{\boldsymbol{\mathscr{h}}^{\prime}}_{\boldsymbol{\xi}_{t}}(f(\boldsymbol{\xi}_{t}))}$ and outer function ${{\ell^{\prime}_{\boldsymbol{\theta}_{t}}\left({\mathbf{g}}_{t+1}\right)}}$ is bounded as

[TABLE]

with $G_{{\mathscr{h}}}$ and $G_{\ell}$ as finite constants with probability $1$ (w.p.$1$).

Assumption 3

The instantaneous derivative of the outer function ${\ell_{\boldsymbol{\theta}}(\cdot)}$ is Lipschitz continuous with respect to its first scalar argument for each realization of $\boldsymbol{\theta}$ so that we may write

[TABLE]

Assumption 4

The functions ${\mathbf{h}}$ is Lipschitz continuous and has bounded variance, i.e.,

[TABLE]

for all $f,f^{\prime}\in\mathcal{H}$. In addition, define ${\boldsymbol{\delta}_{t}:=\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f_{t};\boldsymbol{\xi}_{t})}$ with \bar{\boldsymbol{\delta}}_{t}=\mathbb{E}\left[\boldsymbol{\delta}_{t}{\,\big{|}\,}\boldsymbol{\theta}_{t}\right]. Then $\boldsymbol{\delta}_{t}$ has finite variance as \mathbb{E}[\|\boldsymbol{\delta}_{t}-\bar{\boldsymbol{\delta}}_{t}\|^{2}{\,\big{|}\,}{\mathcal{F}}_{t}]\leq\sigma_{\boldsymbol{\delta}}^{2}\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{ w.p. }1.}

Assumption 5

(Bounded Gradients) The instantaneous gradient of the outer function ${\ell^{\prime}_{\boldsymbol{\theta}}\left({\mathbf{u}}\right)}$ is almost surely bounded $|{\ell^{\prime}_{\boldsymbol{\theta}}\left({\mathbf{u}}\right)}|\leq C_{\ell}$.

Assumption 1 follows from the compactness of the feature space $\boldsymbol{\Xi}\cup\boldsymbol{\Theta}$. Assumption 2 is regarding the second moments of the derivatives which limits the variance of the stochastic approximation error, and is typical in the literature [50]. Assumption 3 and 4 regarding the Lipschitz continuity of the outer and inner function holds for most applications, and holds for most differentiable functions. Note that Assumption 5 is automatically satisfied for most smooth convex functions when data domains are compact. These assumptions are standard and hold for the applications in the next section.

With the technical setting clarified, we are ready to state our theoretical results. That is, first we bound the model order (dictionary size) of the proposed algorithm in terms of the step size $\alpha$ and compression budget $\epsilon$. Specifically, we extend [17][Theorem 3] to compositional settings, and further provide a non-asymptotic dependence of the function complexity on problem constants and the parameter dimension.

Theorem 1

Consider constant step sizes $\alpha$, $\beta$, with finite compression budget $\epsilon$, and regularization parameter as $\lambda=L_{\ell}U^{2}G_{h}^{2}\alpha/\beta+\lambda_{0}=\mathcal{O}(\alpha/\beta+1)$ with $\lambda_{0}<\min\{1,1/\alpha\}$. Let $M_{t}$ denote the model order of the function iterate $f_{t}$ i.e. number of columns in current dictionary $D_{t}$, then under Assumptions 1-5, $M_{t}$ satisfies

[TABLE]

Theorem 1 (see See Appendix A for proof) ensures that the number of data points in the kernel representation of $f_{t}$ generated from Algorithm 1 is finite and dependent upon the ratio of the step-size to the compression parameter. The condition for online sparsification performed by KOMP algorithm boils down to the condition for packing number of the kernelized feature space $\phi(\mathcal{X})$ as described in (50). The packing of kernelized feature space is inversely proportional to the radius $\sqrt{\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}}$. As this radius increases, the packing number reduces, meaning the model order required to cover feature space decreases. A larger radius may be attained by choosing a larger $\epsilon$, meaning that fewer points are required to cover the data domain, and thus yields a lower model order, at the cost of incurring greater sub-optimality.

As presented in Theorem 1, a similar gaurantees exist in the literature but for the non-compositional objective functions. The authors in [35] used a fixed memory size for the function representation but it leads to a constant term as the gradient error in the convergence rate analysis. An online algorithm for streaming data points is proposed in [51] but the memory requirement is linear in $T$. Another functional gradient based technique in [17] achieves a model complexity of $\Omega(T^{p})$. Further, a lower bound on the number of required features for a random feature based kernel methods is established in [52] to be $\Omega(\sqrt{T}\log(d_{\lambda}))$. Here, $d_{\lambda}$ is the effective dimension if the problem and depends upon the kernel matrix. It is difficult to calculate $d_{\lambda}$ in the streaming data applications, and hence deciding the number of features for the problem is difficult. All the existing results are summarized in Table I. A key upshot of our analysis is that it is the first to provide complexity bounds for RKHS function approximations for compositional objectives, both for strongly convex (SC) and non-convex (NC) cases.

We begin establishing a basic property of Algorithm 1, which is inspired by, but distinct from [3]. In particular, the gradient tracking scheme (20) enables us to obtain tighter bounds on the estimation error of the sequence $\{{\mathbf{g}}_{t}\}$ presented as follows.

Lemma 1

Consider ${\boldsymbol{\delta}_{t}:=\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f_{t};\boldsymbol{\xi}_{t})}$ and \bar{\boldsymbol{\delta}}_{t}=\mathbb{E}\left[\boldsymbol{\delta}_{t}{\,\big{|}\,}\boldsymbol{\theta}_{t}\right] as defined in Assumption 4. Then the estimation error of the auxiliary sequence ${\mathbf{g}}_{t}$ [cf. (20)] with respect to $\bar{\boldsymbol{\delta}}_{t}$ satisfies

[TABLE]

where $\beta\in(0,1)$ is a step-size parameter.

The proof of Lemma 1 is provided in Appendix B.

Strongly Convex Objectives Next, we characterize the convergence rate and the memory size requirements when the objective $R(f)$ [cf. (2)] is strongly convex, which in light of Theorem 1, yields a trade-off between the upper-bound on complexity $\mathcal{M}$ and convergence.

Theorem 2

Consider constant step-sizes $\alpha,\beta$, and constant compression budget $\epsilon$ with regularization parameter as $\lambda=L_{\ell}U^{2}G_{h}^{2}\alpha/\beta+\lambda_{0}=\mathcal{O}(\alpha/\beta+1)$ with $\lambda_{0}<\min\{1,1/\alpha\}$. Then when $R(f)$ [cf. (2)] is strongly convex,

• (i)

the iterates of Algorithm 1 converges to a neighborhood of the optimal $f^{\star}$ as

[TABLE]

when run with constant step-sizes $\alpha>0$ and $0<\beta<1$, with problem-dependent constants defined as $Z_{1}=(\sigma_{f}^{2}+2L_{\mathscr{h}}^{2}L_{\ell}U^{2}\sigma_{f}^{2})/\lambda_{0}$, $Z_{2}\!=\!(2L_{\ell}U^{2}\sigma_{\boldsymbol{\delta}}^{2})/\lambda_{0}$, $Z_{3}=(2Z)/\lambda_{0}$, and $Z_{4}\!=\!4L_{\mathscr{h}}^{2}L_{\ell}U^{2}/\lambda_{0}$.

• (ii)

Moreover, to achieve $\delta$ sub-optimality, we require at least $T$ samples and an model complexity upper bound $\mathcal{M}$ as

[TABLE]

for step-sizes $\alpha=\frac{\delta}{5Z_{1}}$, $\beta=\frac{\delta}{5\sqrt{Z_{1}Z_{2}}}$, budget \epsilon\!\leq\!\min\Big{\{}\!\frac{\delta Z_{4}^{1/2}}{Z_{1}^{3/2}},\frac{\delta^{2}}{25Z_{1}Z_{3}}\!\Big{\}}.

See Appendix C for proof. The result in statement (i) of Theorem 2 is comparable to well-known finite sample analysis of stochastic gradient algorithms, which converge to a $\mathcal{O}(\alpha)$ neighborhood of the optimal value for step-size $\alpha$. However, due to the complexity of RKHS parameterizations, we restrict focus to constant compression budget $\epsilon>0$ and study how controlled is the function’s representational complexity. Doing so allows us to obtain specific trade-offs between complexity and convergence, whereas $\epsilon=0$ implies the complexity of the function parameterization grows unbounded with time $T$. A further point of contrast is with approaches that fix the complexity, and may cause divergence [53].

Instead of fixing the complexity a priori, the convergence accuracy parameter $\delta$ determines the number of iterations $T$ the algorithm must run, as well as the lower bound on the compression budget $\epsilon$, to obtain a $\delta$-suboptimal solution. Doing so then yields a specific upper bound on the model order. In doing so, we provide for the first time a non-asymptotic tradeoff between complexity and sub-optimality of solving compositional problems over an RKHS, and thus address a strictly more general setting than existing concentration bounds using random Fourier Features, whose focus is restricted to standard stochastic convex programming [54, 55, 52]. We are able to address the compositional setting in contrast to these works due to the fact that our analysis only depends on function smoothness rather than tying distributional properties of the stochastic approximation error to the eigenvalue decay of the population kernel matrix, as in the aforementioned references.

Non-Convex Objectives As previously noted, many practical risk measures may be non-convex – see Example 1 and [44]. Therefore, we shift focus to analyzing algorithm performance for the setting where $R(f)$ (2) is non-convex.

Theorem 3

For $R(f)$ [cf. (2)] non-convex in $f$,

• (i)

the function sequence $\{f_{t}\}$ generated by Algorithm 1 run for $T$ total iterations with null initialization $f_{0}=0$ converges to an approximate stationary solution as

[TABLE]

under step-sizes $\alpha>0$, $0<\beta<1$, and budget $\epsilon$.

• (ii)

Thus, to achieve an $\delta$-approximate stationary solution, the iteration complexity is $T\geq\mathcal{O}\left(\frac{1}{\delta^{2}}\right)$ and $\mathcal{M}\geq\mathcal{O}\left(\frac{1}{\delta^{p}}\right)$ with step size $\alpha=\mathcal{O}(\delta)$, $\beta=\mathcal{O}(\delta)$, and the compression budget $\epsilon\leq\mathcal{O}(\delta^{3/2})$.

The proof of this theorem is provided in Appendix D. Observe that for $\epsilon=0$, $\mathcal{M}=T$ which grows unbounded with time since we add each new sample to the kernel dictionary, and the right-hand side of (44) reduces to $\mathcal{O}{(T^{-1/2})}$ because $p\geq 1$ and matches existing results [56]. By contrast, for fixed $\epsilon>0$, we can provide a specific tradeoff between convergence accuracy to stationarity and the required complexity of the function representation (Theorem 3(ii), which does not require $\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}}$ to be smooth. Next, we investigate the empirical validity of Algorithm 1 on several practical problems.

V Experiments

To show the efficacy of the proposed algorithm, we consider a problem of nonlinear regression (filtering) over a $p$-dimensional parameter space. We have again have two sets of random variables $({\mathbf{x}},{\mathbf{x}}^{\prime})\in{\mathcal{X}}\subset{\mathbb{R}}^{p}$ but now the target variables are real valued $y,y^{\prime}\in{\mathbb{R}}$. The merit criterion of model fitness for a given training example $({\mathbf{x}}_{n},y_{n})$ is the humble square loss:

[TABLE]

We consider as a surrogate for the approximation error the $p$-th order central moments, as distributions may be completely characterized by their moments [57][Ch. 3]

[TABLE]

which are just deviations raised to the $p$-th power [9]. However, since the computational overhead scales with $P$, we truncate the upper summand index in (46) to $P=4$. Note that the dispersion measure in (46) is non-convex which is used for the experimental purposes which corresponds to the variance, skewness, and kurtosis of the loss distribution. Note that (46) may be convexified through a positive projection of $(\ell(f\!({\mathbf{x}}),\!{\mathbf{y}})\!-\!{\mathbb{E}_{{\mathbf{x}}^{\prime},{\mathbf{y}}^{\prime}}\![\ell\!(f\!({\mathbf{x}}^{\prime}),\!{\mathbf{y}}^{\prime})]})$, in which case the standard deviation becomes a semi-deviation, as do its higher-order analogues [58]. This is the risk function we use for experiments in correspondence to (8) provided for example 1 in Sec.II. However, for simplicity, we omit the positive projection in experiments. Next, we apply the proposed algorithm to solve the nonlinear regression problem which results in the following updates.

[TABLE]

where $m(p,g_{t+1}):=\left[p((f_{t}({\mathbf{x}}_{t})-y_{t})^{2}-g_{t+1})^{p-1}\right]$. Note that the gradient of the outer function is not bounded in the above mentioned problem, but the gradient of the projected version after applying KOMP is bounded.

We evaluate COLK on synthetic and real data sets whose distributions are skewed or heavy-tailed, and compare its test accuracy against existing benchmarks that minimize only bias. Firstly, we evaluate performance on synthetic data regression outliers which has a heavier tailed distribution, i.e., more outliers are present. We inquire as to which methods overfit versus learn successfully: COLK (Algorithm 1), or methods such as BSGD [35], NPBSG [59], POLK [17]. We further consider extensions of BSGD for compositional objectives, as well as a radial basis function (RBF) network of fixed size with $50$ basis functions [3, 56]. We note that [56] is slightly different in that it considers Jacobi rather than Gauss-Seidel-style updates for the auxiliary and main variables, in contrast to SCGD, which incurs fewer calls to a simulation oracle.. These basis functions are $\{\phi_{i}({\mathbf{x}}_{j})+n_{j}\}$ for $j=1$ to $N$ data points ${\mathbf{x}}_{j}$. Our goal in doing so is to elucidate the relationship between a fixed basis representation and the one obtained by the projection step of COLK. Of course, the nonlinearity could be defined by polynomials, sinusoidals, etc. [60]. We consider the RBF for uniformity of experimental comparison, which takes the form $\phi_{i}(x)=\exp(-(x-c_{i})^{2}/r^{2})$ where $c_{i}$ are the uniformly distributed centers points across the feature space ${\mathcal{X}}$ with bandwidth $r=0.06$

We generate $20$ different sets from the same data distribution and then run both POLK and COLK to learn a regression function. To generate the synthetic dataset regression outliers, we used the function $y=2x+3\text{sin}(6x)$ as the original function and target $y$’s observed after adding a zero mean Gaussian noise to $2x+3\text{sin}(6x)$. First we generate $60000$ samples of the data, and then select $20\%$ as the test data set. From the remaining $4800$ samples, we select $50\%$ at random to generate $20$ different training sets. We run COLK over these training set with the following parameter selections: a Gaussian kernel with bandwidth $\sigma=.06$, step-size parameters $\alpha=0.02$, $\beta=0.01$, $\epsilon=K\alpha^{2}$ with parsimony constant $K=5$, variance coefficient $\eta=0.1$, and mini-batch size of $1$. Similarity, for POLK we use $\alpha=0.5$ and $\epsilon=K\alpha^{2}$ with parsimony constant $K=0.09$. We fix the kernel type and bandwidth across the different methods, and the parameters that define comparator algorithms are hand-tuned to optimize performance with the restriction that their model complexity is comparable to each other, with the exception that we choose the number of RBF features in the RBF network to be what is required to obtain comparable performance to COLK. We run these algorithms for different realizations of training data and evaluate their test accuracy as well as its standard deviation.

We present the convergence behavior of Algorithm 1 as well as these alternatives on the regression outliers data in Fig 1. Specifically, Fig. 1(a) shows that the mean plus variance of the loss function is minimized as the number of samples processed increases. Observe that COLK obtains comparable sub-optimality to the RBF network, which is smaller than the comparators. Moreover, the time-series of the test-set error of COLK is given in Fig. 1(b) which converges as the training samples increases. In Fig. 1(c) we plot the model order of the function sequence defined by COLK, and observe it stabilizes over time regardless of the presence of outliers.

Observe that in order to obtain comparable performance to COLK, one requires more than twice as many RBF features (Fig. 1(c)). Moreover, BSGD and other budgeted approaches yield lower accuracy estimators. These preliminary results validate the convergence results established in Section IV. The advantage of minimizing the bias as well as variance is depicted in Fig.  1(d) which plots the learned function for POLK and COLK for two training data sets. It can be observed that how POLK learning varies from one training set to other while COLK is robust to this change. This trend is corroborated in terms of evaluating COLK for $20$ total training runs and reporting the average test error and the standard deviation as the box and whisker plot given in Fig. 2, where we also display the result of the comparators. Observe that COLK yields the lowest error as well as the lowest standard deviation, meaning it yields inferences that are both low bias and low variance.

To check the proposed algorithm for real data, we consider the performance of filtering laser scans to interpolate range to a target via the lidar data [38] (with added outliers) with the results shown in Fig. 3. It is clear from the figure that the proposed algorithm is robust to outliers in the data, whereas alternatives fall short in terms of either representational efficiency or overfitting.

VI Conclusion

In this work, we addressed compositional stochastic programming in Reproducing Kernel Hilbert Space by developing a functional generalization of the stochastic quasi-gradient method operating in parallel with greedy subspace projections. This method, Compositional Online Learning with Kernels (COLK), converges both under attenuating and constant learning rates, and yields memory-efficient parameterizations. Different from classic stochastic quasi-gradient, we introduce a momentum scheme to tighten the estimation error of the faster time-scale.

We experimentally observed both with synthetic and benchmark data that COLK applied to robust supervised learning overcomes the problem of overfitting: by accounting for error variance through coherent risk, we observed consistent performance across training runs. In future work, we hope to investigate the use of hierarchical kernels for larger parameter problems in vision or acoustics, and thus design deep learners that do not overfit.

Appendix A Proof of Theorem 1

The proof of Theorem 1 follows from the proof of [17, Theorem 4] which provides an asymptotic guarantee which establishes finite model order under constant step size $\alpha$ for the standard convex objective (non-compositional). But in this work, we are interested in characterizing the rate at which the model order (or kernel dictionary size) varies for the COLK algorithm with the number of iterations $T$. To achieve that, the initial steps of the proofs are similar to [17, Theorem 4] which we skip here (full derivation is provided in Sec. F of supplementary for completeness) and directly start from the inequality at iteration $t$

[TABLE]

where $\text{dist}(\kappa(\boldsymbol{\xi},\cdot),{\mathcal{H}}_{{\mathbf{D}}})=\|\kappa(\boldsymbol{\xi},\cdot)-[{\mathbf{K}}_{{\mathbf{D}}_{t},{\mathbf{D}}_{t}}^{-1}\boldsymbol{\kappa}_{{\mathbf{D}}_{t}}(\boldsymbol{\xi})]^{T}\boldsymbol{\kappa}_{{\mathbf{D}}_{t}}(\cdot)\|_{{\mathcal{H}}}$ and $V_{t}=\ell^{\prime}_{\boldsymbol{\theta}_{t}}\left({\mathbf{g}}_{t+1}\right){\boldsymbol{\mathscr{h}}^{\prime}_{\boldsymbol{\xi}_{t}}(f_{t}(\boldsymbol{\xi}_{t}))}$. In [17, Theorem 4], the inequality of (48) is obtained for the case when the current sample $\boldsymbol{\xi}_{t}$ is not added to the kernel dictionary $D_{t}$. This implies that when the distance $\text{dist}({\kappa}(\boldsymbol{\xi}_{t},\cdot),{\mathcal{H}}_{{\mathbf{D}}_{t}})$ for the new sample $\boldsymbol{\xi}_{t}$ is less than equal to the term $\frac{\epsilon}{\alpha|V_{t}|}$, then the model order satisfy $M_{t+1}\leq M_{t}$.

Now, consider the contrapositive of the preceding expressions. Observe that the model order growth condition ($M_{t+1}=M_{t}+1$) implies that

[TABLE]

holds. This condition establishes the fact that every time a new data ($\boldsymbol{\xi}$) is appended to kernel dictionary, then the associated product kernel is guaranteed to be at least a distance of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\epsilon}{\alpha|V_{t}|}} from every other kernel function in the current model. Now utilizing the Cauchy Schwartz inequality and Assumption 5 we get

[TABLE]

This upper bound implies that $1/|V_{t}|\geq 1/(C_{\ell}L_{\mathscr{h}})$, therefore we can lower bound the right hand side in (49) as follows

[TABLE]

From (129) (in supplementary), we obtain

[TABLE]

Therefore, the KOMP stopping criterion is violated for the newest point whenever distinct dictionary points ${\mathbf{d}}_{k}$ and ${\mathbf{d}}_{j}$ for $j,k\in\{1,\dots,M_{t}\}$, satisfy the condition $\|\phi({\mathbf{d}}_{j})-\phi({\mathbf{d}}_{k})\|_{2}>{\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}}$. Next, we proceed in a similar manner to that of Theorem 3.1 in [61]. Note that since the space $\mathcal{U}$ is compact and $\kappa$ is continuous, the range $\phi({\mathcal{U}})$ (where $\phi({\mathbf{u}})=\kappa({\mathbf{u}},\cdot)$ for ${\mathbf{u}}\in{\mathcal{U}}$) of the kernel transformation of feature space ${\mathcal{U}}$ is compact. This allows us to conclude that the number of balls of radius $\chi$ (here, $\chi=\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}$) required to completely cover the set $\phi({\mathcal{U}})$ is finite (see, e.g., [62]).

To prove the main result of (40), we consider the result in [61, Proposition 2.2] which states that for a Lipschitz continuous Mercer kernel $\kappa$ on compact set $\mathcal{X}\subseteq\mathbb{R}^{p}$, for any training set $\{{\mathbf{x}}_{t}\}_{t=1}^{\infty}$ and any $\nu>0$, the number of elements in the dictionary is upper bounded as

[TABLE]

where $Y$ is a constant depends upon $\mathcal{X}$ and the kernel function. From the result in (51), we conclude that $\sqrt{\nu}={\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}}$, which we may substitute into (52) to obtain

[TABLE]

as stated in (40). Note that in (53), we have defined a constant $Y^{\prime}=YC_{\ell}^{2p}L_{\mathscr{h}}^{2p}$. We remark that there is a trivial lower bound of $M_{t}\geq 1$. ∎.

Appendix B Proof of Lemma 1

The proof employs a notion of gradient tracking for tightening the rate of convergence [31][Lemma 2], although for an entirely different context: the error in the estimate of the inner expectation, rather than a stochastic gradient estimate error. Whereas in [31], the stochastic gradient sub-sampling error is the primary target hit by the auxiliary use of gradient information, instead here it is the inner-objective evaluation itself. Begin by expanding ${\mathbf{g}}_{t+1}-\bar{\boldsymbol{\delta}}_{t}$ from (20), breaking the terms that multiply $1-\beta$ in two, and then add and subtract $(1-\beta)\bar{\boldsymbol{\delta}}_{t-1}$ to obtain

[TABLE]

where we have gathered like terms. Now, observe that ${\mathbb{E}}[\bar{\boldsymbol{\delta}}_{t-1}-\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f_{t-1}(\boldsymbol{\xi}_{t}))~{}|~{}\mathcal{F}_{t}]={\mathbb{E}}[\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f_{t}(\boldsymbol{\xi}_{t}))-\bar{\boldsymbol{\delta}}_{t}~{}|~{}\mathcal{F}_{t}]=0$ by the definition of $\delta_{t}$ and $\bar{\delta}_{t}$ in Assumption 4. Computing the norm followed by the expectation conditional on ${\mathcal{F}}_{t}$ then yields

[TABLE]

where the cross term vanishes since ${\mathbf{g}}_{t}-\bar{\boldsymbol{\delta}}_{t-1}$ is independent of $\xi_{t}$ and the second term in the cross term is zero mean. The second term in (B) can again be expanded as

[TABLE]

where we have used the inequality $\mathbb{E}\left[\left\|\mathsf{X}-\mathbb{E}\left[\mathsf{X}\right]+\mathsf{Y}\right\|^{2}\right]\leq 2\mathbb{E}\left[\left\|\mathsf{X}\right\|^{2}\right]+2\mathbb{E}\left[\left\|\mathsf{Y}\right\|^{2}\right]$ for any random variables $\mathsf{X}$ and $\mathsf{Y}$ with bounded variances. Take full expectation in (56) and apply (39) to the first term on the right-hand side of the preceding expression. Now, by substituting the resulting expression into the right-hand side of (B), we obtain

[TABLE]

where the last inequality (58) follows from the fact that $\beta\leq 1$ applied to the second term on the right-hand side of (57). $\blacksquare$

Appendix C Proof of Theorem 2

Before analyzing the mean convergence behavior of the distance from optimal $\mathbb{E}\left[\|f_{t+1}-f^{\star}\|_{{\mathcal{H}}}^{2}\right]$, consider the following temporal difference between ${\mathbf{g}}_{t+1}$ and $\bar{\boldsymbol{\delta}}_{t}$ from Lemma 1

[TABLE]

Next, substitute the upper bound on the total expectation of $\|f_{t}-f_{t-1}\|_{{\mathcal{H}}}^{2}$ as described in Lemma 2 (i) (stated in Appendix E of the supplementary) into (C) to obtain

[TABLE]

It is interesting to observe that (C) relates $\mathbb{E}\left[\|{\mathbf{g}}_{t+1}-\bar{\boldsymbol{\delta}}_{t}\|^{2}\right]$ to its previous iterate value. Utilizing this recursion, we can write the following by replacing $t+1$ by $t$

[TABLE]

Substituting (C) into the right-hand side of (C) and repeating the recursion, we can write

[TABLE]

The first term in (62) vanishes due to the initialization $g_{0}=0$ and the convention $\boldsymbol{\delta}_{-1}=0$. Moreover, the second term represents a geometric series and sum can be evaluated using $\sum_{u=0}^{t}(1-\beta)^{u}=[1-(1-\beta)^{t}]/\beta$ provided $\beta<1$. In this geometric sum expression, since the numerator is strictly less than unit, we can further simplify (62) to

[TABLE]

After establishing this auxiliary sequence order in terms of the order of step sizes $\alpha$ and $\beta$, we shift our focus again to the sub-optimality gap $\|f_{t}-f^{\star}\|_{{\mathcal{H}}}$ in expectation. Before proceeding, note that the Hilbert-norm regularizer $(\lambda/2)\|f\|_{{\mathcal{H}}}^{2}$ in (2) makes the objective $R(f)$ strongly convex, i.e.

[TABLE]

Using the inequality in (64) into the expression of Lemma 2(ii) (stated in Appendix E), we obtain

[TABLE]

Now, take the total expectation of (C). Considering that $\left\|f^{\star}\right\|\leq S$ which implies $\|f_{t}-f^{\star}\|_{{\mathcal{H}}}\leq\|f_{t}\|_{{\mathcal{H}}}+\|f^{\star}\|_{{\mathcal{H}}}\leq K+S:=Z$. Now, substitute the regularizer selection $\lambda=L_{\ell}U^{2}G_{h}^{2}\alpha/\beta+\lambda_{0}$ for $\lambda_{0}<1$, and apply (63) to the last term on the right-hand side of the preceding expression:

[TABLE]

Recursively substitute this expression back into itself to obtain:

[TABLE]

Next, since $f_{1}=0$ and using $\sum_{u=0}^{t-1}(1-\alpha\lambda_{0})^{u}=[1-(1-\alpha\lambda_{0})^{t-1}]/\alpha\lambda_{0}\leq 1/\alpha\lambda_{0}$ since $\alpha\lambda_{0}<1$. We get the following result by substituting $t=T$

[TABLE]

where $Z_{1}=(\sigma_{f}^{2}+2L_{\mathscr{h}}^{2}L_{\ell}U^{2}\sigma_{f}^{2})/\lambda_{0}$, $Z_{2}=(2L_{\ell}U^{2}\sigma_{\boldsymbol{\delta}}^{2})/\lambda_{0}$, $Z_{3}=(2Z)/\lambda_{0}$, and $Z_{4}=4L_{\mathscr{h}}^{2}L_{\ell}U^{2}/\lambda_{0}$. From the definition of $\mathcal{M}$ in Theorem 1, we have $\mathcal{M}=Y^{\prime}\left(\frac{\alpha}{\epsilon}\right)^{2p}$, which implies that we can write (C) as

[TABLE]

Next, in order to characterize the rates for iteration complexity as well as model complexity, it is possible to choose the step sizes $\alpha$ and $\beta$ in an optimal manner. Suppose, we require $\mathbb{E}\left[\|f_{T+1}-f^{\star}\|_{{\mathcal{H}}}^{2}\right]$ to be less than some positive value $\delta$ which describes the accuracy for the algorithm. To do so, we ensure that each terms in the right hand side of (C) is bounded by $\delta/5$. In other words, we need to satisfy the following set of inequalities

[TABLE]

The inequalities in (70) can be satisfied if we set $\alpha=\frac{\delta}{5Z_{1}}$, $\beta=\frac{\delta}{5\sqrt{Z_{1}Z_{2}}}$ for small $\delta$ so that $-\log(1-\lambda_{0}\alpha)\approx\lambda_{0}\alpha$ and

[TABLE]

as stated in Theorem 2. From the definition of $\mathcal{M}$, we have $\mathcal{M}=Y^{\prime}(\alpha/\epsilon)^{2p}$ it holds that

[TABLE]

Using the upper bound on $\alpha$ from (71), we substitute $\alpha=\frac{\delta}{5Z_{1}}$ into (72), we get

[TABLE]

After rearranging the terms for $\epsilon$ in (73), we obtain \epsilon\!\leq\!\min\Big{\{}\!\frac{\delta Z_{4}^{1/2}}{Z_{1}^{3/2}},\frac{\delta^{2}}{25Z_{1}Z_{3}}\!\Big{\}}. $\blacksquare$

Appendix D Proof of Theorem 3

Note that the statement Lemma 1 holds independent of the convexity of the objective. Thus, consider the statement of Lemma 1 (stated in Appendix E of the supplementary) and multiply the both sides by $(1+\beta)$, then add to the statement of Lemma 3 to obtain

[TABLE]

by the observation that $1-\beta^{2}\leq 1$, $1+\beta\leq 2$. For the better exposition, let us define a stochastic process $\gamma_{t}=R(f_{t})+\left\|\bar{\boldsymbol{\delta}}_{t-1}-\mathbf{g}_{t}\right\|^{2}$. We can write

[TABLE]

Taking expectation on both sides and noting that $\mathbb{E}[\left\|f_{t}-f_{t-1}\right\|^{2}_{\mathcal{H}}]=\mathbb{E}[\alpha^{2}\mathbb{E}[\left\|\tilde{\nabla}_{f}R(f_{t-1})\right\|^{2}_{\mathcal{H}}\mid\mathcal{F}_{t-1}]]\leq\alpha^{2}\sigma^{2}_{f}$, we obtain

[TABLE]

After rearranging the terms, we get

[TABLE]

Taking sum from $t=1$ to $T$, and then lower bounding the left hand side of (D), we obtain

[TABLE]

Rearranging the terms, we get

[TABLE]

where, $C_{1}=\mathbb{E}[\gamma_{1}]-R(f^{\star})$, $C_{2}=4\sigma^{2}_{\delta}$, $C_{3}=\frac{\sigma^{2}_{f}L_{R}}{2}$, and $C_{4}=4L_{\mathscr{h}}^{2}\sigma^{2}_{f}$. From the definition of $\mathcal{M}$ in Theorem 1, we have $\mathcal{M}=Y^{\prime}\left(\frac{\alpha}{\epsilon}\right)^{2p}$, which implies that we can write (79) as

[TABLE]

To obtain the complexity bounds presented in Theorem 1 statement (ii), we need to consider the case when $\underset{1\leq t\leq T}{\inf}\mathbb{E}[\left\|\nabla_{f}R(f_{t})\right\|^{2}_{\mathcal{H}}]$ is less than the required accuracy $\delta$. From (80), note that $\delta$ accuracy is achieved if it holds that

[TABLE]

Hence for $\alpha=\frac{\delta}{4}$, $\beta\leq\frac{\delta}{4}$, $\mathcal{M}=Y^{\prime}({\alpha}/{\epsilon})^{2p}$, and from the lower bound in (81), it holds that

[TABLE]

Utilizing $\alpha=\frac{\delta}{4}$ and after rearranging the terms in (82), we obtain $\left({\delta}/{4\epsilon}\right)^{2p}\geq\mathcal{O}\left(\left({4}/{\delta}\right)^{p}\right)$. After rearranging the terms for $\epsilon$, we obtain the bound on the compression budget $\epsilon\leq\mathcal{O}\left({\delta\sqrt{\delta}}\right)$. This implies the iteration complexity $T\geq\mathcal{O}\left({1}/{\delta^{2}}\right)$ and the model complexity $\mathcal{M}\geq\mathcal{O}\left({1}/{\delta^{p}}\right)$.

Supplementary Material for

Nonparametric Compositional Stochastic Optimization for Risk-Sensitive Kernel Learning

Appendix E Proofs of Propositions and Technical Lemmas

First, let’s quantify the projection-induced error in terms of the true stochastic quasi-gradient in the following proposition.

Proposition 1

Given independent identical realizations $(\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ of the two associated random variables $(\boldsymbol{\xi},\boldsymbol{\theta})$, for all $t$ it holds that

[TABLE]

where $\alpha>0$ denotes the step-size and $\epsilon>0$ is the compression parameter of Algorithm 2.

Proof: Consider the square-Hilbert-norm difference of $\tilde{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ and $\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ defined in (33) and (34), respectively,

[TABLE]

Multiply and divide $\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$, the last term, by $\alpha$, and reorder terms to write

[TABLE]

where we have taken the common scalar $\frac{1}{\alpha}$ out of the norm. The second equality in (E) holds due to the definition of $\tilde{f}_{t+1}$ and $f_{t+1}$ in (30) and (21) and substituting here. The last inequality in (E) follows from the stopping criteria used for the KOMP algorithm given by $\lVert\tilde{f}_{t+1}-f_{t+1}\rVert_{{\mathcal{H}}}\leq\epsilon$. $\hfill\blacksquare$

This result establishes that norm of the difference between the projected and unprojected stochastic quasi-gradient is the ratio of compression budget $\epsilon$ to the step-size $\alpha$.

Next we present an intermediate lemma which is vital to establishing an approximate stochastic descent relationship, and hence convergence.

Lemma 2

For the given algorithm history ${\mathcal{F}}_{t}\supset\sigma(\xi_{u},\theta_{u}\}_{u\leq t}$ at time $t$, under the Assumptions 1 - 4, consider the sequence of iterates $f_{t}$ generated by Algorithm 1. Then:

1. i.

The conditional expectation of the Hilbert-norm difference of the next function estimate $f_{t+1}$ and current iteration $f_{t}$ $\alpha>0$ satisfies

[TABLE] 2. ii.

The conditional expectation of the Hilbert-norm difference of the next function estimate $f_{t+1}$ and optimal function $f^{\star}$ satisfies that

[TABLE]

Proof of Lemma 2(i): At the current time instant $t$, consider the Hilbert-norm difference between the next function iterate $f_{t+1}$ and current estimate $f_{t}$ using the definition of $f_{t+1}$ in (35), i.e.,

[TABLE]

where we add and subtract the functional stochastic quasi-gradient $\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ on the first line of (88) and apply the inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$ which holds for any $a,b$. Now, we apply the result stated in Proposition 1 to the second term on the right hand side of (88). After performing this, taking the conditional expectation on the filtration ${\mathcal{F}}_{t}$ yields

[TABLE]

Next, utilize the definition of the stochastic quasi functional gradient provided in (33) and again using $(a+b)^{2}\leq 2(a^{2}+b^{2})$, we get

[TABLE]

Applying the Cauchy-Schwartz inequality to the first term on right hand side of (90) yields

[TABLE]

The second inequality in (E) is obtained from the Law of total expectation and using Assumption 1 [cf. (36)] which implies that $\|{\kappa}(\boldsymbol{\xi}_{t},\cdot)\|^{2}_{{\mathcal{H}}}\leq U^{2}$. Next, due to the set projection in (III-B) and Assumption 1, we can conclude that $f_{t}$ has a bounded Hilbert norm, i.e., there exists some $0<K<\infty$ such that $\left\|f_{t}\right\|\leq K$ for all $t$. Using the upper bound of $K$ and the upper bounds as per Assumption 2 [cf. (37)], we get the final equation

[TABLE]

where we define $\sigma_{f}^{2}:=4U^{2}(G_{\mathscr{h}}^{2}G_{\ell}^{2}+\lambda^{2}K^{2})$.

Proof of Lemma 2(ii): This proof is a generalization of Lemma 3 in Appendix G.2 in the Supplementary Material of [3] to a function-valued stochastic quasi-gradient step combined with bias induced by the sparse subspace projections ${\mathcal{P}}_{{\mathcal{H}}_{{\mathbf{D}}_{t+1}}}[\cdot]$ in (III-B). Let us consider the square-Hilbert norm distance of $f_{t+1}$ from the optimal $f^{\star}$, i.e.,

[TABLE]

where we utilized the reformulation of the function update defined in (35) for the first equality, and expand the square in the second. Now, adding and subtracting $\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ (which is the (un-projected) functional stochastic quasi-gradient (33)) to first term in the inner product on right hand side of (E) yields

[TABLE]

Applying the Cauchy-Schwartz inequality to the third inner product term on the right-hand side of (E) and then utilizing the upper bound developed in Proposition 1, we get

[TABLE]

Now, as defined in Lemma 2, let ${\boldsymbol{\delta}_{t}:=\boldsymbol{\mathscr{h}}_{\boldsymbol{\xi}_{t}}(f(\boldsymbol{\xi}_{t}))}$ with \bar{\boldsymbol{\delta}}_{t}=\mathbb{E}\left[\boldsymbol{\delta}_{t}{\,\big{|}\,}\boldsymbol{\theta}_{t}\right]. Add and subtract $\hat{\nabla}_{f}R(f_{t},\bar{\boldsymbol{\delta}}_{t};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$, which is nothing but the stochastic quasi-gradient evaluated at $(f_{t},\bar{\boldsymbol{\delta}}_{t})$ rather than $(f_{t},{\mathbf{g}}_{t+1})$, inside the inner-product term on the right-hand side of (95), to obtain

[TABLE]

where we substitute in the definitions of $\hat{\nabla}_{f}R(f_{t},\bar{\boldsymbol{\delta}}_{t};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ and $\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})$ [cf. (33), (33), respectively] in (E), and cancel out the common regularization term $\lambda f_{t}$. From the derivation in (88) to (E), we can conclude that

[TABLE]

where we define $\sigma_{f}^{2}:=4U^{2}(G_{\mathscr{h}}^{2}G_{\ell}^{2}+\lambda^{2}K^{2}$. To proceed further, let us define the directional error term related to stochastic quasi-gradient and the stochastic gradient as

[TABLE]

From here, compute the conditional expectation on the algorithm history ${\mathcal{F}}_{t}$:

[TABLE]

Utilizing the fact that the compositional objective $R(f)$ defined in (1) is convex with respect to $f$ and utilizing the first order convexity condition, we have

[TABLE]

Using the inequality in (100) and the upper bound in (97) into (E), we get

[TABLE]

It remains to analyze $v_{t}$, the directional error associated with using stochastic quasi-gradients rather than stochastic gradients. Proceed by applying the Cauchy-Schwartz inequality to (98), which allows us to write

[TABLE]

where the second inequality in (102) uses Assumptions 1 [cf. (36)] and 3 [cf. (38)]. Consider Peter-Paul’s inequality $2ab\leq\rho a^{2}+b^{2}/\rho$ for $\rho,a,b>0$, which we apply to (102) with $a=\|\bar{\boldsymbol{\delta}}_{t}-{\mathbf{g}}_{t+1}\|$, $b=\alpha(\|{{\boldsymbol{\mathscr{h}}^{\prime}_{\boldsymbol{\xi}_{t}}(f_{t}(\boldsymbol{\xi}_{t}))}})\|f_{t}-f^{\star}\|_{{\mathcal{H}}}$, and $\rho=\beta_{t}$ so that (102) becomes

[TABLE]

The conditional mean of $v_{t}$ [cf. (98)], using (103), is then

[TABLE]

where we have used the Assumption 1 [cf. (36)]. Now, substitute (104) into the right-hand side of (E) and gather like terms:

[TABLE]

which is as stated in Lemma 2(ii). $\blacksquare$

Lemma 3

Consider the case that $R(f)$ [cf. (2)] is Lipschitz smooth with parameter $L_{R}$. Then under Assumptions 1 - 4, and the condition that the step-size ratio satisfies $\frac{\alpha}{\beta}L_{\ell}^{2}G_{h}^{2}U^{2}\leq\frac{1}{4}$, i.e., $\alpha$ is sufficiently small relative to $\beta$, then the sequence of functions $f_{t}$ generated by the proposed algorithm satisfies

[TABLE]

Proof: For a Lipschitz smooth $R$ with parameter $L_{R}$, Taylor’s expansion implies the quadratic upper-bound:

[TABLE]

where the equality in (106) holds from the definition of the projected gradient iterate (35). Next, we establish the upper bound on the inner product term on the right hand side of (106) as follows

[TABLE]

where we add subtract the term $-\alpha\boldsymbol{\langle}\nabla_{f}R(f_{t}),\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})\boldsymbol{\rangle}_{\mathcal{H}}$. After rearranging the terms, we get

[TABLE]

where we add subtract the term $\alpha\boldsymbol{\langle}\nabla_{f}R(f_{t}),\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})\boldsymbol{\rangle}$. From the inequality $ab\leq(a^{2}+b^{2})/2$, we obtain

[TABLE]

where the inequality in (109) holds from Proposition 1. Next, consider the term first term on the right hand side of (109) and add subtract $\nabla_{f}R(f_{t})$ as follows

[TABLE]

Taking the conditional expectation on the both sides and consider the term $\alpha\mathbb{E}[\boldsymbol{\langle}\nabla_{f}R(f_{t}),\nabla_{f}R(f_{t})-\hat{\nabla}_{f}R(f_{t},{\mathbf{g}}_{t+1};\boldsymbol{\xi}_{t},\boldsymbol{\theta}_{t})\boldsymbol{\rangle}_{\mathcal{H}}\mid\mathcal{F}_{t}]$

[TABLE]

where the equality holds from the definition of the gradient in (33). Using the Cauchy Schwartz inequality and the upper bounds of Assumption 1, we get

[TABLE]

Next, applying the Lipschitz smoothness of the outer function as stated in Assumption 3, we get

[TABLE]

Using the inequality $ab\leq\rho a^{2}+\frac{1}{\rho}b^{2}$ for $0<\rho<1$ and setting $\rho=\beta$, $a=\|\bar{\boldsymbol{\delta}}_{t}-\mathbf{g}_{t+1}\|$, and $b=\alpha UL_{\ell}\|\nabla_{f}R(f_{t})\|_{\mathcal{H}}\|\mathbf{h}^{\prime}_{\xi_{t}}(f(\xi_{t}))\|$, we obtain

[TABLE]

where the second inequality in (113) holds from Assumption 2. Substituting the result of (113) into (110) after taking conditional expectation, we can write

[TABLE]

Similarly, take the conditional expectation of (109) and utilizing the result of (114), we get

[TABLE]

Finally, substituting (115) into (106) and utilizing Assumption 2, we obtain the statement of Lemma 3 as follows

[TABLE]

which is as stated in Lemma 3. $\blacksquare$

Appendix F Proof of Theorem 1

Before discussing the finiteness of the model order, we discuss a lemma which helps us to relate the stopping criterion of specification in KOMP to the Hilbert subspace.

Lemma 4

Let us define the distance of an arbitrary random feature vector $(\boldsymbol{\xi})$ calculated as $\phi(\boldsymbol{\xi})=\kappa(\boldsymbol{\xi},\cdot)$ to, ${\mathcal{H}}_{{\mathbf{D}}}=\text{span}\{\kappa({\mathbf{d}}_{n},\cdot)\}_{n=1}^{M}$, the subspace of the Hilbert space spanned by a dictionary ${\mathbf{D}}$ of size $M$, as

[TABLE]

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

[TABLE]

Proof:

We can write the distance to a Hilbert space ${\mathcal{H}}_{{\mathbf{D}}}$ as follows

[TABLE]

which is obtained since the dictionary $D$ is fixed and the only free parameter left to optimize is ${\mathbf{v}}$. Now similarly to (LABEL:eq:proximal_hilbert_dictionary_polker) - (32), we can obtain an optimal weight vector ${\tilde{\mathbf{v}}}^{*}={\mathbf{K}}_{{\mathbf{D}}_{t},{\mathbf{D}}_{t}}^{-1}\boldsymbol{\kappa}_{{\mathbf{D}}_{t}}(\boldsymbol{\theta},\boldsymbol{\xi})$ and then substitute back into (F), we get

[TABLE]

∎

F-A Proof of Theorem 1

Proof: Consider two arbitrary time instants $t$ and $t+1$, at which $f_{t}$ and $f_{t+1}$ are the iterates generated by Algorithm 1 with corresponding model order denoted by $M_{t}$ and $M_{t+1}$, respectively. We consider a constant step size algorithm with fixed approximation budget $\epsilon=C\alpha^{2}$ for some constant $K>0$. For the sake of analysis, suppose that the model order at $t+1$ is $M_{t+1}\leq M_{t}$ which means that model order does not grow as we go to iterate $t+1$ from $t$. Note that the model order remains the same from $t$ to $t+1$ if the error introduced by the removal of recently appended data points ($\boldsymbol{\xi}_{t}$) to dictionary ${\tilde{\mathbf{D}}}_{t+1}=[{\mathbf{D}}_{t};\boldsymbol{\xi}_{t}]$ [cf. (31)] is less than the approximation budget $\epsilon$. In other words, the model order does not grow if the stopping criteria of KOMP (Algorithm 2), stated as $\min_{j=1,\dots,{M_{t}+1}}\gamma_{j}>\epsilon$, is not satisfied. This leads us to the conclusion that the model order remains the same from $t$ to $t+1$ if

[TABLE]

Further observe that the left hand side of (122) is a lower bound for the approximation error $\gamma_{M_{t}+1}$ at $t+1$ because of the minimization over all $j=1\cdots(M_{t}+1)$. This states that the error $\gamma_{M_{t}+1}$ introduced by removing the recently appended pair $(\boldsymbol{\xi}_{t})$ is such that $\min_{j=1,\dots,{M_{t}+1}}\gamma_{j}\leq\gamma_{M_{t}+1}$. Therefore, the model order does not grow if $\gamma_{M_{t}+1}\leq\epsilon$ holds because then (122) is satisfied. Let use analyze the error term $\gamma_{M_{t}+1}$ as follows. The definition of $\gamma_{M_{t}+1}$ with the substitution of $\tilde{f}_{t+1}$ in (30) with the notation

[TABLE]

allows us to write

[TABLE]

where we denote the $k^{\text{th}}$ data column of ${\mathbf{D}}_{t}$ as ${\mathbf{d}}_{k}$. To obtain the value ${\tilde{\mathbf{u}}}^{*}$ which minimizes the above expression, we calculate the gradient of above expression with respect to ${\mathbf{u}}$ and set it to zero. Following the similar logic to that of (LABEL:eq:proximal_hilbert_dictionary_polker) - (32), we get the following result

[TABLE]

Next, utilizing the optimal value ${\tilde{\mathbf{u}}}^{*}$ from (124) into the expression in (123) along with the short-hand notation ${f}_{t}(\cdot)={\mathbf{w}}^{T}\boldsymbol{\kappa}_{{\mathbf{D}}_{t}}(\cdot)$ and $\sum_{k}u_{k}\kappa({\mathbf{d}}_{k},\cdot)={\mathbf{u}}^{T}\boldsymbol{\kappa}_{{\mathbf{D}}_{t}}(\cdot)$, we get

[TABLE]

Further simplifying the above expression by cancelling the similar terms $(1-\alpha\lambda){\mathbf{w}}^{T}\!\boldsymbol{\kappa}_{{\mathbf{D}}_{t}}(\cdot)$ and taking taking the common term $\alpha|V_{t}|$ outside the norm as

[TABLE]

It is remarked that the norm expression in the right hand side of (126) describes the distance to the subspace ${\mathcal{H}}_{{\mathbf{D}}_{t}}$ as described in (119) and defined in Lemma 4 with a scaling factor of $\alpha|V_{t}|$. The right hand side of (126) can be written as

[TABLE]

where the result in (119) on the right hand side of (F-A) to replace the Hilbert-norm term. Observe that when the stopping criteria of KOMP is violated, (122) holds and thus $\gamma_{M_{t}+1}\leq\epsilon$. Therefore, we have that the right-hand side of (F-A) will be upper-bounded by $\epsilon$, and we can write the following inequality

[TABLE]

The error associated with the model order $M_{t+1}$ is denoted by $\gamma_{M_{t+1}}$. Observe that if (128) holds, then $\gamma_{M_{t+1}}\leq\epsilon$ holds, but since $\gamma_{M_{t+1}}\geq\min_{j}\gamma_{j}$, we may conclude that (122) is satisfied. Consequently the model order at the subsequent step does not grow $M_{t+1}\leq M_{t}$ whenever (128) is valid.

Now, consider the contrapositive of the preceding expressions. Observe that the model order growth condition ($M_{t+1}=M_{t}+1$) implies that

[TABLE]

holds. This condition establishes the fact that every time a new data ($\boldsymbol{\xi}$) is appended to kernel dictionary, then the associated product kernel is guaranteed to be at least a distance of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\epsilon}{\alpha|V_{t}|}} from every other kernel function in the current model.

Now utilizing the Cauchy Schwartz inequality and Assumption 5 we get

[TABLE]

This upper bound implies that $1/|V_{t}|\geq 1/(C_{\ell}L_{\mathscr{h}})$, therefore we can lower bound the right hand side in (129) as follows

[TABLE]

From (129), we obtain

[TABLE]

Therefore, the KOMP stopping criterion is violated for the newest point whenever distinct dictionary points ${\mathbf{d}}_{k}$ and ${\mathbf{d}}_{j}$ for $j,k\in\{1,\dots,M_{t}\}$, satisfy the condition \|\phi({\mathbf{d}}_{j})-\phi({\mathbf{d}}_{k})\|_{2}>{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}}. Next, we proceed in a similar manner to that of Theorem 3.1 in [61]. Note that since the space $\mathcal{U}$ is compact and $\kappa$ is continuous, the range $\phi({\mathcal{U}})$ (where $\phi({\mathbf{u}})=\kappa({\mathbf{u}},\cdot)$ for ${\mathbf{u}}\in{\mathcal{U}}$) of the kernel transformation of feature space ${\mathcal{U}}$ is compact. This allows us to conclude that the number of balls of radius $\chi$ (here, \chi={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}}) required to completely cover the set $\phi({\mathcal{U}})$ is finite (see, e.g., [62]).

To prove the main result of (40), we consider the result in [61, Proposition 2.2] which states that for a Lipschitz continuous Mercer kernel $\kappa$ on compact set $\mathcal{X}\subseteq\mathbb{R}^{p}$, for any training set $\{{\mathbf{x}}_{t}\}_{t=1}^{\infty}$ and any $\nu>0$, the number of elements in the dictionary is upper bounded as

[TABLE]

where $Y$ is a constant depends upon $\mathcal{X}$ and the kernel function. From the result in (131), we conclude that $\sqrt{\nu}={\frac{\epsilon}{\alpha C_{\ell}L_{\mathscr{h}}}}$, which we may substitute into (132) to obtain

[TABLE]

as stated in (40). Note that in (133), we have defined a constant $Y^{\prime}=YC_{\ell}^{2p}L_{\mathscr{h}}^{2p}$. We remark that there is a trivial lower bound of $M_{t}\geq 1$. ∎