# Sparse Representations of Positive Functions via
First and Second-Order Pseudo-Mirror Descent ††thanks: Abhishek Chakraborty is
currently affiliated to NetApp, India. This work was done as part of his
M.Tech thesis at IIT Kanpur, email<EMAIL_ADDRESS>Ketan
Rajawat is with Dept. of EE, IIT Kanpur, India, email<EMAIL_ADDRESS>Alec
Koppel is with Amazon, Supply Chain Optimization Technologies, Bellevue, WA,
USA Email<EMAIL_ADDRESS>A preliminary version of this work appeared as
[1], but under more restrictive assumptions of convexity, and without any
development of second-order methods.

Abhishek Chakraborty, Ketan Rajawat, and Alec Koppel

###### Abstract

We consider expected risk minimization problems when the range of the
estimator is required to be nonnegative, motivated by the settings of maximum
likelihood estimation (MLE) and trajectory optimization. To facilitate
nonlinear interpolation, we hypothesize that the search space is a Reproducing
Kernel Hilbert Space (RKHS). We develop first and second-order variants of
stochastic mirror descent employing (i) _pseudo-gradients_ and (ii)
complexity-reducing projections. Compressive projection in the first-order
scheme is executed via kernel orthogonal matching pursuit (KOMP), which
overcomes the fact that the vanilla RKHS parameterization grows unbounded with
the iteration index in the stochastic setting. Moreover, pseudo-gradients are
needed when gradient estimates for cost are only computable up to some
numerical error, which arise in, e.g., integral approximations. Under constant
step-size and compression budget, we establish tradeoffs between the radius of
convergence of the expected sub-optimality and the projection budget
parameter, as well as non-asymptotic bounds on the model complexity. To refine
the solution’s precision, we develop a second-order extension which employs
recursively averaged pseudo-gradient outer-products to approximate the Hessian
inverse, whose convergence in mean is established under an additional
eigenvalue decay condition on the Hessian of the optimal RKHS element, which
is unique to this work. Experiments demonstrate favorable performance on
inhomogeneous Poisson Process intensity estimation in practice.

##  I Introduction

Function interpolation over nonnegative range arises in numerous settings:
trajectory optimization [2], as well as maximum likelihood estimation [3]. In
this work, we focus on nonnegative function estimation problems defined by an
expected value minimization of a non-convex cost satisfying Polyak–Łojasiewicz
condition. Further, the feasible set is a Reproducing Kernel Hilbert Space
(RKHS) [4], motivated by their nonlinear interpolation attributes. We develop
first and second-order functional variants of pseudo-mirror descent to solve
this problem that, different from prior work, explicitly trade the dynamic
memory or representational complexity of the RKHS estimates with their
accuracy. Our main theoretical results establish trade-offs between sublinear
convergence in expectation and the function complexity. The driving example
throughout is to obtain accurate and efficient representations of the
intensity of an inhomogeneous Poisson Process [3] from incrementally revealed
samples.

To contextualize our approach, let us set aside the issue of range
nonnegativity and the function parameterization for the moment. Finding a
linear statistical model may be formulated as an expected value minimization
over Euclidean space, which, when analytical expressions for the optimizer are
unavailable (i.e., objective is not a quadtratic), one can employ Newton or
gradient method to obtain the globally optimal solution[5]. But, doing so
hinges upon gradient information being efficiently computable, which fails for
expected value objectives over an unknown distribution [6]. For such settings,
stochastic approximations are necessary [7], which use stochastic gradients
(or stochastic Quasi–Newton updates) in lieu of exact updates. Their
performance depends upon the properties of the unknown data distribution, and
one may guarantee their performance only probabilistically [8]. In this work,
we build upon a stochastic variant of proximal gradient method [9], i.e.,
stochastic mirror descent, which substitutes Euclidean distance by a Bregman
divergence in the reformulation of gradient iteration as a local quadratic
minimization. Doing so yields convergence benefits when the feasible set is
structured, e.g., the probability simplex [10]. We further develop
Quasi–Newton iterations based upon tracking a matrix of recursively
constructed gradient outer-products.

To contextualize the RKHS specification, we note that the strong guarantees
for optimization over Euclidean space is belied by the fact that linear models
are outperformed by universal function approximators [11]. Specifically, deep
neural networks (DNNs) [12] and kernel methods [13] have set the standard
across computer vision [14] and natural language processing [15], among other
fields. We focus on RKHS parameterizations since (i) under a suitable choice
of kernel, they, in fact, can be equivalent to certain DNNs [16]; (ii) their
training may explicitly trade-off the accuracy of an iterative update with a
function’s representational complexity [17]; and (iii) one may directly co-
design the kernel and Bregman divergence for likelihood objectives [18].

One challenge associated with the RKHS setting is that the function’s
representational complexity, owing to the Representer Theorem [19], is
proportional to the sample size, which for expected risk minimization,
approaches infinity. A myriad of approaches to approximating RKHS functions
exist: those based upon matrix factorization [20], spectral properties of the
kernel matrix [21], and random feature techniques [22]. Random features may be
combined with random sampling [23] or subspace projections [24] based upon
greedy compression [25] to develop approximated online RKHS updates. Among
kernel approximations, subspace projection methods that fix the per-update
error rather than the complexity have been shown to obtain the strongest
performance in theory and practice [26, 27], and hence we adopt them here to
sparsify stochastic descent directions in RKHS.

Two intertwined challenges emerge when learning nonnegative functions:
enforcing the range constraint and evaluation of the stochastic gradient. Upon
suitable initialization, one may ensure updates remain within the space of
positive functions via specifying Bregman divergence as the KL divergence
[28]. However, the gradient may only be partially evaluable, as in the
numerical evaluation of an integral. This issue arises in point process
likelihoods [29], as well as the case where the true kernel is a latent
convolution of one’s chosen kernels [30]. In these cases, _pseudo-gradients_
[31], i.e., stochastic descent directions positively correlated with the true
gradient, must be used. Doing so in online pseudo-mirror descent has recently
achieved state of the art performance in intensity estimation in point
processes [29]. In this work, we improve upon this approach by designing
projected variants that yield tunable trade-offs between convergence accuracy
and complexity, and in doing so, set a new standard for performance in
practice among first-order schemes.

To refine the numerical precision in practice, we also develop a Quasi–Newton
iteration in mirror space that fixes the RKHS subspace dimension, which
permits using any infinite-dimensional kernel such as a Gaussian, in contrast
to [32, 33]. This iteration constructs Hessian approximations using
recursively averaged pseudo-gradient outer-products. Overall, then, our
technical contributions to:

  1. (i)

formulate the expected risk minimization problem over RKHS functions with
nonnegative range (Sec. II);

  2. (ii)

develop a projected pseudo-mirror descent (Algorithms 1) that uses pseudo-
gradients in the RKHS setting. To maintain a finite model order, the first-
order scheme employs a compressive projection over the dual iterates defined
by the Bregman divergence conjugate (Sec. III);

  3. (iii)

develop a stochastic Quasi-Newton variant as Algorithm 3 that employs employs
a fixed-subspace projection in RKHS. Together, we obtain specific parametric
updates in terms of weights and past training examples; Additionally, we fuse
the merits of the dynamic and fixed-subspace iterations via a hybrid first and
second-order scheme (Algorithm 4 in Sec. IV);

  4. (iv)

establish that the resulting sub-optimality of the first-order scheme
converges linearly in expectation up to a bounded error neighborhood under
constant step-size (Theorem 1), and provide a non-asymptotic complexity
characterization of the learned function’s parameterization (Sec. V). Further,
we establish the convergence of the second-order scheme in Theorem 3 up to a
persistent bias that arises from function approximation error;

  5. (v)

evaluate numerically (Sec. VI) the proposed approach for fitting normalized
intensity of an inhomogeneous Poisson Process on synthetic 1-D Gaussian toy
dataset, real world NBA Stephen Curry dataset [29] and 2-D Chicago crime
dataset [34]. In all instances, we observe a favorable tradeoff relative to
offline and online benchmarks.

In addition, we show the similarity of functional dual-averaging to mirror
descent in Appendix D and extensions of our first order algorithm for
supervised learning problems and experimentation for multi-class kernel
logistic regression on MNIST dataset [35] in Appendices E and F respectively.
We also provide an almost sure convergence proof for first order algorithms
when we restrict the selection of pseudo-gradients to stochastic gradients
(Appendix G).

Additional Context. Most similar stochastic first-order methods in RKHS [26]
mandate the Bregman divergence of mirror descent as the RKHS norm and hence
violate non-negativity. Similarly [29] is an online algorithm that preserves
positivity but requires the dictionary defining the functional estimate to be
static at the outset of training, which incurs bias that is difficult to
quantify. By contrast, our dynamic memory approach leads to a flexible
representation whose error is directly tunable in theory and practice. Offline
reformulations of optimization over RKHS with the non-negative range also have
been encapsulated as a semidefinite program in [36, 37, 38]. Such approaches
are not easily extendable into an incremental setting, as each new sample
introduces an additional constraint.

Quasi–Newton iterations in RKHS, rather than employing classical Broyden
updates [39][Ch. 6], extend gradient outer-product constructions of Hessian
approximations for vector-valued iterates, as in [32, 33]. These approaches
fix RKHS subspace sizes to ensure gradient outer products belong to a static
matrix subspace. In these works, regret bounds are developed in terms of the
RKHS weight sequence, which fail to quantify the function approximation error.
By contrast, our convergence results explicitly establish that the memory
reduction employed to make Hessian estimates computable manifests as a
persistent bias in the algorithm’s radius of convergence, which we quantify in
terms of the decay rate of the eigenvalues of the Hessian of the optimal
objective.

##  II Problem Formulation

Notation. We denote functions by small scalar alphabets such as $f$, $z$ and
$g$. The Hilbert space is denoted by ${\mathcal{H}}$ and its dual space by
${\mathcal{H}}^{\ast}$. All norms applied to elements of ${\mathcal{H}}$ are
Hilbert norms, henceforth denoted simply by $\left\|\cdot\right\|$. The norm
applied on the elements of the ${\mathcal{H}}^{\ast}$ are represented as
$\left\|\cdot\right\|_{\ast}$.

Consider the problem of expected risk minimization in the online setting:
independent identically distributed (i.i.d.) training samples
$\\{{\mathbf{x}}_{t}\\}_{t\geq 1}$ are observed in a sequential manner. Here,
${\mathbf{x}}_{t}\in{\mathcal{X}}\subset{\mathbb{R}}^{d}$ denote the instances
or input data points from an unknown distribution $\mathbb{P}({\mathbf{x}})$.
Based upon this input stream, we seek to fit a (possibly unnormalized) density
$f$ that belongs to a hypothesized function class ${\mathcal{H}}$ according to
its ability to minimize a loss function $\ell(f({\mathbf{x}}))$. This loss is
associated with the negative log-likelihood of a probabilistic model. We are
interested in minimizing its expectation over the unknown data distribution
$\mathbb{P}({\mathbf{x}})$, i.e.,

$\displaystyle
f^{\star}=\operatorname*{arg\,min}_{f\in{\mathcal{H}}_{+}}R(f):=\mathbb{E}[\ell(f({\mathbf{x}}))]\;.$
(1)

Throughout, we occasionally abbreviate the instantaneous cost as
$r_{t}(f):=\ell(f({\mathbf{x}}_{t}))$. We focus on the case that the range of
$f$ is nonnegative, which arises inherently in the estimation of probabilistic
models. An important point to underscore is the fact that the optimization of
function $f$ is over the whole space ${\mathcal{H}}$, and the positivity of
$f$ is not enforced explicitly by hard constraints. Instead, by proper
selection of positivity-preserving Bregman divergence (e.g., the KL
divergence), one may implicitly enforce that constraint, as detailed in Sec.
III.

Next, we specify the function class ${\mathcal{H}}\supset{\mathcal{H}}_{+}$ to
which estimators $f$ belong. Specifically, we hypothesize that ${\mathcal{H}}$
is a Hilbert space associated with symmetric positive definite kernel $\kappa$
that satisfies
${\mathcal{H}}:=\overline{\text{span}(\kappa({\mathbf{x}},\cdot))}$, where
every $f\in{\mathcal{H}}$ can be written as linear combination of the kernel
evaluations. The elements of ${\mathcal{H}}$ also satisfy the reproducing
property which states that $\langle
f,\kappa({\mathbf{x}},\cdot)\rangle=f({\mathbf{x}})$ for all
${\mathbf{x}}\in{\mathcal{X}}$. Hilbert spaces of this type are called RKHS
[4]. Examples range from simple kernel functions such as the Gaussian
$\kappa({\mathbf{x}},{\mathbf{x}}^{\prime})=\exp(-\left\|{\mathbf{x}}-{\mathbf{x}}^{\prime}\right\|^{2}/2\tilde{c})$
and polynomial
$\kappa({\mathbf{x}},{\mathbf{x}}^{\prime})=({\mathbf{x}}^{\top}{\mathbf{x}}+\tilde{b})^{\tilde{c}}$,
to data-dependent convolutional kernels [30].

With the objective and the search space clarified, we present a representative
of Poisson intensity estimation, which is our driving application throughout.

######  Example 1.

Poisson Point Process intensity estimation:  Poisson Processes are a family of
probabilistic models used to count the number of events $N(\mathcal{T})$
within an interval $\mathcal{T}\subset\mathcal{S}\subset\mathbb{R}^{d}$, where
$N(\cdot)$ is a stochastic process. A fundamental question that arises in its
use is the intensity parameter $\lambda(\cdot)$, which determines the rate
$\lambda(s)$ at which new events occur in an infinitesimally small time-
increment, i.e., $\lambda(s)=\lim_{\Delta s\rightarrow 0}\mathbb{E}[N(\Delta
s)]/(\Delta s)$. In inhomogeneous cases, this parameter is a nonlinear
function, such that the likelihood associated with Poisson points
$\\{\mathbf{t}_{n}\\}_{n=1}^{N}$ takes the form:

$\displaystyle
L(f)=\prod_{n=1}^{N}\lambda(\mathbf{t}_{n})\exp\left\\{-\int_{\mathcal{S}}\lambda(\mathbf{t})d\mathbf{t}\right\\}.$
(2)

Then, one may construct an instantiation of (1) by considering the negative
log-likelihood of (2), inspired by [34]:

$\displaystyle
R(f)=-\sum_{n=1}^{N}\log(\lambda(\mathbf{t}_{n}))+\int_{{\mathcal{S}}}\lambda(\mathbf{t})d\mathbf{t},$
(3)

where one may identify that the Poisson points $\mathbf{t}_{n}$ play the role
of ${\mathbf{x}}_{n}$, and $\lambda(\cdot)$ is the unknown function $f(\cdot)$
we seek to estimate which is required to be nonnegative. This instance of
unsupervised learning has been studied both when
$\\{\mathbf{t}_{n}\\}_{n=1}^{N}$ are available all at once or revealed
incrementally (with $N\rightarrow\infty$), respectively, in [34] and [29]. We
focus on online approaches to (3) – see [3] for further details.

In this work, we focus on algorithms to solve (1) by designing search
directions that allow the function to move in the interior of the space of
nonnegative ranged functions. Next, we present a clarifying remark about the
complexity implications of searching over RKHS.

######  Remark 1.

(Empirical Risk Minimization) An important special case of (1) is _empirical
risk minimization_ (ERM) where a fixed collection of data
${\mathcal{D}}:=\\{({\mathbf{x}}_{n})_{n=1}^{N}\\}$ is available, and we seek
to find the empirical predictor

$\displaystyle\hat{f}_{N}=\arg\min_{f\in{\mathcal{H}}}\frac{1}{N}\sum_{n=1}^{N}r_{n}(f).$
(4)

Observe that $\hat{f}_{N}$ in (4) belongs to the Hilbert space, and hence is
infinite-dimensional. However, the Representer Theorem [40] establishes that
$\hat{f}_{N}$ takes the form:

$\displaystyle\hat{f}_{N}(\cdot)=\sum_{n=1}^{N}w_{n}\kappa({\mathbf{x}}_{n},\cdot),$
(5)

where $\\{w_{m}\\}_{m=1}^{N}$ are some real-valued weights. Substituting (5)
into (4) reduces searching over ${\mathcal{H}}$ to real-valued space
$\mathbb{R}^{N}$:

$\displaystyle\hat{{\mathbf{w}}}_{N}$
$\displaystyle=\arg\min_{{\mathbf{w}}\in\mathbb{R}^{N}}\frac{1}{N}\sum_{n=1}^{N}r_{n}({\mathbf{w}}^{\top}{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{n})),$
(6)

where we have collected kernel evaluations
$\\{\kappa({\mathbf{x}}_{n},{\mathbf{x}}_{m})\\}_{n}$ into a vector called the
empirical kernel map
${\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{n})\in{\mathbb{R}}^{N}$ and
$[\kappa({\mathbf{x}}_{n},{\mathbf{x}}_{m})]_{n,m=1}^{N,N}$ into the Gram, or
kernel, matrix
${\mathbf{K}}_{{\mathcal{D}}{\mathcal{D}}}\in{\mathbb{R}}^{N\times N}$.
Observe that as the sample size $N\rightarrow\infty$, due to the curse of
kernelization, one must balance solving (1) to optimality with ensuring that
the memory remains bounded. In this work, we focus on (1) with the
understanding that (4) is a special case when the i.i.d. data
$\\{{\mathbf{x}}_{n}\\}_{n=1}^{N}$ are sampled uniformly at random from a
stationary distribution.

Hereafter, we introduce some technical machinery needed for our algorithmic
development in the following section.

###  II-A Bregman Divergence and Pseudo-Gradients

Bregman Divergence. We proceed by presenting technicalities pertaining to the
mirror map and the Bregman divergence before defining the proposed algorithm.
Let $\psi:{\mathcal{H}}\rightarrow{\mathbb{R}}$ be a proper, closed, smooth,
and strongly convex functional with respect to the RKHS norm
$\left\lVert\cdot\right\rVert$111We denote the RKHS norm
$\|f\|^{2}=\|f\|_{\mathcal{H}}^{2}$ of $f\in\mathcal{H}$ as
$\|f\|_{\mathcal{H}}^{2}=\langle
f,f\rangle=\sum_{i,j}w_{i},w_{j}\kappa(\mathbf{x}_{i},\mathbf{x}_{j})$, with
$f$ admitting the weighted expansion
$f=\sum_{i}w_{i}\kappa(\mathbf{x}_{i},\cdot)$ via the Representer Theorem,
with $w_{i}\in\mathbb{R}$. . The Fenchel conjugate of $\psi$ is denoted as
$\psi^{\ast}:{\mathcal{H}}^{\ast}\rightarrow{\mathbb{R}}$, where
${\mathcal{H}}^{\ast}$ is the dual space of ${\mathcal{H}}$. Define the
shorthand for the objective evaluated at the gradient of the dual
$R_{\psi}(z)=(R\circ\nabla\psi^{\ast})(z)=R(\nabla\psi^{\ast}(z))$ where
$z\in{\mathcal{H}}^{\ast}$. This composition allows one to write $\nabla
R_{\psi}(\nabla\psi(f))=\nabla R(f)$, $f\in{\mathcal{H}}$ since
$\nabla\psi^{\ast}=(\nabla\psi)^{-1}$. This identity will be employed later in
the convergence analysis. The purpose of defining functional $\psi$ is that it
induces a distance-like Bregman divergence
$B_{\psi}:{\mathcal{H}}\times{\mathcal{H}}\Rightarrow{\mathbb{R}}$ [41]:

