Year: 2020

Abstract:
We introduce a framework for Continual Learning (CL) based on Bayesian inference over the function space rather than the parameters of a deep neural network. This method, referred to as functional regularisation for Continual Learning, avoids forgetting a previous task by constructing and memorising an approximate posterior belief over the underlying task-specific function. To achieve this we rely on a Gaussian process obtained by treating the weights of the last layer of a neural network as random and Gaussian distributed. Then, the training algorithm sequentially encounters tasks and constructs posterior beliefs over the task-specific functions by using inducing point sparse Gaussian process methods. At each step a new task is first learnt and then a summary is constructed consisting of (i) inducing inputs – a fixed-size subset of the task inputs selected such that it optimally represents the task – and (ii) a posterior distribution over the function values at these inputs. This summary then regularises learning of future tasks, through Kullback-Leibler regularisation terms. Our method thus unites approaches focused on (pseudo-)rehearsal with those derived from a sequential Bayesian inference perspective in a principled way, leading to strong results on accepted benchmarks.

heading: 1 INTRODUCTION
text: Recent years have seen a resurgence of interest in continual learning, which refers to systems that learn in an online fashion from data associated with possibly an ever-increasing number of tasks (Ring, 1994; Robins, 1995; Schmidhuber, 2013; Goodfellow et al., 2013). A continual learning system must adapt to perform well on all earlier tasks without requiring extensive re-training on previous data. There are two main challenges for continual learning (i) avoiding catastrophic forgetting, i.e. remembering how to solve earlier tasks, and (ii) scalability over the number of tasks. Other possible desiderata may include forward and backward transfer, i.e. learning new tasks faster and retrospectively improving on previously tasks.
Similarly to many recent works on continual learning (Kirkpatrick et al., 2017; Nguyen et al., 2017; Rusu et al., 2016; Li & Hoiem, 2017; Farquhar & Gal, 2018), we focus on the scenario where a sequence of supervised learning tasks are presented to a continual learning system based on a deep neural network. While most methods assume known task boundaries, our approach will be also extended to deal with unknown task boundaries. Among the different techniques proposed to address this problem, we have methods which constrain or regularise the parameters of the network to not deviate significantly from those learnt on previous tasks. This includes methods that frame continual learning as sequential approximate Bayesian inference, including EWC (Kirkpatrick et al., 2017) and VCL (Nguyen et al., 2017). Such approaches suffer from brittleness due to representation drift. That is, as parameters adapt to new tasks the values that other parameters are constrained/regularised towards become obsolete (see Section 2.5 for further discussion on this). On the other hand, we have rehearsal/replay buffer methods, which use a memory store of past observations to remember previous ∗Equal contribution
tasks (Robins, 1995; Robins & McCallum, 1998; Lopez-Paz et al., 2017; Rebuffi et al., 2017). While these methods tend to not suffer from brittleness, uncertainty about the unknown functions is not expressed. Furthermore, they rely on various heuristics to decide which data to store (Rolnick et al., 2018), often requiring large quantities of stored observations to achieving good performance. In this paper we will address the open problem of deriving an optimisation objective to select the best observations for storage.
In this paper, we develop a new approach to continual learning which addresses the shortcomings of both categories. It is based on approximate Bayesian inference, but on the space of functions instead of neural network parameters, so does not suffer from the aforementioned brittleness. Intuitively, while previous approaches constrain the parameters of a neural network to limit deviations from previously learnt parameters, our approach instead constrains the neural network predictions from deviating too far from those that solve previous tasks.
Effectively, our approach avoids forgetting an earlier task by memorising an approximate posterior belief over the underlying task-specific function. To implement this, we consider Gaussian processes (GPs) (Rasmussen & Williams, 2005), and make use of inducing point sparse GP methods (Csato & Opper, 2002; Titsias, 2009; Hensman et al., 2013; Bui et al., 2017b), which summarise posterior distributions over functions using small numbers of so-called inducing points. These inducing points are selected from the training set by optimising a variational objective, providing a principled mechanism to compress the dataset to a meaningful subset of fixed size. They are kept around when moving to the next task and, together with their posterior distributions, are used to regularise the continual learning of future tasks, through Kullback-Leibler regularisation terms within a variational inference framework, thus avoiding catastrophic forgetting of earlier tasks. Our approach bears similarities to replay-based approaches, with inducing points playing the role of the rehersal/replay buffer, but has two important advantages. First the approximate posterior distributions at the inducing points captures the uncertainty of the unknown function as well, rather than providing merely target values. Second, inducing points can be optimised using specialised criteria from the GP literature, achieving better performance than a random selection of observations. An intuitive depiction of our approach is given in Figure 1.
To enable our functional regularisation approach to deal with high-dimensional and complex datasets, we use a linear kernel with features parameterised by neural networks (Wilson et al., 2016). Such GPs can be understood as Bayesian neural networks, where only the weights of the last layer are treated in a Bayesian fashion, while those in earlier layers are optimised. This view allows for a more computationally efficient and accurate training procedure to be carried out in weight space, before the approximation is translated into function space where the inducing points are constructed and then used for regularising learning of future tasks. Finally, note that inducing points are also used to regularize the deep network, even though they were selected to best represent functions given by the GP. See Section 2.3 for further details.

