# Multifractal eigenfunctions
for quantum star graphs

Jonathan P. Keating and Henrik Ueberschär  Sorbonne Université, Université de
Paris, CNRS, IMJ-PRG, F-75006 Paris, France<EMAIL_ADDRESS>Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe
Observatory Quarter, Woodstock Road, Oxford OX2 6GG, UK
<EMAIL_ADDRESS>
###### Abstract.

We prove that the eigenfunctions of quantum star graphs exhibit multifractal
self-similar structure in certain specified circumstances. In the
semiclassical regime, when the spectral parameter and the number of vertices
tend to infinity, we derive an asymptotic condition for the Mellin transform
of a specified function arising from the set of bond lengths which yields an
asymptotic for the Renyi entropy associated with an eigenfunction. We apply
this result to show that one may construct simple quantum star graphs which
satisfy a multifractal scaling law. In the low frequency regime we prove
multifractality by computing the Renyi entropy in terms of a zeta function
associated with the set of bond lengths. In certain arithmetic cases the
fractal exponent $D_{q}$ satisfies a symmetry relation around $q=1/4$ which
arises from the functional equation of the zeta function. Our results are, in
some sense, analogous to the multifractal scaling law that the authors
recently proved for arithmetic S̆eba billiards. However, unlike in that case,
we do not require arithmetic conditions to be satisfied, and nor do we rely on
delicate arithmetic estimates.

##  1\. Introduction

For the most part, the mathematical literature on quantum chaos has focused on
topics such as quantum (unique) ergodicity and scarring of eigenfunctions [2,
13, 38, 24, 36, 14, 15, 16, 17]; however, it is believed that some important
classes of quantum systems possess eigenfunctions whose morphology is far more
complex than being simply ergodic or localized. For example, in systems whose
spectral statistics is intermediate between random-matrix and Poisson there is
a longstanding expectation that the eigenfunctions should exhibit multifractal
self-similar structure.

One such class of intermediate systems are those termed pseudointegrable [28].
Their classical dynamics is close to integrable, yet their quantum dynamics
displays features similar to, but distinct from, that usually observed for
chaotic systems [10]. Examples are rational polygonal billiards, in which the
billiard flow can be studied via the geodesic flow on flat surfaces with
conical singularities. Often included in this class are toy models of such
systems, such as S̆eba billiards [32, 12], quantum star graphs [5, 6, 7, 8]
and parabolic toral automorphisms [25].

Multifractality is ubiquitous in many areas of science (biology, financial
mathematics, quantum physics, etc). It describes self-similarity that is so
complex it cannot be captured by a single fractal exponent. Instead, it can
only be characterised by a spectrum of exponents. While multifractality is an
important concept in the theoretical and experimental study of intermediate
quantum systems, and evidence for its existence has been found in a wide range
of examples [1, 21, 3, 4, 9], there are relatively few rigorous mathematical
results relating to it in that context.

We recently proved [19] the existence of multifractal eigenfunctions for S̆eba
billiards, namely rectangular billiards with a Dirac delta potential. To the
best of our knowledge, this is the first rigorous proof of quantum
multifractality, at least in the context of quantum chaos. In particular, we
obtained a multifractal scaling law for the Renyi entropy in the semiclassical
limit for the case of arithmetic tori, as well as for the ground state in the
case of general rectangular tori. Moreover, we discovered that a symmetry
relation predicted in the physics literature for the fractal exponent
manifests itself, in the case of the ground state, as the functional equation
satisfied by Epstein’s zeta function.

In the present article we deal with another toy model of pseudointegrable
dynamics: quantum star graphs. These systems have been extensively studied in
the mathematical literature on quantum chaos [5, 6, 7, 8]. Star graphs are
graphs with a central vertex that is connected to $v$ outlying vertices. One
can define a self-adjoint operator by introducing the $1d$ Laplacian on each
bond and imposing suitable boundary conditions at the central vertex [20] (cf.
section 3.1). Taking the semiclassical limit, as the eigenvalue tends to
infinity, and at the same time $v\to+\infty$, one can show that the spectral
equation and expressions for the eigenfunctions are of a form which closely
resembles the analogous expressions for S̆eba billiards. Indeed the
statistical properties of the eigenvalues and eigenfunctions in the two
systems are closely related [5, 23].

While our proof of multifractality for S̆eba billiards in the semiclassical
regime was limited to arithmetic billiards and relied heavily on arithmetic
methods, we are able here to obtain analogous results for a generic class of
quantum star graphs. We construct examples which give rise to a multifractal
scaling law for the Renyi entropy. In the case of the ground state regime we
derive a general multifractal scaling law and compute the fractal exponent in
terms of a zeta function defined in terms of the set of bond lengths. In
certain star graphs, the symmetry relation satisfied by the fractal exponent
manifests itself via the functional equation for the associated zeta function.

Acknowledgements. JPK is pleased to acknowledge support from ERC Advanced
Grant 740900 (LogCorRM). H.U. was partially supported by the grant
ANR-17-CE40-0011-01 of the French National Research Agency ANR (project
SpInQS) and by a “délégation CNRS”.

##  2\. Background and main results

We will discuss the spectral properties of quantum graphs in more detail in
section 3. Here we simply introduce what is necessary to formulate our
results. Let us consider a star graph with bond lengths $L_{1},\cdots,L_{v}$.

The spectrum of the self-adjoint operator defined on the star graph forms a
discrete set

$0\leq k_{0}^{2}<k_{1}^{2}\leq k_{2}^{2}\leq\cdots\leq
k_{j}^{2}\leq\cdots\to+\infty$

where the parameters $\pm k_{j}$ are solutions of a spectral equation. are
solutions of the spectral equation

(2.1) $Z(k,L):=\sum_{j=1}^{v}\tan(kL_{j})=\frac{1}{\alpha k}$

where $\alpha$ denotes a coupling parameter. One may consider a renormalized
coupling regime, where $\alpha=\alpha(k)$ may depend on the spectral parameter
$k$.

In particular, the $k_{j}$ interlace with the poles of $Z(k,L)$ which are
given by

$s_{m,j}:=(m+\frac{1}{2})\pi/L_{j},\quad m\in{\mathbb{Z}};\;j=1,\cdots,v.$

Moreover, we have the ordering

$s_{m,1}\geq s_{m,2}\geq\cdots\geq s_{m,v}.$

Let us assume $L_{j}=L+\epsilon(v)\ell_{j}$, where
$\ell_{1}\leq\ell_{2}\leq\cdots\leq\ell_{v}$. Then we have

$s_{m,j}=(m+\frac{1}{2})\frac{\pi}{L}+O(\frac{m\epsilon\ell_{j}}{L^{2}})=(m+\frac{1}{2})\frac{\pi}{L}+O(\frac{m\epsilon\ell_{v}}{L^{2}}).$

In particular, one has

$|s_{m,j}-s_{m,k}|\asymp\frac{m\epsilon}{L^{2}}|\ell_{j}-\ell_{k}|.$

We will consider the semiclassical limit $k\to+\infty$ and at the same time we
will take $v\to+\infty$ in such a way that $k_{n}\epsilon(v)\ell_{v}\to 0$.
This ensures $s_{m,j}=(m+1/2)\pi/L+o(1)$, and because the $k_{j}$ interlace
with the poles of the spectral function it follows that they form clusters
around the values $(m+1/2)\pi/L$ which are well separated. So with each
$m\in{\mathbb{Z}}$ we may associate a cluster ${\mathcal{C}}_{m}$.

Moreover, there are precise formulae for the eigenfunctions, and for the
maximum amplitude squared on each bond. We, thus, obtain a probability measure
on the finite set $\Omega_{v}=\\{1,\cdots,v\\}$ given by (cf. subsection 3.2)

(2.2) $\mu_{k}(i)=\frac{\sec(kL_{i})^{2}}{\sum_{j=1}^{v}\sec(kL_{j})^{2}}.$

The Renyi entropy associated with the probability measure $\mu$ on
$\Omega_{v}$ is then given by

(2.3)
$H_{q}(\mu_{k})=\frac{1}{1-q}\log\left(\sum_{j=1}^{v}\mu_{k}(j)^{q}\right).$

Let us define the Mellin transform of a compactly supported piecewise
continuous function $f$ as

$({\mathcal{M}}f)(s)=\int_{0}^{\infty}t^{-s}f(t)dt.$

Let $\mathfrak{j}\in C^{1}({\mathbb{R}}_{+})$ be an invertible function which
satisfies $\mathfrak{j}(\ell_{j})=j$, $\forall j=1,\dots,v$. We introduce the
length function

$F(v,\mathfrak{j})(t)=(\mathbbm{1}_{[\ell_{i+2},\ell_{v}]}\mathfrak{j}^{\prime})(t+\ell_{i+1})$

We have the following result which provides a condition in terms of the Mellin
transform of the length function to ensure that an asympotic of the Renyi
entropy may be obtained.

######  Theorem 2.1.

If

$\lim_{v\to+\infty}\frac{\log({\mathcal{M}}F(v+1,\mathfrak{j}))(s)}{\log({\mathcal{M}}F(v,\mathfrak{j}))(s)}=1,$

then the asymptotics of the Renyi entropy may be computed as

(2.4)
$H_{q}(\mu_{k})\sim\frac{\log({\mathcal{M}}F(v,\mathfrak{j}))(2q)-q\log({\mathcal{M}}F(v,\mathfrak{j}))(2)}{1-q}.$

We may apply the theorem to deduce the following corollary which shows that
one may construct simple quantum star graphs with a multifractal scaling law.

######  Corollary 2.2.

Let $\mathfrak{j}(t)=t^{n}$, $n\in{\mathbb{N}}$. Denote by $k_{m}$ the first
eigenvalue inside the cluster ${\mathcal{C}}_{m}$ for $m\in{\mathbb{N}}$ and
by $\sigma_{m}$ the distance to the nearest pole. Assume that

$\liminf_{m\to\infty}\sigma_{m}>0.$

For $q>n/2$ and $n\geq 2$, we may compute the fractal exponent as

(2.5) $D_{q}=\lim_{k\to+\infty}\frac{H_{q}(\mu_{k})}{\log
v}=\begin{cases}1-\frac{2q-n}{n(q-1)},\quad\text{if}\;n\geq 3,\\\ \\\
0,\quad\text{if}\;n=2.\end{cases}$

Next we state our results in the ground state regime. These are in some sense
analogous to the results we obtained for S̆eba’s billiard, however, a wide
variety of zeta functions may appear. Note that in the regime of small
eigenvalues the width of the clusters does not blow up with $k$, and therefore
we do not have to take $v\to+\infty$ at the same time, as we vary the spectral
parameter. Rather, let us suppose that $\epsilon(v)\ell_{v}$ is small but
fixed. Then we are in the clustering regime and we may compute the Renyi
entropy for any fixed $v$.

######  Theorem 2.3.

Let us denote by $k_{0}^{2}$ the lowest eigenvalue associated with the cluster
${\mathcal{C}}_{0}$ and we write $k_{0}=\pi/(2L)+\sigma_{0}$. For $v$ fixed,
we have

(2.6) $D_{q}(v)=\lim_{\sigma_{0}\to
0}H_{q}(\mu_{k_{0}})=\frac{\log\zeta_{G,v}(2q)-q\log\zeta_{G,v}(2)}{1-q}$

where

(2.7)
$\zeta_{G,v}(s)=\sum_{j=1}^{v}\left(\frac{\ell_{j}}{L+\epsilon(v)\ell_{j}}\right)^{-s}$

If $\operatorname{Re}q>1/2$, and $\epsilon(v)\ell_{v}\to 0$, as $v\to\infty$,
then

(2.8)
$D_{q}=\lim_{v\to\infty}D_{q}(v)=\frac{\log\zeta_{G}(2q)-q\log\zeta_{G}(2)}{1-q},$

where

(2.9) $\zeta_{G,v}(s)=\sum_{j=1}^{\infty}\ell_{j}^{-s}.$

Depending on the choice of the set $S=\\{\ell_{j}\\}_{j=1}^{\infty}$, which
corresponds to a star graph $G(S)$, different types of zeta functions may
arise which may satisfy a functional equation. For example $S={\mathbb{N}}$
gives rise to Riemann’s zeta function

$\zeta_{G({\mathbb{N}})}(s)=\zeta(s)=\sum_{n=1}^{\infty}n^{-s}$

whereas $S={\mathcal{N}}=\\{|\xi|^{2}\mid\xi\in{\mathcal{L}}\\}$, where
${\mathcal{L}}={\mathbb{Z}}(1,0)\oplus{\mathbb{Z}}(0,\sqrt{D})$ denotes a
rectangular lattice, gives rise to Epstein’s zeta function associated with the
quadratic form $Q(x,y)=x^{2}+Dy^{2}$,

$\zeta_{G({\mathcal{N}})}(s)=\zeta_{Q}(s)=\sum_{(x,y)\in{\mathbb{Z}}^{2}}Q(x,y)^{-s}.$

In these cases, the functional equation for the respective zeta function