$\displaystyle
B_{\psi}(f,{\tilde{f}}):=\psi(f)-\psi({\tilde{f}})-\langle\nabla\psi({\tilde{f}}),f-{\tilde{f}}\rangle_{{\mathcal{H}}}.$
(7)

The functional Bregman divergence in (7) behaves mostly like its vector-valued
counterpart, i.e., it satisfies non-negativity, strong-convexity in the first
argument, as well as a generalized Pythagorean theorem. We present a few
examples next.

  1. i)

KL-divergence or I-divergence: Let ${\mathcal{H}}$ be the space of probability
density functions, constructed from a radial kernel such as Gaussian,
Student-t, or Laplacian, such that $f$ takes positive values and integrates to
$1$. Defining $\psi(f)=\langle f,\log(f)-1\rangle_{{\mathcal{H}}}$ similar to
[29] yields $B_{\psi}(f,\tilde{f})=\langle
f,\log(f/{\tilde{f}})\rangle_{{\mathcal{H}}}$ if $f$ is absolutely continuous
with respect to ${\tilde{f}}$. For this case, $B_{\psi}$ is the KL-divergence
or I-divergence associated with the convex map $\psi(f)$ [28].

  2. ii)

Squared RKHS norm difference: The choice
$\psi(f)=\frac{1}{2}\left\|f\right\|^{2}$ leads to
$B_{\psi}(f,{\tilde{f}})=\frac{1}{2}\left\|f-{\tilde{f}}\right\|^{2}$.

  3. iii)

Squared Mahalanobis difference: Let ${\mathcal{L}}$ be a compact self-adjoint
operator, then the choice
$\psi(f)=\frac{1}{2}\langle{\mathcal{L}}f,f\rangle_{{\mathcal{H}}}$ results in
the Bregman divergence
$B_{\psi}(f,{\tilde{f}}))=\frac{1}{2}\langle{\mathcal{L}}(f-{\tilde{f}}),f-{\tilde{f}}\rangle_{{\mathcal{H}}}^{2}$,
which generalizes the squared difference. For the simple case when
${\mathcal{L}}$ is a diagonal operator, it amounts to pointwise multiplication
with a weight function.

Pseudo-Gradient. With Bregman divergence clarified, we shift to defining
pseudo-gradients, which are those directions $g_{t}$ with nonnegative expected
(unnormalized) cosine similarity with the gradient $\nabla_{f}R(f)$:

$\displaystyle\langle\nabla
R(f),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle\geq 0,$ (8)

where $g_{t}$ is the pseudo gradient associated with sample ${\mathbf{x}}_{t}$
and ${\mathcal{F}}_{t}$ denotes the past sigma algebra which contains all the
past data points one iteration back, i.e.
$\\{{\mathbf{x}}_{i}\\}_{i=1}^{t-1}$. Equivalently, a pseudo-gradient $g_{t}$
is any search direction that forms an acute angle with the original gradient
$\nabla R(f)$ in the dual space, as stated in [31], and may be used in lieu of
the true gradient when its evaluation is costly or intractable [29]. Next we
present a few instances to build intuition.

  1. i)

Stochastic Gradients: If $g_{t}$ is the stochastic gradient of $R(\cdot)$ at
$f_{t}$, then ${\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]=\nabla R(f_{t})$. So the
inner product $\langle\nabla
R(f_{t}),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle=\left\|\nabla
R(f_{t})\right\|^{2}_{\ast}\geq 0$.

  2. ii)

Kernel embedding:
$\kappa:{\mathcal{X}}\times{\mathcal{X}}\rightarrow{\mathbb{R}}_{+}$ is a
symmetric and positive definite kernel. So $\langle
G,\langle\kappa,G\rangle\rangle\geq 0$ for all vector $G\in{\mathcal{H}}$. Now
$g_{t}$ defined as $g_{t}=\langle\kappa({\mathbf{x}}_{t},\cdot),\nabla
R(f_{t})\rangle$ is a pseudo-gradient since (8) is satisfied by the symmetric
and positive definiteness property as given above.

  3. iii)

Gradient sign: Define the gradient as $g_{t}=\text{sgn}(\nabla R(f_{t}))$
where sgn denotes the sign operator. Then for all data points
${\mathbf{x}}\in{\mathcal{X}}$, the inner product of (8) is given as
$\langle\nabla R(f_{t}),\text{sgn}(\nabla
R(f_{t}))\rangle=\int_{{\mathcal{X}}}|\nabla
R(f_{t}({\mathbf{x}}))|d{\mathbf{x}}\geq 0$.

The Bregman divergence and pseudo-gradient are employed in the following
section for our main algorithm development.

##  III Projected Pseudo-Mirror Descent

Now we shift to deriving an iterative approach to solving (1) via a functional
extension of mirror descent [42]. We select a Bregman divergence that ensures
positivity of the range of the function $f$ during optimization, such as the
KL divergence, and present a generalization of a gradient called a “pseudo-
gradient” [31]. The merit of employing pseudo-gradients is the ability to
define approximate search directions that arise in point process intensity
estimation as stated in Example 1, as well as a broader ability to incorporate
samples into the functional representation through a kernel embedding.

To define our iterative approach to solving (1) based upon streams of samples,
we first define a first-order functional variant of stochastic mirror descent
using pseudo-gradients. Next, we develop for a Quasi–Newton kernelized
algorithm in mirror space. While theoretically we can only show it matches the
convergence rate of its first-order counterparts, experimentally it achieves
state of the art accuracy for inhomogeneous intensity estimation in Poisson
Processes (Sec. VI).

###  III-A First-Order Method for Sparse Positive Functions

We build upon first-order functional variant of stochastic mirror descent:

$\displaystyle f_{t+1}=\arg\min\limits_{f\in{\mathcal{H}}}$
$\displaystyle\Big{(}\langle\nabla
r_{t}(f_{t}),f\rangle_{{\mathcal{H}}}+\frac{1}{\eta_{t}}B_{\psi}(f,f_{t})\Big{)},$
(9)

where $\eta_{t}>0$ is a nonnegative step-size that may be constant or
diminishing, but we subsequently fix as constant $\eta_{t}=\eta$. Note that
$B_{\psi}(f,\tilde{f})=\frac{1}{2}\|f-\tilde{f}\|^{2}$ reduces (9) to
functional SGD. Moreover, for (9) to be implementable, we require closed-form
expressions for Bregman divergence. Note that for the squared Mahalanobis
distance, the update takes the form
$f_{t+1}=f_{t}-\eta{\mathcal{L}}^{-1}\nabla r_{t}(f_{t})$, where the inverse
${\mathcal{L}}^{-1}$ must be computable in closed-form. Worth mentioning is
also the fact that finding functions with non-negative range (1) necessitates
specifying the Bregman divergence to preserve positivity.

We prioritize the use of pseudo-gradients $g_{t}$ in lieu of stochastic
gradients $\nabla r_{t}(f_{t})$ in (9) since they inherently arise in
evaluating gradient estimates of likelihood models arising in Poisson
Processes. The resultant iteration then takes the form starting from function
$f_{t}$ using samples $\\{{\mathbf{x}}_{t}\\}$:

$\displaystyle{\tilde{f}}_{t+1}=\arg\min\limits_{f\in{\mathcal{H}}}\Big{(}\langle
g_{t},f\rangle_{{\mathcal{H}}}+\frac{1}{\eta}B_{\psi}(f,f_{t})\Big{)},$ (10)

where $g_{t}$ is any pseudo-gradient (8) at sample ${\mathbf{x}}_{t}$.

Mirror Descent with Pseudo-Gradients. To clarify how (10) is executable in
practice, suppose that true stochastic gradient:

$\displaystyle\nabla
r_{t}(f_{t})=\ell^{\prime}(f_{t}({\mathbf{x}}_{t}))\kappa({\mathbf{x}}_{t},\cdot)$
(11)

is used in the mirror descent update [cf. (9)] instead of the pseudo-gradient,
i.e., $g_{t}=\nabla r_{t}(f_{t})$. For this case, the first-order optimality
condition for (10) may be written as

$\displaystyle
g_{t}+\frac{1}{\eta}\nabla_{\tilde{f}_{t+1}}B_{\psi}(\tilde{f}_{t+1},f_{t})=0.$
(12)

Noting the fact
$\nabla_{\tilde{f}_{t+1}}B_{\psi}(\tilde{f}_{t+1},f_{t})=\nabla\psi(\tilde{f}_{t+1})-\nabla\psi(f_{t})$,
(12) can be expressed as

$\displaystyle
g_{t}+\frac{1}{\eta}(\nabla\psi({\tilde{f}}_{t+1})-\nabla\psi(f_{t}))=0.$ (13)

Define ${\tilde{z}}_{t+1}=\nabla\psi({\tilde{f}}_{t+1})$ and
$z_{t}=\nabla\psi(f_{t})$222Subsequently, we note that
$\\{z_{t}\\}\subset{\mathcal{H}}_{\ast}$, i.e., $z_{t}$ is an element of the
dual space ${\mathcal{H}}_{\ast}$ of the RKHS ${\mathcal{H}}$ defined at the
outset of Sec. II., which after rearranging (13), yields

$\displaystyle{\tilde{z}}_{t+1}=z_{t}-\eta g_{t}.$ (14)

Observe that $z_{t}=\nabla\psi(f_{t})$, which implies that the function f can
be recovered through the gradient of the conjugate
$f_{t}=\nabla\psi^{\ast}(z_{t})$. Specifically, $\psi^{\ast}$ is the Fenchel
conjugate of $\psi$ and $\nabla\psi^{\ast}=(\nabla\psi)^{-1}$. Hence the
update (14) becomes

$\displaystyle{\tilde{z}}_{t+1}$
$\displaystyle=z_{t}-\eta\ell^{\prime}(\nabla\psi^{\ast}(z_{t}),{\mathbf{x}}_{t})\kappa({\mathbf{x}}_{t},\cdot),$
(15)

where (15) is derived via applying the chain rule and the reproducing property
of the kernel to $\nabla r_{t}(\nabla\psi^{\ast}(z_{t}))$ analogous to (11).
Thus, the functional update is executable parametrically via a data set
$\mathcal{D}_{t}$ consisting of $t-1$ points that grows by one per time step
with an associated coefficient vector $\mathbf{w}_{t}$:

$\displaystyle\\!\\!\\!{\mathcal{D}}_{t+1}\\!=\\!{\mathcal{D}}_{\\!t}\\!\cup\\{{\mathbf{x}}_{t}\\},\
$
$\displaystyle[\mathbf{w}_{t+1}]_{n}\\!=\\!\\!\begin{cases}\\![\mathbf{w}_{t}]_{n}&\\!\\!n<t\\!\\\
\\!-\eta\ell^{\prime}(\nabla\psi^{\ast}(z_{t}\\!),{\mathbf{x}}_{t})&\\!\\!n=t,\\!\end{cases}\\!\\!\\!$
(16)

where $[\mathbf{w}_{t}]_{n}$ denotes the $n$-th coordinate of the vector
$\mathbf{w}_{t}$. Different from classical online learning with kernels [17],
here the weights $[\mathbf{w}_{t}]_{n}$ and basis points ${\mathcal{D}}_{t}$
represent the function $z_{t}$ associated with the gradient of the Fenchel
dual in (15).

We note that the function evaluation $f$ at point
${\mathbf{x}}\in{\mathcal{X}}$ then takes the form at time $t+1$

$\displaystyle
f_{t+1}({\mathbf{x}})=\nabla\psi^{\ast}(z_{t+1}({\mathbf{x}}))=\nabla\psi^{\ast}({\mathbf{w}}_{t+1}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{t+1}}({\mathbf{x}})).$
(17)

A technical challenge emerges due to the fact that for general pseudo-
gradients $g_{t}\neq\nabla r_{t}(f_{t})$, the expression (13) cannot be
evaluated in closed form. However, for the pseudo-gradient defined by the
kernel embedding (ii), (16) takes the form

$\displaystyle{\mathcal{D}}_{t+1}\\!=\\!{\mathcal{D}}_{t}\\!\cup\\!\\{\\!{\mathbf{x}}_{t}\\!\\}\;,\
$
$\displaystyle[\mathbf{w}_{t+1}]_{n}\\!=\\!\\!\begin{cases}\\![\mathbf{w}_{t}]_{n}&\\!{\mathbf{x}}_{n}\\!\in\\!{\mathcal{D}}_{t}\\\
\\!-\eta
g_{t}^{\prime}&\\!{\mathbf{x}}_{n}\\!=\\!{\mathbf{x}}_{t},\end{cases}$ (18)

where the pseudo-gradient admits an expansion via the chain rule
$g_{t}=g_{t}^{\prime}\kappa({\mathbf{x}}_{t},\cdot)$, where $g_{t}^{\prime}$
is the scalar component of the derivative of the cost $r_{t}(f)$. The kernel
embedding, or any other pseudo-gradient, gives rise to the specific form of
$g_{t}^{\prime}$. This condition is valid when the pseudo-gradient of the
objective with respect to RKHS element $f$ respects the chain rule.

Subspace Projections. Owing to the RKHS parameterization in terms of weights
and feature vectors ${\mathbf{x}}_{t}$ [cf. (18)], the complexity of the
function grows unbounded with time $t$. To surmount this issue, we project
functions onto low-dimensional subspaces near the current search direction.
Doing so may be achieved via Kernel Orthogonal Matching Pursuit (KOMP)[25],
that takes an input dictionary ${\tilde{{\mathcal{D}}}}_{t+1}$ and a weight
vector ${\tilde{\mathbf{w}}}_{t+1}$, returning a lower-dimensional
(compressed) function $z_{t+1}$ with dictionary ${\mathcal{D}}_{t+1}$ and
weights ${\mathbf{w}}_{t+1}$ that are $\epsilon$-away in the RKHS norm, where
$\epsilon$ denotes the compression budget. Overall, then, Sparse Positive
Projected Pseudo-Mirror Descent (SPPPOT), takes the form

$\\{z_{t+1},{\mathcal{D}}_{t+1},{\mathbf{w}}_{t+1}\\}=\text{KOMP}(\tilde{z}_{t+1},{\tilde{{\mathcal{D}}}}_{t+1},{\tilde{\mathbf{w}}}_{t+1},\epsilon),$
(19)

where $\tilde{z}_{t+1}$ defined in (14) employs general pseudo-gradients. The
iterative process is summarized as Algorithm 1. Note that different from [26]:
the RKHS-norm approximation criterion is in terms of the dual sequence
$\\{z_{t}\\}$ (15) – see Algorithm 2.

Algorithm 1 SPPPOT: Sparse Positive Functions via Projected Pseudo-Mirror
Descent

0: kernel $\kappa$, step-size $\eta$, compression parameter $\epsilon$

1: Initialize  Arbitrary small $z_{0}$

2: for $t=1,2,\ldots$ do

3: Read:  data realization ${\mathbf{x}}_{t}$

4: Evaluate:  Pseudo Gradient $g_{t}$

5: Update:  ${\tilde{z}}_{t+1}$ as per (14)

6: Compress:
$\\{{\mathcal{D}}_{t+1},{\mathbf{w}}_{t+1}\\}=\text{KOMP}({\tilde{{\mathcal{D}}}}_{t+1},{\tilde{\mathbf{w}}}_{t+1},\epsilon)$

7: Broadcast:  $z_{t+1}$

8: end for

9: Evaluation of actual function $f_{t+1}$ at ${\mathbf{x}}$: As per (17)

Examples. With Algorithm 1 defined, we next specify the updates in terms of
two Bregman divergences defined in Sec. III, one that ensures positivity, (i)
KL-divergence, and one that does not: (ii) squared RKHS norm.

  1. i)

KL-divergence or I-divergence: As in (i), the update for the auxiliary
variables $\tilde{z}_{t+1}$ and $z_{t}$ (14) have the explicit forms
${\tilde{z}}_{t+1}=\log({\tilde{f}}_{t+1})$ and $z_{t}=\log(f_{t})$, which may
be substituted into (14) to yield an expression for $\tilde{f}_{t+1}$ as

$\displaystyle\log({\tilde{f}}_{t+1})=\log(f_{t})-\eta g_{t},$ (20)

which, following exponentiation, permits us to write

$\displaystyle{\tilde{f}}_{t+1}=f_{t}\exp(-\eta g_{t}).$ (21)

Observe from (21) that the update for $f_{t}$ is non-linear and hence a
vanilla application of the Representer Theorem (5) is not possible. However,
by focusing on the auxiliary sequences
${\tilde{z}}_{t+1}=\log({\tilde{f}}_{t+1})$ and $z_{t}=\log(f_{t})$, (21)
permits linearization of the form (14), and hence allows a parametric weight
update and dictionary update with respect to the auxiliary function
${\tilde{z}}_{t+1}$ akin to (16). Note that from (21), the function
${\tilde{f}}_{t+1}$ can be expressed recursively as a product of functions
from $f_{0}$ to $f_{t}$. So if the initialization $f_{0}$ is nonnegative, and
the projection in (19) is onto the space of functions with nonnegative range,
then function positivity is ensured throughout training. This is salient for
Example 1, and MLE more broadly.

Algorithm 2 Destructive Kernel Orthogonal Matching Pursuit

Require: ${\tilde{z}}$ in form of
$({\tilde{{\mathcal{D}}}},{\tilde{\mathbf{w}}})$, budget $\epsilon$

Initialize: $z={\tilde{z}}$ so that
$({\mathcal{D}},{\mathbf{w}})=({\tilde{{\mathcal{D}}}},{\tilde{\mathbf{w}}})$

while ${\mathcal{D}}\neq\emptyset$ do

for ${\mathbf{x}}_{j}\in{\mathcal{D}}$ do