heading: 2 FUNCTIONAL REGULARISATION FOR CONTINUAL LEARNING
text: We consider supervised learning of multiple tasks, with known task boundaries, that are processed sequentially one at a time. At each step we receive a set of examples (Xi,yi) where Xi = {xi,j}Nij=1
are input vectors and yi = {yi,j}Nij=1 are output targets so that each yi,j is assigned to the input xi,j ∈ RD. We assume the most extreme case (and challenging in terms of avoiding forgetting) where each dataset (Xi,yi) introduces a new task, while less extreme cases can be treated similarly. We wish to sequentially train a shared model or representation from all tasks so that catastrophic forgetting is avoided, i.e. when the model is trained on the i-th task it should still provide accurate predictions for all tasks j < i seen in the past. As a model we consider a deep neural network with its final hidden layer providing the feature vector φ(x; θ) ∈ RK where x is the input vector and θ are the model parameters. This representation is shared across tasks and θ is a task-shared parameter. To solve a specific task i we additionally construct an output layer
fi(x;wi) ≡ fi(x;wi, θ) = w>i φ(x; θ), (1) where for simplicity we assume that fi(x;wi) is a scalar function and wi is the vector of task-specific weights. Dealing with vector-valued functions is straightforward and is discussed in the Appendix. By placing a Gaussian prior on the output weights, wi ∼ N (wi|0, σ2wI), we obtain a distribution over functions. While each task has its own independent/private weight vector wi the whole distribution refers to the full infinite set of tasks that can be tackled by the same feature vector φ(x; θ). We can marginalise out wi and obtain the equivalent function space view of the model, where each task-specific function is an independent draw from a GP (Rasmussen & Williams, 2005), i.e.
fi(x) ∼ GP(0, k(x, x′)), k(x, x′) = σ2wφ(x; θ)>φ(x′; θ), where the kernel function is defined by the dot product of the neural network feature vector. By assuming for now that all possible tasks are simultaneously present similarly to multi-task GPs (Bonilla et al., 2008; Álvarez et al., 2012), the joint distribution over function values and output data for all tasks is written as ∏ i p(yi|fi)p(fi) = ∏ i p(yi|fi)N (fi|0,KXi),where the vector fi stores all function values for the input dataset Xi, i.e. fi,j = f(xi,j), j = 1, . . . , Ni. Also the kernel matrix KXi is obtained by evaluating the kernel function on Xi, i.e. each element [KXi ]j,k = σ 2 wφ(xi,j ; θ)
>φ(xi,k; θ) where xi,j , xi,k ∈ Xi. The form of each likelihood function p(yi|fi) depends on the task, for example if the i-th task involves binary classification then p(yi|fi) = ∏Ni j=1 p(yi,j |fi,j) = ∏Ni j=1 1 1+e−yi,jfi,j where yi,j ∈ {−1, 1} indicates the binary class label. Inference in this model requires estimating each posterior distribution p(fi|yi, Xi), which can be approximated by a multivariate Gaussian N (fi|µi,Σi). Given this Gaussian we can express our posterior belief over any function value fi,∗ at some test input xi,∗ using the posterior GP (Rasmussen & Williams, 2005),
p(fi,∗|Xi, yi) = ∫ pθ(fi,∗|fi)N (fi|µi,Σi)dfi.
Given that the tasks arrive one at a time, the above suggests that one way to avoid forgetting the i-th task is to memorise the corresponding posterior belief N (fi|µi,Σi). While this can regularise continual learning of subsequent tasks (similarly to the more general variational framework in the next section), it can be prohibitively expensive since the non-parametric nature of the model means that for each N (fi|µi,Σi) we need to store O(N2i ) parameters and additionally we need to keep in memory the full set of input vectors Xi.
Therefore, in order to reduce the time and memory requirements we would like to apply data distillation and approximate each full posterior by applying sparse GP methods. As shown next, by applying variational sparse GP inference (Titsias, 2009) in a sequential fashion we obtain a new algorithm for function space regularisation in continual learning.

heading: 2.1 LEARNING THE FIRST TASK
text: Suppose we encounter the first task with data (X1,y1). We introduce a small set Z1 = {z1,j}M1j=1 of inducing inputs where each z1,j lives in the same space as each training input x1,j . The inducing set Z1 can be a subset of X1 or it can contain pseudo inputs (Snelson & Ghahramani, 2006), i.e. points lying between the training inputs. For simplicity next we consider Z1 as pseudo points, although in practice for continual learning it can be more suitable to select them from the training inputs (see Section 2.4). By evaluating the function output at each z1,j we obtain a vector of auxiliary function values u1 = {u1,j}M1j=1, where each u1,j = f(z1,j). Hence, we obtain the joint distribution
p(y1, f1,u1) = p(y1|f1)pθ(f1|u1)pθ(u1). (2)
The exact posterior distribution pθ(f1|u1,y1)pθ(u1|y1) is approximated by a distribution of the form, q(f1,u1) = pθ(f1|u1)q(u1), where q(ui) is a variational distribution and pθ(f1|u1) is the GP prior conditional, pθ(f1|u1) = N (f1|KX1Z1K−1Z1 u1,KX1 − KX1Z1K −1 Z1 KZ1X1). Here, KX1Z1 is the cross kernel matrix between the sets X1 and Z1, KZ1X1 = K > X1Z1
and KZ1 is the kernel matrix on Z1. The method learns (q(u1), Z1) by minimising the KL divergence KL(pθ(f1|u1)q(u1)||pθ(f1|u1,y1)pθ(u1|y1)). The ELBO is also maximised over the neural network feature vector parameters θ that determine the kernel matrices. This ELBO is generally written in the form (Hensman et al., 2013; Lloyd et al., 2015; Dezfouli & Bonilla, 2015; Hensman et al., 2015; Sheth et al., 2015),
F(θ, q(u1)) = N1∑ j=1 Eq(f1,j)[log p(y1,j |f1,j)]− KL(q(u1)||pθ(u1)), (3)
where q(f1,j) = ∫ p(f1,j |u1)q(u1)du1 is an univariate Gaussian distribution with analytic mean and variance that depend on (θ, Z1, q(u1), x1,j). Each expectation Eq(f1,j)[log p(y1,j |f1,j)] is a one-dimensional integral and can be estimated by Gaussian quadrature. The variational distribution q(u1) is chosen to be a Gaussian, parameterised as q(u1) = N (u1|µu1 , Lu1L>u1), where Lu1 is a square root matrix such as a lower triangular Cholesky factor. Then, based on the above we can jointly apply stochastic variational inference (Hensman et al., 2013) to maximise the ELBO over (θ, µu1 , Lu1) and optionally over the inducing inputs Z1.