$\zeta_{G}(1-s)=\varphi_{G}(s)\zeta_{G}(s)$

yields a symmetry relation around the critical point $q=1/4$:

(2.10)
$D_{1/2-q}=D_{q}\times\frac{1-q}{1/2+q}+\frac{\log\varphi(2q)+(2q-1/2)\log\zeta(2)}{1/2+q}$

##  3\. Spectral theory of quantum star graphs

In this section we briefly recall the definition and spectral properties of
quantum star graphs which are necessary to introduce quantities such as the
Renyi entropy in this context. For a more detailed discussion we refer the
reader to the introduction of the paper [18].

###  3.1. Schrödinger operator on star graphs

A star graph $G_{v}$ consists of a central vertex and $v$ outlying vertices.
There are, thus, $v$ bonds and we denote the associated bond lengths by
$L_{1},\cdots,L_{v}$. We denote the vector of bond lengths by
$L=(L_{1},\cdots,L_{v})$. Moreover we assume these to be ordered according to
increasing length:

$L_{1}\leq L_{2}\leq\cdots\leq L_{v}.$

Let us introduce a space of $C^{2}$ functions $\psi$ on $G_{v}$, where we
denote by $\psi_{j}$ the component of the function on the $j$th bond. We
identify each bond with an interval $[0,L_{j}]$ and introduce a coordinate $x$
on the interval. Then $\psi_{j}:[0,L_{j}]\mapsto{\mathbb{R}}$. For each bond
the point $x=0$ corresponds to the central vertex, whereas $x=L_{j}$
corresponds to the outlying vertex respectively.

The Schrödinger operator on $G_{v}$ is given by the Euclidean Laplacian
$-d^{2}/dx^{2}$ and the components of the function $\psi$ must satisfy the
following matching conditions:

(3.1) $\psi_{i}(0)=\psi_{j}(0)=C,\quad\text{for all}\;i,j=1,\cdots,v,$

where $C$ is a real constant, and

(3.2) $\sum_{j=1}^{v}\psi_{j}^{\prime}(0)=\frac{C}{\alpha},$

as well as

(3.3) $\psi_{j}^{\prime}(0)=0,\quad\text{for all}\;j=1,\cdots,v,$

where $\alpha\in{\mathbb{R}}\setminus\\{0\\}$ denotes a coupling parameter
[20].

###  3.2. Renyi entropy for star graphs

As explained in [18], it is natural to study the properties of the
eigenfunctions of a quantum star graph in terms of the maximum amplitude
squared on each bond length. This can be calculated explicitly for an $L^{2}$
normalized eigenfunction with spectral parameter $k$ on the $i$th bond as

(3.4)
$A_{i}(k,L;v)=\frac{2\sec(kL_{i})^{2}}{\sum_{j=1}^{v}L_{j}\sec(kL_{j})^{2}}.$

This gives rise to a probability measure on the finite set
$\Omega_{v}=\\{1,\cdots,v\\}$, given by

(3.5) $\mu_{k}(i)=\frac{\sec(kL_{i})^{2}}{\sum_{j=1}^{v}\sec(kL_{j})^{2}}.$

The Renyi entropy associated with $\mu$ on $\Omega_{v}$ is then given by

(3.6) $H_{q}(\mu_{k})=\frac{1}{1-q}\log\left(\sum_{i=1}^{v}\mu(i)^{q}\right).$

We may rewrite this as

(3.7) $H_{q}(\mu_{k})=\frac{\log M_{q}(k,L;v)-q\log M_{1}(k,L;v)}{1-q}$

where we introduce the moment sum

(3.8) $M_{q}(k,L;v)=\sum_{j=1}^{v}|\sec(kL_{j})|^{2q}.$

##  4\. Semiclassical Multifractality

###  4.1. Semiclassical asymptotics of the Renyi entropy

Let us suppose that the eigenvalue $k$ is in the $m$th cluster and that
$s_{m,i+1}<k<s_{m,i}$.

We have

(4.1) $\begin{split}\cos(kL_{j})&=(-1)^{m+1}\sin((k-s_{m,j})L_{j})\\\
&=(-1)^{m+1}(k-s_{m,j})L_{j}+O(\frac{m^{3}\epsilon^{3}\ell_{v}^{3}}{L^{3}}).\end{split}$

Therefore, we have, as $v\to\infty$,

$M_{q}(n,L;v)\asymp\sum_{j=1}^{v}(k-s_{m,j})^{-2q}L_{j}^{-2q}.$

Moreover, the spacing is of order

$|s_{m,j}-s_{m,k}|\asymp\frac{m\epsilon}{L^{2}}|\ell_{j}-\ell_{k}|.$

Hence, we renormalize: $k^{\prime}=kL^{2}/(m\epsilon)$,
$s^{\prime}_{m,j}=s_{m,j}L^{2}/(m\epsilon)$, so that we have
$|s^{\prime}_{m,j}-s^{\prime}_{m,k}|\asymp|\ell_{j}-\ell_{k}|$. And, since
$L_{j}\asymp L$, we have

$M_{q}(n,L;v)\asymp(\frac{m\epsilon}{L^{2}})^{-2q}\sum_{j=1}^{v}|k^{\prime}-s^{\prime}_{m,j}|^{-2q}.$

We have
$k^{\prime}-s^{\prime}_{m,j}=(k^{\prime}-s^{\prime}_{m,i})+(s^{\prime}_{m,i}-s^{\prime}_{m,j})$
and assume we have a subsequence where
$\sigma_{m}=k^{\prime}-s^{\prime}_{m,i}\to\sigma>0$ as $m\to+\infty$ along
this subsequence. So, it is sufficient to study the sum

$\sum_{j=1}^{v}|\sigma_{m}+\ell_{i}-\ell_{j}|^{-2q}=\sum_{j\leq
i+2}|\sigma_{m}+\ell_{i}-\ell_{j}|^{-2q}+\sum_{j>i+2}(\ell_{j}-\ell_{i}-\sigma_{m})^{-2q}.$

We have

$\sum_{j>i+2}(\ell_{j}-\ell_{i})^{-2q}\leq\sum_{j>i+2}(\ell_{j}-\ell_{i}-\sigma_{m})^{-2q}\leq\sum_{j>i+2}(\ell_{j}-\ell_{i+1})^{-2q}.$

Let us introduce an invertible function $\mathfrak{j}\in
C^{1}({\mathbb{R}}_{+})$ which satisfies $\mathfrak{j}(\ell_{j})=j$, $\forall
j=1,\dots,v$. We then have

$\sum_{j=i+3}^{v}(\ell_{j}-\ell_{i+1})^{-2q}\leq\int_{i+2}^{v}(\mathfrak{j}^{-1}(t)-\ell_{i+1})^{-2q}dt$

and

$\sum_{j=i+3}^{v}(\ell_{j}-\ell_{i})^{-2q}\geq\int_{i+3}^{v+1}(\mathfrak{j}^{-1}(t)-\ell_{i})^{-2q}dt.$

Define the Mellin transform of a compactly supported piecewise continuous
function $f$ as

$({\mathcal{M}}f)(s)=\int_{0}^{\infty}t^{-s}f(t)dt.$

We may rewrite

$\int_{i+2}^{v}(\mathfrak{j}^{-1}(t)-\ell_{i+1})^{-2q}dt=\int_{\ell{i+2}-\ell_{i+1}}^{\ell_{v}-\ell_{i+1}}t^{-2q}\mathfrak{j}^{\prime}(t+\ell_{i+1})dt=({\mathcal{M}}F(v,\mathfrak{j}))(2q)$

where

$F(v,\mathfrak{j})(t)=(\mathbbm{1}_{[\ell_{i+2},\ell_{v}]}\mathfrak{j}^{\prime})(t+\ell_{i+1}).$

If

$\lim_{v\to+\infty}\frac{\log({\mathcal{M}}F(v+1,\mathfrak{j}))(s)}{\log({\mathcal{M}}F(v,\mathfrak{j}))(s)}=1,$

then the asymptotics of the Renyi entropy may be calculated in terms of the
Mellin transform as

(4.2) $\begin{split}H_{q}(\mu_{k})=&\frac{\log M_{q}(k,L;v)-q\log
M_{1}(k,L;v)}{1-q}\\\
\sim&\frac{\log({\mathcal{M}}F(v,\mathfrak{j}))(2q)-q\log({\mathcal{M}}F(v,\mathfrak{j}))(2)}{1-q}.\end{split}$

###  4.2. An example of multifractality

Let us consider $\ell_{j}=j^{1/n}$ for some integer $n>0$. Then
$\mathfrak{j}(t)=t^{n}$. Let us suppose $q>n/2$. Set $\sigma_{1}=1+\sigma$.

Asymptotically, the Mellin transform can be evaluated as

(4.3)
$\begin{split}({\mathcal{M}}F(v,\mathfrak{j}))(2q)=&\int_{2^{1/n}}^{v^{1/n}}(s-\sigma_{1})^{-2q}ns^{n-1}ds\\\
=&n\int_{2^{1/n}-\sigma_{1}}^{v^{1/n}-\sigma_{1}}s^{-2q}(s+\sigma_{1})^{n-1}ds\\\
=&n\sum_{k=0}^{n-1}c_{n,k}\sigma_{1}^{n-1-k}\int_{2^{1/n}-\sigma_{1}}^{v^{1/n}-\sigma_{1}}s^{-2q+k}ds\\\
=&n\sum_{k=0}^{n-1}c_{n,k}\sigma_{1}^{n-1-k}\frac{(2^{1/n}-\sigma_{1})^{-2q+k+1}}{-2q+k+1}+O(v^{1-2q/n})\end{split}$

For $n\geq 2$, we have

(4.4)
$({\mathcal{M}}F(v,\mathfrak{j}))(2)=\int_{2}^{v}(t^{1/n}-1-\sigma)^{-2}dt\sim\begin{cases}\frac{n}{n-2}v^{1-2/n},\quad\text{if}\;n\geq
3,\\\ \\\ \log v,\quad\text{if}\;n=2.\end{cases}$

This yields the following asymptotics for the Renyi entropy:

(4.5) $H_{q}(k_{n})=\frac{\log M_{q}-q\log
M_{1}}{1-q}\sim\begin{cases}\frac{q}{q-1}(1-\frac{2}{n})\log
v,\quad\text{if}\;n\geq 3,\\\ \\\ \frac{q}{q-1}\log\log
v,\quad\text{if}\;n=2.\\\ \end{cases}$

We note that the dependence on $\sigma$ only appears in lower order terms. In
particular, for $q>n/2$ and $n\geq 2$, we have the fractal exponent

(4.6) $D_{q}=\lim_{k\to+\infty}\frac{H_{q}(k)}{\log
v}=\begin{cases}1-\frac{2q-n}{n(q-1)},\quad\text{if}\;n\geq 3,\\\ \\\
0,\quad\text{if}\;n=2,\end{cases}$

where we used $v\to+\infty$, as $k\to+\infty$.

For an example where $D_{q}=1$, consider $\ell_{j}=(v+j)^{\delta}$,
$\delta>0$, which leads to moment sums of the form

$\sum_{j=1}^{v}(\sigma^{\prime}+(v+j)^{\delta})^{-2q}\asymp v^{1-2\delta q}.$

So $\log M_{q}\sim(1-2q)\delta\log v$, which yields the asymptotic

$H_{q}(k)=\frac{\log M_{q}-q\log M_{1}}{1-q}\sim\frac{(1-2\delta
q-q(1-2\delta))\log v}{1-q}=\log v,$

and, hence, $D_{q}=\lim_{m\to+\infty}H_{q}(k)/\log v=1$.

##  5\. Ground state multifractality

Let us now fix $m=0$ and take $\epsilon=\epsilon(v)$ to be a decaying function
such that $\epsilon\ell_{v}\to 0$, which ensures separate clusters.

Note that this regime is markedly different from the semiclassical regime
above. In the semiclassical regime we considered a sequence of eigenvalues
which tends to infinity ($m\to+\infty$) and the number of vertices
$v=v(m)\to+\infty$ is coupled to the semiclassical parameter to ensure
separation of the clusters.

In the ground state regime we consider the lowest eigenvalue in the bottom
cluster in a weak coupling regime such that $k_{0}=\pi/(2L)+\sigma_{0}$ and
$\sigma_{0}\to 0$. In this regime we have $m=0$ and, therefore, $m$ is not
coupled to $v$. Analogously to our definition for the ground state of S̆eba’s
billiard we, therefore, define the fractal exponent as

(5.1) $D_{q,v}=\lim_{\sigma_{0}\to 0}\frac{\log M_{q}(k_{0})-q\log
M_{1}(k_{0})}{1-q}.$

The moment sum for ${\mathcal{C}}_{0}$ is of the form

$M_{q}(k_{0})=\sum_{j=1}^{v}(k_{0}-s_{0,j})^{-2q}=\sum_{j=1}^{v}\left(\sigma_{0}+\frac{\pi\epsilon(v)\ell_{j}}{2L(L+\epsilon\ell_{j})}\right)^{-2q}$

so that