Evaluate
${\mathbf{w}}_{j}^{\star}\\!=\\!\arg\min_{{\mathbf{w}}}\\!\gamma_{j}({\mathbf{w}})\\!\\!:=\\!\|{\tilde{z}}-\\!\\!\\!\\!\\!\sum\limits_{{\mathbf{x}}_{n}\in{\mathcal{D}}\setminus\\{{\mathbf{x}}_{j}\\}\\!\\!\\!\\!\\!\\!}\\!\\!\\!\\!\\!\\!\\!\\!w_{n}\kappa(\\!{\mathbf{x}}_{n},\\!\cdot\\!)\\!\|_{\ast}$

end for

if $\gamma_{j}({\mathbf{w}}_{j}^{\star})>\epsilon$ for all
${\mathbf{x}}_{j}\in{\mathcal{D}}$ then

break

else

Prune
${\mathcal{D}}\leftarrow{\mathcal{D}}\setminus\\{{\mathbf{x}}_{j^{\star}}\\}$
where $j^{\star}=\arg\min\gamma_{j}$

Update weights ${\mathbf{w}}\leftarrow{\mathbf{w}}_{j^{\star}}^{\star}$

end if

end while

return $z$, such that $\left\|z-{\tilde{z}}\right\|_{\ast}\leq\epsilon$

  2. ii)

Contrasting example: Squared RKHS-norm difference: With (ii), the quantities
in (14) take the form
${\tilde{z}}_{t+1}=\nabla\psi({\tilde{f}}_{t+1})={\tilde{f}}_{t+1}$ and
$z_{t}=\nabla\psi(f_{t})=f_{t}$, which yields

$\displaystyle{\tilde{f}}_{t+1}=f_{t}-\eta g_{t},$ (22)

where the dictionary stacks past points [cf. (16)], and coefficients are
updated as

$\displaystyle{\tilde{\mathbf{w}}}_{t+1}=[{\mathbf{w}}_{t},\quad-\eta g_{t}].$
(23)

With KOMP-based subspace projections, this update is identical to POLK [26] if
one additionally equates the pseudo-gradient with the stochastic gradient
$g_{t}=\nabla r_{t}(f_{t})$ [cf. (11)]. For this choice of Bregman divergence,
under positive initialization, positivity will be violated.

A key observation from the case of KL-Divergence is that since the update for
$f$ does not yield a linear basis expansion, it cannot be computed through
weighted combinations of kernel evaluations (5). However, it can be recovered
through inverting the logarithmic transformation, i.e.,
$f_{t+1}(\cdot)=\exp(z_{t+1}(\cdot))=\exp({\mathbf{w}}_{z,t+1}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{z,t+1}}(\cdot))$.
Consequently, as each new datum ${\mathbf{x}}$ arrives, we can obtain its
evaluation under $f$ as
$f_{t+1}({\mathbf{x}})=\exp({\mathbf{w}}_{z,t+1}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{z,t+1}}({\mathbf{x}}))$.
In addition, $f_{t+1}({\mathbf{x}})$ will be positive since it is obtained by
positivity preserving exponential transformation on $z_{t+1}({\mathbf{x}})$.
Observe that this phenomenon holds more broadly for any positivity preserving
Bregman divergence when its update is nonlinear in $f$ and the gradient of the
Fenchel dual $\nabla\psi^{*}(f)$ of its inducing map $\psi(f)$ is efficiently
computable.

Next, we expand upon the form of Algorithm 1[Steps 4-5] for Poisson point
process intensity estimation (Example 1), which are employed experimentally in
Sec. VI.

######  Example 2.

Observe that for (3), the instantaneous loss for single data for the function
$f_{t}$ takes the form

$\displaystyle\ell(f_{t}({\mathbf{x}}_{t}))=-\log(f_{t}({\mathbf{x}}_{t}))+\int_{{\mathcal{X}}}f_{t}({\mathbf{x}})d{\mathbf{x}}.$
(24)

When no special form of $f_{t}$ is known for analytically, one must evaluate
the integral numerically using, e.g., Bayesian quadrature or kernel smoothing.
We adopt the later approach, inspired by [34]: the quantity
$\int_{{\mathcal{X}}}f_{t}({\mathbf{x}})d{\mathbf{x}}$ is approximated by
$h\sum_{j\in{\mathcal{U}}}f_{t}({\mathbf{u}}_{j})$, where ${\mathbf{u}}_{j}$
are uniform grid points over the sample space ${\mathcal{X}}$ whose indices
lie in set ${\mathcal{U}}:=\\{{\mathbf{u}}_{j}\\}$, and $h$ is the
infinitesimally small uniform grid area, similar to Trapezoidal rule for
discrete approximations of integrals. Then differentiating with respect to
$f_{t}$ yields the pseudo-gradient

$\displaystyle
g_{t}=-\frac{1}{f_{t}({\mathbf{x}}_{t})}\kappa({\mathbf{x}}_{t},\cdot)+h\sum_{j=1}^{|\mathcal{U}|}\kappa({\mathbf{u}}_{j},\cdot).$
(25)

(25) is a pseudo-gradient as the second term of (25) can be seen as an
approximated version of the integral using kernel embedding. Note that for the
KL divergence/I-divergence, the transformation $f_{t}=\exp(z_{t})$ yields the
pseudo-gradient $g_{t}$:

$\displaystyle
g_{t}=-\frac{1}{\exp(z_{t}({\mathbf{x}}_{t}))}\kappa({\mathbf{x}}_{t},\cdot)+h\sum_{j=1}^{|\mathcal{U}|}\kappa({\mathbf{u}}_{j},\cdot).$
(26)

The subtlety of grid points versus Poisson points causes slight differences in
the updates relative to (16). Specifically, in (26), the point
${\mathbf{x}}_{t}$ are samples from the unknown Poisson process, while the
points ${\mathbf{u}}_{j}$ are fixed grid points across updates. These
components come together to specify Algorithm 1 as follows:

  * •

Initialize the dictionary ${\tilde{{\mathcal{D}}}}_{1}$ with uniform grid
points ${\mathbf{u}}_{j}$ for all ${\mathbf{u}}_{j}\in\mathcal{U}$, with
corresponding weight vector elements as $[{\tilde{\mathbf{w}}}_{1}]_{j}=-\eta
h$ for each ${\mathbf{u}}_{j}\in{\mathcal{U}}$.

  * •

Receive Poisson samples ${\mathbf{x}}_{t}$, compute the pseudo-gradient (26),
and update dictionary
${\tilde{{\mathcal{D}}}}_{t+1}={\mathcal{D}}_{t}\cup\\!\\{\\!{\mathbf{x}}_{t}\\!\\}$.
Given that the set of grid points satisfies
$\mathcal{U}\subset\mathcal{D}_{t}$, the corresponding weight update is

$\displaystyle[\tilde{\mathbf{w}}_{t+1}]_{n}=\begin{cases}[\mathbf{w}_{t}]_{n}-\eta
h&{\mathbf{x}}_{n}\in\mathcal{U}\\\
[\mathbf{w}_{t}]_{n}&{\mathbf{x}}_{n}\notin\mathcal{U}\\\
\frac{\eta}{\exp(z_{t}({\mathbf{x}}_{t}))}&{\mathbf{x}}_{n}={\mathbf{x}}_{t}\;.\end{cases}$
(27)

KOMP is then applied to the dictionary ${\tilde{{\mathcal{D}}}}_{t+1}$ and
weights $\tilde{\mathbf{w}}_{t+1}$ to yield ${\mathcal{D}}_{t+1}$ and
$\mathbf{w}_{t+1}$ . Since grid points $\\{{\mathbf{u}}_{j}\\}$ are required
to approximate the integral in the pseudo-gradient, their presence in the
dictionary is fixed. Therefore, our selection scheme discerns which Poisson
samples ${\mathbf{x}}_{t}$ and associated weights $w_{t}$ are statistically
significant for estimating the inhomogeneous intensity parameter
$f(\cdot)=\lambda(\cdot)$ as in Example 1.

Next, we verify that (26) is indeed a pseudo-gradient. Following similar
reasoning to[29], the inner product between $\nabla R(f_{t})$ and
$\mathbb{E}[g_{t}|\mathcal{F}_{t}]$ can be written as $\langle\nabla
R(f_{t}),\mathbb{E}[g_{t}|\mathcal{F}_{t}]\rangle=\int_{\mathcal{X}}\int_{\mathcal{X}}\left(-\frac{1}{f_{t}(\mathbf{x})}+1\right)\kappa(\mathbf{x},\mathbf{y})\left(-\frac{1}{f_{t}(\mathbf{y})}+1\right)d\mathbf{x}d\mathbf{y}$,
where we have taken a continuous domain data for algebraic simplicity and
hence writing in terms of empirical sum is avoided. The above equation is
lower bounded by
$\zeta_{\min}\left\lVert-\frac{1}{f_{t}}+\mathbf{1}\right\rVert^{2}$ using the
positive definiteness property of the kernel, where $\mathbf{1}$ is vector of
$1$’s matching the dimension of $f_{t}$ and $\zeta_{\min}$ is the minimum
eigenvalue of the integral operator associated with $\kappa(\cdot,\cdot)$.
Thus, the inner-product will be non-negative, meaning the pseudo-gradient
property [cf. (8)] holds.

With Poisson intensity estimation example detailed, we next shift to
incorporating approximate Hessian information.

##  IV Second-Order Online Pseudo-Mirror Descent

Iterative updates that incorporate Hessian information, i.e., variants of
Newton’s Method, often surpass schemes that only use gradient information in
convergence rate and numerical precision. However, they often are inapplicable
to online settings due to the quadratic computational cost in the sample size
of Hessian inversion [39][Ch. 6]. Recently, memory-efficient Quasi–Newton
schemes in the online setting based on incremental Hessian updates have been
developed [43, 44].

However, these techniques do not generalize to the RKHS setting due to their
dependence on finite-dimensional parametric forms for the search space.
Efforts to bridge this gap have been addressed by introducing Nyström
approximations to the kernel matrix [20], resulting in memory-efficient
schemes [32]. Noticeably, however, Nyström approximation does not permit the
points of approximation to be well-calibrated to the data domain. By contrast,
methods which optimize the location from which kernel matrices are sampled via
_inducing_ inputs can refine their statistical error [45, 24]. Hence,
hereafter, we develop a Quasi–Newton scheme via inducing input approximations
to the kernel matrix, which additionally employs pseudo-gradients rather than
true stochastic gradients.

To continue, we denote as ${\tilde{g}}_{t}$ [cf. (8)] the pseudo-gradient
defined over a fixed subspace $|\mathcal{D}|$ associated with weight vector
$\mathbf{w}_{t}\in\mathbb{R}^{|\mathcal{D}|}$. Associated with it is a a
Hessian approximation
$\mathbf{A}_{t+1}\in\mathbb{R}^{|\mathcal{D}|\times|\mathcal{D}|}$ that is
approximated using rank-1 outer products of pseudo-gradients observed so far
as

$\displaystyle\mathbf{A}_{t+1}=\mathbf{A}_{t}+{\tilde{g}}_{t}{\tilde{g}}_{t}^{\top}=\delta\mathbf{I}+\sum_{j=1}^{t}{\tilde{g}}_{j}{\tilde{g}}_{j}^{\top},$
(28)

where, $\delta\mathbf{I}$ is the initialized Hessian matrix $\mathbf{A}_{0}$.
Observe that without any subspace approximation, i.e., if
$\mathcal{D}=\mathcal{D}_{t}$ is allowed to grow unbounded as (14), then the
dimensionality of $\mathbf{A}_{t}$ will be $(t-1)\times(t-1)$ at step $t$. To
address this issue, we fix our subspace of inputs defined by data matrix
${\mathcal{D}}_{t}={\mathcal{D}}$ for all $t$, in contrast to the Nyström
approximations in [32]. Doing so ensures sufficient coverage of data matrix
${\mathcal{D}}$ by assigning grid points over the feature space $\mathcal{X}$.
We note that to more closely adhere to the concept of inducing inputs [45,
24], one could optimize the locations by searching over a space of
dimensionality $|\mathcal{D}|\times{\color[rgb]{0,0,1}d}$ after training
completes. With this approximation, then, (28) may be implemented online when
${\tilde{g}}_{t}$ is computable, which incurs storage complexity
$|\mathcal{D}|\times|\mathcal{D}|$.

Observe that with this dimensionality fixed, optimization over the dual
sequence of functions $z_{t}$ [cf. (14)] is transformed to the search over
fixed dimension of weights
${\mathbf{w}}_{t}\in{\mathbb{R}}^{\lvert{\mathcal{D}}\rvert}$ for all $t$,
where $\lvert{\mathcal{D}}\rvert$ is the cardinality of the set
${\mathcal{D}}$. Here, the pseudo-gradient, computed by differentiating the
loss with respect to the fixed dimensional weights, is generically expressed
as
${\tilde{g}}_{t}={\tilde{g}}_{t}^{\prime}{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{t})$.
For the case of a stochastic gradient, this simplifies to
${\tilde{g}}_{t}^{\prime}=\nabla_{{\mathbf{w}}_{t}}\ell(f_{t}({\mathbf{x}}_{t}))=\ell^{\prime}(f_{t}({\mathbf{x}}_{t}))$.
Note the fact that ${\tilde{g}}_{t}\in\mathbb{R}^{|\mathcal{D}|}$. Now, the
gradient outer product on the right hand side of (28) can be expressed as
${\tilde{g}}_{t}{\tilde{g}}_{t}^{\top}=[{\tilde{g}}_{t}^{\prime}]^{2}{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{t}){\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{t})^{\top}$,
which is always positive. Hence the Hessian matrix $\mathbf{A}_{t+1}$ is
always positive definite and also symmetric. The Hessian matrix estimate is
positive definite via the initialization $\delta\mathbf{I}$.

Overall, the Hessian estimate (28) and the feature space approximation
together yield the online Quasi–Newton update

$\displaystyle{\mathbf{w}}_{t+1}={\mathbf{w}}_{t}-\eta_{t}\mathbf{A}_{t+1}^{-1}{\tilde{g}}_{t}.$
(29)

Observe that since the matrix
$\mathbf{A}_{t}\in\mathbb{R}^{|\mathcal{D}|\times|\mathcal{D}|}$ belongs to a
fixed subspace of $|\mathcal{D}|\times|\mathcal{D}|$, its inversion may be
executed with complexity $|\mathcal{D}|^{3}$ in the worst case with a naive
implementation, or with $|\mathcal{D}|^{2}$ through application of matrix
inversion lemmas. While this is less favorable than a near-linear dependence
on the parameter dimension achievable by the best memory-limited Quasi–Newton
methodologies, it is far better than the computational cost of a full Newton
step. In particular, computing a full Newton step for this setting would
require evaluating the Hessian with respect to the $|\mathcal{D}|$-dimensional
kernel truncation, which depends on doubly infinitely many realizations
$\mathbf{x}_{n}$ of the random variable $\mathbf{x}$. Via (29), the auxiliary
functional iterates may be expressed as

$\displaystyle
z_{t+1}={\mathbf{w}}_{t+1}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)=z_{t}-\eta_{t}{\tilde{g}}_{t}^{\top}\mathbf{A}_{t+1}^{-1}{\mathbf{k}}_{{\mathcal{D}}}(\cdot).$
(30)

(30) defines a Quasi–Newton update on the dual iterate $z_{t}$ over the fixed-
dimensional feature subspace $|{\mathcal{D}}|$. One may recover the original
function estimate via computation of the Fenchel dual (17) when the Bregman
divergence admits a closed-form expression for the evaluation of the gradient
of its conjugate (as is the case for KL divergence, and others). Next we
specify the form of these updates for Poisson process.

Algorithm 3 Online Quasi–Newton Mirror Descent in RKHS

0: kernel $\kappa$, step-size $\eta_{t}$, regularizer $\delta$

1: Initialize  Arbitrary small negative
${\mathbf{w}}_{0}\in{\mathbb{R}}^{\lvert{\mathcal{U}}\rvert}$,
$\mathbf{A}_{0}=\delta\mathbf{I}\in{\mathbb{R}}^{\lvert{\mathcal{U}}\rvert\times\lvert{\mathcal{U}}\rvert}$,
dictionary as grid points $\\{{\mathcal{D}}_{t}\\}=\\{{\mathcal{U}}\\}\
\forall t$.

2: for $t=1,2,\ldots$ do

3: Read:  sample ${\mathbf{x}}_{t}$

4: Evaluate:  Pseudo Gradient ${\tilde{g}}_{t}$ as (33)

5: Update:
$\mathbf{A}_{t+1}=\mathbf{A}_{t}+{\tilde{g}}_{t}{\tilde{g}}_{t}^{\top}$

6: Update:
${\mathbf{w}}_{t+1}={\mathbf{w}}_{t}-\eta_{t}\mathbf{A}_{t+1}^{-1}{\tilde{g}}_{t}$

7: To evaluate functional estimate, invert mirror map:
$f_{t+1}({\mathbf{x}})=(\nabla\psi)^{-1}(z_{t}({\mathbf{x}}))=\exp({\mathbf{w}}_{t+1}{\mathbf{k}}_{{\mathcal{D}}_{t+1}}({\mathbf{x}}))$
for KL divergence.

8: end for

######  Example 3.

Quasi–Newton Poisson Intensity Estimation. Consider the objective in (24) that
integral approximations with respect to the uniform grid points
${\mathbf{u}}_{j}\in{\mathcal{U}}$. In the context of the second-order scheme,
we fix the dictionary with this specification, i.e.
${\mathcal{D}}_{t}={\mathcal{D}}={\mathcal{U}}$ for all $t$, meaning that no
compression is required. Then, the first-order component of the update
manifests by updating the auxiliary variable $z_{t}$ according to the gradient
of the Fenchel dual (14), as in mirror descent (natural gradient) [46]. In
particular, we rewrite the loss in (24) with respect to the auxiliary function
$z_{t}$ as

$\displaystyle\ell(f_{t}({\mathbf{x}}_{t}))=\ell((\nabla\psi)^{-1}(z_{t}({\mathbf{x}}_{t})))=-\log((\nabla\psi)^{-1}(z_{t}({\mathbf{x}}_{t})))+h\sum_{{\mathbf{u}}_{j}\in{\mathcal{U}}}(\nabla\psi)^{-1}(z_{t}({\mathbf{u}}_{j})).$
(31)

Specialized to the case of KL divergence, the expression
$(\nabla\psi)^{-1}(z_{t}({\mathbf{x}}_{t}))=\exp(z_{t}({\mathbf{x}}_{t}))$
substituted into (31) yields