heading: 2.2 LEARNING THE SECOND AND SUBSEQUENT TASKS
text: The functional regularisation framework for continual learning arises from the variational sparse GP inference method as we encounter the second and subsequent tasks.
Once we have learned the first task we throw away the dataset (X1,y1) and we keep in memory only a task summary consisting of the inducing inputs Z1 and the variational Gaussian distribution q(u1) (i.e. its parameters µu1 and Lu1 ). Note also that θ (that determines the neural network feature vector φ(x; θ)) has a current value obtained by learning the first task. When the dataset (X2,y2) for the second task arrives, a suitable ELBO to continue learning θ and also estimate the second task summary (Z2, q(u2)) is
N1∑ j=1 Eq(f1,j)[log p(y1j |f1,j)] + N2∑ j=1 Eq(f2,j)[log p(y2,j |f2,j)]− ∑ i=1,2 KL(q(ui)||pθ(ui)),
which is just the sum of the corresponding ELBOs for the two tasks. We need to approximate this ideal objective by making use of the fixed summary (Z1, q(u1)) that we have kept in memory for the first task. By considering Z1 ⊂ X1 as our replay buffer with outputs ỹ1 ⊂ y and ũ1 ⊂ f1 the above can be approximated by
N1 M1 M1∑ j=1 Eq(u1,j)[log p(ỹ1,j |u1,j)] + N2∑ j=1 Eq(f2,j)[log p(y2,j |f2,j) − ∑ i=1,2 KL(q(ui)|pθ(ui)),
where each q(u1,j) is a univariate marginal of q(u1). However, since q(u1) is kept fixed the whole expected log-likelihood term N1M1 ∑M1 j=1Eq(u1,j)[log p(ỹ1,j |u1,j)] is just a constant that does not depend on the parameters θ any more. Thus, the objective function when learning the second task reduces to maximising,
F(θ, q(u2)) = N2∑ j=1 Eq(f2,j)[log p(y2,j |f2,j)]− ∑ i=1,2 KL(q(ui)||pθ(ui)).
The only term associated with the first task is KL(q(u1)||pθ(u1)). While q(u1) is fixed (i.e. its parameters are constant), the GP prior pθ(u1) = N (u1|0,KZ1) is still a function of the feature vector parameters θ, since KZ1 depends on θ. Thus, this KL term regularises the parameters θ so that, while learning the second task, the feature vector still needs to explain the posterior distribution over the function values u1 at input locations Z1. Notice that −KL(q(ui)||pθ(ui) is further simplified as∫ q(u1) log pθ(u1)du1 + const, which shows that the regularisation is such that pθ(u1) needs to be consistent with all infinite draws from q(u1) in a moment-matching or maximum likelihood sense.
Similarly for the subsequent tasks we can conclude that for any new task k the objective will be
F(θ, q(uk)) = Nk∑ j=1
Eq(fk,j) log p(yk,j |fk,j)− KL(q(uk)||pθ(uk))︸ ︷︷ ︸ objective for the current task
− k−1∑ i=1
KL(q(ui)||pθ(ui))︸ ︷︷ ︸ regularisation from previous tasks . (4)
Thus, functional regularisation when learning a new task is achieved through the sum of the KL divergences ∑k−1 i=1 KL(q(ui)||pθ(ui)) of all previous tasks, where each q(ui) is the fixed posterior distribution which encodes our previously obtained knowledge about task i < k. Furthermore, in order to keep the optimisation scalable over tasks, we can form unbiased approximations of this latter sum by sub-sampling the KL terms, i.e. by performing minibatch-based stochastic approximation over the regularisation terms associated with these previous tasks.

heading: 2.3 ACCURATE WEIGHT SPACE INFERENCE FOR THE CURRENT TASK
text: While the above framework arises by applying sparse GP inference, it can still be limited. When the budget of inducing variables is small, the sparse GP approximation may lead to inaccurate estimates of the posterior belief q(uk), which will degrade the quality of regularisation when learning new tasks. This is worrisome as in continual learning it is desirable to keep the size of the inducing set as small as possible.
One way to deal with this issue is to use a much larger set of inducing points for the current task or even maximise the full GP ELBO ∑Nk j=1 Eq(fk,j) log p(yk,j |fk,j) − KL(q(fk)||pθ(fk)) (i.e. by using as many inducing points as training examples), and once training is completed to distill the small subset Zk,uk ⊂ Xk, fk, and the corresponding marginal distribution q(uk) from q(fk), for subsequently regularising continual learning. However, carrying out this maximisation in the function space can be extremely slow since it scales as O(N3k ) per optimisation step. To our rescue, there is an alternative computationally efficient way to achieve this, by relying on the linear form of the kernel, that performs inference over the current task in the weight space. While this inference does not immediately provides us with the summary (induced points) for building the functional regularisation term, we can distill this term afterwards as discussed next. This allows us to address the continual learning aspect of the problem. Given that the current k-th task is represented in the weight space as fk(x;wk) = w > k φ(x; θ), wk ∼ N (0, σ2wI), we introduce a full Gaussian variational approximation q(wk) = N (wk|µwk ,Σwk), where µk is a K dimensional mean vector and Σwk is the corresponding K ×K full covariance matrix parameterised as Σwk = LwkL>wk . Learning the k-th task is carried out by maximising the objective in equation 4, with the only difference that the ELBO for the current task is now in the weight space. The objective becomes
F(θ, q(wk)) = Nk∑ j=1 Eq(wk)[log p(yk,j |w > k φ(xk,j ; θ))]− KL(q(wk)||p(wk))− k−1∑ i=1 KL(q(ui)||pθ(ui)),
where Eq(wk)[log p(yk,j |w>k φ(xk,j , θ))] can be re-written as one-dimensional integral and estimated using Gaussian quadrature. Once the variational distribution q(wk) has been optimised, together with the constantly updated feature parameters θ, we can rely on this solution to select inducing points Zk. See Section 2.4 for more detail. We also compute the posterior distribution over their function values uk according to q(uk) = N (uk|µuk , LukL>uk), where
µuk = ΦZkµwk , Luk = ΦZkLwk (5)
and the matrix ΦZk stores as rows the feature vectors evaluated at Zk. Subsequently, we store the k-th task summary (Zk, µuk , Luk) and use it for regularising continual learning of subsequent tasks, by always maximising the objective F(θ, q(wk)). Pseudo-code of the procedure is given in Algorithm 1.

heading: 2.4 SELECTION OF THE INDUCING POINTS
text: After having seen the k-th task, and given that it is straightforward to compute the posterior distribution q(uk) for any set of function values, the only issue remaining is to select the inducing inputs Zk. A simple choice that works well in practice is to select Zk as a random subset of the training inputs Xk. The question is whether we can do better with some more structured criterion.
Algorithm 1 Functional Regularised Continual Learning (FRCL) with task boundary detection Input: Feature vector φ(x; θ) with initial value of θ, task k = 0, starting_time(k) = 10. Construct output weights w0 and initialise variational parameters µw0 (e.g. around zero) and Lwk = I . for t = 1, 2, . . . , do
Receive next data minibatch (Xt,yt). Compute KL values `t = {`t,i}, for any xt,i ∈ Xt. if t - starting_time(k) > min_time_in(k) and StatisticalTest(`t, `old) is significant then
# A new task has been detected. Select inducing inputs Zk for current task. Compute the parameters of q(uk) from equation 5 and store them. Construct new output weights wk+1 and variational parameters (µwk+1 , Lwk+1). k = k + 1; starting_time(k) = t.
else `old = `t. end if Gradient step to update (θ, µwk , Lwk) by maximising F(θ, q(wk)).
end for
In our experiments we will investigate several criteria where the most effective one will be an unsupervised criterion that only depends on the training inputs, while the other supervised criteria are described in the Appendix. This unsupervised criterion quantifies how well we reconstruct the full kernel matrix KXk from the inducing set Zk and it can be expressed as the trace of the covariance matrix of the prior GP conditional p(fk|uk), i.e.
T (Zk) = tr ( KXk −KXkZKK −1 Zk KZkXk ) = Nk∑ j=1 [ k(xk,j , xk,j)− k>ZK ,xk,jK −1 Zk kZk,xk,j ] , (6)
where each k(xk,j , xk,j) − k>ZK ,xk,jK −1 Zk kZk,xk,j ≥ 0 is a reconstruction error for an individual training point. The above quantity appears in the ELBO in (Titsias, 2009), is also used in (Csato & Opper, 2002) and it has deep connections with the Nystróm approximation (Williams & Seeger, 2001) and principal component analysis. The criterion in equation 6 promotes finding inducing points Zk that are repulsive with one another and are spread evenly in the input space under a similarity/distance implied by the dot product of the feature vector φ(x; θk) (with θk being the current parameter values after having trained with the k-th task). An illustration of this repulsive property is given in Section 4.
To select Zk, we minimise T (Zk) by applying discrete optimisation where we select points from the training inputs Xk. The specific optimisation strategy we use in the experiments is to start with an initial random set Zk ⊂ Xk and then further refine it by doing local moves where random points in Zk are proposed to be changed with random points of Xk.

heading: 2.5 PREDICTION AND DIFFERENCES WITH WEIGHT SPACE METHODS
text: Prediction at any i-th task that has been encountered in the past follows the standard sparse GP predictive equations. Given a test data point xi,∗ the predictive density of its output value yi,∗ takes the form p(yi,∗) = ∫ p(yi,∗|fi,∗)pθ(fi,∗|ui)q(ui)duidfi,∗ = ∫ p(yi,∗|fi,∗)qθ(fi,∗)dfi,∗ where qθ(fi,∗) = N (fi,∗|µi,∗, σ2i,∗) is an univariate posterior Gaussian with mean and variance,
µi,∗ = µ > uiK −1 Zi ΦZ1φ(xi,∗; θ), σ 2 i,∗ = k(xi,∗, xi,∗) + k > Zixi,∗K −1 Zi [ LuiL > ui −KZi ] K−1Zi kZixi,∗ ,
where in µi,∗ we have explicitly written the cross kernel vector kZkxi,∗ = ΦZiφ(xi,∗; θ) (assuming σ2w = 1 for simplicity) to reveal a crucial property of this prediction. Specifically, given that fi,∗ = w > i φ(xi,∗; θ) the vector µ > uiK −1 Zi
ΦZi acts as a mean prediction for the task-specific parameter row vector w>i . As we learn future tasks and the parameter θ changes, this mean parameter prediction automatically adapts (recall that KZi = ΦZiΦ > Zi
and ΦZi vary with θ and only µui is constant) in order to counteract changes in the feature vector φ(xi,∗; θ), so that the overall prediction for the function value, i.e. µi,∗ = E[fi,∗], does not become obsolete. For instance, the prediction of the function values at the inducing inputs Zi always remains constant to our fixed/stored mean belief µui
since by setting xi,∗ = Zi the formula gives µ>uiK −1 Zi ΦZiΦ > Zi = µ>ui . Similar observations can be made for the predictive variances.
The above analysis reveals an important difference between continual learning in function space and in weight space, where in the latter framework task-specific parameters such as wi might not automatically adapt to counteract changes in the feature vector φ(x; θ), as we learn new tasks and θ changes. For instance, if as a summary of the task, instead of the function space posterior distribution q(ui), we had kept in memory the weight space posterior q(wi) (see Section 2.3), then the corresponding mean prediction on the function value, E[fi,∗] = µ>wiφ(xi,∗; θ), can become obsolete as φ(xi,∗; θ) changes and µwi remains constant.

heading: 3 DETECTING TASK BOUNDARIES USING BAYESIAN UNCERTAINTIES
text: So far we have made the assumption that task switches are known, which may not always be a realistic setting. Instead, we now introduce a novel approach for detecting task boundaries in continual learning arising naturally from our method by a simple observation: The GP predictive uncertainty grows as the model is queried far away from observed data, eventually falling back to the prior. When a minibatch of data {xi, yi}bi=1 from a new task arrives, we thus expect the distance between prior and posterior to be small (see Figure 2). Thus, a simple way to detect a change in the input distribution is to compare the GP univariate posterior density
q(f(xi)) = N (f(xi)|µi, σ2i ) ≈ ∫ p(f(xi)|f)p(f |y, X)df(xi),
where (µi, σ2i ) are predictive mean and variance, with the prior GP density p(f(xi)) = N (f(xi)|0, k(xi, xi)). This can be achieved by using a divergence measure between distributions such as the symmetrised KL divergence,
`i = 0.5 (∫ q(f(xi)) log q(f(xi))
p(f(xi)) df(xi) +
∫ p(f(xi)) log p(f(xi))
q(f(xi)) df(xi)
) , i = 1, . . . , b,
computed separately for any xi in the minibatch. Given that all distributions are univariate Gaussians the above can be obtained analytically. When each score `i is close to zero this indicate that the input distribution has changed so that a task switch can be detected. Each `i ≥ 0 can be thought of as expressing a degree of surprise about xi, i.e. the smaller is `i the more surprising is xi. Thus our idea has close links to Bayesian surprise (Itti & Baldi, 2006).
In order to use this intuition to detect task switches, we can perform a statistical test between the values {`i}bi=1 for the current batch and those from the previous batch {`oldi }bi=1 before making any updates to the parameters using the current batch. A suitable choice is Welch’s t-test (due to unequal variances), demanding that with high statistical significance the mean of {`i}bi=1 is smaller than the mean of {`oldi }bi=1. The ability to detect changes based on the above procedure arises from our framework as we construct posterior distributions over function values f(xi) that depend on inputs xi (while in contrast a posterior over weights alone does not depend on any input). Subsequently, these predictive densities contain information about the distribution of these inputs in the sense that when an xi is close to the
training inputs from the same task we expect reduced uncertainty, while for inputs of a different task we expect high uncertainty that falls back to the prior uncertainty.

heading: 4 EXPERIMENTS
text: We now test the scalability and competitiveness of our method on various continual learning problems, referring to the proposed approach as Functional Regularised Continual Learning (FRCL). Throughout this section, we will aim to answer three key questions:
(i) How does FRCL compare to state-of-the-art algorithms for Continual Learning? (ii) How important is a principled criterion for inducing point selection? How do varying
numbers of inducing points/task affect overall performance? (iii) If ground truth task boundaries are not given, does the detection method outlined in Section
3 succeed in detecting task changes?
In order to answer these questions, we consider experiments on three established Continual Learning classification problems: Split-MNIST, Permuted-MNIST and sequential Omniglot (Goodfellow et al., 2013; Zenke et al., 2017; Schwarz et al., 2018), described in the Appendix. FRCL methods have been implemented using GPflow (Matthews et al., 2017).
In addition to comparing our method with other approaches in the literature by quoting published results, we also show results for an additional baseline (BASELINE) corresponding to a simple replay-buffer method for Continual Learning (explained in the Appendix).

heading: 4.1 IS FRCL A COMPETITIVE MODEL FOR CONTINUAL LEARNING?
text: Addressing first question (i), we show results on the MNIST-variations in Table 1 and on the more challenging Omniglot benchmark in Table 2. Note that we also specify the inducing point optimisation criterion in brackets, i.e. FRCL (TRACE-TERM) corresponds to the loss in equation 6. We observe strong results on all benchmarks, setting a new state-of-the-art results on Permuted-MNIST & Omniglot while coming close to existing results on Split-MNIST. The improvement on the BASELINE shows that approximate posterior distributions over functions values can lead to more effective regularisation for CL compared to just having a replay buffer of input-output pairs. Furthermore, despite its simplicity, the simple BASELINE strategy performs competitively. In line with other results, we conclude that rehearsal-based ideas continue to provide a strong baseline. This also gives justification to a main motivation of this work:
To unite the two previously separate lines of CL work on rehearsal-based and Bayesian methods. Nevertheless, other methods may be more suitable when storing data is not feasible.

heading: 4.2 INDUCING POINT OPTIMISATION AND TASK SWITCH DETECTION
text: An appealing theoretical property of our method is the principled selection of inducing points through optimisation. Answering question (ii), we now proceed to investigate the importance of the criterion used as well as the dependence on the number of inducing points. These results are shown in Figure 3. Note that the definition of objectives Class. Error, ELBO & Log pred. density are given in the Appendix. In accordance with our intuition, we observe that optimisation becomes increasingly important as the number of inducing points is reduced. The results also give strong statistical motivation to use the trace-term motivated before. Further, as can be seen looking at the results for Log pred. density, a poorly chosen criterion may behave worse than random.
To provide an insight into the solutions obtained by the trace-term criterion, we provide a visualisation of inducing points in Figure 4. Remarkably, even though the objective is unsupervised, it results in a consistent allocation of one example per class. Furthermore, the optimised inducing points are spread across the input space as shown by the TNSE (Maaten & Hinton, 2008) visualisation, which is in line with the intuition that the objective encourages repulsive inducing points.
Finally, we answer question (iii) by first showing both the mean of the terms {`i}bi=1 (top) as well as the result of Welch’s t-test (bottom) between terms {`i}bi=1, {`oldi }bi=1 in Figure 5, using only a small number Omniglot alphabets and 1000 training iterations per task for illustrative purposes. We note that the intuition built up in Section 3 holds, with clear spikes being shown whenever the t-test returns a positive result.
Furthermore, we provide a quantitative comparison in Table 3. On the positive side, we note very strong results for Split- & Perm-MNIST and further observe that we find a similar t-test threshold value to apply to all dataset, making this an easy hyper-parameter to set. While the task boundary detection results for Omniglot are less strong, which may due to the smaller batch size (32 for Omniglot, ≥ 100 for the MNIST-versions), resulting in a noisier test result. Note that this could be easily mitigated by using a larger set {`oldi }, e.g. the last 10 minibatches, which would make this test more robust.
0.0
0.5
1.0
1.5
2.0 Mean symmetric. KL divergence between GP prior and posterior.
0 1000 2000 3000 4000 5000 10
0
10
20
30
40
Task 1: Korean Task 2: Ge'ez Task 3: Oriya Task 4: Tibetan Task 5: Jawi Arabic
Welch's t-test statistic between subsequent mini-batches.
Figure 5: Visualising KL terms and test statistics on multiple Omniglot tasks.
Table 3: Task boundary detection evaluated as a binary classification task (Positive labels corresponds to task switches).
Precision/Recall/F1 @ Threshold
Split-MNIST Perm-MNIST Omniglot
4.0 0.94/1.0/0.97 0.43/0.98/0.60 0.13/0.90/0.23 5.0 1.0/1.0/1.0 0.95/0.94/0.94 0.71/0.75/0.73 6.0 1.0/1.0/1.0 1.0/0.91/0.95 0.89/0.59/0.71 7.0 1.0/1.0/1.0 1.0/0.91/0.95 0.89/0.46/0.60 8.0 1.0/1.0/1.0 1.0/0.86/0.92 0.88/0.32/0.46

heading: 5 DISCUSSION
text: We introduced a functional regularisation approach for supervised continual learning that combines inducing point GP inference with deep neural networks. Our method constructs task-specific posterior beliefs or summaries on inducing inputs. Subsequently, the task-specific summaries allow us to regularise continual learning and avoid catastrophic forgetting. Our approach unifies the two extant approaches to continual learning, of parameter regularisation and replay/rehersal. Viewed from the regularisation perspective, our approach regularises the functional outputs of the neural network, thus avoid the brittleness due to representation drift. Viewed from a rehearsal method perspective, we provide a principled way of compressing data from previous task, by means of optimizing the selection of inducing points. By investigating the behaviour of the posterior beliefs, we also proposed a method for detecting task boundaries. All these improvements lead to strong empirical gains compared to state-of-the-art continual learning methods.
Regarding related work on online learning using GPs, notice that previous algorithms (Bui et al., 2017a; Csato & Opper, 2002; Csato, 2002) learn in an online fashion a single task where data from this task arrive sequentially. In contrast in this paper we developed a continual learning method for dealing with a sequence of different tasks.
A direction for future research is to enforce a fixed memory buffer (or a buffer that grows sub-linearly w.r.t. the number of tasks), in which case one would need to compress the summaries of all previous seen tasks into a single summary. Finally, while in this paper we applied the method to supervised classification tasks, it will be interesting to consider also applications in other domains such as reinforcement learning.

heading: ACKNOWLEDGEMENTS
text: We would like to thank Hyunjik Kim and Raia Hadsell for useful discussions and feedback.

heading: A TASK DESCRIPTIONS
text: Split-MNIST and Permuted-MNIST. Among a large number of diverse experiments in continual learning publications, two versions of the popular MNIST dataset have recently started to become increasingly popular benchmarks: Permuted- and Split-MNIST. In Permuted-MNIST (e.g. Goodfellow et al., 2013; Kirkpatrick et al., 2017; Zenke et al., 2017), each task is a variant of the initial 10-class MNIST classification task where all input pixels have undergone a fixed (random) permutation. The Split-MNIST experiment was introduced by Zenke et al. (2017): Five binary classification tasks are constructed from the classes in the following order: 0/1, 2/3, 4/5, 6/7, and 8/9.
Omniglot. To assess our method under more challenging conditions, we consider the sequential Omniglot task proposed for continual learning in Schwarz et al. (2018). Omniglot Lake et al. (2011) is a dataset of 50 alphabets, each with a varying number of classes/characters which we treat of as a sequence of distinct classification problems. As suggested in Schwarz et al. (2018), for Omniglot, we apply data-augmentation and use the same train/validation/test split. Following the same experimental setup proposed, we used an identical convolutional network to construct the feature vector φ(x; θ). Results reported are obtained by training on the union of training and validation set after choosing any hyper-parameters based on the validation set. Note that all experiments were run with data processing and neural network construction code kindly provided by the authors in Schwarz et al. (2018), ensuring directly comparable results.
Given that Permuted-MNIST and Omniglot are multi-class classification problems, where each k-th task involves classification over Ck classes, we need to generalise the model and the variational method to deal with multiple GP functions per task. This is outlined in the next section.

heading: B EXTENSION TO MULTI-CLASS (OR MULTIPLE-OUTPUTS) TASKS
text: For multi-class classification problems, such as permuted MNIST and Omniglot considered in our experiments, where in general each k-th task involves classification over varying number of classes, we need to extend the method to deal with multiple functions per task. For instance, assume that the k-th task is a multi-class classification problem that inolves Ck classes. To model this we need Ck independent draws from the GP, such that f ck(x) ∼ GP(0, k(x, x′)) with c = 1, . . . , Ck, which are combined based on a multi-class likelihood such as softmax. While in the main text we presented the method assuming a single GP function per task, the generalisation to multiple functions is straightforward by assuming that all variational approximations factorise across different functions/classes. For example, the variational distribution over all inducing variables Uk = {uck} Ck c=1
takes the form q(Uk) = ∏Ck c=1 q(u c k) and similarly the variational approximation over the task weights Wk = {wck}, needed in the ELBO in Section 3.2, also factorises across classes. Notice also that all inducing variables Uk are evaluated on the same inputs Zk. Furthermroe, the KL regularization term for each task takes the form of a sum over the different functions, i.e. KL(q(Uk)||pθ(Uk)) =∑Ck c=1 KL(q(u c k)||pθ(uck)).

heading: C BASELINE MODEL
text: The BASELINE model (see main text) is based on storing an explicit replay buffer (ỹi, X̃i), i.e. a subset of the training data where ỹi ⊂ yi and X̃i ⊂ Xi, for each past task. Then, at each step when we encounter the k-th task training is performed by optimising an unbiased estimate of the full loss
(i.e. if we had all k tasks at once), given by
L(θ, w1:k) = `k(yk, Xk;wk, θ) + k−1∑ i=1 Ni Mi `i(ỹi, X̃i;wi, θ),
where each `i(·) is a task-specific loss function, such as cross entropy for classification, and each scalar NiMi corrects for the bias on the loss value caused by approximating the initial full loss by a random replay buffer of size Mi. All output weights wi:k of the current and old tasks, in the multi-head architecture, are constantly updated together with the feature parameter vector θ. Also at each step a fresh set of output weights is constructed in order to deal with the current task.

heading: D SELECTION OF THE INDUCING POINTS AND OPTIMISATION CRITERIA
text: After having seen the k-th task, and given that it is straightforward to compute the posterior distribution q(uk) for any set of function values, the only issue remaining is to select the inducing inputs Zk. A simple choice that works well in practice is to select Zk as a random subset of the training inputs Xk. The question is whether we can do better with some more structured criterion. In our experiments we will investigate three supervised criteria, that make use of class labels, and one unsupervised that only depends on inputs and the neural network feature vector. All criteria below optimise over Zk using discrete search.
The first supervised criterion is to minimize the negative average log predictive density,
Log pred. density(Zk) = − 1 |Xk \ Zk| ∑
xk,j∈Xk\Zk
log p(yk,j |Zk), (7)
computed at all remaining training inputs by excluding the selected inducing inputs Zk. Each predictive density p(yk,j |Zk) is obtained through the inducing inputs Zk and it takes the form
p(yk,j |Zk) = ∫ p(yk,j |fk,j)q(uk|Zk)duk = ∫ p(yk,j |fk,j)N (uk|µuk , LukL>uk)duk, (8)
which for the classification problems, where the likelihood p(yk,j |fk,j) is not Gaussian, is computed numerically by one-dimensional Gaussian quadrature. Notice that for the multi-class case the predictive density is obtained by integrating over q(Uk) ≡ q(Uk|Zk) (see previous Appendix B). Modern GP packages, such as GPflow, provide efficient implementation of the above predictive densities.
The second supervised criterion is similar, but − log p(yk,j |ZK) is replaced by classification error score I(yk,j 6= y∗k,j),
Class. Error(Zk) = 1 |Xk \ Zk| ∑
xk,j∈Xk\Zk
I(yk,j 6= y∗k,j), (9)
where y∗k,j is the predicted label, which is obtained by first computing the predictive density p(yk,j |Zk) conditional on a given set of inducing inputsZk and then choosing y∗k,j by taking argmax (i.e. selecting the class with the largest predictive probability).
Both criteria above essentially capture how well we predict the remaining points from the inducing points Zk, thus good choices for Zk are presumably those that lead to good prediction for the remaining points.
The third supevised criterion is the standard sparse GP variational lower bound,
ELBO(Zk) = Nk∑ j=1 Eq(fk,j)[log p(yk,j |fk,j)]− KL(q(uk|Zk)||pθ(uk|Zk)), (10)
which is viewed purely as a function of Zk while all remaining quantities are constant. This criterion tries to optimise over Zk in order to approximate as best as possible the marginal likelihood on the training data and it is the one used by all variational sparse GP training methods (although there the ELBO is maximized jointly over both Zk with the remaining parameters).
The last unsupervised criterion (also discussed in the main paper) expresses how well we reconstruct the full kernel matrix KXk from the inducing set Zk, which can be described by
Trace term(Zk) = ( KXk −KXkZkK −1 Zk KZkXk ) = Nk∑ j=1 [ k(xk,j , xk,j)− k>ZK ,xj,kK −1 Zk kZk,xi,k ] ,
(11) where each k(xk,j , xk,j) − k>ZK ,xj,kK −1 Zk kZk,xi,k ≥ 0 is a reconstruction error for an individual training point. The above quantity appears in the ELBO in (Titsias, 2009), is also used in (Csato & Opper, 2002) and it has deep connections with the Nystróm approximation (Williams & Seeger, 2001) and principal component analysis. The criterion in equation 6 promotes finding inducing points Zk that are repulsive with one another and are spread evenly in the input space under a similarity/distance implied by the dot product of the feature vector φ(x; θk) (with θk being the current parameter values after having trained with the k-th task). An illustration of this is given in the Figure 4.
Figure 4 in the Experiments Section illustrates the optimisation of the inducing inputs in PermutedMNIST. Also, Figure 6(a) shows randomly chosen (3 examples per class) for the Greek alphabet/task in the Omniglot experiments, while Figure 6(b) shows the corresponding points after have been optimised by using the trace criterion above.
Figure 3 shows the evolution of the performance for all selection criteria on individual tasks for Split-MNIST and Permuted-MNIST.

heading: E TASK BOUNDARY DETECTION
text: In order to apply the task-boundary detection for multiple classes, we perform Welch’s t-test separately for each function. We consider mean, median and the maximum of those t-test results to make a decision. Furthermore, we also found that conducting the test in log-space significantly improved robustness of the test results and thus find that tests using the max over functions in log-space allow for higher thresholds and more robust results. In order to give justification to those conclusions we show our experiments on Permuted-MNIST in Figure 7, once again using 10 random task permutations.
Note also that after a task boundary was detected, we do not consider a new test for the next 10 iterations.
At this point it should be said that while the boundary detection method is simple, cheap to compute, unsupervised and general, it is not applicable when the Continual Learning faces a continuum of tasks without clear switches or when the input distribution is constant and different tasks merely correspond to varying labels.

heading: F EXPERIMENTAL DETAILS
text: Experimental details for all experiments are shown in Tables 4, 5 and 6. Note that for the MNIST results, we obtain final results after optimising hyperparameters on the validation set and using those values to train on the union of training & validation set.
For Omniglot on the other hand, we report final test-set results only after training on the training set in order to remain consistent with the results in Schwarz et al. (2018)

heading: G COMPARISON TO VCL ON OMNIGLOT
text: To provide a further comparison to VCL (Nguyen et al., 2017), we show results on Omniglot using Multi-Layer Perceptrons (MLP). A comparison with MLPs is due to the fact that reliable variational inference methods for CNNs (which are usually used for Omniglot) are yet to be developed. All our results for VCL are obtained using code provided by the authors.
For all experiments, we used an MLP with 4 hidden layers of 256 units each and ReLU activations, a batch size of 100 and the Adam Optimiser (Step size of 0.001 for VCL and 0.0001 for FRCL). We optimised both types of algorithms independently and found the following parameters:
1. VCL: 100 training epochs per task, 50 adaptation epochs to coreset, Multi-Head 2. FRCL (TRACE): 1500 training steps per task, 2000 discrete optimization steps, inducing
points initialised as a uniform distribution over classes.

Review summaries:
The paper develops a continual learning method based on Gaussian Processes (GPs) applied in the way introduced by prior work as Deep Kernel Learning (DKL). The proposed method summarizes tasks as sparse GPs and use them as regularizers for the subsequent tasks in order to avoid catastrophic forgetting. Salleviating the instability resulting from the representation drift.