$\lim_{\sigma_{0}\to
0}M_{q}(k_{0})=\left(\frac{\pi\epsilon(v)}{2L}\right)^{-2q}\sum_{j=1}^{v}\left(\frac{\ell_{j}}{L+\epsilon\ell_{j}}\right)^{-2q}$

and, in turn, we compute

(5.2) $D_{q,v}=\frac{\log\zeta_{G,v}(2q)-q\log\zeta_{G,v}(2)}{1-q}$

where

(5.3)
$\zeta_{G,v}(s)=\sum_{j=1}^{v}\left(\frac{\ell_{j}}{L+\epsilon\ell_{j}}\right)^{-s}.$

Now, since $m$ is fixed, we may take $v\to\infty$ which gives rise to the zeta
function associated with the length spectrum of the graph

(5.4) $\zeta_{G,\infty}(2q)=\sum_{j=1}^{\infty}\ell_{j}^{-2q}$

For example take $\ell_{j}=j$, then the arithmetic graph $G_{\mathbb{N}}$
which arises in the limit $v\to+\infty$ has length spectrum ${\mathbb{N}}$ and
gives rise to the Riemann zeta function

(5.5)
$\zeta_{G_{\mathbb{N}},\infty}(2q)=\zeta(2q)=\sum_{n=1}^{\infty}n^{-2q}.$

If one takes the set of lengths $\ell_{j}$ to be the set of lengths of vectors
in a rectangular lattice ${\mathcal{L}}$, then the associated graph
$G_{\mathcal{L}}$ gives rise to the Epstein zeta function

(5.6)
$\zeta_{\mathcal{L}}(2q)=\sum_{\xi\in{\mathcal{L}}\setminus\\{0\\}}\|\xi\|^{-2q}=\zeta_{Q}(q)$

where $Q$ is the quadratic form induced by the Euclidean norm $\|\cdot\|$.

If we recall the functional equation for Riemann’s zeta function

$\zeta(1-s)=\varphi(s)\zeta(s)$

then we have the following symmetry relation of the fractal exponent $D_{q}$
around the value $q=1/4$ under the transformation $q\to 1/2-q$:

(5.7)
$\begin{split}D_{1/2-q}=&\frac{\log\zeta(2(1/2-q))-(1/2-q)\log\zeta(2)}{1-(1/2-q)}\\\
=&\frac{\log\varphi(2q)+\log\zeta(2q)+(q-1/2)\log\zeta(2)}{1/2+q}\\\
=&D_{q}\times\frac{1-q}{1/2+q}+\frac{\log\varphi(2q)+(2q-1/2)\log\zeta(2)}{1/2+q}\end{split}$

## References

  * [1] B. L. Altshuler, V. E. Kravtsov, I. V. Lerner, Statistical properties of mesoscopic fluctuations and similarity theory, JEPT Lett. 43 (1986), no. 7, 441–44.
  * [2] N. Anantharaman, Entropy and the localization of eigenfunctions, Ann. Math. (2) 168 (2008), no. 2, 435–475.
  * [3] Y. Y. Atas and E. Bogomolny, Multifractality of eigenfunctions in spin chains, Phys. Rev. E 86 (2012), 021104.
  * [4] A. Bäcker, M. Haque, and I. M. Khaymovich, Multifractal dimensions for random matrices, chaotic quantum maps, and many-body systems, Phys. Rev. E 100 (2019), 032117.
  * [5] G. Berkolaiko, E. B. Bogomolny, J. P. Keating, Star graphs and Seba billiards, J. Phys. A 34 (2001), 335–350.
  * [6] G. Berkolaiko, J. P. Keating, Two-point spectral correlations for star graphs, J. Phys. A 32 (1999), 7827–7841.
  * [7] G. Berkolaiko, J. P. Keating, B. Winn, Intermediate wave function statistics, Phys. Rev. Lett. 91 (2003), 134103.
  * [8] G. Berkolaiko, J. P. Keating, B. Winn, No quantum ergodicity for star graphs, Commun. Math. Phys. 250 (2004), 259–285.
  * [9] A. M. Bilen, I. Garcia-Mata, B. Georgot, O. Giraud, Multifractality of open quantum systems, Phys. Rev. E 100 (2020), 032223.
  * [10] E. Bogomolny, U. Gerland, C. Schmit, Singular Statistics, Phys. Rev. E 63 (2001), 036206.
  * [11] E. Bogomolny, O. Giraud, Perturbation approach to multifractal dimensions for certain critical random matrix ensembles, Phys. Rev. E 84 (2011), 036212.
  * [12] E. Bogomolny, P. Leboeuf, C. Schmit, Spectral statistics of chaotic systems with a pointlike scatterer, Phys. Rev. Lett. 85 (2000), 2486–2489.
  * [13] Y. Colin de Verdière, Ergodicité et fonctions propres du laplacien, Comm. Math. Phys. 102 (1985), no. 3, 497–502.
  * [14] F. Faure, S. Nonnenmacher, S. de Bièvre, Scarred eigenstates for quantum cat maps of minimal periods, Comm. Math. Phys. 239 (2003), no. 3, 449–492.
  * [15] A. Hassell, Ergodic billiards that are not quantum unique ergodic (with an appendix by Andrew Hassell and Luc Hillairet) Ann. Math. (2) 171 (2010), no. 1, 605–618.
  * [16] R. Holowinsky, Sieving for mass equidistribution, Ann. Math. (2) 172 (2010), no. 2, 1499–1516.
  * [17] R. Holowinsky, K. Soundararajan, Mass equidistribution for Hecke eigenforms, Ann. Math. (2) 172 (2010), no. 2, 1517–28.
  * [18] J. P. Keating, J. Marklof and B. Winn, Value distribution of the eigenfunctions and spectral determinants of quantum star graphs, Comm. Math. Phys. 241 (2003), 421–452.
  * [19] J. P. Keating, H. Ueberschär, Multifractal eigenfunctions for a singular quantum billiard, Comm. Math. Phys., to appear.
  * [20] T. Kottos. U. Smilansky, Periodic orbit theory and spectral statistics for quantum graphs, Ann. Phys. 274 (1999), 76–124.
  * [21] V. E. Kravtsov, K. A. Muttalib, New Class of Random Matrix Ensembles with Multifractal Eigenvectors, Phys. Rev. Lett. 79 (1997), 1913–1916.
  * [22] P. Kurlberg, H. Ueberschär, Quantum Ergodicity for Point Scatterers on Arithmetic Tori, Geom. Funct. Anal. 24 (2014), no. 5, 1565–1590.
  * [23] P. Kurlberg, H. Ueberschär, Superscars in the S̆eba billiard, J. Eur. Math. Soc. 19 (2017), no. 10, pp. 2947–2964.
  * [24] E. Lindenstrauss, Invariant measures and arithmetic quantum unique ergodicity, Ann. Math. (2) 163 (2006), no. 1, 165–219.
  * [25] J. Marklof, Z. Rudnick, Quantum unique ergodicity for parabolic maps, Geom. Funct. Anal. 10 (2000) 1554–1578
  * [26] J. Martin, I. Garcia-Mata, O. Giraud, and B. Georgeot, Multifractal wave functions of simple quantum maps, Phys. Rev. E 82 (2010), 046206.
  * [27] A. D. Mirlin, Statistics of energy levels and eigenfunctions in disordered systems, Physics Reports 326 (2000), 259–382.
  * [28] P. J. Richens, M. V. Berry, Pseudointegrable systems in classical and quantum mechanics, Physica D 2 (1981), 495–512.
  * [29] Z. Rudnick, P. Sarnak, The behaviour of eigenstates of arithmetic hyperbolic manifolds, Comm. Math. Phys 161 (1994), no. 1, 195–213.
  * [30] Z. Rudnick, H. Ueberschär, Statistics of Wave Functions for a Point Scatterer on the Torus, Comm. Math. Phys. 316 (2012), no. 3, 763–782.
  * [31] M. Schreiber, H. Grussbach, Multifractal wave functions at the Anderson transition Phys. Rev. Lett. 67 (1991), 607–610.
  * [32] P. S̆eba, Wave chaos in singular quantum billiard, Phys. Rev. Lett. 64 (1990), 1855–1858.
  * [33] C. E. Shannon, A mathematical theory of communication, The Bell System Technical Journal 27 (1948), no. 3, 379–423.
  * [34] C. E. Shannon, A mathematical theory of communication, The Bell System Technical Journal 27 (1948), no. 4, 623–666.
  * [35] A. I. Snirelman, Ergodic properties of eigenfunctions, Uspehi Mat. Nauk. 29 (1974), no. 6 (180), 181–182.
  * [36] K. Soundararajan, Quantum Unique Ergodicity for SL(2,Z)\H, Ann. Math. (2) 172 (2010), no. 2, pp. 1529–1538.
  * [37] H. Ueberschär, Quantum chaos for point scatterers on flat tori, Phil. Trans. R. Soc. A.372 (2014), 20120509.
  * [38] S. Zelditch, Uniform distribution of eigenfunctions on compact hyperbolic surfaces, Duke Math. J. 55 (1987), no. 4, 919–941.


# Asymptotic methods for transaction costs

Eberhard Mayerhofer11endnote: 1University of Limerick, Department of
Mathematics & Statistics, V94 T9PX, Limerick, Ireland,
<EMAIL_ADDRESS>


# Adaptive Sparse Gaussian Process

Vanessa Gómez-Verdejo1, Manel Martínez-Ramón2, .

Corresponding author: Vanessa Gómez-Verdejo (email: vanessag@ing.uc3m.es).
This work is partially supported by grants PID2020-115363RB-I00 funded by
MCIN/AEI/10.13039/50110001103, Spain, EPSCoR OIA-1757207, funded by NSF, USA,
and the King Felipe VI Endowed Chair of the UNM. 1Signal Processing and
Communications Department, University Carlos III de Madrid, Leganés, 28911,
Madrid, Spain.  2Electrical and Computer Engineering Department, The
University of New Mexico, Albuquerque, 87106, NM, USA.

###### Abstract

Adaptive learning is necessary for non-stationary environments where the
learning machine needs to forget past data distribution. Efficient algorithms
require a compact model update to not grow in computational burden with the
incoming data and with the lowest possible computational cost for online
parameter updating. Existing solutions only partially cover these needs. Here,
we propose the first adaptive sparse Gaussian Process (GP) able to address all
these issues. We first reformulate a variational sparse GP algorithm to make
it adaptive through a forgetting factor. Next, to make the model inference as
simple as possible, we propose updating a single inducing point of the sparse
GP model together with the remaining model parameters every time a new sample
arrives. As a result, the algorithm presents a fast convergence of the
inference process, which allows an efficient model update (with a single
inference iteration) even in highly non-stationary environments. Experimental
results demonstrate the capabilities of the proposed algorithm and its good
performance in modeling the predictive posterior in mean and confidence
interval estimation compared to state-of-the-art approaches.

###### Index Terms:

Sparse Gaussian Process, Variational learning, Online learning, Adaptive
learning

##  I Introduction

Signal processing applications [1] usually require online learning methods,
i.e., algorithms capable of self-updating as new data arrive in the system.
Moreover, in applications such as spectral estimation, signal modeling,
adaptive filtering, or array processing, the signals to be processed are not
stationary (their distributions change over time), so we need algorithms
capable of adapting to the new data distributions and forgetting the past.
This is known as adaptive learning [2, 3].

In this sense, Bayesian models [4] and in particular Gaussian Processes (GP)
[5] seem to be an ideal solution since each time a new sample arrives the
predictive posterior can be updated using the previous posterior as new prior
of the model and multiplying it by the likelihood of the new data. This idea
can be further extended to models which include hidden latent variables and
require learning algorithms such as Expectation Maximization, providing what
is known in the literature as online Variational Bayes [6]. The main problem
of these approaches relies on the fact that the predictive posterior
parameters usually depend on all the training data and therefore, as every a
new sample is available, its complexity tends to grow with
$\mathcal{O}(N^{3})$, $N$ being the number of data, which makes them
intractable if not combined with pruning schemes [7].

A possible solution to this problem is to use low computational cost GP
versions [8] or, more specifically, to resort to sparse GP [9]. In these
approaches, the model solution only depends on $M<<N$ points in the
observation space, called inducing points. This avoids the complexity of the
model growing as each a new sample arrives. However, this strategy complicates
the model inference since, in each iteration, the algorithm has to learn and
update the kernel and noise parameters and, additionally, the positions of
these inducing points and their variational distribution if we resort to
variational versions of these models [10].

Despite the number of online applications, there are few proposals for online
GP able to provide, on the one hand, a compact solution or, in other words,
whose complexity does not grow with the incoming training data, while, on the
other hand, keep a low computational cost and capable of efficiently updating
the parameters of the model with every new data. For example, [11] proposed a
compact Bayesian online algorithm based on an approximation to the real
predictive posterior with a computational cost of $\mathcal{O}(NM^{2})$. The
drawback of this approach is that, despite sequentially updating the subset of
$M$ relevant data, the remaining model hyperparameters, as the kernel
parameters and noise variance, keep fixed.