$\displaystyle\ell(f_{t}({\mathbf{x}}_{t}))=-z_{t}({\mathbf{x}}_{t})+h\sum_{{\mathbf{u}}_{j}\in{\mathcal{U}}}\exp(z_{t}({\mathbf{u}}_{j}))=-{\mathbf{w}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{t})+h\sum_{{\mathbf{u}}_{j}\in{\mathcal{U}}}\exp({\mathbf{w}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{u}}_{j})).$
(32)

Then, differentiating (32) with respect to weights
${\mathbf{w}}_{t}\in\mathbb{R}^{|\mathcal{U}|}$ (since
$\mathcal{U}=\mathcal{D}$) yields the pseudo gradient
${\tilde{g}}_{t}\in{\mathbb{R}}^{\lvert{\mathcal{U}}\rvert}$:

$\displaystyle{\tilde{g}}_{t}\\!=\\!\nabla_{\\!{\mathbf{w}}_{t}}\ell(f_{t}({\mathbf{x}}_{t}))\\!=\\!-{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}}_{t})\\!+\\!h\\!\sum_{j\in{\mathcal{U}}}\\!\exp({\mathbf{w}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}({\mathbf{u}}_{j})){\mathbf{k}}_{{\mathcal{D}}}({\mathbf{u}}_{j}).$
(33)

We proceed to formalizing a special case of update direction in (30) and
ensure it defines pseudo-gradient in the sense of (8). Begin by noting that we
may write the pseudo-gradient as an RKHS element through appropriate
multiplication with the feature map ${\mathbf{k}}_{{\mathcal{D}}}(\cdot)$
associated with fixed dictionary ${\mathcal{D}}$. Then, the functional
gradient representation
${\tilde{g}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)$ represents the
derivative of the loss in (32) with respect to auxiliary variable $z_{t}$,
i.e.,
$\nabla_{z_{t}}(\ell(f_{t}(x_{t})))={\tilde{g}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)$.
Further note that $\nabla R_{\psi}(z_{t})=\nabla R(f_{t})$ where
$f_{t}\in{\mathcal{H}}$ and $z_{t}\in{\mathcal{H}}_{\ast}$. Also denote the
Gram matrix of kernel evaluations as
${\mathbf{K}}_{{\mathcal{D}}{\mathcal{D}}}$, whose entries are
$[\kappa(\mathbf{d}_{n},\mathbf{d}_{m})]_{n\times
m=1}^{|{\mathcal{D}}|,|{\mathcal{D}}|}$. We compute the conditional
expectation of ${\tilde{g}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)$ and
$\nabla R_{\psi}(z_{t})$ as

$\displaystyle\langle\nabla
R_{\psi}(z_{t}),\mathbb{E}[{\tilde{g}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)|\mathcal{F}_{t}]\rangle=\int_{\mathcal{X}}\int_{\mathcal{X}}\left(-1+f_{t}(\mathbf{x})\right){\mathbf{k}}_{{\mathcal{D}}}({\mathbf{x}})^{\top}{\mathbf{K}}_{{\mathcal{D}}{\mathcal{D}}}{\mathbf{k}}_{{\mathcal{D}}}(\mathbf{y})\left(-1+f_{t}(\mathbf{y})\right)d\mathbf{x}d\mathbf{y}\;,$
(34)

which can be lower bounded by
$\zeta_{\min}\left\lVert-\mathbf{1}+f_{t}\right\lVert^{2}\geq 0$ by the same
logic as detailed in Example 2. This provides a basis for showing our second-
order step also defines a pseudo-gradient. In particular, the functional
gradient and expected value of Quasi–Newton pseudo-gradient [cf. (30)] given
the past sigma algebra which can be expressed as

$\displaystyle\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}[{\tilde{g}}_{t}^{\top}\mathbf{A}_{t+1}^{-1}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)|{\mathcal{F}}_{t}]=\left\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}\left[\frac{{\tilde{g}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)}{\delta+\sum_{j=1}^{t}{\tilde{g}}_{j}^{\top}{\tilde{g}}_{j}}|{\mathcal{F}}_{t}\right]\right\rangle,$
(35)

where the matrix simplification is similar to [33] since $\mathbf{A}_{t+1}$ is
positive definite. Now since the inner product between $\nabla
R_{\psi}(z_{t})$ and
$\mathbb{E}[{\tilde{g}}_{t}^{\top}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)|\mathcal{F}_{t}]$
is non-negative and
$\delta+\sum_{j=1}^{t}{\tilde{g}}_{j}^{\top}{\tilde{g}}_{j}$ is positive, we
may conclude the quantity in (35) is non-negative. Therefore, it defines a
pseudo-gradient.

Algorithm 4 Hybrid First and Second-order Online Iteration

0: kernel $\kappa$, constant step-size $\eta$, regularizer $\delta$

1: Initialize  Arbitrary small negative
${\mathbf{w}}_{0}\in{\mathbb{R}}^{\lvert{\mathcal{U}}\rvert}$ and
$\\{{\mathcal{D}}_{0}\\}=\\{{\mathcal{U}}\\}$.

2: Run:  Alg. 1 until stable model order $M^{\infty}$ at step $t=T$.

3: Denote $(\mathbf{w}_{T},\mathcal{D}_{T})$ as output of Alg. 1 on step $T$.

4: Initialize  Step size $\eta_{T}$ and
$\mathbf{A}_{0}=\delta\mathbf{I}\in{\mathbb{R}}^{M\times M}$.

5: for $k=T,T+1,\ldots$ do

6: Read:  sample ${\mathbf{x}}_{k}$

7: Run:  Algorithm 3

8: To obtain function estimate, invert mirror map:
$f_{t+1}({\mathbf{x}})=(\nabla\psi)^{-1}(z_{t}({\mathbf{x}}))=\exp({\mathbf{w}}_{t+1}{\mathbf{k}}_{{\mathcal{D}}_{t+1}}({\mathbf{x}}))$
for KL divergence.

9: end for

The overall procedure for incorporating Quasi–Newton updates with pseudo-
gradients is presented in Algorithm 3. It is essentially similar to Algorithm
1 with the additional outer product of gradients defining a recursively
constructed pre-conditioner for the pseudo-gradient updates that are executed
in the dual space. We also present a hybrid algorithm as a mixture of
Algorithms 1 and 3, which is presented in Algorithm 4. This combination of
SPPPOT and the Quasi–Newton Mirror Step allows one to flexibly learn the
appropriate dictionary to obtain a functional estimate of reasonable
precision, and then employ the learned dictionary to define the approximation
subspace for executing Newton steps to obtain a function estimate of maximum
precision. In this case, SPPPOT employs updates based upon functional
differentiation of (24), whereas the gradient for kernel online Newton mirror
step is with respect to the weights of the auxiliary function [cf. (33)] over
a fixed dictionary. Next, we establish convergence and complexity tradeoffs of
Algorithms 1 and 3.

##  V Convergence analysis

We shift towards analyzing the convergence behavior of Algorithms 1 and 3 in
terms of iteratively solving (1), when the compression budget $\epsilon$ and
the step-size parameter $\eta$ are held constant for Algorithm 1. For
Algorithm 3, we prove the convergence under constant step-size in terms of the
fixed sub-space approximation for the Hessian. A key point of departure in the
analysis of Algorithm 3 from prior works (e.g. [33, 32]) is the need to
quantify the subspace approximation error associated with a fixed dictionary
as compared to employing a fully infinite dimensional feature map
${\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)$ in the Quasi-Newton direction.
Denote as ${\tilde{g}}_{\infty,t}$ and $\mathbf{A}_{\infty,t+1}$,
respectively, as corresponding infinite dimensional gradient and doubly
infinite dimensional Hessian at iterate $t$. Moreover, for notational
simplicity, we denote the associated Quasi-Newton directions using a fixed
subspace of points ${\mathcal{D}}$ its infinite-dimensional counterpart as

$\displaystyle
G_{t}={\tilde{g}}_{t}^{\top}\mathbf{A}_{t+1}^{-1}{\mathbf{k}}_{{\mathcal{D}}}(\cdot)\;,\;\;G_{t}^{(\infty)}={\tilde{g}}_{\infty,t}^{\top}\mathbf{A}_{\infty,t+1}^{-1}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot).$
(36)

To proceed with our analysis, some technical conditions are required on the
objective function and its gradients, as well as stochastic approximation
errors, which we state next.

######  Assumption 1.

The inner product between the gradient and the expectation of the pseudo-
gradient (13) and its fixed-subspace Quasi–Newton variant (29) conditioned on
the filtration $\mathcal{F}_{t}=\sigma(\\{{\mathbf{x}}_{i}\\}_{i=1}^{t-1})$,
are nonnegative.

$\displaystyle\\!\\!\\!\\!\langle\nabla
R(f_{t}),\\!{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle\\!\geq
0\;,\langle\nabla
R_{\psi}(z_{t}),\\!{\mathbb{E}}[{{\tilde{g}}_{\infty,t}}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)|{\mathcal{F}}_{t}]\rangle\\!\geq\\!0.$
(37)

Moreover for $f_{t}\neq f^{\ast}$, the product-moment is lower bounded by the
second-moment of the gradient in the dual norm:

$\displaystyle\\!\\!\\!\mathbb{E}[\langle\nabla
R(f_{t}),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle]\geq
D{\mathbb{E}}[||\nabla R(f_{t})||_{\ast}^{2}]\;,$ (38)
$\displaystyle\\!\\!\\!{\mathbb{E}}[\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}[{{\tilde{g}}_{\infty,t}}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)|{\mathcal{F}}_{t}]\rangle]\geq
D{\mathbb{E}}[||\nabla R_{\psi}(z_{t})||_{\ast}^{2}],$ (39)

where $D$ is a positive constant.

######  Assumption 2.

The optimizer of (1) is finite $\left\|f^{\star}\right\|^{2}\leq B$.

######  Assumption 3.

The function $R(\cdot)$ satisfies the Polyak–Łojasiewicz (P–Ł) condition

$\displaystyle\left\|\nabla R(f_{t})\right\|^{2}_{\ast}\geq
2\lambda[R(f_{t})-R(f^{\ast})]\;,$ (40)

where $\lambda$ is a positive constant.

######  Assumption 4.

The function $R_{\psi}(\cdot)$ which takes as inputs the dual functions
$z=\nabla\psi(f)$ is $L_{1}$-smooth.

######  Assumption 5.

The total expectation of the pseudo-gradient $g_{t}$ [cf. (13)] is upper-
bounded in dual norm as

$\displaystyle{\mathbb{E}}[\left\|g_{t}\right\|^{2}_{\ast}]$
$\displaystyle\leq b^{2}\\!+\\!c^{2}{\mathbb{E}}[\langle\nabla
R(f),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle]\;,$ (41)

whereas we assume the Quasi–Newton direction satisfies

