
\subsection{Gaussian process regression}

Consider a set of $N$ input-output pairs ${\cal D} = \{\vect{x}_i, y_i\}_{i=1}^N$, where $\vect{x}_i \in \Real^D$ are $D$-dimensional input vectors and $y_i\in\Real$ are scalar outputs. We assume that the observed data can be described by the following model,
\begin{equation}
y_i = f(\mathbf{x}_i) + \varepsilon_i,
\label{eq:regmodel}
\end{equation}
in which $f$ represents an unobservable \emph{latent function} and $\varepsilon_i \sim \mathcal{N}(0,\sigma_n^2)$ is zero-mean Gaussian noise.

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution \cite{rasmussen2006gaussian}. 
To indicate that a random function $f(\mathbf{x})$ follows a Gaussian process we write it as
\begin{equation*}
f(\mathbf{x}) \sim \GP(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}')).
\end{equation*}
All values of $f$ at any locations $\mathbf{x}$ are jointly normally distributed, with $m(\mathbf{x})$ and $k(\mathbf{x},\mathbf{x}')$ representing the mean function and covariance function, respectively.

In a Bayesian regression setting, we are interested in inferring the predictive distribution of a new, unseen output $y_{\ast}$ given the corresponding input $\mathbf{x}_{\ast}$ and the data ${\cal D}$.
In particular, we take a Gaussian process as the prior over the latent function, and the vector of observations $[y_1,\dots,y_n]^\top$ is related to the latent function through the likelihood function $p(\vect{y}|f)$. When the observations are contaminated with zero-mean Gaussian noise, as in Eq.~\eqref{eq:regmodel},
there exists a closed-form solution for the \emph{posterior distribution} over functions, i.e. the distribution over the unknown function $f(x)$ after incorporating all the observed data. Specifically, the posterior  of the function at any new location $\mathbf{x}_{\ast}$ is described by
\begin{equation*}
p(f_{\ast}|\mathbf{x}_{\ast},{\cal D}) = \Normal(\hat f_\ast,\hat \sigma_\ast^2).
\end{equation*}
When $m(\mathbf{x})=0$, which is a very common assumption, we obtain the following expressions \cite{rasmussen2006gaussian}
\begin{subequations}
\begin{align}
\hat f_\ast &= \vect{k}_\ast^\top (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \vect{y} \label{eq:pred_mean_gpml}\\
\hat \sigma_\ast^2 & = k_{\ast\ast} - \vect{k}_\ast^\top (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \vect{k}_\ast^\top \label{eq:pred_var_gpml}
\end{align}
\label{eq:pred_distr_gpml}%
\end{subequations}
where the covariances (or \emph{kernel}) matrices $\mathbf{K}$ contain the elements $[\mathbf{K}]_{ij} = k(\mathbf{x}_i,\mathbf{x}_j)$ and we have introduced the shorthand notations $\vect{k}_\ast = [k(\mathbf{x}_1,\mathbf{x}_\ast),\dots,k(\mathbf{x}_N,\mathbf{x}_\ast)]^\top$ and $k_{\ast\ast} = k(\mathbf{x}_\ast,\mathbf{x}_\ast)$.
The matrix inversion involved in Eqs.~\eqref{eq:pred_mean_gpml} and \eqref{eq:pred_var_gpml} leads to $\mathcal{O}(N^3)$ complexity.

\subsection{Incremental GP updates}

In an \emph{online} scenario, the data pairs are made available on a one-at-a-time basis, i.e. $(\mathbf{x}_t,y_t)$ arrives at time $t$. Instead of recalculating the predictive distribution entirely once a new data pair $(\mathbf{x}_{t+1},y_{t+1})$ arrives, i.e. by solving \eqref{eq:pred_distr_gpml}, it is more interesting to perform an incremental update.
In this section we briefly review the sequential updates for online GP regression as presented in \cite{vanvaerenbergh2012kernel}. In order to avoid the unbounded growth of the involved matrices, the online learning process is typically coupled with a \emph{sparsification} procedure.

At the $t$-th iteration of the online GP, the model contains the variables
\begin{equation*}
\mathcal{M}_t = \{{\cal D}_t,\vect{\mu}_t,\mat{\Sigma}_t,\mathbf{Q}_t\},
\end{equation*}
where ${\cal D}_t$ is the observed data set that contains the data pairs $\{\mathbf{x}_i,y_i\}_{i=1}^t$; $\vect{\mu}_t$ and $\mat{\Sigma}_t$ are the mean and covariance matrices of the posterior $p(\mathbf{f}_t|{\cal D}_t) = \Normal(\vect{\mu}_t,\mat{\Sigma}_t)$; and $\mathbf{Q}_t = \mathbf{K}_t^{-1}$ is the inverse covariance matrix corresponding to ${\cal D}_t$.

When a new data pair $(\mathbf{x}_{t+1},y_{t+1})$ is obtained, the posterior distribution
\begin{equation*}
p(\mathbf{f}_{t+1}|{\cal D}_{t+1}) = \Normal(\mathbf{f}_{t+1}| \vect{\mu}_{t+1}, \mat{\Sigma}_{t+1})
\end{equation*}
is updated as
\begin{subequations}
\begin{align}
\vect{\mu}_{t+1} & =
\begin{bmatrix} \vect{\mu}_t \\ \hat{y}_{t+1} \end{bmatrix}
+ \frac{y_{t+1} - \hat{y}_{t+1}}{\hat{\sigma}_{yt+1}^2}
\begin{bmatrix} \mathbf{h}_{t+1} \\ \hat{\sigma}_{ft+1}^2, \end{bmatrix}, \label{eq:pmean} \\
\mat{\Sigma}_{t+1} & = 
\begin{bmatrix} \mat{\Sigma}_t & \mathbf{h}_{t+1}\\ \mathbf{h}_{t+1}^\top & \hat{\sigma}_{ft+1}^2 \end{bmatrix}
- \frac{1}{\hat{\sigma}_{yt+1}^2}
\begin{bmatrix} \mathbf{h}_{t+1} \\ \hat{\sigma}_{ft+1}^2 \end{bmatrix}
\begin{bmatrix} \mathbf{h}_{t+1} \\ \hat{\sigma}_{ft+1}^2 \end{bmatrix}^\top, \label{eq:pcov}
\end{align}
\end{subequations}
where
\begin{align*}
\mathbf{q}_{t+1} & = \mathbf{Q}_t \vect{k}_{t+1}, \\
\mathbf{h}_{t+1} & = \mat{\Sigma}_t \mathbf{q}_{t+1},
\end{align*}
and the vector $\vect{k}_{t+1}$ has elements $[\vect{k}_{t+1}]_i = k(\mathbf{x}_i,\mathbf{x}_{t+1})$. Furthermore, the output variance is obtained as
\begin{equation*}
\hat{\sigma}_{yt+1}^2 = \sigma_n^2 + \hat{\sigma}_{ft+1}^2,
\label{eq:sigma_y}
\end{equation*}
and the variance of the latent function evaluations is
\begin{align*}
\hat{\sigma}_{ft+1}^2 & = k_{t+1} + \vect{k}_{t+1}^\top (\mathbf{Q}_t\mat{\Sigma}_t\mathbf{Q}_t-\mathbf{Q}_t)\vect{k}_{t+1}\\
& = \gamma_{t+1}^2 + \mathbf{q}_{t+1}^T \mathbf{h}_{t+1}.
\end{align*}
In order to further simplify equations, the variable
\begin{equation*}
\gamma_{t+1}^2 = k_{t+1} - \vect{k}_{t+1}^T \mathbf{Q}_t \vect{k}_{t+1}
\end{equation*}
is introduced.
The inverse kernel matrix $\mathbf{Q}_t$ can be updated efficiently through
\begin{equation}
\mathbf{Q}_{t+1} = 
\begin{bmatrix}
\mathbf{Q}_t & \mathbf{0} \\ \mathbf{0}^\top & 0
\end{bmatrix}
+ \frac{1}{\gamma_{t+1}^2}
\begin{bmatrix}
\mathbf{q}_{t+1} \\ -1
\end{bmatrix}
\begin{bmatrix}
\mathbf{q}_{t+1} \\ -1
\end{bmatrix}^\top.
\label{eq:q_update}
\end{equation}
After applying Eqs.~\eqref{eq:pmean}, \eqref{eq:pcov}, and \eqref{eq:q_update}, we obtain the updated model $\mathcal{M}_{t+1} = \{{\cal D}_{t+1},\vect{\mu}_{t+1},\mat{\Sigma}_{t+1},\mathbf{Q}_{t+1}\}$.

At each step $t$ of the learning process, the predictive distribution of a new observation $y_{t+1}$ given all past data is a Gaussian
$
p(y_{t+1}|{\cal D}_t) = \Normal(y_{t+1}|\hat y_{t+1},\hat \sigma_{yt+1}^2)
$
with
\begin{subequations}
\begin{align}
\hat y_{t+1} & = \mathbf{q}_{t+1}^\top \vect{\mu}_t = \vect{k}_{t+1}^\top \mathbf{Q}_t \vect{\mu}_t \label{eq:pred_mean_ogp} \\
\hat \sigma_{yt+1}^2 & = \sigma_n^2 + k_{t+1} + 
\vect{k}_{t+1}^\top (\mathbf{Q}_t\mat{\Sigma}_t\mathbf{Q}_t-\mathbf{Q}_t)\vect{k}_{t+1}. \label{eq:pred_cov_ogp}
\end{align}
\label{eq:pred_distr_ogp}
\end{subequations}
For more details we refer the reader to \cite{vanvaerenbergh2012kernel,lazaro2011bayesian}.

\subsection{Kernel recursive least-squares}
The kernel-methods counterpart of online GP regression is kernel recursive least-squares (KRLS, see for instance \cite{engel2004kernel}), which obtains the KRR solution recursively. After receiving $t$ data points, the kernel weights obtained by KRLS are those that solve the batch problem
\begin{equation}
\vect{\alpha}_t = (\mathbf{K}_t + \sigma_n^2\mat{I})^{-1}\vect{y}_t.
\label{eq:krls_alpha}
\end{equation}
The same weights can be obtained in online GP regression by computing
\begin{equation}
\vect{\alpha}_t = \mathbf{K}_t^{-1}\vect{\mu}_t = \mathbf{Q}_t\vect{\mu}_t,
\end{equation}
which follows from Eq.~\eqref{eq:pred_mean_ogp}.

\subsection{KLMS algorithms} 

Similar to online GP regression, updating the KRLS estimate in Eq.~\eqref{eq:krls_alpha} with a new data point requires quadratic complexity. Kernel least-mean-squares (KLMS) algorithms alleviate this computational burden by performing stochastic gradient descent of the mean square error, resulting in linear complexity per time step \cite{vanvaerenbergh2014online}.

KLMS algorithms can be categorized into two classes, depending on which kernel weights they update in each iteration to account for the prediction error\footnote{Both approaches are different approximations of the same update formulation in the \emph{kernel feature space}; see \cite{vanvaerenbergh2014online} for details.}. We outline both approaches briefly in the remainder of this section.

\subsubsection{Type-I KLMS: concentrating the novelty}

The first type of KLMS algorithm updates only one coefficient in order to compensate for the prediction error, at each time step. When the weight vector is allowed to grow, this update takes the form
\begin{equation}
\vect{\alpha}_{t+1}^{(I)} = \begin{bmatrix}
\vect{\alpha}_t\\ \eta e_{t+1}
\end{bmatrix},
\label{eq:klms_update}
\end{equation}
in which $e_{t+1} = y_{t+1} - \hat y_{t+1}$ is the instantaneous error.

Eq.~\eqref{eq:klms_update} represents the basic update of what is known as the KLMS algorithm, proposed in \cite{liu2008kernel}. In order to avoid the infinite growth of $\vect{\alpha}_t$, a more sophisticated version of this algorithm was presented in \cite{chen2012quantized}, known as Quantized Kernel Least Mean Square (QKLMS). When QKLMS receives a datum similar to a previously seen datum, for instance the $i$-th base it has stored, it does not expand $\vect{\alpha}_t$ but instead updates the corresponding weight $\alpha_i$.

\subsubsection{Type-II KLMS: spreading the novelty}

A different strategy consists in updating all coefficients of $\vect{\alpha}_t$ in each iteration. This approach is followed for instance by the Kernel Normalized Least Mean Square (KNLMS) algorithm \cite{richard2009online}, whose update reads
\begin{equation}
\vect{\alpha}_{t+1}^{(II)} = \begin{bmatrix}
\vect{\alpha}_t\\ 0
\end{bmatrix} + \eta \frac{e_{t+1}}{\epsilon + k_{t+1}^2 + \|\vect{k}_{t+1}\|^2} \begin{bmatrix}
\vect{k}_{t+1}\\ k_{t+1}
\end{bmatrix}
\label{eq:knlms_update}
\end{equation}
when $\vect{\alpha}_t$ is allowed to grow. Note that the rule \eqref{eq:knlms_update} updates all coefficients in each iteration. In order to avoid unbounded growth, KNLMS follows a \emph{coherence} criterion that promotes sparsity.

\subsection{Specific cases: $\beta=0$ and $\beta=1$}

Setting $\beta=0$ in Eq.~\eqref{eq:update_alpha_subs} yields the update rule
\begin{equation}
\vect{\alpha}_{t+1}^{\beta=0} = \begin{bmatrix}
\vect{\alpha}_t\\ \frac{1}{\sigma_n^2 + k_{t+1}} e_{t+1}
\end{bmatrix}.
\label{eq:update_alpha_subs1}
\end{equation}
This rule is very similar to the KLMS update \eqref{eq:klms_update}, and even identical when the learning rate is set to $\eta = 1/(\sigma_n^2 + k_{t+1})$. Note, furthermore, that $\beta=0$ indicates that the posterior covariance from Eq.~\eqref{eq:pcov_simple} reduces to
\begin{equation}
\mat{\Sigma}_t^{\beta=0} = \mathbf{K}_t.
\end{equation}
In other words, for $\beta=0$ we obtain a KLMS model that implies a fixed posterior covariance, equal to the prior covariance.
The predictive variance, Eq.~\eqref{eq:beta_sigma}, simplifies to
\begin{equation}
\hat\sigma^2_{ft+1}|^{\beta=0} = k_{t+1}.
\label{eq:pred_cov_beta0}
\end{equation}
Under the adopted framework, this is the model that underlies algorithms such as KLMS from \cite{liu2008kernel} and QKLMS from \cite{chen2012quantized}.

By setting $\beta=1$ in Eq.~\eqref{eq:update_alpha_subs}, the update rule becomes
\begin{equation}
\vect{\alpha}_{t+1}^{\beta=1} = \begin{bmatrix}
\vect{\alpha}_t\\ 0
\end{bmatrix} + \frac{e_{t+1}}{\sigma_n^2 + k_{t+1} + \|\vect{k}_{t+1}\|^2} \begin{bmatrix}
\vect{k}_{t+1}\\ 1
\end{bmatrix}
\label{eq:update_alpha_subs2}
\end{equation}
which is very similar to the update of KNLMS, shown in Eq.~\eqref{eq:knlms_update}. The implied posterior covariance then reads
\begin{equation}
\mat{\Sigma}_t^{\beta=1} = \mathbf{K}_t \mathbf{K}_t + \mathbf{K}_t,
\end{equation}
and the predictive covariance from Eq.~\eqref{eq:beta_sigma} now becomes
\begin{equation}
\hat\sigma^2_{ft+1}|^{\beta=1} = k_{t+1} + \vect{k}_{t+1}^\top \vect{k}_{t+1}.
\label{eq:pred_cov_beta1}
\end{equation}
Interestingly, this KLMS model implies a predictive covariance, for each point in space, that is larger than the prior covariance $k_{t+1}$. Furthermore, a closer inspection of the second term in Eq.~\eqref{eq:pred_cov_beta1} shows that the predictive variance grows as more training data is processed.
According to the GP framework, this is the model that underlies type-II KLMS algorithms, in particular KNLMS from \cite{richard2009online}, for which $\beta = 1$, and in general any algorithm that corresponds to $\beta>0$.

As seen through the adopted probabilistic perspective, the capability of Type-II KLMS algorithms to update all coefficients with each new data point is related to an increase in prediction uncertainty.
These effects are illustrated in Fig. \ref{fig:uncertainty}, which compares the interpreted predictive uncertainty for the different algorithms. GP regression exhibits the expected behavior, i.e. uncertainty shrinks around the observed data points. Type-I KLMS ($\beta=0$) assumes a fixed predictive covariance, regardless of the number of observed data and the distances between them. The behavior of type-II KLMS (for instance with $\beta=1$) is rather counterintuitive: its update procedure increases the predictive variance as more data points are observed, and this happens in the neighborhoods of those points. While an in-depth study is needed to extract solid conclusions from this observation, it is interesting to note that a similar behavior has been observed in the linear recursive least squares (RLS) algorithm with forgetting factor, see \cite{vanvaerenbergh2012kernel}.

\begin{figure}[t!]
\centering
\includegraphics[width=\linewidth]{uncertainty_plot_3}
\includegraphics[width=\linewidth]{uncertainty_plot_8}
\includegraphics[width=\linewidth]{uncertainty_plot_25}
\caption{Comparison of the predictive variance of three algorithms, for $3$ data points (top plot, data is marked as blue crosses), $8$ data points (middle), and $25$ data points (bottom). The predictive mean of GP regression is indicated as the black curve that passes through the observations. The grey zone marks the GP mean plus/minus two standard deviations $\hat \sigma_{y}$, corresponding to the GP's $95\%$ confidence interval. The dashed line marks the confidence interval for type-I KLMS ($\beta=0$), and the dash-dot line marks the confidence interval for type-II KLMS with $\beta=1$.
}
\label{fig:uncertainty}
\end{figure}




\subsection{Online regression on stationary data}

In the first experiment, we wish to study the effect of different values of $\beta$ in the KLMS algorithm from Eq.~\eqref{eq:update_alpha_subs}, which we will denote by $\beta$-KLMS. We evaluate this algorithm and three established kernel adaptive filtering algorithms on the stationary KIN40K benchmark.\footnote{Available at \url{http://www.cs.toronto.edu/~delve/data/datasets.html}} This data set is obtained from the forward kinematics of an $8$-link all-revolute robot arm, and it represents a very difficult regression problem. We randomly select $5000$ data points for online training, and $5000$ points for testing the regression.

The algorithms are considered in their evergrowing version here, in order to highlight only the influence of $\beta$. The KRLS-T algorithm, however, which implements the full online GP regression, is given a limited memory of $500$ bases, for computational reasons. A Gaussian kernel was used for all algorithms, with parameters determined offline by standard GP regression, as detailed in \cite{vanvaerenbergh2012kernel}.

The results are shown in Fig.~\ref{fig:kin40k}. Each point of the learning curves corresponds to the test error on the entire test set. We observe that for $\beta=1$, the $\beta$-KLMS algorithm obtains similar performance to KNLMS from \cite{richard2009online}. For $\beta=0$, the performance is almost identical to that of KLMS from \cite{liu2008kernel}.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{kin40k_nmse}
\caption{Performance comparison for online prediction on the KIN40K benchmark regression problem.} \label{fig:kin40k}
\end{figure}

\subsection{Online prediction of time series with a model switch}

In Fig. \ref{fig:reconvergence} we repeat a standard experiment from the kernel adaptive filtering literature. In particular, we consider a nonlinear system comprised of a linear channel followed by a nonlinearity, and the online regression task consists in predicting the next sample of the time series produced at its output. An abrupt change in the linear channel is triggered after $500$ time steps, in order to test the algorithms' ability to reconverge. The experiment is repeated $5$ times with random channel coefficients, and the average results are shown in Fig. \ref{fig:reconvergence}. The same kernel was used for each algorithm. 

In this experiment we observe that, as is usual in the LMS and KLMS literature, the $\beta$-KLMS model exhibits tracking behavior although it is based on a static data model. Furthermore, it does so without the need of including an additional parameter to control the update step size.



\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{reconvergence}
\caption{Performance comparison for online prediction on a nonlinear channel with an abrupt change.} \label{fig:reconvergence}
\end{figure}



\section{Introduction}
\label{sec:intro}
\input{1_introduction.tex}


\section{Online GP regression}
\label{sec:ogp}
\input{2_ogp.tex}


\section{Kernel adaptive filtering and KLMS}
\label{sec:klms}
\input{3_klms.tex}


\section{From online GP to KLMS}
\label{sec:approx}
\input{4_approximation.tex}




\section{Experiments}
\label{sec:exp}
\input{6_experiments.tex}


\section{Conclusions}
\label{sec:concl}
\input{7_conclusions.tex}


\bibliographystyle{IEEEbib}