Other approaches directly propose to use sparse GP models and exploit the
structure of the resulting kernel matrices. For example, [12] approximates the
kernel matrix using $K$ eigenfunctions to obtain a fast computation of the
likelihood derivatives. A more recent work [13] introduces a structured kernel
interpolation approach of order $K$ to reduce the computational cost of
updating the predictive distribution to $\mathcal{O}(MK^{2})$. This model has
the limitation of having to predefine (and keep fixed) the set of inducing
points. Besides, it has the additional drawback of being based on sparse GPs,
so the obtained posterior is an approximation to the exact one, and its
performance tends to be worse than the batch approach.

To overcome this last limitation, it is preferable to rely on variational
sparse GPs algorithms, where the fact of including a variational distribution
over the inducing points avoids the error of approximating the posterior. In
this line, we can find the works [14] or [15] where the authors propose an
incremental learning algorithm for variational sparse GPs, although their
hyperparameters are fixed during the training procedure.

Another interesting approach is the model presented in [16], where the authors
propose to update the variational bound with an online variational Bayesian
scheme. This way the model is efficiently updated at each iteration with
complexity similar to that of a standard sparse GP, but with the advantage of
having a variational bound expression to be optimized with respect to all the
model hyperparameters, including the inducing point locations.

Within the models based on variational SGP formulations, we can find other
methods [17, 16] that use Variational Stochastic Inference (VSI) to obtain
online versions, since the VSI strategy allows them to work with mini-batches
of data, providing a straightforward scheme to add new data to the minibatch
as they arrive. However, these approaches present two problems: first, the
stochastic optimization assumes that the data subsampling process is uniformly
random, an assumption that does not necessarily hold in online or non-
stationary environments; second, this optimization usually requires that each
mini-batch is processed in several iterations, which is often not compatible
with the computational constraints of an online approach.

Finally, the major shortcoming that we have found in all the algorithms above-
mentioned for its use in real signal processing applications is that they are
designed to work in an online environment, but not in nonstationary scenarios.
To our knowledge, there are no adaptive versions of compact or sparse GP
formulations.

So, to provide a solution to these needs, in this paper, we propose an
adaptive GP that includes into a single model:

  * •

Adaptive learning, ability to forget the information of the remote samples
and, thus, capability to self-adapting to non-stationary environments.

  * •

Efficiently update of the predictive distribution with a cost of
$\mathcal{O}(M^{3})$ each time a new sample arrives.

  * •

If it is desired, at each iteration, it can update the set of inducing points
and the rest of the hyperparameters of the model with a cost of
$\mathcal{O}(NM^{2})$.

  * •

To update the inducing point set at each iteration, only a single inducing
point is updated. This significantly reduces the number of parameters to be
inferred, and a single inference iteration is usually sufficient to obtain an
accurate value of the model parameters.

  * •

Finally, unlike reference approaches, the proposed model is able to recover
the solution of the batch variational sparse GP formulation. That is, if at a
given time instant we consider the same model parameters and eliminate the
forgetting factor, its predictive distribution converges to that of the
original model. As we will see in the experimental section, in practice, this
translates into an improvement in the model performance.

The rest of the paper is structured as follows. Section 2 reviews GP models,
paying special attention to sparse GPs and the variational version used as
starting point of our proposal. Then, Section 3 introduces the proposed
algorithm, starting by developing its adaptive formulation and then explaining
how to do efficient online learning on this model. Section 4 analyzes the
performance of the model showing over a load forecasting application the
advantages of having a compact, efficient, and adaptive GP model able to track
the signal changes over time. Finally, Section 5 presents the main conclusions
of this work.

The implementation of all proposed models is available at
https://github.com/vgverdejo/AdaptiveSparseGP.

##  II A review of Sparse Gaussian Processes

###  II-A Gaussian Processes

Assume a regression problem defined by a training dataset $\\{{\bf
x}_{n},y_{n}\\}_{n=1}^{N}$, with ${\bf x}_{n}\in\mathbb{R}^{D}$ and
$y_{n}\in\mathbb{R}$, where each $y_{n}$ is generated by a noisy version of a
latent function $f({\bf x}_{n})=f_{n}$, i.e.

$y_{n}=f_{n}+\varepsilon_{n},\leavevmode\nobreak\ \leavevmode\nobreak\
\varepsilon_{n}\sim\mathcal{N}(0,\sigma^{2}).$ (1)

Assume further that function $f_{n}$ can be expressed as

$f_{n}=\langle{\bf w},\varphi({\bf x}_{n})\rangle={\bf w}^{\top}\varphi({\bf
x}_{n}),$ (2)

where ${\bf w}$ are the primal parameters, which are modeled as a multivariate
Gaussian with zero mean and covariance matrix ${\boldsymbol{\Sigma}}_{p}$.
Function $\varphi(\cdot)$ is a nonlinear mapping of the input data ${\bf
x}_{n}$ into a Hilbert space $\mathcal{H}$ endowed with a dot product
$\langle\varphi({\bf x}),\varphi({\bf x}^{\prime})\rangle$, usually called a
Mercer’s kernel function, which, by virtue of the Mercer’s Theorem [18], is a
dot product if and only if such function is positive semidefinite. From
expression (1) the covariance between instances of the latent function can be
computed as

$\begin{split}\mathbb{E}\left[f_{n}\cdot
f_{m}\right]=&\mathbb{E}\left[\varphi({\bf x}_{n})^{\top}{\bf w}{\bf
w}^{\top}\varphi({\bf x}_{m})\right]\\\ =&\varphi({\bf
x}_{n})^{\top}\mathbb{E}\left[{\bf w}{\bf w}^{\top}\right]\varphi({\bf
x}_{m})\\\ =&\varphi({\bf x}_{n})^{\top}{\boldsymbol{\Sigma}}_{p}\varphi({\bf
x}_{m})\\\ =&\langle{\boldsymbol{\Sigma}}_{p}^{1/2}\varphi({\bf
x}_{n}),{\boldsymbol{\Sigma}}_{p}^{1/2}\varphi({\bf x}_{m})\rangle=k({\bf
x}_{n},{\bf x}_{m}).\end{split}$ (3)

Since ${\boldsymbol{\Sigma}}_{p}$ is positive semidefinite, then function
$k({\bf x}_{n},{\bf x}_{m})$ is positive semidefinite and hence it is a valid
kernel dot product.

The GP model for regression [19] considers that the joint process of the
latent functions ${\bf f}=[f({\bf x}_{1}),\ldots,f({\bf x}_{N})]$ is drawn
from a multivariate Gaussian prior $p({\bf f})$ with with zero mean and
covariance ${\bf K}_{xx}$, where this matrix contains all the kernel dot
products $k({\bf x}_{n},{\bf x}_{m})$ between input data ${\bf x}_{n}$ and
${\bf x}_{m}$ for $1\leq n,m\leq N$. From this model, we can obtain the
predictive posterior for a test sample ${\bf x}_{*}$ as

$\begin{split}p(f_{*}|{\bf y})&=\int p(f_{*}|{\bf f})p({\bf f}|{\bf y})d{\bf
f}=\\\ &=\mathcal{N}\left(f_{*}|{\bf k}_{*x}\left({\bf K}_{xx}+\sigma^{2}{\bf
I}\right)^{-1}{\bf y},\right.\\\ &\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \left.k_{**}-{\bf k}_{*x}\left({\bf
K}_{xx}+\sigma^{2}{\bf I}\right)^{-1}{\bf k}_{x*}\right),\end{split}$ (4)

where ${\bf k}_{*x}$ is a row vector with the kernel products $k({\bf
x}_{*},{\bf x}_{n})$, for $1\leq n\leq N$, $k_{**}=k({\bf x}^{*},{\bf x}^{*})$
is the kernel product of ${\bf x}_{*}$ with itself, and ${\bf
y}=[y_{1},\ldots,y_{N}]\in\mathbb{R}^{N}$ is the vector representation of all
the training labels. The computation of this posterior requires the inversion
of an $N\times N$ matrix, which takes a computational burden of
$\mathcal{O}(N^{3})$.

###  II-B Sparse Gaussian Processes

The motivation of the Sparse GP (SGP) is to reduce this computational cost by
including a set of $M$ inducing variables or inducing points ${\bf u}_{m}$,
$1\leq m\leq M$, that approximate the dual space to produce a compact GP model
with computational burden $\mathcal{O}(NM^{2})$.

The first versions of the SGP can be unified in a common framework [9]
considering that the covariance of the GP prior is modified by an approximate
matrix. For example, in [20] ${\bf K}_{xx}$ is replaced by a Nyströn
approximation ${\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf K}_{ux}$, ${\bf K}_{uu}$ and
${\bf K}_{xu}$ being the kernel matrix of the inducing points with the
training data and with themselves, respectively. An evolution of this method,
introduced in [21], corrects the approximated matrix with the term
$\rm{diag}({\bf K}_{xx})-\rm{diag}({\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf K}_{ux})$
so that the diagonal of ${\bf K}_{xx}$ is exact.

The main drawback of these approximations is that they do not tend to the
exact GP, since they start from an approximation to the real prior, so their
posterior will be an approximation to the real one. Moreover, the inducing
inputs constitute additional parameters to be inferred, increasing the
overfitting risk.

###  II-C Variational Sparse Gaussian Processes

The Variational Sparse Gaussian Process (VSGP) introduced by Titsias in [10]
is intended to overcome the abovementioned limitations by minimizing the
divergence between the exact GP posterior and a variational approximation
where the inducing points are modeled as variational parameters. This way,
minimizing the divergence with respect to these inducing points, we get the
SGP to tend to the original one and, introducing a variational prior over the
inducing variables, we get to reduce the overfitting risk.

To introduce this formulation, let’s consider a set of inducing inputs ${\bf
U}=[{\bf u}_{1},\ldots,{\bf u}_{M}]$, with their corresponding latent
variables ${\bf f}_{u}=[f({\bf u}_{1}),\ldots,f({\bf u}_{M})]^{\top}$. Then,
the predictive posterior of the VSGP is given by

$p(f_{*}|{\bf y})=\int\int p(f_{*}|{\bf f}_{u},{\bf f})p({\bf f}|{\bf
f}_{u},{\bf y})p({\bf f}_{u}|{\bf y})d{\bf f}d{\bf f}_{u}.$ (5)

Assuming that ${\bf f}_{u}$ is a sufficient statistic for ${\bf f}$,
$p(f_{*}|{\bf f}_{u},{\bf f})$ can be approximated by $p(f_{*}|{\bf f}_{u})$.
Moreover, approximating the posterior distribution $p({\bf f}_{u}|{\bf y})$ by
a variational distribution $q({\bf f}_{u})$, the following approximate
posterior can be used

$\begin{split}q(f_{*})=\int\int p(f_{*}|{\bf f}_{u})p({\bf f}|{\bf
f}_{u})q({\bf f}_{u})d{\bf f}d{\bf f}_{u}=\int p(f_{*}|{\bf f}_{u})q({\bf
f}_{u})d{\bf f}_{u}.\end{split}$ (6)

Considering that $q({\bf f}_{u})$ is a Gaussian distribution with mean
${\boldsymbol{\mu}}$ and covariance ${\bf A}$, the following approximate
predictive posterior is obtained

$\begin{split}q(f_{*})&=\mathcal{N}\left(f_{*}|m_{*},v_{*}\right)\\\
m_{*}&={\bf k}_{*u}{\bf K}_{uu}^{-1}{\boldsymbol{\mu}}\\\ v_{*}&=k_{**}-{\bf
k}_{*u}{\bf K}^{-1}_{uu}{\bf k}_{u*}+{\bf k}_{*u}{\bf K}_{uu}^{-1}{\bf A}{\bf
K}_{uu}^{-1}{\bf k}_{u*},\end{split}$ (7)

where ${\bf k}_{*u}={\bf k}_{u*}^{\top}$ is a row vector with the kernel
products between ${\bf x}_{*}$ and ${\bf U}$. The main advantage of this
approach relies on the fact that, once ${\boldsymbol{\mu}}$, ${\bf A}$, and
${\bf U}$ are obtained, the prediction of a new sample does not depend on
training data, but only on these parameters. So the complexity of predicting a
new sample is $\mathcal{O}(M^{3})$ (it does not depend on $N$).

To complete this model, one needs to find the variational parameters
${\boldsymbol{\mu}}$, ${\bf A}$, and ${\bf U}$. For this purpose, [10]
introduces the following variational lower bound of the true log marginal
likelihood

$\log p({\bf y})\geq F_{V}({\bf U},q({\bf f}_{u}))=\int q({\bf
f}_{u})\log\frac{G({\bf f}_{u},{\bf y})p({\bf f}_{u})}{q({\bf f}_{u})}d{\bf
f}_{u}$ (8)

where

$\begin{split}\log G({\bf f}_{u},{\bf y})&=\log\mathcal{N}\left({\bf y}|{\bf
K}_{xu}{\bf K}_{uu}^{-1}{\bf f}_{u},\sigma^{2}{\bf I}\right)\\\
&-\frac{1}{2\sigma^{2}}\text{tr}\left\\{{\bf K}_{xx}-{\bf K}_{xu}{\bf
K}_{uu}^{-1}{\bf K}_{ux}\right\\},\end{split}$ (9)