$\displaystyle{\mathbb{E}}[\left\|{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{\\!{\mathcal{D}}_{\\!\infty}}\\!(\cdot\\!)\right\|^{2}_{\ast}]\\!$
$\displaystyle\leq\\!\\![b^{\prime}]^{2}\\!\\!+\\![c^{\prime}]^{2}{\mathbb{E}}[\langle\nabla
R_{\psi}\\!(z_{t}\\!),{\mathbb{E}}[{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{\\!{\mathcal{D}}_{\\!\infty}}\\!(\cdot\\!)|{\mathcal{F}}_{t}]\rangle],$
(42)

for all $f\in{\mathcal{H}}$, $t\in{\mathbb{N}}$, and some real constants $b,c$
and $b^{\prime},c^{\prime}$.

Assumption 1 asserts that for $f\neq f^{\ast}$, an acute angle exists between
the pseudo-gradient $g_{t}$ and actual gradient $\nabla R(f)$ which is no
larger than 90 degrees, and the constant $D$ determines the minimal degree of
collinearity. Gradients for SPPPOT (end of Sec. III) and Quasi–Newton Mirror
Step [cf. (34)] satisfies this assumption. Note that any stochastic gradient
is also a pseudo gradient. This condition is used in establishing decrement-
like relationships in the second term on the right hand side of (52).
Assumption 2 is employed to ensure that $R(f_{0})-R(f^{\ast})$ is finite. The
P–Ł condition in Assumption 3, often presented as a consequence of strong
convexity, holds for many non-convex functions (such as matrix completion
[47]), and connects decrease in the gradient norm to progress towards the
optimal objective $R(f^{*})$. Assumption 4 is standard in the analysis of
mirror descent and proximal methods [9], and is used in the Bregman “three-
point inequality” (49). Estimating the exact parameters of this condition is
difficult, but may be conducted numerically. Doing so is beyond the scope of
this work. Assumption 5 (similar to [29][Theorem 3]) is weaker than the
standard second-moment boundedness condition on the stochastic gradient, and
instead permits possibly _unbounded_ gradients under certain growth conditions
(see [48][Sec. 4], [49], and [50] eqn. (4.9) and Assumptions 4.3.). It holds
in the RKHS setting when the dual norm of the universal kernel is bounded, and
the optimizer belongs to a Sobolev space. The exact values of these constants
are difficult to estimate, but we note that experimentally all data
distributions under consideration had bounded noise, which translates to no
arbitrarily large spikes in the pseudo-gradient. Assumption 5 (eq. (41)) is
used in analyzing the direction associated with projected pseudo-gradients
(Lemma 1). Next, we shift to presenting our convergence results regarding the
first-order projected mirror descent in RKHS.

###  V-A Accuracy and Complexity Tradeoffs of Algorithm 1

We proceed to establish the convergence of Algorithm 1 under constant
compression budget. To do so, we define the projected pseudo-gradient through
the difference of the optimality condition reformulation of the pseudo-
gradient (13) before and after projection. Specifically, KOMP-based projection
is applied to the auxiliary sequence ${\tilde{z}}_{t+1}$ to obtain the dual
function iterates $z_{t+1}$. Therefore, define the _projected_ pseudo-gradient
${\hat{g}}_{t}$ as

$\displaystyle{\hat{g}}_{t}:=\frac{1}{\eta}(z_{t}-z_{t+1})=\frac{1}{\eta}\left(\nabla\psi(f_{t})-\nabla\psi(f_{t+1})\right).$
(43)

which simplifies the evolution of the projected iterates $z_{t}$ (14):

$\displaystyle z_{t+1}=z_{t}-\eta{\hat{g}}_{t}.$ (44)

Next we define some operations employed in the analysis: ${\mathbb{E}}$
denotes the total expectation taken with respect to the unknown joint
distribution $P({\mathbf{x}})$, and ${\mathbb{E}}[g|{\mathcal{F}}_{t}]$
denotes the expectation of $g$ conditioned on the filtration
${\mathcal{F}}_{t}$. For brevity, we also define:

$\displaystyle{\Gamma}_{t}:={\mathbb{E}}{\left\|\nabla
R(f_{t})\right\|^{2}_{\ast}}.$ (45)

where $f^{\star}$ is defined in (1), $B_{\psi}$ is given in (7), and
$\left\|\cdot\right\|_{\ast}$ denotes the RKHS dual-norm as defined in Sec.
II. We begin by bounding the directional error of the projected pseudo-
gradient (43) with respect to its un-projected variant (13).

######  Lemma 1.

The directional error of the projected pseudo-gradient (43) relative to the
pseudo-gradient [cf. (13)] as quantified by the RKHS dual-norm is bounded by
the ratio of the compression budget to the step-size $\epsilon/\eta$.
Moreover, it has bounded second-moment in the dual-norm. That is,

$\displaystyle\left\|{\hat{g}}_{t}-g_{t}\right\|_{\ast}\leq\frac{\epsilon}{\eta},$
(46) $\displaystyle{\mathbb{E}}\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}\leq
2\left(\\!\left(\\!\frac{\epsilon}{\eta}\right)^{\\!\\!2}+b^{2}+c^{2}{\mathbb{E}}[\langle\nabla
R(f_{t}),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle]\right).$ (47)

With Lemma 1 (proof in Appendix A) established, we now focus on establishing
the mean convergence of Algorithm 1 under suitable selections of $\eta$ and
$\epsilon$.

######  Theorem 1.

Under Assumptions 1-5, upon running Algorithm 1 for $t+1$ iterations, the
objective sub-optimality attenuates linearly up to a bounded neighborhood when
run with constant step-size $\eta<\min(\frac{1}{q_{1}},\frac{q_{1}}{q_{2}})$
and compression budget $\epsilon=\alpha\eta$,

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq(1-\rho)^{t}{\mathbb{E}}[R(f_{0})-R(f^{\ast})]+\frac{1}{\rho}\left[L_{1}\eta^{2}b^{2}+\left(\frac{\eta\omega_{1}}{2}+L_{1}\eta^{2}\right)\alpha^{2}\right],$
(48)

where $\rho=q_{1}\eta-q_{2}\eta^{2}$ with
$q_{1}=2\lambda\left(D-\frac{1}{2\omega_{1}}\right)$ and $q_{2}=2\lambda
DL_{1}c^{2}$, $D$ are positive constants: $D$ is the correlation constant in
Assumption 1, $\omega_{1}$ is a constant of Young’s inequality (aka Peter–Paul
inequality) [51] such that $\omega_{1}>\frac{1}{2D}$, $\alpha$ is the
parsimony constant, $\lambda>0$ is the P–Ł inequality constant of Assumption
3, and positive constant $c^{2}$ coming from Assumption 5.

###### Proof.

This proof generalizes [52][cf. eqn. (2.38)] to the case of function sequences
in RKHS with pseudo-gradients that are observed with deterministic errors. We
begin by recalling the definition of $R_{\psi}$ for the convex function $\psi$
defined by Bregman divergence (7), i.e.,
$R_{\psi}(z)=(R\circ\nabla\psi^{\ast})(z)=R(\nabla\psi^{\ast}(z))$ for
$z\in{\mathcal{H}}_{\ast}$. Consider this quantity evaluated at auxiliary
functions $z_{t+1}=\nabla\psi(f_{t+1})$ and $z_{t}=\nabla\psi(f_{t})$, and
apply Assumption 4 regarding its Lipschitz continuity:

$\displaystyle
R_{\psi}(\nabla\psi(f_{t+1}))-R_{\psi}(\nabla\psi(f_{t}))-\langle\nabla
R_{\psi}(\nabla\psi(f_{t})),\nabla\psi(f_{t+1})-\nabla\psi(f_{t})\rangle\leq\frac{L_{1}}{2}||\nabla\psi(f_{t+1})-\nabla\psi(f_{t})||_{\ast}^{2}.$
(49)

Now, consider the expression for the projected pseudo-gradient in (43), which
may be rearranged to obtain
$\nabla\psi(f_{t+1})-\nabla\psi(f_{t})=-\eta{\hat{g}}_{t}$. Taken together
with the fact that $R_{\psi}(\nabla\psi(f_{t}))=R(f_{t})$ since
$\nabla\psi^{\ast}=(\nabla\psi)^{-1}$ one may rewrite (49) as an approximate
descent relation as

$\displaystyle R(f_{t+1})\leq R(f_{t})-\langle\nabla
R(f_{t}),\eta{\hat{g}}_{t}\rangle+\frac{L_{1}\eta^{2}}{2}||{\hat{g}}_{t}||_{\ast}^{2}$
$\displaystyle\qquad\quad\ =R(f_{t})-\eta\langle\nabla
R(f_{t}),g_{t}\rangle+\eta\langle\nabla
R(f_{t}),g_{t}-{\hat{g}}_{t}\rangle+\frac{L_{1}\eta^{2}}{2}||{\hat{g}}_{t}||_{\ast}^{2}.$
(50)

The right-hand side of this expression decomposes into four terms: (i) the
objective at the previous step, (ii) the product-moment between the gradient
and the pseudo-gradient, (iii) the directional error associated with pseudo-
gradient projections, (iv) a second-moment of the projected pseudo-gradient.

We proceed to the third term of the right-hand side of (50). Via Young’s
inequality with constant $\omega_{1}$ and the KOMP stopping criteria (46), its
expected value is upper-estimated as

$\displaystyle\eta{\mathbb{E}}\left[\nabla
R(f_{t}),g_{t}-{\hat{g}}_{t}\right]\leq\frac{\eta}{2\omega_{1}}{\Gamma}_{t}+\frac{\eta\omega_{1}}{2}\left(\frac{\epsilon}{\eta}\right)^{2},$
(51)

where we employ the short-hand notation (45) for ${\Gamma}_{t}$. (iv) Observe
that the expected value of the fourth term of the right-hand side of (50) can
be bounded as (47) (Lemma 1). Now, we substitute the right-hand sides of (51),
(47) and gather terms:

$\displaystyle\\!\\!{\mathbb{E}}[R(f_{t+1})]\leq{\mathbb{E}}[R(f_{t})]-(\eta\\!-\\!L_{1}c^{2}\eta^{2}){\mathbb{E}}[\langle\nabla
R(f_{t}),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle]+\left(\frac{\eta}{2\omega_{1}}-D\eta\right){\Gamma}_{t}+\left(\frac{\eta\omega_{1}}{2}+L_{1}\eta^{2}\right)\left(\frac{\epsilon}{\eta}\right)^{2}+L_{1}\eta^{2}b^{2}.$
(52)

Proceed by subtracting ${\mathbb{E}}[R(f^{\ast})]$ from both sides, and then
apply (38) (Assumption 1) to the gradient-pseudo-gradient inner-product term
to obtain the expected decrement relation:

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq{\mathbb{E}}[R(f_{t})-R(f^{\ast})]-\left(D\eta-\frac{\eta}{2\omega_{1}}-DL_{1}c^{2}\eta^{2}\right){\Gamma}_{t}+\left(\frac{\eta\omega_{1}}{2}+L_{1}\eta^{2}\right)\alpha^{2}+L_{1}\eta^{2}b^{2},$
(53)

where we have substituted in the choice of compression budget
$\epsilon=\alpha\eta$ for some scalar $\alpha>0$ in order to simplify
fractions of the compression budget to the step-size that appear in the
preceding expressions. Now we use the P–Ł condition as stated in Assumption 3
to the term ${\Gamma}_{t}$ [cf. (45)] in (53) as follows:

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq(1-\rho){\mathbb{E}}[R(f_{t})-R(f^{\ast})]+L_{1}\eta^{2}b^{2}+\left(\frac{\eta\omega_{1}}{2}+L_{1}\eta^{2}\right)\alpha^{2},$
(54)

where constant $\rho=q_{1}\eta-q_{2}\eta^{2}$ determines the transient rate of
convergence, with $q_{1}=2\lambda\left(D-\frac{1}{2\omega_{1}}\right)$ and
$q_{2}=2\lambda DL_{1}c^{2}$. Next, we obtain conditions on constants $D$ and
$\omega_{1}$ such that:

$\displaystyle 0\leq\rho\leq 1.$ (55)

The constant $q_{1}$ is required to be positive to satisfy (55). Imposing this
constraint implies that the Peter Paul inequality constant $\omega_{1}$ in
(51) satisfies $\omega_{1}>\frac{1}{2D}$, where $D$ is the correlation
constant in (38). These conditions together imply

$\displaystyle\eta<\frac{1}{q_{1}\left(1-\frac{q_{2}}{q_{1}}\eta\right)}.$
(56)

Now if $\rho\leq 1$, and if $\left(1-\frac{q_{2}}{q_{1}}\eta\right)\leq 1$,
then $\eta<\frac{1}{q_{1}}$ implies (56) holds. On the contrary, $\rho\geq 0$
implies $\eta<\frac{q_{1}}{q_{2}}$. Overall, then, we obtain the valid step-
size range as $\eta<\min(\frac{1}{q_{1}},\frac{q_{1}}{q_{2}})$.

Returning to (54), we iteratively break down the right-hand together with the
fact that
$\sum\limits_{i=0}^{t}(1-\rho)^{i}\leq\sum\limits_{i=0}^{\infty}(1-\rho)^{i}=\frac{1}{\rho}$
to obtain (48), which concludes the proof. ∎

Theorem 1 characterizes the trade-off between the rate of the convergence and
the radius of the ball to which the iterates converge at steady state. First
note that regardless of the choice of $\eta$ and $t$, the mean distance from
the optimum will always be ${\mathcal{O}}(\alpha^{2})$ in the worst case. The
bound in (48) is for $\epsilon>0$, which causes the additional $\alpha^{2}$ to
appear. For $\epsilon=0$, the $\alpha^{2}$ term of (48) vanishes and hence
simplifies to ${\mathcal{O}}(\eta)$ asymptotically since $\rho$ is
approximately of order $\eta$ for $\eta<1$.

Further observe that the second term on the right-hand side of (48) simplifies
to ${\mathcal{O}}(\eta b^{2}+(1+\eta)\alpha^{2})$, which is in accordance with
the convergence rates of the iterates of stochastic mirror descent for vector-
valued problems [52, 53]. Theorem 1 is a generalization to the RKHS setting,
where we additionally require the range of functions to be nonnegative.
Moreover, we explicitly characterize the error in the convergence behavior
incurred due to subspace projections of the auxilliary sequence $z_{t}$ [cf.
(19)]. Relative to [29][Theorem 6], our convergence result holds under
comparable conditions, but incorporates the additional aspect of trading off
parameterization efficiency and convergence accuracy due to projections. Thus,
for $\epsilon=0$ our result simplifies to the aforementioned result, but
requires a parameterization that grows unbounded with $t$ [cf. (18)].

Parameterization efficiency. We analyze the complexity of the function
parameterization associated with Algorithm 1 when employing sparse projections
defined by KOMP. To do so, we require additional two conditions.

######  Assumption 6.

The pseudo-gradient may be written in the form
$g_{t}=g_{t}^{\prime}\kappa({\mathbf{x}}_{t},\cdot)$ with scalar
$g_{t}^{\prime}$ bounded by constant $C$

$\displaystyle|g_{t}^{\prime}|\leq C.$ (57)

######  Assumption 7.

The feature space $\mathcal{X}$ is compact.

Assumption 6 implies that either the objective is differentiable or there
exists a suitable kernel embedding such that the chain rule is applicable.
Moreover, Assumption 7 ensures that the data domain has finite covering number
[12]. Under these conditions, we establish that the model order of the
function parameterization defined by Algorithm 1 is finite via analogous logic
to [26][Theorem 4], stated here as a corollary.

######  Corollary 2.

Denote as $M_{t}$ the model order of the function $z_{t}$ obtained from
running Algorithm 1 with fixed compression budget $\epsilon>0$. Then we have
that $M_{t}\leq M^{\infty}$, where $M^{\infty}$ is the maximum model order
upper-estimated as

$\displaystyle
M^{\infty}\leq{\mathcal{O}}\left(\frac{1}{\epsilon}\right)^{d}.$ (58)

A detailed proof is given in Appendix B. We note that in contrast to [26],
compactness is only required to establish the finite model order property in
Corollary 2, but not for Theorem 1. By contrast, Assumption 7 is directly
required for convergence in [26]. Moreover, [26] necessitates the use of exact
stochastic gradients, and therefore is inapplicable to the primary focus of
this work, which is the interpolation of the intensity parameter in
inhomogeneous Poisson Processes.

###  V-B Convergence in Mean of Algorithm 3

To analyze Algorithm 3, we require similar conditions employed in the
previously conducted analysis for the first-order scheme. We impose
Assumptions 1 [cf. (39)], 2, 3, 4 and 5 [cf. (42)] and also introduce
assumptions on maximum eigen-value $\mu_{t+1}^{\max}$ of the Hessian matrix
$\mathbf{A}_{\infty,t+1}$ and on the second moment of norm difference between
finite and infinite dictionary Quasi–Newton gradients.

######  Assumption 8.

For finite $t$, the maximum eigenvalue $\mu_{t+1}^{\max}$ of the Hessian
estimate $\mathbf{A}_{\infty,t+1}$ [cf. (28)] is finite.

######  Assumption 9.

The mean-square error between finite and infinite dictionary Quasi–Newton
directions satisfies the growth condition:

$\displaystyle{\mathbb{E}}[\left\|G_{t}^{(\infty)}-G_{t}\right\|^{2}_{\ast}]\leq\nu_{\mathcal{D}}^{2}{\mathbb{E}}[\left\|G_{t}^{(\infty)}\right\|^{2}_{\ast}],$
(59)

for $t\in{\mathbb{N}}$ and some real positive constant $\nu_{{\mathcal{D}}}$
that exponentially decays as the dictionary subspace size $|\mathcal{D}|$
increases.

Observe that the analysis we conduct is for the Quasi–Newton algorithm
developed in Section IV, which employs finite-dimensional truncations. Its un-
truncated counterpart defines an update direction, which we use as a barometer
of the quality of our truncation through Assumptions 8 and 9. In particular,
these conditions quantify the error due to truncation, and allow us to express
algorithm performance in terms of the objective at the true optimizer in the
RKHS, rather than any truncated optimizer.

Assumption 8 enforces that for any finite $t$, the doubly infinite dimensional
Hessian matrix formed by the iterative update (28) with
${\tilde{g}}_{\infty,t}$ does not blow up to infinity in terms of its maximum
eigenvalue $\mu_{t+1}^{\max}$, which determines its operator norm. As noted in
[29] for the case of diminishing step-size, i.e. when $f_{t}\rightarrow
f^{\ast}$, the inner product between the true gradient and conditional
expectation of the first order pseudo-gradient becomes orthogonal. Analogous
reasoning is applicable to the approximate Newton step. Assumption 9
generalizes the eigenvalue decay conditios that are typically imposed on the
solution path of kernel ridge regression, i.e., that the eigenvalues of the
kernel matrix form a summable series [54, 55], which holds when one restricts
the RKHS to a Sobolov space. Here we impose a related condition, which is that
the tails of this series are a finite quantifiable number satisfying
$\nu_{{\mathcal{D}}}\propto\exp(-\tilde{a}|\mathcal{D}|)$, where $\tilde{a}$
is a positive constant.

We begin by employing Assumption 8 to characterize the upper and lower-
boundedness of Hessian approximations (28).

######  Lemma 2.

The inverse of the Hessian matrix generated through the Quasi–Newton updates
is bounded as

$\displaystyle\frac{1}{\mu_{t+1}^{\max}}\mathbf{I}\preccurlyeq\mathbf{A}_{\infty,t+1}^{-1}\preccurlyeq\frac{1}{\delta}\mathbf{I}.$
(60)

Proof of Lemma 2 is given in Appendix C. We now focus on the convergence
analysis of Algorithm 3. For brevity, here also we use
${\Gamma}_{t}:={\mathbb{E}}{\left\|\nabla R(f_{t})\right\|^{2}_{\ast}}$ as
detailed in (45).

######  Theorem 3.

(i) Under Assumptions 1-5, 8 and 9 and for finite $t$, upon selecting a
_constant step-size_ $\eta<\min\left(\frac{1}{2\lambda
q_{3}},\frac{q_{3}}{q_{4}}\right)$ and regularizer
$\delta>\frac{c^{\prime}\nu_{{\mathcal{D}}}\mu_{t+1}^{\max}}{D}$, the average
sub-optimality satisfies

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq(1-\rho)^{t}{\mathbb{E}}[R(f_{t})-R(f^{\ast})]+\frac{{b^{\prime}}^{2}L_{1}(1+\nu_{{\mathcal{D}}}^{2})}{\rho\delta^{2}}\eta^{2}+\frac{{b^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2\rho\delta^{2}}\eta,$
(61)

where $\rho_{t}=2\lambda q_{3}(\eta_{t}-\frac{q_{4}}{q_{3}}\eta_{t}^{2})$ with
$q_{3}=\frac{D}{\mu_{t+1}^{\max}}-\frac{D{c^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2\delta^{2}}-\frac{1}{2\omega_{2}}$
and
$q_{4}=\frac{D{c^{\prime}}^{2}}{\delta^{2}}L_{1}(1+\nu_{{\mathcal{D}}}^{2})$.
The positive constants $D$, $\lambda$, $L_{1}$, ${b^{\prime}}^{2}$,
${c^{\prime}}^{2}$, $\mu_{t+1}^{\max}$ and $\nu_{{\mathcal{D}}}^{2}$ are
defined in Assumptions 1, 3, 4, 5, 8 and 9 respectively, and the choice of
Young’s inequality constant $\omega_{2}$ is given in (72).
(ii) Under the same conditions, if we select a _diminishing step-size_
$\eta_{t}=\min\left(\frac{q_{3}}{2q_{4}},\frac{(2t+1)}{\lambda
q_{3}(t+1)^{2}}\right)$, the average sub-optimality converges as

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq\frac{4L_{1}{b^{\prime}}^{2}(1+\nu_{{\mathcal{D}}}^{2})}{\delta^{2}\lambda^{2}q_{3}^{2}(t+1)}+\frac{{b^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2\lambda\delta^{2}q_{3}}.$
(62)

###### Proof.

The proof follows analogous logic to that which yields Theorem 1. Begin by
employing Assumption 4 to write

$\displaystyle R_{\psi}(z_{t+1})\\!-\\!R_{\psi}(z_{t})\\!-\\!\langle\nabla
R_{\psi}(z_{t}),z_{t+1}\\!-\\!z_{t}\rangle\\!\leq\\!\frac{L_{1}}{2}\\!\left\|z_{t+1}\\!-\\!z_{t}\right\|^{2}_{\ast}.$
(63)

Now, substitute (30) in the definition of $z_{t+1}$ associated with the
Quasi–Newton update and add and subtract the infinite dimensional second order
gradient $G_{t}^{(\infty)}$ [cf (36)] with $z_{t+1}$ term to write

$\displaystyle R_{\psi}(z_{t+1})-R_{\psi}(z_{t})\leq-\eta_{t}\langle\nabla
R_{\psi}(z_{t}),G_{t}^{(\infty)}\rangle+L_{1}\eta_{t}^{2}\left\|G_{t}^{(\infty)}\right\|^{2}_{\ast}+\eta_{t}\langle\nabla
R_{\psi}(z_{t}),G_{t}^{(\infty)}-G_{t}\rangle+L_{1}\eta_{t}^{2}\left\|G_{t}^{(\infty)}-G_{t}\right\|^{2}_{\ast},$
(64)

where we have used the inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$ to split the
last term. We first focus on the expectation of the third term on the right
hand side of (64) which can be bounded using Young’s inequality of constant
$\omega_{2}$ as

$\displaystyle\eta_{t}{\mathbb{E}}\langle\nabla
R_{\psi}(z_{t}),G_{t}^{(\infty)}\\!-\\!G_{t}\rangle\\!\leq\\!\frac{\eta_{t}}{2\omega_{2}}{\Gamma}_{t}\\!+\\!\frac{\eta_{t}\omega_{2}}{2}{\mathbb{E}}\left\|G_{t}^{(\infty)}\\!-\\!G_{t}\right\|^{2}_{\ast},$
(65)

where ${\Gamma}_{t}$ is given in (45) and noting that $\nabla
R_{\psi}(z_{t})=\nabla R(f_{t})$. The last terms of (65) and (64) can be
grouped together and use of assumption 9 further simplifies expectation of
(64) as

$\displaystyle{\mathbb{E}}[R_{\psi}(z_{t+1})-R_{\psi}(z_{t})]\leq-\eta_{t}{\mathbb{E}}\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}[G_{t}^{(\infty)}|{\mathcal{F}}_{t}]\rangle+\frac{\eta_{t}}{2\omega_{2}}{\Gamma}_{t}+\left(\frac{\eta_{t}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2}+L_{1}(1+\nu_{{\mathcal{D}}}^{2})\eta_{t}^{2}\right){\mathbb{E}}\left\|G_{t}^{(\infty)}\right\|^{2}_{\ast}.$
(66)

Next focusing on the third term on the right hand side of (66), the term
$\left\|G_{t}^{(\infty)}\right\|^{2}_{\ast}$ [cf (36)] can be factorized as

$\displaystyle\\!\\!\\!\left\|G_{t}^{(\infty)}\right\|^{2}_{\ast}={\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)^{\top}\mathbf{A}_{\infty,t+1}^{-1}{\tilde{g}}_{\infty,t}{\tilde{g}}_{\infty,t}^{\top}\mathbf{A}_{\infty,t+1}^{-1}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot).$
(67)

Using the definition
${\tilde{g}}_{\infty,t}{\tilde{g}}_{\infty,t}^{\top}=[{\tilde{g}}_{\infty,t}^{\prime}]^{2}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}({\mathbf{x}}_{t}){\mathbf{k}}_{{\mathcal{D}}_{\infty}}({\mathbf{x}}_{t})^{\top}$,
which is non-negative, along with using the upper-bound on the Hessian inverse
estimate $\mathbf{A}_{t+1}^{-1}$ in Lemma 2, we may simplify the right-hand
side of (67) as
$\|G_{t}^{(\infty)}\|_{\ast}^{2}\leq\frac{1}{\delta^{2}}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)^{\top}{\tilde{g}}_{\infty,t}{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)=\frac{1}{\delta^{2}}\left\|{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)\right\|^{2}_{\ast}$.
Using Assumption 5, its expected value can be bounded as

$\displaystyle\\!\\!{\mathbb{E}}\\!\left\|G_{t}^{(\infty)}\right\|^{2}_{\ast}\\!\leq\frac{{b^{\prime}}^{2}}{\delta^{2}}\\!+\\!\frac{{c^{\prime}}^{2}}{\delta^{2}}{\mathbb{E}}[\langle\nabla
R_{\psi}(z_{t}\\!),{\mathbb{E}}[{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}\\!(\cdot\\!)|{\mathcal{F}}_{t}]\rangle].\\!$
(68)