with $\text{tr}\left\\{\cdot\right\\}$ being the trace operator. To optimize
(8) w.r.t the variational parameters, [10] proposes to apply the Jensen’s
inequality to move the logarithm out of the integral and, thus, to cancel
distribution $q({\bf f}_{u})$ to obtain a collapsed version of this bound,
i.e.,

$\begin{split}F_{V}({\bf U})&=\log\mathcal{N}\left({\bf y}|{\bf
0},\sigma^{2}{\bf I}+{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf K}_{ux}\right)\\\
&-\frac{1}{2\sigma^{2}}\text{tr}\left\\{{\bf K}_{xx}-{\bf K}_{xu}{\bf
K}_{uu}^{-1}{\bf K}_{ux}\right\\}.\end{split}$ (10)

This new bound can be maximized to obtain the optimal values of ${\bf U}$, as
well as other model hyperparameters, such as, the kernel parameters and the
noise variance.

Finally, to obtain the optimal value of $q({\bf f}_{u})$, we can differentiate
(8) respect to $q({\bf f}_{u})$ and set it to zero. This leads to the fact
that the optimal distribution is proportional to $\mathcal{N}\left({\bf
y}|{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf f}_{u},\sigma^{2}{\bf I}\right)p({\bf
f}_{u})$ and identifying terms, we see that the optimal expression for $q({\bf
f}_{u})$ is

$\begin{split}q({\bf f}_{u})&=\mathcal{N}\left({\bf
f}_{u}|{\boldsymbol{\mu}},{\bf A}\right)\\\
{\boldsymbol{\mu}}&=\sigma^{-2}{\bf K}_{uu}{\bf B}{\bf K}_{ux}{\bf y}\\\ {\bf
A}&={\bf K}_{uu}{\bf B}{\bf K}_{uu},\end{split}$ (11)

where ${\bf B}=\left({\bf K}_{uu}+\sigma^{-2}{\bf K}_{ux}{\bf
K}_{xu}\right)^{-1}$. Eq. (11) allows to obtain ${\boldsymbol{\mu}}$ and ${\bf
A}$ to completely define the predictive posterior (7). Note that both the
optimization of (10) and the computation of (11) have a computational burden
of $\mathcal{O}(NM^{2})$. This provides fundamental advantages that allow
novel and efficient adaptive GP variants, which will become apparent next.

##  III Adaptive Online VSGPs

In this section we consider that we are working in a time-dependent framework,
so in each time instant $t$ a new data pair $({\bf x}_{t},y_{t})$ arrives. To
be able to efficiently deal with this new scenario, we next reformulate the
previous VSGP model to make it able to: (1) forget older samples (_adaptive
capabilities_); and, (2) be updated with the information of each new sample
reusing the previous time instant model solution (_online learning_).

###  III-A Adaptive VSGP

To endow the VSGP formulation with adaptive properties, we can modify the
marginal likelihood bound (8) by including a forgetting factor $\lambda$
(similar to that of the Recursive Least Squares filter [2]) so that older data
weigh less in the bound. Then, we can maximize this $\lambda$ dependent bound
to obtain the model hyperparameters, as well as the new parameters of the
$q({\bf f}_{u})$ distribution.

For this purpose, as the only data dependent term in $F_{V}({\bf U},q({\bf
f}_{u}))$ is $\log G({\bf f}_{u},{\bf y})$, we include $\lambda$ into this
term and we start considering time instant $t$, so only data from initial
sample at $t^{\prime}=1$ to $t$ are available, then,

$\begin{split}&\log G^{\lambda}({\bf f}_{u},{\bf
y})=\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\log G({\bf
f}_{u},y_{t^{\prime}})\\\
&=\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\log\mathcal{N}\left(y_{t^{\prime}}|{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf f}_{u},\sigma^{2}\right)\\\
&-\frac{1}{2\sigma^{2}}\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right),\end{split}$
(12)

where ${\bf y}=[y_{1},\ldots,y_{t}]^{\top}$, ${\bf k}_{t^{\prime}u}={\bf
k}_{{\bf x}_{t^{\prime}},u}$ is a vector constructed with the kernel products
between sample ${\bf x}_{t^{\prime}}$ and ${\bf U}$, and
$k_{t^{\prime}t^{\prime}}$ is the kernel of ${\bf x}_{t^{\prime}}$ with
itself. After some manipulations (see Appendix A), we can transform (12) into

$\begin{split}\log G^{\lambda}({\bf f}_{u},{\bf
y})&\propto\log\mathcal{N}\left({\bf y}|{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf
f}_{u},\sigma^{2}{{\boldsymbol{\Lambda}}}^{-1}\right)\\\
&-\frac{1}{2}\sum_{t^{\prime}=1}^{t}(\lambda^{t-t^{\prime}}-1)\log(2\pi\sigma^{2})\\\
&+\frac{1}{2\sigma^{2}}\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right),\end{split}$
(13)

where ${\boldsymbol{\Lambda}}$ is a $t\times t$ diagonal matrix with
${\boldsymbol{\Lambda}}_{t^{\prime}t^{\prime}}=\lambda^{t-t^{\prime}}$ and now
${\bf K}_{xu}$ is constructed with the kernel products between samples ${\bf
x}_{1},\ldots,{\bf x}_{t}$ and the inducing points. The $\lambda$ dependent
$F_{V}({\bf U})$ bound is thus given by:

$\begin{split}F_{V}^{\lambda}({\bf U})&=\log\mathcal{N}\left({\bf y}|{\bf
0},\sigma^{2}{\boldsymbol{\Lambda}}^{-1}+{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf
K}_{ux}\right)\\\
&-\frac{1}{2\sigma^{2}}\sum_{t^{\prime}=1}^{t}(\lambda^{t-t^{\prime}}-1)\log(2\pi\sigma^{2})\\\
&+\frac{1}{2\sigma^{2}}\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right).\end{split}$
(14)

The effect of the introduction of parameter $\lambda$ can be interpreted in
Eq. (14) as an enhancement of the noise term that increases when $t^{\prime}$
decreases. That is, an increasing uncertainty is attributed to the older
training samples.

To complete this model, we obtain the adaptive version of the variational
distribution of the inducing points, $q_{\lambda}({\bf f}_{u})$. For this
purpose, its optimum value has to be proportional to $\mathcal{N}\left({\bf
y}|{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf
f}_{u},\sigma^{2}{\boldsymbol{\Lambda}}^{-1}\right)p({\bf f}_{u})$ and
therefore its mean ${\boldsymbol{\mu}}_{\lambda}$ and covariance ${\bf
A}_{\lambda}$ are given by:

$\begin{split}{\boldsymbol{\mu}}_{\lambda}&=\sigma^{-2}{\bf K}_{uu}{{\bf
B}_{\lambda}}{\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf y}\\\ {\bf
A}_{\lambda}&={\bf K}_{uu}{{\bf B}_{\lambda}}{\bf K}_{uu},\end{split}$ (15)

with ${\bf B}_{\lambda}=\left({\bf K}_{uu}+\sigma^{-2}{\bf
K}_{ux}{\boldsymbol{\Lambda}}{\bf K}_{xu}\right)^{-1}$. Finally, replacing the
optimum values of ${\boldsymbol{\mu}}_{\lambda}$ and ${\bf A}_{\lambda}$ into
the predictive distribution of the VSGP (see (7)), we can obtain the mean and
variance of adaptive predictive distribution as

$\begin{split}m_{\lambda,*}&=\sigma^{-2}{\bf k}_{*u}{\bf B}_{\lambda}{\bf
K}_{ux}{\boldsymbol{\Lambda}}{\bf y}\\\ v_{\lambda,*}&=k_{**}+{\bf
k}_{*u}\left({\bf B}_{\lambda}-{\bf K}^{-1}_{uu}\right){\bf
k}_{u*}.\end{split}$ (16)

###  III-B Online update of the Adaptive VSGP

With the above model, we are able to train a GP for non-stationary
environments so that the model will learn the distribution of the most recent
samples. However, to make this model work efficiently in such scenarios, we
need to be able to efficiently update its solution as new samples arrive.

To carry out this online learning there are two possible working modes.
Firstly, a fast implementation where noise and kernel parameters are
considered fixed and, for every new data, we update the predictive
distribution and, if needed, the inducing point set. Secondly, a more accurate
solution where all model parameters are updated to efficiently track the data
distribution changes. Next, these approaches are explained in detail.

####  III-B1 Model update without inference over the model parameters

If we consider both the noise and kernel parameters are fixed, when new
training sample arrives, $({\bf x}_{t+1},y_{t+1})$, we only have to update the
mean and variance of the adaptive predictive distribution (see Eq. (16)).
Considering these values were $(m_{\lambda,*})_{t}$ and $(v_{\lambda,*})_{t}$,
the new values at $t+1$ will be given by:

$\begin{split}(m_{\lambda,*})_{t+1}&=\sigma^{-2}{\bf k}_{*u}\left({\bf
B}_{\lambda}\right)_{t+1}\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
y}\right)_{t+1}\\\ (v_{\lambda,*})_{t+1}&=k_{**}+{\bf k}_{*u}\left(\left({\bf
B}_{\lambda}\right)_{t+1}-{\bf K}^{-1}_{uu}\right){\bf k}_{u*},\end{split}$
(17)

where

$\begin{split}\left({\bf B}_{\lambda}\right)_{t+1}&=\left({\bf
K}_{uu}+\sigma^{-2}\left(\lambda\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
K}_{xu}\right)_{t}+{\bf k}_{t+1,u}^{\top}{\bf k}_{t+1,u}\right)\right)^{-1}\\\
\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
y}\right)_{t+1}&=\lambda\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
y}\right)_{t}+{\bf k}_{t+1,u}^{\top}y_{t+1},\end{split}$ (18)

where ${\bf k}_{t+1,u}$ is a row vector constructed with the kernel dot
products between sample ${\bf x}_{t+1}$ and ${\bf U}$.

Note that as we are working with an adaptive algorithm, the model can forget
the past without explicitly removing the old data. Besides, the data dependent
terms, $\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf y}\right)_{t}$ and
$\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf K}_{xu}\right)_{t}$, do not
increase its size when new data are added, so the predictive distribution
complexity does not increase with each new data. However, if we want to limit
the memory needed to store past samples (necessary, for example, to recompute
${\bf K}_{ux}$ when the inducing points are updated), it is advisable to
remove old data and keep a window of length $T$. This can be easily done
during the addition of the data using these update rules, which leads to an
expression alternative to those in Eq. (18):

$\begin{split}\left({\bf B}_{\lambda}\right)_{t+1}&=\left({\bf
K}_{uu}+\sigma^{-2}\left(\lambda\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
K}_{xu}\right)_{t}+{\bf k}_{t+1,u}^{\top}{\bf k}_{t+1,u}\right)\right.\\\
&-\left.\lambda^{T}{\bf k}_{t-T,u}^{\top}{\bf k}_{t-T,u}\right)^{-1}\\\
\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
y}\right)_{t+1}&=\lambda\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
y}\right)_{t}+{\bf k}_{t+1,u}^{\top}y_{t+1}-\lambda{\bf
k}_{t-T,u}^{\top}y_{t-T}.\end{split}$ (19)

If the values of $\left({\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
K}_{xu}\right)_{t}$ from the previous iteration are saved, the complexity of
these operations is $O(M^{3})$ and $O(M)$, respectively.

In nonstationary scenarios, it is expected that as new data arrive, their
distribution changes and we also need to update the positions of the inducing
points to be able to represent the support of the new distribution. To avoid
applying a computationally expensive inference process to update the position
of all inducing points, when the inducing set is not representative enough, we
propose to directly add the new data as a new inducing point. Besides, to
limit the length of this set, we can remove the less representative inducing
points.

To analyze the representativeness of the inducing set, we can use as criterion
the regularization term of the adaptive collapsed marginal likelihood (14),
i.e.,

$R_{\rm
tot}=\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right).$ (20)

This quantity measures the $\lambda$ weighed error generated by predicting the
training latent values from the inducing points. In fact, the non-adaptive
version of this criterion has been already used by [22, 23] for the selection
of inducing points from the training data. So, we can analyze this quantity in
each iteration, and in the case is larger than a given threshold
($R_{th,tot}$), we can decide to include the observation ${\bf x}_{t+1}$ as a
new inducing point. Besides, according to this criterion, we can consider that
the relevance of each inducing point, $R_{m}$ for $m=1,\ldots,M$, is given by

$R_{m}=\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}{\bf k}_{t^{\prime}m}{\bf
K}_{mm}^{-1}{\bf k}_{mt^{\prime}},$ (21)

such that $R_{\rm
tot}=\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}k_{t^{\prime}t^{\prime}}-\sum_{m=1}^{M}R_{m}$.
So, we can now set a relevance threshold ($R_{th}$) over $R_{m}$ to detect the
useless inducing points and remove them.

Finally, to complete this algorithm, we again need to update the predictive
distribution when the inducing point set is modified. For this purpose, given
that there are $k$ elements in the inducing point set, we add a new inducing
point ${\bf u}_{k+1}$ to this set, so the new predictive mean and variance
are:

$\begin{split}(m_{\lambda,*})_{k+1}&=\sigma^{-2}\left({\bf
k}_{*u}\right)_{k+1}\left({\bf B}_{\lambda}\right)_{k+1}\left({\bf
K}_{ux}\right)_{k+1}{\boldsymbol{\Lambda}}{\bf y}\\\
(v_{\lambda,*})_{k+1}&=k_{**}+\left({\bf k}_{*u}\right)_{k+1}\left(\left({\bf
B}_{\lambda}\right)_{k+1}-\left({\bf
K}^{-1}_{uu}\right)_{k+1}\right)\left({\bf k}_{u*}\right)_{k+1}.\end{split}$
(22)

Updating $\left({\bf k}_{*u}\right)_{k+1}$ and $\left({\bf
K}_{ux}\right)_{k+1}$ is straightforward since it only implies adding a new
element and row to this vector and matrix, respectively. The update of
$\left({\bf K}^{-1}_{uu}\right)_{k+1}$ can be efficiently computed from
$\left({\bf K}^{-1}_{uu}\right)_{k}$ using the properties of the block matrix
inversion111https://en.wikipedia.org/wiki/Block_matrix#Block_matrix_inversion
, if we consider that $\left({\bf K}_{uu}\right)_{k}$ is an update of
$\left({\bf K}_{uu}\right)_{k+1}$ with a new row and column. And a similar
procedure can be applied to $\left({\bf B}_{\lambda}\right)_{k+1}$ taking into
account that

$\begin{split}\left({\bf B}_{\lambda}\right)_{k+1}&=\left({\bf
K}_{uu}+\sigma^{-2}{\bf K}_{ux}{\boldsymbol{\Lambda}}{\bf
K}_{xu}\right)^{-1}\\\ &=\left(\begin{array}[]{ll}\left({\bf
B}_{\lambda}^{-1}\right)_{k}&\left({\bf b}_{\lambda}\right)_{k+1}\\\
\left({\bf
b}_{\lambda}^{\top}\right)_{k+1}&\left(b_{\lambda}\right)_{k+1}\end{array}\right)^{-1},\end{split}$
(23)

where

$\begin{split}\left({\bf b}_{\lambda}\right)_{k+1}&={\bf
k}_{u,k+1}+\sigma^{-2}\left({\bf K}_{ux}\right)_{k}{\boldsymbol{\Lambda}}{\bf
k}_{x,k+1}\\\ \left(b_{\lambda}\right)_{k+1}&=k_{k+1,k+1}+\sigma^{-2}{\bf
k}_{k+1,x}{\boldsymbol{\Lambda}}{\bf k}_{x,k+1}.\end{split}$ (24)

This way, the addition of a new inducing point has a computational cost of
$\mathcal{O}(M^{2})$.

Note that for the computation of $\left({\bf K}_{ux}\right)_{k+1}$, we need to
have stored the matrix ${\bf X}$ with all previous training data, although due
to $\lambda$ factor, only the most recent data have influence over the
solution; so, as we have argued before, we can keep a window of $T$ data to
reduce the computational requirements and without any performance degradation.

When we need to remove an irrelevant inducing point for being useless or to
keep the inducing point set with $M$ elements, we will have to incur a higher
computational burden. This is due to the fact that the inducing point to be
removed occupies an arbitrary position in the matrix; therefore, a rank 1
update is not possible and we have to compute the inversions of $\left({\bf
B}_{\lambda}\right)_{k}$ and $\left({\bf K}_{uu}\right)_{k}^{-1}$ from scratch
with a computational cost of $\mathcal{O}(M^{3})$ in the worst case that
$k=M$.

A summary of this fast model update, called fast-AGP, is included in Algorithm
1, where the model at the previous instant or a model initialized with some
data history is considered as input.

$\rm{(fast-AGP)}_{t+1}$ = Update_fast-AGP $\left(\rm{(fast-AGP)}_{t},\right.$

$\left.({\bf x}_{t+1},y_{t+1}),T,R_{th},R_{th,tot}\right)$

Inputs: Previous model (fast-AGP)t, new data sample ${\bf x}_{t+1}$,
$y_{t+1}$, length of data window (optional) $T$, and relevance thresholds to
update the inducing point set ($R_{th}$ and $R_{th,tot}$).

// Update predictive distribution with new data

Compute new $(m_{\lambda,*})_{t+1}$ and $(v_{\lambda,*})_{t+1}$ with Eq. (17)
and (18).

If we keep a data window of length $T$, use Eq. (17) and (19).

// Add a new inducing point if it is needed

If $R_{\rm tot}>R_{th,tot}$ (see Eq. (20)), set ${\bf u}_{k+1}={\bf x}_{t+1}$
and update the predictive distribution with Eq. (22) -(24).

// Remove useless inducing points

For $m=1,\ldots,k+1$, compute (21) and remove inducing points with
$R_{m}<R_{th}$ or those with the lowest $R_{m}$ so that the set length is no
larger than $M$.

Then, if needed, update $\left({\bf B}_{\lambda}\right)$ and ${\bf
K}_{uu}^{-1}$and, later, $(m_{\lambda,*})_{k+1}$ and $(v_{\lambda,*})_{k+1}$.

// Return model fast-AGP model at $t+1$

Return $\rm{(fast-AGP)}_{t+1}$

Algorithm 1 Fast Adaptive VSGP (fast-AGP)

####  III-B2 Model update with inference over the model parameters

We can improve the performance of the previous algorithm by updating all the
model parameters (kernel parameters, noise variance, and the position of the
new inducing point) in each iteration.

To optimize these parameters, we need to maximize Eq. (14), which leads to a
computational burden of $\mathcal{O}(TM^{2})$, considering we limit the
training data to a window of length $T$. In this case, since we will have to
recalculate the kernel matrices every time we update the kernel parameters, we
cannot save computation in the updates with rank one updates, so the
predictive distribution will have to be calculated from scratch with a cost of
$\mathcal{O}(TM^{2})$. A summary of this approach is included in Algorithm 2.

$\rm{(AGP)}_{t+1}$ = Update_AGP $\left(\rm{(AGP)}_{t},({\bf
x}_{t+1},y_{t+1}),T,R_{th}\right)$

Inputs: previous model (AGP)t, new data sample ${\bf x}_{t+1}$, $y_{t+1}$,
length of data window (optional) $T$ and relevance threshold for inducing
points ($R_{th}$).

// Update the data window

Include $({\bf x}_{t+1},y_{t+1})$ to the model data set and remove $({\bf
x}_{t-T},y_{t-T})$.

// Remove useless inducing points and limit set size to $M-1$

For $m=1,\ldots,k+1$, compute (21) and remove inducing points with
$R_{m}<R_{th}$ or those with the lowest $R_{m}$ so that the set length is no
larger than $M-1$.

// Add a new inducing point

Set ${\bf u}_{k+1}={\bf x}_{t+1}$

// Update model parameters

Compute new values of $(\sigma^{2})_{t+1}$, $({\bf u}_{k+1})_{t+1}$ and kernel
parameters $(\theta)_{t+1}$ by maximizing Eq. (14).

// Update predictive distribution with new parameters

Compute the predictive distribution at $t+1$ with (16).

// Return model AGP model at $t+1$

Return $\rm{(AGP)}_{t+1}$

Algorithm 2 Adaptive VSGP (AGP)

Note that even though we need to compute the predictive function from scratch,
the computational burden of this approach is significantly lower than one of
the updates of a standard VSGP starting from the parameters optimized in the
last time instant. This is because the number of parameters to update is
significantly lower (we only update the position of one inducing point, not
the whole set), so the convergence is likely to be achieved in a lower number
of inference iterations; even, as we will show in the experimental section, a
single iteration per time instant. Moreover, this procedure is subject to
further simplifications. For example, the full parameter update may not be
needed every time and, if we need to add a new inducing point, its point
position can be set to that of a new data point, i.e., ${\bf u}_{k+1}={\bf
x}_{t+1}$; this way, no inference is needed in these iterations and the
computational cost of the parameter update is lowered further.

##  IV Experiments

In this section, we evaluate the performance of the proposed adaptive models
against several reference methods. For this purpose, we first consider a
synthetic unidimensional problem where we can easily control the non-
stationary of the environment to evaluate the advantages and drawbacks of the
different methods. Second, we test these methods on a load forecasting
application.

###  IV-A Experimental setup

To analyze the performance of the proposed model, we have considered the two
versions introduced:

  * •

Adaptive VSGP without inference (fast-AGP) (see Algorithm 1), where only the
set of inducing points is updated to represent the new data space, but the
kernel parameters and the model noise are considered constant. Besides, the
inducing point set is updated without any inference, since new data are
included as new inducing points, and the useless ones are efficiently removed.

  * •

Adaptive VSGP with inference (AGP) (see Algorithm 2), where all model
parameters are updated using a single inference iteration every time new data
arrives.

In both cases, for automatic update of the set of inducing points,
$R_{th,tot}$ at time instant $t$ is set as
$\frac{1}{T}\sum_{t^{\prime}=t-T}^{t}\lambda^{t-t^{\prime}}k_{t^{\prime}t^{\prime}}$
and $R_{th}=r_{th}\max_{m}(R_{m})$ where $r_{th}$ is a value which has been
arbitrarily set to $10^{-4}$ in all the experiments below.

To analyze the performance of the proposed methods, we use the following
reference approaches:

  * •

An adaptive VSGP where we use VSI (Variational Stochastic Inference) (AGP-VSI)
to update the model parameters in each iteration (see Appendix B for details
of this formulation). For this model, we also consider a window of length $T$
so that only these $T$ samples are used for each data chunk. To make this
algorithm converge, in each time instant, we have needed to use $50$
iterations.

  * •

Windowed VSGP (w-VSGP): This is the standard variational SGP of [10] using the
implementation of [24], where we use a sliding window of length $T$ to
simulate the adaptive scenario. As in the previous model, for each new data,
we need to use $50$ inference iterations to adjust the model parameters.

  * •

Online SGP (OSGP): this is the online approach proposed in [16]. In this case,
instead of using the original implementation of the
authors222https://github.com/thangbui/streaming_sparse_gp, for comparison
purposes, we have preferred to use its PyTorch
implementation333https://github.com/wjmaddox/online_gp/blob/main/online_gp/models/streaming_sgpr.py.
Here, we have used the learning rate parameters recommended by the authors
($lr=0.1$ for the noise and kernel parameters and $lr=0.2$ for the inducing
points) with a Cosine Annealing optimizer. Here, as for the proposed methods,
the online inference of the model uses a single iteration.

  * •

The Kernel Interpolation for Scalable Online GP (WISKI) approach [13] using
the authors implementation444https://github.com/wjmaddox/online_gp. For this
model, since the inducing points are not updated by the model, we have had to
predefine a grid over the input space to set them. Besides, as the definition
of this grid raises computational problems in high-dimensional spaces, when
the input dimension is larger than two, we have projected the data into a two-
dimensional space with a one-layer neural network (as the authors suggest);
the weights of this net are learned with the remaining parameters of the
model.

For all the models and experiments, we have used a square exponential kernel
with the width and variance as parameters. To initialize the different models,
we have considered a window of length $T$ and we have used the first $T$
samples to train a first model (in batch) during $200$ iterations. For the
inference of all these models (either during their initialization or, if used,
during online learning), except for OSGP which has used the optimization
proposed by the authors, we have used an Adam algorithm with a learning rate
of $0.05$. For adaptive algorithms, the value of the forgetting factor
($\lambda$) was set as a function of the window length so that
$\lambda^{T}=0.1$.

###  IV-B Toy example

As a first adaptive scenario, we have generated a synthetic but challenging
data set. In particular, we have generated a sinusoidal signal:

$y(t)=A_{t}{\rm sin}(2\pi ft)+\epsilon_{t}$

with $500$ samples uniformly distributed in the interval from $0$ to $5$. The
first $300$ samples (in the interval $[0,3]$) use $f=4/2\pi$ with an
increasing amplitude, $A_{t}$, from $0.5$ to $2$. The remaining $200$ samples
consider $A_{t}=2$ and a frequency of $8/2\pi$. In both cases, these
sinusoidal signals have been contaminated with a gaussian noise of zero mean
and standard deviation of $0.2$. In Figure 1 we include an example of this
data set.

Figure 1: Example of the synthetic problem. Black points have been used for
the model initialization and blue samples for online learning.

For all models under study, we have used a data window of $T=100$ samples, so
we have trained a first model with first $100$, and we have used $M=10$
inducing points; for WISKI model as we have to predefine an inducing points
grid we have set $40$ inducing points in the interval from $0$ to $5$.
Finally, for the adaptive models, $\lambda$ has been set to $\lambda=0.97724$,
so $\lambda^{100}\approx 0.1$.