Now, let us return to the first-term on the right hand side of (66). Apply
Lemma 2 to replace $\mathbf{A}_{t+1}^{-1}$ by its upper-bound in terms of
identity and $\mu_{t+1}^{\max}$:

$\displaystyle-\eta_{t}{\mathbb{E}}\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}[G_{t}^{(\infty)}|{\mathcal{F}}_{t}]\rangle\leq-\frac{\eta_{t}}{\mu_{t+1}^{\max}}{\mathbb{E}}\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}[{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)|{\mathcal{F}}_{t}]\rangle.$
(69)

So now, via (68) and (69), the expected value of the objective difference
across time in (66) is bounded as

$\displaystyle{\mathbb{E}}[R_{\psi}(z_{t+1})-R_{\psi}(z_{t})]\leq$
$\displaystyle\frac{\eta_{t}}{2\omega_{2}}{\Gamma}_{t}+\frac{{b^{\prime}}^{2}}{\delta^{2}}\left[\frac{\eta_{t}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2}+L_{1}(1+\nu_{{\mathcal{D}}}^{2})\eta_{t}^{2}\right]$
$\displaystyle-\left[\frac{\eta_{t}}{\mu_{t+1}^{\max}}-\frac{{c^{\prime}}^{2}}{\delta^{2}}\left(\frac{\eta_{t}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2}+L_{1}(1+\nu_{{\mathcal{D}}}^{2})\eta_{t}^{2}\right)\right]\times{\mathbb{E}}\langle\nabla
R_{\psi}(z_{t}),{\mathbb{E}}[{\tilde{g}}_{\infty,t}^{\top}{\mathbf{k}}_{{\mathcal{D}}_{\infty}}(\cdot)|{\mathcal{F}}_{t}]\rangle$
(70)

The third term on the right-hand side of (70) can be bounded via the pseudo-
gradient inequality [cf. (39)] of Assumption 1, after which we gather like
terms. Then, by substituting $\nabla R_{\psi}(z_{t})=\nabla R(f_{t})$ and
subtracting $R(f^{\ast})$ on both sides of (70), and using (45), we may write

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq{\mathbb{E}}[R(f_{t})-R(f^{\ast})]-(q_{3}\eta_{t}-q_{4}\eta_{t}^{2}){\Gamma}_{t}+\frac{{b^{\prime}}^{2}}{\delta^{2}}\left[\frac{\eta_{t}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2}+L_{1}(1+\nu_{{\mathcal{D}}}^{2})\eta_{t}^{2}\right],$
(71)

where
$q_{3}=\frac{D}{\mu_{t+1}^{\max}}-\frac{D{c^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2\delta^{2}}-\frac{1}{2\omega_{2}}$
and
$q_{4}=\frac{D{c^{\prime}}^{2}}{\delta^{2}}L_{1}(1+\nu_{{\mathcal{D}}}^{2})$.
The quantity $q_{3}$ has to be positive, which yields a quadratic in terms of
$\omega_{2}$. Via the quadratic formula, we have:

$\displaystyle\frac{\frac{D}{\mu_{t+1}^{\max}}\\!-\\!\sqrt{\frac{D^{2}}{(\mu_{t+1}^{\max})^{2}}\\!-\\!\frac{{c^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}}{\delta^{2}}}}{\frac{{c^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}}{\delta^{2}}}<\omega_{2}<\frac{\frac{D}{\mu_{t+1}^{\max}}\\!+\\!\sqrt{\frac{D^{2}}{(\mu_{t+1}^{\max})^{2}}\\!-\\!\frac{{c^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}}{\delta^{2}}}}{\frac{{c^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}}{\delta^{2}}}.$
(72)

The lower bound on $\omega_{2}$ is always positive and from the discriminant
term we get the regularizer selection as
$\delta\geq\frac{c^{\prime}\nu_{{\mathcal{D}}}\mu_{t+1}^{\max}}{D}$.
Assumption 3 applied to the second term of the right-hand side of (71)
followed by grouping terms yields

$\displaystyle{\mathbb{E}}[R(f_{t+1})$
$\displaystyle-R(f^{\ast})]\leq(1-\rho_{t}){\mathbb{E}}[R(f_{t})-R(f^{\ast})]+\frac{{b^{\prime}}^{2}L_{1}(1+\nu_{{\mathcal{D}}}^{2})}{\delta^{2}}\eta_{t}^{2}+\frac{{b^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}}{2\delta^{2}}\eta_{t},$
(73)

where $\rho_{t}=2\lambda q_{3}(\eta_{t}-\frac{q_{4}}{q_{3}}\eta_{t}^{2})$ with
$q_{3}$ and $q_{4}$ defined after (71).
Constant step-size: Set the step-size as constant $\eta_{t}=\eta$ and hence
$\rho_{t}$ is time-invariant: $\rho_{t}=\rho$. To ensure (73) defines a
contraction up to an error bound same condition as (55) shall be satisfied for
$\rho$. Following the similar analysis for Algorithm 1 step-size selection
criterion in (54) - (56), here we get the condition on step-size as

$\displaystyle\eta<\min\left(\frac{1}{2\lambda
q_{3}},\frac{q_{3}}{q_{4}}\right).$ (74)

Note that the maximum eigen value $\mu_{t+1}^{\max}$ (inside $q_{3}$) of
$\mathbf{A}_{\infty,t+1}$, obtained by addition of rank-1 positive outer
products of pseudo-gradient vectors $\tilde{g}_{\infty,t}$, is the maximum
eigenvalue of Hessian estimates up to the present time $t$, and increases
slowly as the iterate progresses. Hence our convergence results pertain to
finite $t$. Returning to (73) with the step-size selection (74), we may
simplify via recursive substitution of the first term on the right-hand side
of (73) yielding (61). We know that $\rho={\mathcal{O}}(\eta)$. Hence the mean
distance from optima depends on
${\mathcal{O}}(\eta)+{\mathcal{O}}(\frac{\nu_{{\mathcal{D}}}^{2}\omega_{2}}{\delta^{2}})$,
where the last term with fixed bias came as a result of approximating the
infinite dimensional dictionary in terms of finite dimension ${\mathcal{D}}$.

Diminishing step-size: For this case, we choose
$\eta_{t}=\min\left(\frac{q_{3}}{2q_{4}},\frac{(2t+1)}{\lambda
q_{3}(t+1)^{2}}\right)$. For such a step-size choice, the factor
$\eta_{t}-\frac{q_{4}}{q_{3}}\eta_{t}^{2}$ in $\rho_{t}$ of (73) can be lower
bounded by

$\displaystyle\eta_{t}-\frac{q_{4}}{q_{3}}\eta_{t}^{2}\geq\eta_{t}-\frac{q_{4}}{q_{3}}\times\frac{q_{3}}{2q_{4}}\eta_{t}=\frac{\eta_{t}}{2}.$
(75)

Hence the lower bound of (75) allows us to lower bound $\rho_{t}$ with
$\lambda q_{3}\eta_{t}$. Using this lower bound on $\rho_{t}$ finally in (73)
and using the fact that $\eta_{t}=\frac{(2t+1)}{\lambda q_{3}(t+1)^{2}}$ for
large $t$, we get finally

$\displaystyle{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq\frac{t^{2}}{(t+1)^{2}}{\mathbb{E}}[R(f_{t})-R(f^{\ast})]+\frac{L_{1}{b^{\prime}}^{2}(1+\nu_{{\mathcal{D}}}^{2})(2t+1)^{2}}{\delta^{2}\lambda^{2}q_{3}^{2}(t+1)^{4}}+\frac{{b^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}(2t+1)}{2\lambda\delta^{2}q_{3}(t+1)^{2}}.$
(76)

Multiplying both sides of (76) with $(t+1)^{2}$ and taking
$\beta(t)=t^{2}{\mathbb{E}}[R(f_{t})-R(f^{\ast})]$ and using the fact that
$\frac{2t+1}{t+1}<2$ for all $t$, (76) simplifies as

$\displaystyle\beta(t+1)-\beta(t)\leq\frac{4L_{1}{b^{\prime}}^{2}(1+\nu_{{\mathcal{D}}}^{2})}{\delta^{2}\lambda^{2}q_{3}^{2}}+\frac{{b^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}(2t+1)}{2\lambda\delta^{2}q_{3}}.$
(77)

Doing a telescopic sum of (77) from $t=0$ to $t$ and noting that $\beta(0)=0$
yields

$\displaystyle(t+1)^{2}{\mathbb{E}}[R(f_{t+1})-R(f^{\ast})]\leq\frac{4L_{1}{b^{\prime}}^{2}(1+\nu_{{\mathcal{D}}}^{2})(t+1)}{\delta^{2}\lambda^{2}q_{3}^{2}}+\frac{{b^{\prime}}^{2}\nu_{{\mathcal{D}}}^{2}\omega_{2}(t+1)^{2}}{2\lambda\delta^{2}q_{3}}.$
(78)

Dividing both sides of (78) by $(t+1)^{2}$ yields (62). ∎

The convergence result of (61) regarding Algorithm 3 is similar to its first-
order counterpart, but deals with the additional technicalities associated
with the second-order information incorporated into the pseudo-gradients as
preconditioning. Moreover, due to the use of a fixed subspace size
$|\mathcal{D}|$, there is no need to discuss projection-induced errors in
update directions. In addition, for Algorithm 3, we achieve a diminishing rate
of convergence of order
${\mathcal{O}}(\frac{1+\nu_{{\mathcal{D}}}^{2}}{\delta^{2}q_{3}^{2}(t+1)})$
plus an additional bias
${\mathcal{O}}(\frac{\nu_{{\mathcal{D}}}^{2}\omega_{2}}{\delta^{2}q_{3}})$ due
to approximating the subspace with finite dictionary ${\mathcal{D}}$, for
diminishing step-size selection as given in (62). We acknowledge that faster
local rate analyses of Quasi–Newton methods exist (such as superlinear [56]
and quadratic [39][Ch. 6]), but our priority was global rate analysis and
analyzing the evolution in the RKHS. Thus, we defer sharpening the rate
analysis to future work.

##  VI Experiments

We shift to experimental validation of Algorithms 1, 3 and 4 as compared with
state of the art second order offline and first order online benchmarks for
normalized intensity estimation in inhomogeneous Poisson Point Processes (PPP,
Example 1). Throughout, we select as our pseudo-gradient the kernel embedding
– see Secs. III and IV.

We employ KOMP in SPPPOT with both constant (same as [26]) and adaptive budget
selection. Constant budgeting requires trial and error to yield a stabilizing
model order, whereas for adaptive budget, we tune $\alpha$ defined by budget
specification $\epsilon=\alpha\eta$ at time $t$ as

$\displaystyle\alpha(t+1)\\!=\\!\alpha(t)\\!\left[1\\!+\\!\max\left(\min\left((M_{t}\\!-\\!d_{mo}\\!)\\!\times\\!0.1\%,0.1\\!\right)\\!,\\!-\\!0.1\\!\right)\right],$

where $\alpha(t+1)$ has an enforced upper and lower bound to ensure it remains
nonzero even when near the target dictionary size $d_{mo}$. The intuition is
that this constant increases when the model order $M_{t}$ starts increasing
from $d_{mo}$ so that we allow a bit higher error tolerance budget
$\epsilon_{t+1}=\alpha(t+1)\eta$, that results in reduction of $M_{t+1}$. The
contrary is true when $M_{t}<d_{mo}$.

Algorithm |  Epochs |  a |  Kernel B.W. |  $\eta$ |  Mini-Batch Size |  Budget |  $\delta$ |  $\lvert{\mathcal{U}}\rvert$ |  Model order |  Time
---|---|---|---|---|---|---|---|---|---|---
1-D Gaussian toy example
BFGS |  1 |  0.15 |  0.0065 |  – |  – |  – |  – |  100 |  10211 |  $\sim$ 8 hrs
PMD |  8 |  – |  0.0065 |  0.05 |  30 |  – |  – |  100 |  100 |  $\sim$ 5 mins
POLK |  8 |  – |  0.0065 |  0.006 |  30 |  $\epsilon=3.5\times 10^{-8}$ |  – |  100 |  179 |  $\sim$ 5 mins
Algorithm 1 ($\epsilon$ fixed) |  8 |  – |  0.0065 |  0.012 |  30 |  $\epsilon=6.6\times 10^{-6}$ |  – |  100 |  100 |  $\sim$ 5 mins
Algorithm 1 ($\epsilon_{t}$ adaptive) |  8 |  – |  0.0065 |  0.012 |  30 |  $\alpha(0)=2\times 10^{-6}$, $d_{mo}=105$ |  – |  100 |  105 |  $\sim$ 5 mins
Algorithm 3 |  8 |  – |  0.0065 |  1.25 |  1 |  – |  1 |  100 |  100 |  $\sim$ 20 mins
Algorithm 3 |  8 |  – |  0.0065 |  1 |  1 |  – |  1 |  50 |  50 |  $\sim$ 20 mins
Algorithm 4 |  8 |  – |  0.0065 |  SPPPOT: 0.1, Quasi–Newton: 1.8182 |  SPPPOT: 30, Quasi–Newton: 1 |  $\alpha(0)=1.5\times 10^{-6}$, $d_{mo}=90$ |  1 |  50 |  90 |  $\sim$ 20 mins
Stephen Curry Data
BFGS |  1 |  1 |  0.0025 |  – |  – |  – |  – |  100 |  8298 |  $\sim$ 8 hrs
PMD |  10 |  – |  0.0025 |  0.1 |  30 |  – |  – |  100 |  100 |  $\sim$ 5 mins
Algorithm 1 |  10 |  – |  0.0025 |  0.03 |  30 |  $\epsilon=10^{-5}$ |  – |  100 |  100 |  $\sim$ 5 mins
Algorithm 3 |  10 |  – |  0.0025 |  0.2381 |  1 |  – |  1 |  100 |  100 |  $\sim$ 20 mins
Algorithm 3 |  10 |  – |  0.0025 |  0.2381 |  1 |  – |  1 |  50 |  50 |  $\sim$ 20 mins
Algorithm 4 |  10 |  – |  0.0025 |  SPPPOT: 0.03, Quasi–Newton: 0.1408 |  SPPPOT: 30, Quasi–Newton: 1 |  $\epsilon=10^{-5}$ |  1 |  100 |  100 |  $\sim$ 20 mins
Chicago Crime Spatial Data
BFGS |  1 |  1 |  0.01 |  – |  – |  – |  – |  441 |  10000 |  $\sim$ 12 hrs
PMD |  4 |  – |  0.01 |  0.8 |  30 |  – |  – |  441 |  441 |  $\sim$ 5 mins
Algorithm 1 |  4 |  – |  0.01 |  0.15 |  30 |  $\epsilon=3\times 10^{-4}$ |  – |  441 |  441 |  $\sim$ 5 mins
Algorithm 3 |  4 |  – |  0.01 |  6.6667 |  1 |  – |  1 |  441 |  441 |  $\sim$ 20 mins
Algorithm 4 |  4 |  – |  0.01 |  SPPPOT: 0.4, Quasi–Newton: 5 |  SPPPOT: 30, Quasi–Newton: 1 |  $\epsilon=3\times 10^{-4}$ |  1 |  441 |  441 |  $\sim$ 20 mins
TABLE I: Parameters selected for the various experiments.

For PPP pdf estimation (Example 1) of function $f$, the negative log-
likelihood (24) (identical to [34, Eq. (3.1)]) defines the loss, whose
minimization over time may be executed via Algorithms 1, 3 and 4. For
comparison, batch offline BFGS (Quasi–Newton) minimization method [34] and
also online Pseudo Mirror Descent (PMD) [29, Eq. (8)] are implemented, where
the later uses only grid points with no dictionary adaptation. We also note
that the PMD loss function [29, Eq. (8)] requires the optimal intensity
$f^{\ast}$ term, which is unknown for the real world data. BFGS is included in
’fminunc’ package of MATLAB with intensity given by
$f(\cdot)=af^{\prime}(\cdot)^{2}$ where $a$ is a positive scalar used for
tuning. $f(\cdot)$ is the actual intensity function expressed via
$f^{\prime}(\cdot)^{2}$ to preserve positivity. We also compare with the
online benchmark of POLK [26], i.e., functional stochastic gradient method
with the RKHS norm as the Bregman divergence, combined with KOMP in the primal
RKHS, not the dual.

For all online schemes that use mini-batching, dividing the cumulative
gradient by the mini-batch size is needed to remove scale ambiguity of
integral approximations in the pseudo-gradient (25). The function learnt in
BFGS is scaled by the number of data points. Thus, the function learnt for all
the schemes is the normalized intensity. To obtain the intensity, the density
must be multiplied by the sample size. For all the experiments we use Gaussian
kernel, whose bandwidth is chosen using via Silverman’s rule of thumb [57] for
synthetic 1-D Gaussian toy example, and by trial and error for the real world
data sets. The rest of the system parameters for different experiments are
chosen by trial and error – see Table I.

1-D Gaussian toy example. We consider 10211 training and 1001 test data points
generated from a non-homogeneous PPP model having Gaussian density
$f(x)=\frac{10}{\sqrt{2\pi}}\exp\left(50(x-0.5)^{2}\right)$. All data points
$x\in(0,1)$ are univariate and have restricted range.

Results. The normalized intensity estimate obtained through different
algorithms and their corresponding model order, log of RMSE and the test loss
evaluation with respect to the number of samples processed are given in Figs.
1(a), 1(d), 1(b) and 1(c) respectively. From the intensity plot in Fig. 1(a),
we observe that BFGS, PMD, Algorithms 1, 3 and 4 yield densities near the
ground truth, whereas POLK, with high test loss: $2.0148$ and test RMSE:
$5.8317$, yields density estimates that become negative at a few places. This
underscores the importance of using the KL-Divergence in Algorithms 1, 3 and 4
to preserve positivity. Algorithm 4 with model order $90$ ($|\mathcal{U}|=50$
uniform grid points and $40$ Poisson points) outperforms alternatives
including Algorithm 3 since application of SPPPOT in Algorithm 4 allows it to
discern the optimal subspace for the Quasi–Newton updates. Overall, then,
SPPPOT yields a favorable accuracy/complexity tradeoff among first order
algorithms whereas Algorithms 3 and 4 achieve superior performance.

Stephen Curry Data. We experiment with the NBA dataset of Stephen Curry which
has shot distances from the basket up to $40$ feet normalized as samples
$x\in\mathbb{R}$. These samples are split into $8298/1000$ points for
training/testing.

Results. The training loss, the estimated normalized intensity, and the model
order are given in Figs. 2(a), 2(b), and 2(c), respectively. The runtime
difference between BFGS and the other online methods (given in Table I) is
reflected in the model complexity difference in 2(c). Among first order
algorithms, SPPPOT outperforms the batch solver after a few training epochs
(Fig. 2(a)), and yields an intensity much closer to the offline baseline BFGS
compared to PMD. Algorithms 3 and 4 outperforms the existing state of the art
schemes and the offline BFGS by a healthy margin.

Chicago Crime Spatial Data. Chicago 2018 crime dataset has 2D features (X-Y
Coordinates) determining the crime location in the city. $10,000$ training
points and $2000$ test points are randomly selected and normalized to lie
within a unit square.

Results. Fig. 3(a) shows the histogram of the $10,000$ training data
(normalized) with $30\times 30$ bins in the square grid. The histogram scale
is not to be confused with the learned density since histogram is dependent on
the number of data points and the size of each bin whereas the later is not.
Fig. 3(b) is negative at someplaces since density aka pdf is used in the
negative log-likelihood. Observe that the precision and accuracy of Algorithms
3 and 4 in terms of test loss evaluation (Fig. 3(b)) and determination of the
peaks/troughs of the normalized intensity with respect to the ground truth
(Figs. 3(g) and 3(h)) helps it achieve the state of the art. The run time
complexity (given in Table I) is justified from Fig. 3(c).

(a) Fitted Normalized intensity

(b) log of RMSE vs. samples

(c) Test Loss vs. samples

(d) Model Complexity vs. samples

Figure 1: We showed the performance of Algorithms 1, 3 and 4 in comparison to
PMD and offline BFGS algorithm on the 1-D synthetic Gaussian toy example. POLK
is also implemented with constant budget to emphasize that it violates the
range constraint. Online second order Algorithms 3 and 4 achieves the state of
the art.

(a) Test Loss vs. samples processed

(b) Fitted normalized Poisson Intensity

(c) Model Complexity vs. samples processed

Figure 2: We compare our Algorithms 1, 3 and 4 with BFGS solver for an offline
MLE problem [34] as well as with online Pseudo-mirror descent (PMD) [29],
i.e., which do not incorporate complexity-reducing projections on the NBA
dataset of Stephen Curry. Online second order Algorithms 3 and 4 is much more
efficient and achieves the state of the art.

(a) Histogram of training data

(b) Test Loss vs. samples

(c) Model Complexity vs. samples

(d) Pdf using offline BFGS

(e) Pdf using PMD

(f) Pdf using Algorithm 1

(g) Pdf using Algorithm 3

(h) Pdf using Algorithm 4

Figure 3: We compare the performance of our first order Algorithm 1, second
order Algorithm 3 and their combination Algorithm 4 with PMD and offline BFGS
on Chicago crime data set of 2018. Second order Algorithms 3 and 4 beats the
state of the art algorithms in normalized intensity estimation.

##  VII Conclusion

We studied expected value minimization when the decision variable belongs to a
Reproducing Kernel Hilbert Space (RKHS) whose target domain is required to be
nonnegative. We developed efficient algorithms using first and second order
gradient information. For the first order algorithm, we put forth a variant of
stochastic mirror descent that employs (i) _pseudo-gradients_ and (ii)
projections. Compressive projections are executed via KOMP, which mitigate the
complexity issues of RKHS parameterizations. For the second order algorithm,
we fixed our underlying subspace and used second order _pseudo-gradients_ to
make a descent in the weights of the auxiliary function. We established
accuracy/complexity tradeoffs between sub-optimality convergence in mean and
bounds on the model complexity under standard assumptions. Experiments
demonstrated competitive performance with the state of the art for Poisson
intensity estimation. Future directions include scaling the intensity
estimation approach to higher dimensions through convolutional kernels, and
employing the resulting intensity estimates for event triggers in queueing
systems.

##  Appendix A Proof of Lemma 1

###### Proof.

This result mirrors the proof of [26][Prop. 8] but instead studies the
directional error of pseudo-gradients with respect to the dual sequence
$\\{z_{t}\\}$ (14). First, consider the dual-norm difference between
${\hat{g}}_{t}$ [cf. (43)] and $g_{t}$ [cf. (14)], to get

$\displaystyle\left\|{\hat{g}}_{t}-g_{t}\right\|_{\ast}$
$\displaystyle=\|\frac{1}{\eta}(z_{t}-z_{t+1})-g_{t}\|_{\ast}$
$\displaystyle=\frac{1}{\eta}\left\|z_{t+1}-{\tilde{z}}_{t+1}\right\|_{\ast}\leq\frac{\epsilon}{\eta},$
(79)

where we have grouped terms and used the definition of $\tilde{z}_{t+1}$ in
(14), together with the KOMP stopping criterion
$\left\|z_{t+1}-{\tilde{z}}_{t+1}\right\|_{\ast}\leq\epsilon$. To obtain (47),
consider ${\mathbb{E}}\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}$, and add and
subtract $g_{t}$ inside the expectation to write

$\displaystyle{\mathbb{E}}\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}\leq
2{\mathbb{E}}\left\|g_{t}\right\|^{2}_{\ast}+2{\mathbb{E}}\left\|{\hat{g}}_{t}-g_{t}\right\|^{2}_{\ast}$
(80) $\displaystyle\qquad\leq 2({b^{2}}+c^{2}{\mathbb{E}}[\langle\nabla
R(f_{t}),{\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]\rangle])+2{\mathbb{E}}\left\|{\hat{g}}_{t}-g_{t}\right\|^{2}_{\ast},$