To analyze the performance of the different methods, we have computed the Mean
Square Error (MSE) and the percentage of samples whose prediction error is
less than $\pm 2\left(v_{\lambda,*}+\sigma^{2}\right)$ (see Eqs. (17) and
(22)), known as $95\%$ confidence interval (95% CI) when the model at the time
instant $t$ (i.e, trained with samples from $t-T+1$ to $t$) predicts the
output for sample at $t+1$. Besides, to obtain representative results, we have
generated $1000$ different runs of the dataset and averaged the results for
all the runs. In Table I we include these results for all the methods
including their Training Times (TT). Besides, in Figures 2 and 3 we include,
respectively, the evolution of the MSE and the predictive distribution during
the online learning of the algorithms.

| w-SVGP | OSGP | WISKI | AGP-VSI | fast-AGP | AGP
---|---|---|---|---|---|---
MSE | 1.0989 | 0.1333 | 0.4507 | 0.5211 | 0.1960 | 0.0630
$95\%$ CI | 93.10 | 94.21 | 75.28 | 93.71 | 71.68 | 95.02
Tr. Time | 3266.2 | 299.5 | 120.4 | 3022.0 | 23.1 | 108.6
Table I: Comparison of the different models under study in the synthetic
dataset.
Figure 2: Evolution of the MSE during the online learning in the synthetic
dataset. Figure 3: Predictive distribution for the different methods under
study along the online learning in the synthetic dataset.
Figure 4: Evolution of the AGP MSE for different $\lambda$ values in the
synthetic dataset.

Table I clearly shows the algorithm presenting the best performance is AGP,
both in terms of MSE and estimation of the $95\%$ CI followed by OSGP and
fast-AGP, being the latter a good choice if the computational burden is
critical. Analyzing Figures 2 and 3, we can check as all the models work
properly in the first part of the signal (when the frequency is low and only
have to follow a change of amplitude), but when the signal frequency changes
most all the methods under study, but AGP, fail. For example, the fast-AGP is
not able to track the signal because it does not update the kernel and noise
parameters, so it is not able to follow the signal variations. w-SVGP and AGP-
VSI (mainly the former) fail because they need more than $50$ iterations, at
each time stamp, to track the signal changes. WISKI tries to track the signal
variations, but it struggles to capture the amplitude, and their estimations
are delayed (see Fig. 3 bottom right); the OSGP is the only reference
algorithm able to track the signal changes, but the confidence interval
estimation grows along the time being every time more pessimistic.

To complete this analysis, Figure 4 shows the evolution of the MSE in AGP for
different values of $\lambda$ focusing this analysis on the region around the
sinusoidal frequency change (from $2.5$ to $4$). We set up parameter
$\lambda\leq 1$ so the most remote sample of the window has an attenuation
$\lambda^{T}=10^{-v}$, where $v$ takes the values $0,0.25,0.5,0.75,1,2$ and
$3$.

As we can see, the introduction of adaptive learning is critical to fastly
track the signal transition and to match the high frequency of the new
sinusoidal. If $\lambda=1$ (no forgetting is applied), the AGP model has a
performance similar to the baselines. Setting $\lambda$ to any of the other
values, the algorithm adapts the model to the new signal frequency. The
tracking error decreases as $\lambda$ is reduced from $1$ to $0.933$ (see the
MSE differences in the interval from $t=3.2$ to $t=3.4$).

###  IV-C Load forecast

The dataset used in this article is the archived data ISO-New England555The
dataset is publicly available at https://www.iso-
ne.com/isoexpress/web/reports/load-and-demand/-/tree/zone-info., from 2011 to
2020, including real-time demand of each one of the eigth load zones of New
England. ISO-New England Inc. is a Regional Transmission Organization (RTO),
an Independent System Operator that operates the bulk power electric system
and transmission lines serving the territory of the six New England states:
Connecticut (CT), Maine (ME), Northeast Massachusetts and Boston (NEMA), New
Hampshire (NH), Rhode Island (RI), Southeast Massachusetts (SEMA), Vermont
(VT), and Western/Central Massachusetts (WCMA).

The input patterns consist of 24 hourly power measures of a single day to
predict every hour of the next day, and we have used an observation window of
90 days ($T=90$). Assuming that the most remote sample of the window has to be
attenuated by a factor of 0.1, this leads to $\lambda=0.9747$. The initial
model has also been trained with $90$ data, and we have used $M=20$ inducing
points, but the fast-AGP and AGP algorithms select between 8 and 10 inducing
points in all scenarios.

As the model is using $24$ input variables, the WISKI algorithm cannot work in
such a high space, so we include a preprocessing step (as suggested in [13])
to project the data into a two-dimensional space with a one layer neural
network consisting in a set of linear weights, a batch normalization and an
output hyperbolic tangent to transform the data into the $[-1,1]$ scale. In
this way, the inducing point grid is predefined in this dimensional space with
$400$ points equally spaced on the grid $[-1,1]\times[-1,1]$. Using less dense
grids leads to poor results.

The OSGP reference approach has to be removed from this experiment since it
has shown convergence errors in most regions, providing poor performance.
These converge problems are due to the expression of its inference bound,
which has to invert a matrix defined by the subtraction of two terms of a
similar order, causing its cholesky decomposition, despite being computed with
a robust implementation, to fail in many realizations.

Figure 5: MAPE comparison of all algorithms for the power load experiment in
the regions of Connecticut and West/Central Massachussets.

Figure 5 shows the mean absolute percentage of error (MAPE) of the 24 hour
prediction achieved by all methods in the regions of Connecticut (CT) and
West/Central Massachusetts. In these results, we also include the performance
of the _persistence_ approach, that is, an estimation given by the load
consumption on the previous day at the same hour. The best MAPEs correspond to
the fast-AGP and AGP algorithms, the WISKI algorithm shows a performance
similar to that of the simple persistence, and AGP-VSI and w-SVGP show
performances close to those of the AGP and fast-AGP ones. The results for all
areas are summarized in Table II.

|  | Persist. | w-SVGP | WISKI | AGP-VSI | fast-AGP | AGP
---|---|---|---|---|---|---|---
| MSE | 0.1282 | 0.0982 | 0.1202 | 0.0876 | 0.0876 | 0.0933
CT |  $95\%$ CI | – | 89.91 | 90.22 | 89.82 | 90.70 | 96.16
| Tr. Time | – | 5262.4 | 280.4 | 5459.9 | 36.6 | 196.0
| MSE | 0.0092 | 0.0062 | 0.0090 | 0.0112 | 0.0058 | 0.0066
ME |  $95\%$ CI | – | 91.96 | 91.44 | 92.78 | 93.88 | 96.70
| Tr. Time | – | 5535.3 | 278.8 | 5519.2 | 60.5 | 198.5
| MSE | 0.0802 | 0.0593 | 0.0706 | 0.0535 | 0.0490 | 0.0533
NEMA |  $95\%$ CI | – | 89.74 | 89.14 | 90.01 | 89.73 | 96.20
| Tr. Time | – | 5178.2 | 280.2 | 5162.3 | 38.7 | 197.7
| MSE | 0.0160 | 0.0107 | 0.0148 | 0.0121 | 0.0097 | 0.0107
NH |  $95\%$ CI | – | 91.50 | 90.46 | 91.21 | 92.38 | 96.63
| Tr. Time | – | 5209.2 | 284.9 | 5107.9 | 54.6 | 196.8
| MSE | 0.0088 | 0.0099 | 0.0100 | 0.0134 | 0.0061 | 0.0069
RI |  $95\%$ CI | – | 90.75 | 89.42 | 92.26 | 88.00 | 96.35
| Tr. Time | – | 5180.2 | 289.9 | 5154.2 | 65.1 | 196.3
| MSE | 0.0347 | 0.0280 | 0.0375 | 0.0281 | 0.0247 | 0.0263
SEMA |  $95\%$ CI | – | 90.19 | 91.33 | 89.65 | 89.15 | 95.84
| Tr. Time | – | 5240.5 | 286.6 | 5105.7 | 45.6 | 196.2
| MSE | 0.0056 | 0.0042 | 0.0049 | 0.0044 | 0.0038 | 0.0040
VT |  $95\%$ CI | – | 92.58 | 92.58 | 93.77 | 97.14 | 97.34
| Tr. Time | – | 5128.1 | 284.7 | 5186.4 | 54.2 | 196.1
| MSE | 0.0360 | 0.0255 | 0.0309 | 0.0240 | 0.0237 | 0.0248
WCMA |  $95\%$ CI | – | 89.84 | 92.00 | 89.90 | 93.15 | 96.84
| Tr. Time | – | 5042.0 | 287.6 | 5290.0 | 37.9 | 194.4
| MSE | 0.0398 | 0.0303 | 0.0372 | 0.0293 | 0.0263 | 0.0282
Average |  $95\%$ CI | – | 90.81 | 90.82 | 91.17 | 91.77 | 96.51
| Tr. Time | – | 5222.0 | 284.1 | 5248.2 | 49.2 | 196.5
Table II: Comparison of the different models under study for the different
load forecasting regions.

Analyzing the detailed results of Table II, we can corroborate that fast-AGP
and AGP show the best performance. Although fast-AGP gets the lowest MSE, AGP
is the only one able to obtain an accurate $95\%$ CI estimation. This
advantage of AGP in terms of the mean estimation can be due to a slight
overfitting of fast-AGP.

Regarding reference methods w-SVGP, AGP-VSI and WISKI present an acceptable
performance, but if we analyze the MAPE results (see Figure 5), we see as
WISKI is quite unstable and in some time hours present an error higher than
the Persistence approach. We think this punctual performance degradation is
probably due to the fact that the predefined grid of inducing points is
insufficient or is badly placed for some hour estimations.

However, when we analyze the confidence interval estimation, we observe that
AGP is the only one that systematically provides an accurate estimate (perhaps
slightly conservative, being 1% above the exact estimate), but the remaining
methods are around 90%.

Finally, if we analyze the computational cost, we corroborate as w-SGP and
AGP-VSI are not affordable, since their computational burden is $25$ times
higher than AGP and 100 times higher than fast-AGP. WISKI presents a
computational burden similar to the AGP one, but fast-AGP is the fastest with
$5.7$ times lower computational time than WISKI.

In an environment such as energy load forecast, where kernel and noise
parameters do not need to be recalculated quickly, fast-AGP can be a good
choice. However, if an accurate confidence interval estimation is needed, the
AGP can be a better option. Even for these semi-stationary environments both
models can be combined (for instance running AGP every ten days stamps to
track the model parameters) and obtain better variance estimation and higher
speed.

##  V Conclusions

Adaptive algorithms are used in signal processing and in these applications
where the environment is non-stationary. We introduced an algorithm that
implements an adaptive SGP for regression based on the VSGP introduced in
[10], where a variational approach is used to maximize a variational lower
bound of the marginal likelihood of the training regressors.

To endow the VSGP with adaptive properties, the components of the log-
likelihood are weighted with a factor that exponentially decays with time, so
the information included in the log-likelihood is conveniently forgotten in a
way similar to what is done in the RLS algorithm. To make this model
efficient, the algorithm used for the inference is formulated in two
alternative ways that make it adequate for online learning. In one of the
approaches, the noise and kernel parameters are frozen after initial training,
and only the inducing points are adapted, to obtain the fast-AGP, which
combines good performance and low computational burden. The second one, called
AGP, updates all parameters in an online way, thus improving its capabilities
by increasing the computational burden. In all cases, the computational burden
is lower than the state-of-the-art approaches. To limit the number of inducing
points, we propose a relevance criterion, based on the change of the marginal
likelihood, to automatically add a new sample as an inducing point or remove
an existing inducing point.

The performance of the algorithm has been tested in two problems. The first
one is an artificial problem where a sinusoidal signal has to be tracked (with
one-step ahead prediction). This signal shows a sudden frequency change that
is properly tracked by the proposed algorithms but not by the state-of-the-art
algorithm used in the comparison. The fast-AGP shows the second best
performance, while the AGP shows the best performance. A test with different
values of the forgetting factor $\lambda$ shows the importance of the
forgetting factor, and it also shows the robustness of the algorithm regarding
the parameter choice.

A second experiment consists of electric load forecasting with the ISO New
England load databases. The results show that both AGP and fast-AGP show the
best performance in mean. Besides, the AGP is the only one that achieves an
accurate 95% CI, while the rest of the approaches show a 91% of samples inside
the estimated 95% CI. This is important in these applications, where it is
useful to know what is the confidence interval of the prediction to decide
whether to trust the prediction or not.

##  Appendix A Derivation of the $\lambda$-dependent log likelihood

In this appendix, we derive the different expressions of Section III-A to
demonstrate the solution for the adaptive version of the VSGP formulation.

Let’s start including $\lambda$ into $\log G({\bf f}_{u},{\bf y})$

$\begin{split}&\log G^{\lambda}({\bf f}_{u},{\bf
y})=\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\log G({\bf
f}_{u},y_{t^{\prime}})\\\
=&\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\log\mathcal{N}\left(y_{t^{\prime}}|{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf f}_{u},\sigma^{2}\right)\\\
&-\sum_{t^{\prime}=1}^{t}\frac{\lambda^{t-t^{\prime}}}{2\sigma^{2}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right)\end{split}$ (25)

and, now, focusing on the term of the first sum and defining
$\alpha_{t^{\prime}}={\bf k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf f}_{u}$, we
can be reformulated this term as:

$\begin{split}&\lambda^{t-t^{\prime}}\log\mathcal{N}\left(y_{t^{\prime}}|\alpha_{t^{\prime}},\sigma^{2}\right)\\\
&=-\frac{\lambda^{t-t^{\prime}}}{2}\left(\log\left(2\pi\sigma^{2}\right)+\frac{1}{\sigma}\left(y_{t^{\prime}}-\alpha_{t^{\prime}}\right)\right)\\\
&=-\frac{1}{2}\lambda^{t-t^{\prime}}\log\left(2\pi\sigma^{2}\right)+\frac{1}{2}\frac{\lambda^{t-t^{\prime}}}{\sigma}\left(y_{t^{\prime}}-\alpha_{t^{\prime}}\right)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{1}{2}\log\left(\frac{2\pi\sigma^{2}}{\lambda^{2(t-t^{\prime})}}\right)-\frac{1}{2}\log\left(\frac{2\pi\sigma^{2}}{\lambda^{2(t-t^{\prime})}}\right)\\\
&=\log\mathcal{N}\left(y_{t^{\prime}}|\alpha_{t^{\prime}},\frac{\sigma^{2}}{\lambda^{t-t^{\prime}}}\right)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{1}{2}\lambda^{t-t^{\prime}}\log\left(2\pi\sigma^{2}\right)+\frac{1}{2}\log\left(\frac{2\pi\sigma^{2}}{\lambda^{2(t-t^{\prime})}}\right)\end{split}$
(26)

Now, we including this term into (25), we get that

$\begin{split}&\log G^{\lambda}({\bf f}_{u},{\bf y})=\log\mathcal{N}\left({\bf
y}|{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf
f}_{u},\sigma^{2}{{\boldsymbol{\Lambda}}}^{-1}\right)\\\
&-\frac{1}{2}\sum_{t^{\prime}=1}^{t}\left(\log\lambda^{2(t-t^{\prime})}+(\lambda^{t-t^{\prime}}-1)\log(2\pi\sigma^{2})\right)\\\
&-\frac{1}{2}\sum_{t^{\prime}=1}^{t}\frac{\lambda^{t-t^{\prime}}}{\sigma^{2}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right).\end{split}$
(27)

where ${\boldsymbol{\Lambda}}$ is a $t\times t$ diagonal matrix with
${\boldsymbol{\Lambda}}_{t^{\prime}t^{\prime}}=\lambda^{t-t^{\prime}}$.

Once $\log G^{\lambda}({\bf f}_{u},{\bf y})$ is defined, we can obtain the
$\lambda$ dependent marginal likelihood or variational bound as

$\begin{split}&F_{V}^{\lambda}({\bf U})=\log\mathcal{N}\left({\bf y}|{\bf
0},\sigma^{2}{\boldsymbol{\Lambda}}^{-1}+{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf
K}_{ux}\right)\\\
&-\frac{1}{2}\sum_{t^{\prime}=1}^{t}\left(\log\lambda^{2(t-t^{\prime})}+(\lambda^{t-t^{\prime}}-1)\log(2\pi\sigma^{2})\right)\\\
&-\frac{1}{2}\sum_{t^{\prime}=1}^{t}\frac{\lambda^{t-t^{\prime}}}{\sigma^{2}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right)\end{split}$ (28)

and simplifying it, by removing the constant terms, we get

$\begin{split}&F_{V}^{\lambda}({\bf U})=\log\mathcal{N}\left({\bf y}|{\bf
0},\sigma^{2}{\boldsymbol{\Lambda}}^{-1}+{\bf K}_{xu}{\bf K}_{uu}^{-1}{\bf
K}_{ux}\right)\\\
&-\frac{1}{2}\sum_{t^{\prime}=1}^{t}(\lambda^{t-t^{\prime}}-1)\log(2\pi\sigma^{2})\\\
&+\frac{1}{2}\sum_{t^{\prime}=1}^{t}\frac{\lambda^{t-t^{\prime}}}{\sigma^{2}}\left(k_{t^{\prime}t^{\prime}}-{\bf
k}_{t^{\prime}u}{\bf K}_{uu}^{-1}{\bf k}_{ut^{\prime}}\right)\end{split}$ (29)

##  Appendix B An Adaptive SGP with Stochastic Variational Inference

An alternative optimization of (8) can be carried out by means of Stochastic
Variational Inference (SVI) [25]. For this purpose, we firstly rewrite the
true marginal bound (8) as an Evidence Lower Bound (ELBO)

$F_{V}({\bf U},q({\bf f}_{u}))=\mathbb{E}_{q({\bf f}_{u})}\left[\log{G({\bf
f}_{u},{\bf y})p({\bf f}_{u})}-\log{q({\bf f}_{u})}\right]$ (30)

Secondly, we take into account that $G({\bf f}_{u},{\bf y})$ factorizes with
respect to the data,

$\begin{split}\log G({\bf f}_{u},{\bf y})=&\sum_{n=1}^{N}\log G({\bf
f}_{u},y_{n})\\\ =&\sum_{n=1}^{N}\log\mathcal{N}\left(y_{n}|{\bf k}_{nu}{\bf
K}_{uu}^{-1}{\bf f}_{u},\sigma^{2}\right)\\\
&-\sum_{n=1}^{N}\frac{1}{2\sigma^{2}}\left(k_{nn}-{\bf k}_{nu}{\bf
K}_{uu}^{-1}{\bf k}_{un}\right)\end{split}$ (31)

where ${\bf k}_{nu}$ is a vector with the kernel products of ${\bf x}_{n}$
with the inducing points and $k_{nn}$ is the kernel of ${\bf x}_{n}$ with
itself.

Finally, a stochastic expression of the marginal bound can be obtained as

$F_{SVI}({\bf U},q({\bf f}_{u}))=\mathbb{E}_{q({\bf
f}_{u})}\left[\sum_{n=1}^{N}\log{G({\bf f}_{u},y_{n})p({\bf
f}_{u})}-\log{q({\bf f}_{u})}\right].$ (32)

An interesting point of expression (32) relies on the fact that the bound over
the marginal likelihood is now expressed as a sum over the training data,
which facilitates endowing this formulation with adaptive capabilities by
weighting each data term with forgetting factor $0<\lambda<1$. This way,
taking into account that we are in a time-dependent scenario at time $t$ and
including $\lambda$ into (32) new bound becomes

$\begin{split}&F_{SVI}^{\lambda}({\bf U},q({\bf f}_{u}))\\\
&=\mathbb{E}_{q({\bf
f}_{u})}\left[\sum_{t^{\prime}=1}^{t}\lambda^{t-t^{\prime}}\log{G({\bf
f}_{u},y_{t^{\prime}})p({\bf f}_{u})}-\log{q({\bf f}_{u})}\right].\end{split}$
(33)

By maximizing this lower bound, the inducing points ${\bf U}$ can be updated
at every instant, as well as the variational parameters ${\boldsymbol{\mu}}$
and ${\bf A}$, the kernel parameters, and the noise variance. Besides, this
can be easily carried out with probabilistic programming packages, such as
Pyro [24], which samples from the variational distribution $q({\bf f}_{u})$ to
approximate this bound with a computational cost of $\mathcal{O}(NM^{2})$. If
we want to reduce this computational burden or the required memory, we can
work with a data window of length $T$ covering training data from instant
$t-T+1$ to $t$. Thus, the sum over the data in Eq. (33) can be replaced by a
sum from $t-T+1$ to $t$ to obtain the values of ${\bf U}$ ${\boldsymbol{\mu}}$
and ${\bf A}$ at instant $t$. During the next instant, $t+1$, the sliding
window advances one position and the parameters are updated in a stochastic
fashion using the previous solution as the starting point. This way, the
computational burden is $\mathcal{O}(TM^{2})$.

So, including the forgetting factor $\lambda$, this formulation provides an
adaptive SGP formulation, as experimental results show (see Section IV).
However, the final computational burden of this model is going to be similar
to retrain an SGP in each iteration since we are going to need many iterations
of the stochastic optimization to make the algorithm to converge to an
adequate solution. Altough this number of iterations can be reduced by using
the optimum values of the parameters at instant $t-1$ to compute its new value
at time $t$, it is not low enough since the stochastic nature of the model
forces us to sample in several iterations over the variational distribution to
obtain a workable solution.

## References

  * [1] L. Tan and J. Jiang, _Digital signal processing: fundamentals and applications_. Academic Press, 2018.
  * [2] S. Haykin, _Adaptive filter theory_ , 4th ed. Upper Saddle River, NJ: Prentice Hall, 2002.
  * [3] V. Ingle, S. Kogon, and D. Manolakis, _Statisical and adaptive signal processing_. Artech, 2005.
  * [4] K. P. Murphy, _Machine learning: a probabilistic perspective_. MIT press, 2012.
  * [5] C. K. Williams and C. E. Rasmussen, _Gaussian processes for machine learning_. MIT press Cambridge, MA, 2006\.
  * [6] Z. Ghahramani and H. Attias, “Online variational bayesian learning.” NIPS Workshop on Online Learning, 2000. [Online]. Available: http://mlg.eng.cam.ac.uk/zoubin/papers/nips00w.pdf
  * [7] S. Van Vaerenbergh, M. Lázaro-Gredilla, and I. Santamaria, “Kernel recursive least-squares tracker for time-varying regression,” _Neural Networks and Learning Systems, IEEE Transactions on_ , vol. 23, pp. 1313–1326, 08 2012\.
  * [8] H. Liu, Y.-S. Ong, X. Shen, and J. Cai, “When gaussian process meets big data: A review of scalable GPs,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 31, no. 11, pp. 4405–4423, 2020.
  * [9] J. Quinonero-Candela and C. E. Rasmussen, “A unifying view of sparse approximate gaussian process regression,” _The Journal of Machine Learning Research_ , vol. 6, pp. 1939–1959, 2005.
  * [10] M. Titsias, “Variational learning of inducing variables in sparse gaussian processes,” in _Artificial intelligence and statistics_. PMLR, 2009, pp. 567–574.
  * [11] L. Csató and M. Opper, “Sparse on-line gaussian processes,” _Neural computation_ , vol. 14, no. 3, pp. 641–668, 2002.
  * [12] T. Evans and P. Nair, “Scalable gaussian processes with grid-structured eigenfunctions (GP-GRIEF),” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1417–1426.
  * [13] S. Stanton, W. Maddox, I. Delbridge, and A. G. Wilson, “Kernel interpolation for scalable online gaussian processes,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2021, pp. 3133–3141.
  * [14] T. N. Hoang, Q. M. Hoang, and B. K. H. Low, “A unifying framework of anytime sparse gaussian process regression models with stochastic variational inference for big data,” in _International Conference on Machine Learning_. PMLR, 2015, pp. 569–578.
  * [15] C.-A. Cheng and B. Boots, “Incremental variational sparse gaussian process regression,” _Advances in Neural Information Processing Systems_ , vol. 29, 2016.
  * [16] T. D. Bui, C. Nguyen, and R. E. Turner, “Streaming sparse gaussian process approximations,” _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * [17] J. Hensman, N. Fusi, and N. D. Lawrence, “Gaussian processes for big data,” _arXiv preprint arXiv:1309.6835_ , 2013.
  * [18] M. A. Aizerman, E. M. Braverman, and L. Rozoner, “Theoretical foundations of the potential function method in pattern recognition learning,” _Automation and remote Control_ , vol. 25, pp. 821–837, 1964.
  * [19] C. E. Rasmussen, “Gaussian processes in machine learning,” in _Summer school on machine learning_. Springer, 2003, pp. 63–71.
  * [20] M. W. Seeger, C. K. Williams, and N. D. Lawrence, “Fast forward selection to speed up sparse gaussian process regression,” in _International Workshop on Artificial Intelligence and Statistics_. PMLR, 2003, pp. 254–261.
  * [21] E. Snelson and Z. Ghahramani, “Sparse gaussian processes using pseudo-inputs,” _Advances in neural information processing systems_ , vol. 18, 2005.
  * [22] A. J. Smola and B. Schölkopf, “Sparse greedy matrix approximation for machine learning.” Morgan Kaufmann, 2000, pp. 911–918.
  * [23] N. Lawrence, M. Seeger, and R. Herbrich, “Fast sparse gaussian process methods: The informative vector machine,” _Advances in neural information processing systems_ , vol. 15, 2002.
  * [24] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. Szerlip, P. Horsfall, and N. D. Goodman, “Pyro: Deep universal probabilistic programming,” _The Journal of Machine Learning Research_ , vol. 20, no. 1, pp. 973–978, 2019.
  * [25] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley, “Stochastic variational inference,” _Journal of Machine Learning Research_ , 2013.