where we have applied $(a+b)^{2}\leq 2a^{2}+2b^{2}$ to the result of this
subtraction, followed by Assumption 5. ∎

##  Appendix B Proof of Corollary 2

###### Proof.

Begin by recalling the definition of the function approximation error of KOMP
(Algorithm 2)

$\gamma_{j}:=\left\|{\tilde{z}}-\sum\limits_{{\mathbf{x}}_{n}\in{\mathcal{D}}\setminus\\{{\mathbf{x}}_{j}\\}}w_{n}\kappa({\mathbf{x}}_{n},.)\right\|_{\ast}$

used to determine the compression stopping criteria. $M_{t}$ is the model
order of $z_{t}$ and the model order of ${\tilde{z}}_{t+1}$ is $M_{t}+1$ for
addition of single data point. Then, $\gamma_{M_{t}+1}$ may be written as

$\displaystyle\gamma_{M_{t}+1}$
$\displaystyle=\left\|{\tilde{z}}_{t+1}-\sum\limits_{{\mathbf{x}}_{n}\in{\mathcal{D}}_{t+1}\setminus\\{{\mathbf{x}}_{M_{t}+1}\\}}w_{n}\kappa({\mathbf{x}}_{n},.)\right\|_{\ast}.$
(81)

Substitute in the definition ${\tilde{z}}_{t+1}=z_{t}-\eta g_{t}$ [cf. (14)],
and note that dictionary elements are deleted if
$\min\limits_{j=1,\dots,M_{t+1}}\gamma_{j}\leq\epsilon$. We note that
$\gamma_{j}$ for ${j=1,\dots,M_{t+1}}$ is lower-bounded by $\gamma_{M_{t}+1}$,
so it suffices to consider $\gamma_{M_{t}+1}$. Since the pseudo-gradient
admits an expansion via the chain rule as
$g_{t}=g_{t}^{\prime}\kappa({\mathbf{x}}_{t},\cdot)$ (Assumption 6), one may
reduce this condition to the set distance between the dual of the subspace
defined by the current dictionary and the kernel evaluated at the latest
sample $\mathbf{x}_{t}$ by applying similar logic to [26, Appendix D.1 Eq.
(70)-(75) and Lemma 10], that is, the following holds

$\displaystyle\text{dist}\left(\kappa(\mathbf{x}_{t},\cdot),\mathcal{H}_{{\mathcal{D}}_{t}}^{\ast}\right)\leq\frac{\epsilon}{\eta|g_{t}^{\prime}|}=\frac{\alpha}{|g_{t}^{\prime}|}.$
(82)

There are two possibilities: either a new point is added or it is not. If it
is not, then $M_{t+1}\leq M_{t}+1$ and the model order does not change. Let’s
consider when this condition is violated, such that the model order increases,
i.e., $M_{t+1}=M_{t}+1$:

$\displaystyle\text{dist}\left(\kappa(\mathbf{x}_{t},\cdot),\mathcal{H}_{{\mathcal{D}}_{t}}^{\ast}\right)>\frac{\epsilon}{\eta|g_{t}^{\prime}|}=\frac{\alpha}{|g_{t}^{\prime}|}.$
(83)

Assumption 6 implies $\frac{1}{|g_{t}^{\prime}|}\geq\frac{1}{C}$, which allows
us to lower bound (83) with $\frac{\alpha}{C}$, where
$\alpha=\frac{\epsilon}{\eta}$. KOMP stopping criteria is violated when
distinct dictionary points $\mathbf{d}_{n}$ and $\mathbf{d}_{j}$ for
$j,n\in\\{1,...,M_{t}\\}$ satisfy the condition
$\left\|\phi(\mathbf{d}_{n})-\phi(\mathbf{d}_{j})\right\|_{\ast}>\frac{\epsilon}{\eta
C}=\frac{\alpha}{C}$. By Assumption 7, the number of Euclidean balls of radius
$\frac{\epsilon}{\eta C}$ required to cover the feature kernelized space
$\phi({\mathbf{x}})=\kappa({\mathbf{x}},\cdot)$ is always finite by the
continuity of the kernel. Therefore, there exists a finite $t$ such that (82)
holds which implies that equation
$\min\limits_{j=1,\dots,M_{t+1}}\gamma_{j}\leq\epsilon$ is always valid,
meaning the model order does not grow. More specifically, there exists a
finite maximum model order $M^{\infty}<\infty$ such that $M_{t}<M^{\infty}$,
holds for all $t$.

Now if the kernel taken is Lipschitz continuous on the set ${\mathcal{X}}$
which is compact by Assumption 7, then similar to [58, Proposition 2.2], the
maximum model order $M^{\infty}$ is bounded as (58), which yields an upper-
bound on $M^{\infty}$ for all $t$. ∎

##  Appendix C Proof of Lemma 2

###### Proof.

The Hessian matrix $\mathbf{A}_{\infty,t+1}$ has finite maximum eigenvalue
$\mu_{t+1}^{\max}$ via Assumption 8, which implies
$\frac{1}{\mu_{t+1}^{\max}}$ is the minimum eigenvalue of
$\mathbf{A}_{\infty,t+1}^{-1}$. Thus we may write the lower bound of (60). The
upper-bound employs the update (28) with ${\tilde{g}}_{\infty,t}$ recursively
backwards in time as

$\displaystyle\mathbf{A}_{\infty,t+1}^{-1}=\left(\delta\mathbf{I}+\sum_{j=1}^{t}{\tilde{g}}_{\infty,j}{{\tilde{g}}_{\infty,j}}^{\top}\right)^{-1}\preccurlyeq\frac{1}{\delta}\mathbf{I}.$
(84)

∎

##  Appendix D Functional Pseudo-Dual Averaging

Dual averaging (see [59]) is closely related to mirror descent. Here, we
expand upon these connections. Specifically, dual averaging with pseudo-
gradients in the functional RKHS setting takes the form

$\displaystyle{\tilde{h}}_{t+1}$ $\displaystyle=h_{t}+g_{t}$ (85)
$\displaystyle f_{t+1}$
$\displaystyle=\operatorname*{arg\,min}\limits_{f\in{\mathcal{H}}}\Big{(}\langle
h_{t+1},f\rangle_{\mathcal{H}}+\frac{1}{\eta}\psi(f)\Big{)}$
$\displaystyle=\nabla\psi^{\star}(-\eta h_{t+1})$ (86)

and then one may apply KOMP on ${\tilde{h}}_{t+1}$ with budget parameter
$\epsilon$ in order to ensure an efficient RKHS parameterization in terms of
weights and dictionary elements. Here, $\psi^{\star}$ denotes the conjugate of
$\psi$. In order for (86) to be valid, one requires
$\nabla\psi^{\star}:{\mathcal{H}}^{\ast}\rightarrow{\mathcal{H}}$ so that
$f_{t+1}\in{\mathcal{H}}$. For the stochastic case (a special case of pseudo
gradient), the update can be written recursively as

$\displaystyle{\tilde{h}}_{t+1}=h_{t}+\nabla r_{t}(\nabla\psi^{\star}(-\eta
h_{t}))=h_{t}-\eta{\tilde{g}}_{t}$ (87)

where ${\tilde{g}}_{t}:=-\frac{1}{\eta}\nabla r_{t}(\nabla\psi^{\star}(-\eta
h_{t}))$. To see the link between the (85)-(86) and Algorithm 1, multiply (85)
by $-\eta$ and use the second update equation to obtain

$\displaystyle\nabla\psi({\tilde{f}}_{t+1})=\nabla\psi(f_{t})-\eta g_{t}$ (88)

where ${\tilde{f}}_{t+1}:=\nabla\psi^{\star}(-\eta{\tilde{h}}_{t+1})$. This
means that under an additional hypothesis that the gradient of the Fenchel
conjugate of $\psi$ does not take the functions $h_{t}$ outside the RKHS,
update expressions for dual averaging in updating $h_{t}$ exactly coincides
with $z_{t}$ in (14). This link provides alternative potential ways of
approaching the convergence theory of Algorithm 1, which are the subject of
future work.

##  Appendix E Extensions of SPPPOT to supervised learning

The first order algorithm SPPPOT can be easily extended to high dimensional
supervised learning problems owing to its data adaptive dictionary and can
give excellent performances compared to the current state of the art
techniques. Let us first elucidate here the problem of supervised learning and
then shall demonstrate the example of multi-class kernel logistic regression
as an use case.

Supervised Learning. For supervised learning, target variables
$y_{t}\in{\mathcal{Y}}$ are also received corresponding to each
${\mathbf{x}}_{t}$, where $y_{t}\in{\mathcal{Y}}$ denote the labels or real
values, i.e., ${\mathcal{Y}}=\\{1,\dots,C\\}$ in the case of classification
and ${\mathcal{Y}}\subset{\mathbb{R}}$ in the case of regression. The examples
$({\mathbf{x}}_{t},y_{t})$ are realizations from an unknown distribution
$\mathbb{P}({\mathbf{x}},y)$. The goal is to fit a predictive model
$f:{\mathcal{X}}\rightarrow{\mathcal{Y}}$ that belongs to a particular
hypothesized class of functions ${\mathcal{H}}$, i.e., to form estimates of
the target of the form $\hat{y}_{t}=f({\mathbf{x}}_{t})$. The merit of a given
estimator $f$ is quantified by a convex loss
$\ell:\mathcal{H}\times{\mathcal{X}}\rightarrow{\mathcal{Y}}$ which is small
when $f({\mathbf{x}})$ and $y$ are close. We seek to minimize the statistical
loss $R(f):=\mathbb{E}[\ell(f({\mathbf{x}}),y)]$ in expectation over the data.
Example include the squared loss $(y-f({\mathbf{x}}))^{2}$, the hinge loss
$\max\\{0,1-yf({\mathbf{x}})\\}$, and the logistic loss derived from
minimizing the class-conditional misclassification probability – see [60]. For
this case also, we solve the same optimization problem as (1), where the
instantaneous cost $r_{t}(f):=\ell(f({\mathbf{x}}_{t}),y_{t})$.

Next we present the demonstrative example of multi-class Kernel Logistic
Regression (KLR) as a representative of supervised learning.

######  Example 4.

Kernel Logistic regression for multi-class classification:  In this case, the
target domain ${\mathcal{Y}}=\\{1,\dots,C\\}$ is the set of classes, and the
loss is defined by the negative log-likelihood of a generalized odds ratio.
Specifically, define a vector-function associated with a one-hot encoding,
i.e., the target domain is $\\{0,1\\}^{C}$, such that each ${\mathbf{x}}$ is
assigned to a binary vector of length $C$. Then, one hypothesizes that the
probability of the predicted label $y$ to lie in the class $c$ is given by the
logit (softmax) function
${\mathbb{P}}(y=c|{\mathbf{x}}):=\frac{\exp(f_{c}({\mathbf{x}}))}{\sum_{\tilde{c}}f_{\tilde{c}}({\mathbf{x}})}$,
which gives rise to the negative log-likelihood: [60]

$\displaystyle\ell(\textbf{f},{\mathbf{x}}_{t},y_{t})$
$\displaystyle=-\log({\mathbb{P}}(y=y_{t}|{\mathbf{x}}_{t}))$ (89)
$\displaystyle=\log\left(\sum_{\tilde{c}}f_{\tilde{c}}({\mathbf{x}}_{t})\right)-f_{y_{t}}({\mathbf{x}}_{t})$
(90)

where f denotes vector of functions $f_{\tilde{c}}$ where $\tilde{c}$ denotes
each of the class labels, and the classification of a point ${\mathbf{x}}$ is
defined by the maximum a posteriori assignment:
$\hat{c}=\operatorname*{arg\,max}_{c\in\\{1,\dots,C\\}}f_{c}({\mathbf{x}})$.
Observe that the space of binary sequences $\\{0,1\\}^{C}$ is nonnegative,
meaning that minimizing (90) in expectation over $\mathbb{P}({\mathbf{x}},y)$
is an instance of (1).

First order gradient calculation for KLR.  Consider the setting of multi-class
classification using a probabilistic generalization of the odds ratio which
defines the loss (90). We differentiate it with respect to each
$f_{\tilde{c}}$ to obtain the functional stochastic gradient for each class
where $\tilde{c}\in\\{1,\ldots,C\\}$. The gradient expression takes the form:

$\displaystyle
g_{t}=\begin{cases}\frac{\exp\left(f_{c^{\prime}}({\mathbf{x}}_{t})\right)}{\sum_{\tilde{c}}\exp\left(f_{\tilde{c}}({\mathbf{x}}_{t})\right)}\kappa({\mathbf{x}}_{t},\cdot)&c^{\prime}\neq
y_{t}\\\
\left(\frac{\exp\left(f_{y_{t}}({\mathbf{x}}_{t})\right)}{\sum_{\tilde{c}}\exp\left(f_{\tilde{c}}({\mathbf{x}}_{t})\right)}-1\right)\kappa({\mathbf{x}}_{t},\cdot)&c^{\prime}=y_{t}\end{cases}$
(91)

The parametric updates for the entire function vector f are carried out via
(91), which operates upon the basis of a stream of samples
$\\{{\mathbf{x}}_{t},y_{t}\\}$ in tandem with subspace projections. Here also
pseudo-gradient property holds since simple stochastic gradients are used.

In Appendix F, we will numerically experiment on multi-class MNIST handwritten
dataset using the KLR formulation. There, we will also show the merit of
SPPPOT by using the Bregman divergence as KL-divergence compared to its
squared RKHS norm difference variant that is nothing but POLK [26].

##  Appendix F First-order Kernel logistic regression on multi-class MNIST

(a) Test Loss vs. samples processed

(b) Test Error vs. samples processed

(c) Model Complexity vs. samples processed

Figure 4: We compared SPPPOT algorithm (both constant and adaptive budgets)
with POLK [26] and offline FGD over MNIST dataset for solving KLR. SPPPOT
obtains competitive accuracy with LIBSVM package, although it is computable
online.

We consider the problem of Kernel Logistic Regression (KLR) (Example 4), whose
loss is (90) and stochastic gradient takes the form (91). For SPPPOT, we have
taken KL divergence as the Bregman divergence and considered both constant and
adaptive parsimony constants. We have compared SPPPOT with the online
algorithm POLK [26] (without the regularizer term) and with the offline batch
algorithm Functional Gradient Descent (FGD) without any dictionary
compression. We ran FGD for 1000 iterations. We have also compared with the
multi-class kernel SVM (KSVM) solver in LIBSVM package [61], which is also a
batch algorithm. Even though it uses the KSVM algorithm, not KLR and the loss
function is also different, still it is a classification algorithm and is used
as a baseline to plot the classification error at the end which is shown as a
dotted baseline in Fig. 4(b).

Dataset. We consider algorithm performance for the problem of classifying
MNIST handwritten digits [35], which consists of 60000 data points for
training and 10000 data points for testing. Out of 10000 training points, we
randomly sampled 2000 points for our training purpose. Each data point is an
image of a handwritten digit, and the goal is to label each images into the
corresponding digit. We simply read the images into a $784$-dimensional
feature vector and scaled the pixel values ranging from 0 to 255 to the
intensities in $[0,1]$.

Parameter selection. We select the Gaussian kernel with bandwidth selected by
cross-validation over $3.33\%$ of the training data (randomly chosen $2000$
points) and the KOMP budget is kept as zero during cross-validation for both
POLK and SPPPOT, which yielded bandwidth is taken as $16$ (standard deviation
$4$) for all experiments. A single epoch of training data is fed to both FGD
and LIBSVM (C-SVC option) whereas $2$ epochs are being run for the online
algorithms POLK and SPPPOT (both constant and adaptive budget). The constant
learning rate is taken by trial and error to be $24$ for FGD, POLK and SPPPOT
for both budget combinations. The KOMP budget is taken as $\epsilon=1.4579$
for SPPPOT with constant budget and $\epsilon=1.8812$ for POLK, which yield
converged model orders, respectively, of $401$ and $404$. Matching the model
orders for both the experiments is hard and here adaptive budgeting becomes
handy. For SPPPOT with adaptive budget, we selected the desired model order as
$400$ and the constant $\alpha(0)=0.1$. It is seen in Fig. 4(c) that the mean
model order is $400$ and the model order fluctuations with respect to the
samples processed are all bounded and eventually converge. The storage
requirement is finite and comparable to POLK with model order $401$.

Results. The test loss, error and the model order plots with respect to
samples processed has been presented in Figs. 4(a), 4(b) and 4(c)
respectively. Observe that SPPPOT with both constant and adaptive budgets
outperforms POLK and obtains performance competitive with the offline
benchmarks FGD and LIBSVM (C-SVC) whose parameterization complexity is $60000$
and $16122$, respectively, which are not depicted in the same figure due to
scale. SPPPOT, after the loss and error settles down, performs comparably to
FGD with 1000 iterates, whereas for POLK, there is still a considerable gap.
FGD took about $12$ hours to compute a single epoch with $1000$ iterations. As
seen in Fig. 4(b), KVSM has better accuracy than the asymptotic performance of
KLR algorithms. Still we can say the test error for SPPPOT is within $2\%$ of
the LIBSVM solver and yields a more efficient parameterization.

##  Appendix G Almost sure convergence of SPPPOT for stochastic gradients

Stochastic gradients are special case of pseudo-gradients which are available
whenever the instantaneous cost is evaluable in closed form as in, for
instance, kernel logistic regression [cf. (91)]. Note that stochastic
gradients satisfy the stochasticity rule, i.e.
${\mathbb{E}}[g_{t}|{\mathcal{F}}_{t}]=\nabla R(f_{t})$. To prove the almost
sure convergence of RKHS-valued stochastic mirror descent with projections,
only Assumption 2 and the next three conditions are required.

######  Assumption 10.

The function $R_{\psi}(\cdot)$ is $\lambda_{1}$-strongly convex.

######  Assumption 11.

The function $\psi$ is $1$-strongly convex.

######  Assumption 12.

The projected stochastic gradient has bounded second moment given the
filtration $\mathcal{F}_{t}=\sigma(\\{{\mathbf{x}}_{i}\\}_{i=1}^{t-1})$.

$\displaystyle\mathbb{E}[\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}|{\mathcal{F}}_{t}]\leq
G^{2}$ (92)

Assumption 11 is required to lower bound the norm difference between the
auxiliary iterates from the optima, i.e.
$\left\|z_{t}-z^{\ast}\right\|_{\ast}$ by the the norm difference between the
actual function iterates and the optima, i.e. $\left\|f_{t}-f^{\ast}\right\|$.
Note that Assumption 10 can be enforced explicitly by adding a regularizer
term in the dual space with respect to the auxiliary function $z_{t}$.
Assumptions 10 and 12 are standard in the literature and is required to show
the almost sure convergence theorem for the stochastic case.

######  Theorem 4.

Under Assumption 2 and Assumptions 10-12, upon running Algorithm 1 with
constant step-size $\eta$ and constant compression budget $\epsilon$, the
iterates $f_{t}$ converges almost surely to a neighbourhood in limit infrimum
as $t\rightarrow\infty$.

$\displaystyle\liminf_{t\rightarrow\infty}\left\|f_{t}-f^{\ast}\right\|\leq\xi:=\frac{\epsilon+\sqrt{\epsilon^{2}+\eta^{3}\lambda_{1}G^{2}}}{2\eta\lambda_{1}}$
(93)

where $G$ is defined in (92) as the bound on the conditional second-moment of
the stochastic gradient in the dual norm and $\lambda_{1}$ is the strong-
convexity parameter.

###### Proof.

The proof starts by taking the dual Hilbert norm difference between $z_{t+1}$
and $z^{\ast}=\nabla\psi(f^{\ast})$, where $z^{\ast}$ is taken to be the
transformed function of the optimum $f^{\ast}$ in the dual space, and then
expanding $z_{t+1}$ same as (44) as shown below

$\displaystyle\left\|z_{t+1}-z^{\ast}\right\|^{2}_{\ast}$
$\displaystyle=\left\|z_{t}-\eta{\hat{g}}_{t}-z^{\ast}\right\|^{2}_{\ast}$
(94) $\displaystyle=\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}-2\eta\langle
z_{t}-z^{\ast},{\hat{g}}_{t}\rangle+\eta^{2}\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}$
(95) $\displaystyle=\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}-2\eta\langle
z_{t}-z^{\ast},g_{t}\rangle+2\eta\langle
z_{t}-z^{\ast},g_{t}-{\hat{g}}_{t}\rangle+\eta^{2}\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}\;,$
(96)

where the last equation is obtained by addition and subtraction of the actual
stochastic gradient $g_{t}$. The third term on the right hand side of (96) can
be upper bounded using Cauchy-Schwartz inequality and then application of (46)
on it yields the bound

$\displaystyle\left\|z_{t+1}-z^{\ast}\right\|^{2}_{\ast}\leq$
$\displaystyle\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}-2\eta\langle
z_{t}-z^{\ast},g_{t}\rangle+2\epsilon\left\|z_{t}-z^{\ast}\right\|_{\ast}+\eta^{2}\left\|{\hat{g}}_{t}\right\|^{2}_{\ast}\;.$
(97)

So now taking conditional expectation given the past sigma algebra
${\mathcal{F}}_{t}$ and using Assumption 12, we get

$\displaystyle\mathbb{E}[\left\|z_{t+1}-z^{\ast}\right\|^{2}_{\ast}|{\mathcal{F}}_{t}]$
$\displaystyle\leq\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}-2\eta\langle
z_{t}-z^{\ast},\nabla
R(f_{t})\rangle+2\epsilon\left\|z_{t}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}\;.$
(98)

We assumed $R_{\psi}(\cdot)$ to be the risk function which takes the auxiliary
function $z_{t}$ as an argument such that
$R_{\psi}(z_{t})=R_{\psi}(\nabla\psi(f_{t}))=R(f_{t})$, where $R(\cdot)$ is
our actual risk function. Hence, using this property and convexity of
$R_{\psi}$ and then Assumption 10, the second term on the right hand side of
(98) can be simplified as

$\displaystyle\langle z_{t}-z^{\ast},\nabla R(f_{t})\rangle$
$\displaystyle=\langle z_{t}-z^{\ast},\nabla R_{\psi}(z_{t})\rangle$ (99)
$\displaystyle\geq R_{\psi}(z_{t})-R_{\psi}(z^{\ast})$ (100)
$\displaystyle\geq\frac{\lambda_{1}}{2}\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}\;.$
(101)

So now using (101) in (98), the inequality further simplifies to

$\displaystyle\mathbb{E}[\left\|z_{t+1}-z^{\ast}\right\|^{2}_{\ast}|{\mathcal{F}}_{t}]$
$\displaystyle\leq(1-\eta\lambda_{1})\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{t}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}\;.$
(102)

A stopping stochastic process can be constructed using inequality (102), that
tracks when $\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}$ reaches the sub-optimal
threshold. When the threshold is achieved, the function converges to the
neighbourhood of the optimum. In order to establish this fact, we construct a
stochastic process $\delta_{t}$ such that the supermartingale condition is
satisfied, i.e. $\mathbb{E}[\delta_{t+1}|{\mathcal{F}}_{t}]\leq\delta_{t}$. We
would like to find out the threshold from (102) by finding out the roots for
$\left\|z_{t}-z^{\ast}\right\|$ when the following holds true

$\displaystyle\mathbb{E}[\left\|z_{t+1}-z^{\ast}\right\|^{2}_{\ast}|{\mathcal{F}}_{t}]\leq\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}$
(103)
$\displaystyle(1-\eta\lambda_{1})\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{t}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}\leq\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}\;.$
(104)

(104) simplifies to the following quadratic inequality

$\displaystyle-\eta\lambda_{1}\left\|z_{t}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{t}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}\leq
0\;.$ (105)

The roots of the polynomial are

$\displaystyle\left\|z_{t}-z^{\ast}\right\|_{\ast}$
$\displaystyle=\frac{-2\epsilon\pm\sqrt{4\epsilon^{2}+4(\eta\lambda_{1})(\eta^{2}G^{2})}}{-4\eta\lambda_{1}}$
(106)
$\displaystyle=\frac{\epsilon\mp\sqrt{\epsilon^{2}+\eta^{3}\lambda_{1}G^{2}}}{2\eta\lambda_{1}}\;.$
(107)

The inequality (105) opens downwards and since
$\left\|z_{t}-z^{\ast}\right\|_{\ast}$ is always positive, hence only positive
root will be taken. So the radius of the convergence can be represented as

$\displaystyle\xi:=\frac{\epsilon+\sqrt{\epsilon^{2}+\eta^{3}\lambda_{1}G^{2}}}{2\eta\lambda_{1}}\;.$
(108)

So using (108), let’s now define the stopping stochastic process $\delta_{t}$
as

$\displaystyle\delta_{t}=\left\|z_{t}-z^{\ast}\right\|_{\ast}{\mathbb{I}}\Big{\\{}\min_{u\leq
t}-\eta\lambda_{1}$
$\displaystyle\left\|z_{u}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{u}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}>\xi\Big{\\}}\;,$
(109)

where ${\mathbb{I}}$ denotes the indicator function. Note $\delta_{t}\geq 0$
since both the norm and the indicator are non-negative. So when the indicator
condition $\min_{u\leq
t}-\eta\lambda_{1}\left\|z_{u}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{u}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}>\xi$
is satisfied, we may compute the square root of (103) yielding

$\displaystyle\mathbb{E}[\delta_{t+1}|{\mathcal{F}}_{t}]\leq\delta_{t}.$ (110)

Now taking the contrary condition, i.e., $\min_{u\leq
t}-\eta\lambda_{1}\left\|z_{u}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{u}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}\leq\xi$,
makes the indicator to be $0$ for the subsequent iterates, since min is used
inside the indicator. So, in either case, (110) holds, which implies that
$\delta_{t}\xrightarrow{a.s.}0$ as $t\rightarrow\infty$. This further implies
the fact that either
$\lim_{t\rightarrow\infty}\left\|z_{t}-z^{\ast}\right\|_{\ast}-\xi=0$ or the
indicator goes to null for large $t$ almost surely, i.e.
$\lim_{t\rightarrow\infty}{\mathbb{I}}\Big{\\{}\min_{u\leq
t}-\eta\lambda_{1}\left\|z_{u}-z^{\ast}\right\|^{2}_{\ast}+2\epsilon\left\|z_{u}-z^{\ast}\right\|_{\ast}+\eta^{2}G^{2}>\xi\Big{\\}}=0$.
Hence any of these conditions terminating to null translates to our
convergence in the neighbourhood

$\displaystyle\liminf_{t\rightarrow\infty}\left\|z_{t}-z^{\ast}\right\|_{\ast}\leq\xi\;.$
(111)

Now using 11, we can write

$\displaystyle\left\|z_{t}-z^{\ast}\right\|_{\ast}=\left\|\nabla\psi(f_{t})-\nabla\psi(f^{\ast})\right\|_{\ast}\geq\left\|f_{t}-f^{\ast}\right\|\;.$
(112)

Using (112) in (111) concludes the proof of theorem 4. ∎

Note that the almost sure convergence established for stochastic gradients in
theorem 4 is a much stronger bound than that of the pseudo-gradient
convergence result given in theorem 1. Hence when stochastic gradient
computation is feasible, as in Kernel logistic regression for multi-class
classification, on running SPPPOT we can almost surely reach to the
neighborhood of the optima.

## References

  * [1] A. Chakraborty, K. Rajawat, and A. Koppel, “Projected pseudo-mirror descent in reproducing kernel hilbert space,” in _2021 55th Asilomar Conference on Signals, Systems, and Computers_ , 2021, pp. 1008–1012.
  * [2] G. Francis, L. Ott, and F. Ramos, “Stochastic functional gradient for motion planning in continuous occupancy maps,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 3778–3785.
  * [3] L. C. Drazek, “Intensity estimation for poisson processes,” _The University of Leeds, School of Mathematics_ , 2013.
  * [4] A. Berlinet and C. Thomas-Agnan, _Reproducing kernel Hilbert spaces in probability and statistics_. Springer Science & Business Media, 2011.
  * [5] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [6] E. Dall’Anese, A. Simonetto, S. Becker, and L. Madden, “Optimization and learning with information streams: Time-varying algorithms and applications,” _IEEE Signal Processing Magazine_ , vol. 37, no. 3, pp. 71–83, 2020.
  * [7] K. Slavakis, G. B. Giannakis, and G. Mateos, “Modeling and optimization for big data analytics:(statistical) learning tools for our era of data deluge,” _IEEE Signal Processing Magazine_ , vol. 31, no. 5, pp. 18–31, 2014.
  * [8] A. Shapiro, D. Dentcheva, and A. Ruszczyński, _Lectures on stochastic programming: modeling and theory_. SIAM, 2014.
  * [9] R. T. Rockafellar and R. J.-B. Wets, _Variational analysis_. Springer Science & Business Media, 2009, vol. 317\.
  * [10] A. Beck and M. Teboulle, “Mirror descent and nonlinear projected subgradient methods for convex optimization,” _Operations Research Letters_ , vol. 31, no. 3, pp. 167–175, 2003.
  * [11] V. Kůrková, “Kolmogorov’s theorem and multilayer neural networks,” _Neural networks_ , vol. 5, no. 3, pp. 501–506, 1992.
  * [12] M. Anthony and P. L. Bartlett, _Neural network learning: Theoretical foundations_. cambridge university press, 2009.
  * [13] T. Hofmann, B. Schölkopf, and A. J. Smola, “Kernel methods in machine learning,” _The annals of statistics_ , pp. 1171–1220, 2008.
  * [14] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [15] D. Yu, G. Hinton, N. Morgan, J.-T. Chien, and S. Sagayama, “Introduction to the special section on deep learning for speech and language processing,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 20, no. 1, pp. 4–6, 2011.
  * [16] J. Mairal, “End-to-end kernel learning with supervised convolutional kernel networks,” in _Advances in neural information processing systems_ , 2016, pp. 1399–1407.
  * [17] J. Kivinen, A. J. Smola, and R. C. Williamson, “Online learning with kernels,” _IEEE transactions on signal processing_ , vol. 52, no. 8, pp. 2165–2176, 2004.
  * [18] Y. Lei and D.-X. Zhou, “Convergence of online mirror descent,” _Applied and Computational Harmonic Analysis_ , vol. 48, no. 1, pp. 343–373, 2020.
  * [19] G. Kimeldorf and G. Wahba, “Some results on tchebycheffian spline functions,” _Journal of mathematical analysis and applications_ , vol. 33, no. 1, pp. 82–95, 1971.
  * [20] C. Williams and M. Seeger, “Using the nyström method to speed up kernel machines,” in _Proceedings of the 14th annual conference on neural information processing systems_ , no. CONF, 2001, pp. 682–688.
  * [21] C. Richard, J. Carlos, M. Bermudez, and P. Honeine, “Online prediction of time series data with kernels,” _IEEE Trans. Signal Process._ , vol. 57, no. 3, pp. 1058–1067, 2009.
  * [22] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” in _Advances in neural information processing systems_ , 2008, pp. 1177–1184.
  * [23] Y. Yang, M. Pilanci, M. J. Wainwright _et al._ , “Randomized sketches for kernels: Fast and optimal nonparametric regression,” _The Annals of Statistics_ , vol. 45, no. 3, pp. 991–1023, 2017.
  * [24] Z. Wang, K. Crammer, and S. Vucetic, “Breaking the curse of kernelization: Budgeted stochastic gradient descent for large-scale svm training,” _Journal of Machine Learning Research_ , vol. 13, no. Oct, pp. 3103–3131, 2012.
  * [25] P. Vincent and Y. Bengio, “Kernel matching pursuit,” _Machine Learning_ , vol. 48, pp. 165––187, 2002.