{'fraction_non_alphanumeric': 0.044689889549702634, 'fraction_numerical': 0.013593882752761258, 'mean_word_length': 4.23979818964238, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SOLA: Continual Learning with Second-Order Loss Approximation

Neural networks have achieved remarkable success in many cognitive tasks. However, when they are trained sequentially on multiple tasks without access to old data, it is observed that their performance on old tasks tend to drop significantly after the model is trained on new tasks. Continual learning aims to tackle this problem often referred to as catastrophic forgetting and to ensure sequential learning capability. We study continual learning from the perspective of loss landscapes and propose to construct a second-order Taylor approximation of the loss functions in previous tasks. Our proposed method does not require any memorization of raw data or their gradients, and therefore, offers better privacy protection. We theoretically analyze our algorithm from an optimization viewpoint and provide a sufficient and worst-case necessary condition for the gradient updates on the approximate loss function to be descent directions for the true loss function. Experiments on multiple continual learning benchmarks suggest that our method is effective in avoiding catastrophic forgetting and in many scenarios, outperforms several baseline algorithms that do not explicitly store the data samples.


Introduction
Neural networks are achieving human-level performance on many cognitive tasks including image classification [23] and speech recognition [16]. However, as opposed to humans, their acquired knowledge is comparably volatile and can be easily dismissed. Especially, the catastrophic forgetting phenomenon refers to the case when a neural network forgets the past tasks if it is not allowed to retrain or reiterate on them again [13,28].
Continual learning is a research direction that aims to solve the catastrophic forgetting problem. Recent works tried to tackle this issue from a variety of perspectives. Regularization methods (e.g., [21,47]) aim to consolidate the weights that are important to previous tasks while expansion based methods (e.g., [40,46]) typically increase the model capacity to cope with the new tasks. Repetition based methods (e.g., [26,5]) usually do not require additional and complex modules, however, they have to maintain a small memory of previous data and use them to preserve knowledge. Unfortunately, the performance boost of repetition based methods comes at the cost of storing previous data which may be undesirable whenever privacy is important. To address this issue, authors in [8] proposed a method to work with the gradients of the previous data to constrain the weight updates; however, this may still be subject to privacy issues as the gradient associated with each individual data point may disclose information about the raw data.
In this paper, we study the continual learning problem from the perspective of loss landscapes. We explicitly target minimizing an average over all tasks' loss functions. The proposed method stores neither the data samples nor the individual gradients on the previous tasks. Instead, we propose to construct an approximation to the loss surface of previous tasks. More specifically, we approximate the loss function by estimating its second-order Taylor expansion. The approximation is used as a surrogate added to the loss function of the current task. Our method only stores information based on the statistics of the entire training dataset, such as full gradient and full Hessian matrix (or its low rank approximation), and thus better protects privacy. In addition, since we do not expand the model capacity, the neural network structure is less complex than that of expansion based methods.
We study our algorithm from an optimization perspective, and make the following theoretical contributions: • We prove a sufficient and worst-case necessary condition under which by conducting gradient descent on the approximate loss function, we can still minimize the actual loss function.
• We further provide convergence analysis of our algorithm for both non-convex and convex loss functions. Our results imply that early stopping can be helpful in continual learning.
• We make connections between our method and elastic weight consolidation (EWC) [21].
In addition, we make the following experimental contributions: • We conduct a comprehensive comparison among our algorithm and several baseline algorithms [21,5,8] on a variety of combinations of datasets and models. We observe that in many scenarios, especially when the learner is not allowed to store the raw data samples, our proposed algorithm outperforms them. We also discuss the conditions under which the proposed method or any of the alternatives are effective.
• We provide experimental evidence validating the importance of accurate approximation of the Hessian matrix and discuss scenarios in which early stopping is helpful for our algorithm.
The expansion based methods allocate new neurons or layers or modules to accommodate new tasks while utilizing the shared representation learned from previous ones [40,45,46,25,17]. Although being a very natural approach the mechanism of dynamic expansion can be quite complex and can add considerable overhead to the training process.
The repetition and memory based methods store previous data or, alternatively, train a generative model of them and replay samples from them interleaved with samples drawn from the current task [41,19,48,38,27,26,8]. They achieve promising performance however at the cost of higher risk of users' privacy by storing or learning a generative model of their data.
The regularization based approaches impose limiting constraints on the weight updates of the neural network according to some relevance score for previous knowledge [21,33,42,39,29,47,35]. These methods provide a better privacy guarantee as they do not explicitly store the data samples.
In general, SOLA also belongs to this category as we use the second-order Taylor expansion as the regularization term in new tasks. Many of the regularization methods are derived from a Bayesian perspective of estimating the posterior distribution of the model parameters given the data from a sequence of tasks [21,33,42,39]; some of these methods use other heuristics to either estimate the importance of the weights of the neural network [47,35] or implicitly limit the capacity of the network [29]. Similar to our approach, several regularization based methods use quadratic functions as the regularization term, and many of them use the diagonal form of quadratic functions [21,47,35]. In Section 5.3, we demonstrate that in some cases, the EWC algorithm [21] can be considered as the diagonal approximation of our approach. Here, we note that the diagonal form of quadratic regularization has the drawback that it does not take the interaction between the weights into account.
Among the regularization based methods, the online Laplace approximation algorithm [39] is the most similar one to our proposed method. Despite the similarity in the implementations, the two algorithms are derived from very different perspectives: the online Laplace approximation algorithm uses a Bayesian approach that approximates the posterior distribution of the weights with a Gaussian distribution, whereas our algorithm is derived from an optimization viewpoint using Taylor approximation of loss functions. More importantly, the Gaussian approximation in [39] is proposed as a heuristic; whereas in this paper, we provide rigorous theoretical analysis on how the approximation error affects the optimization procedure. We believe that our analysis provides deeper insights to the loss landscape of continual learning problems, and explains some important implementation details such as early stopping.
We also note that continual learning is broader than just solving the catastrophic forgetting and is connected to many other areas such as meta learning [37], few-shot Learning [44,12], learning without explicit task identifiers [36,2], to name a few.

Problem formulation
We consider a sequence of K supervised learning tasks T k , k ∈ [K]. 1 For task T k , there is an unknown distribution D k over the space of feature-label pairs X × Y. Let W ⊆ R d be a model parameter space, 2 and for the k-th task, let k (w; x, y) : W → R be the loss function of w associated with data point (x, y). The population loss function of task T k is defined as L k (w) := E (x,y)∼D k k (w; x, y). Our general objective is to learn a parametric model with minimized population loss over all the K tasks. More specifically, in continual learning, the learner follows the following protocol: When learning on the k-th task, the learner obtains access to n k data points (x k,i , y k,i ), i ∈ [n k ] sampled i.i.d. according to D k and we define L k (w) := 1 n k n k i=1 k (w; x k,i , y k,i ) as the empirical loss function; the learner then updates the model parameter w using these n k training data, and after the training procedure is finished, the learner loses access to the training data, but can store some side information about the task. Our goal is to avoid forgetting previous tasks when trained on new tasks by utilizing the side information. We provide details of our algorithm design in the next section.

Our approach
To measure the effectiveness of a continual learning algorithm, we use a simple criterion that after each task, we hope the average population loss over all the tasks that have been trained on to be small, i.e., for every k ∈ [K], after training on T k , we hope to solve min w∈W 1 k k k =1 L k (w). Since minimizing the loss function is the key to training a good model, we propose a straightforward method for continual learning: storing the second-order Taylor expansion of the empirical loss function, and using it as a surrogate of the loss function for an old task when training on new tasks. We start with a simple setting. Suppose that there are two tasks, and at the end of T 1 , the we obtain a model w 1 . Then we compute the gradient and Hessian matrix of L 1 (w) at w 1 , and construct the second-order Taylor expansion of L 1 (w) at w 1 : When training on T 2 , we try to minimize 1 2 ( L 1 (w) + L 2 (w)). The basic idea of this design is that, we hope in a neighborhood around w 1 , the quadratic function L 1 (w) stays as a good approximation of L 1 (w), and thus approximately we still minimize the average of the empirical loss functions 1 2 ( L 1 (w) + L 2 (w)), which in the limit generalizes to the population loss function 1 2 (L 1 (w) + L 2 (w)). We rely on the assumption that the second-order Taylor approximation of loss function can capture their local geometry well. For a general nonlinear function and arbitrary displacement, this approximation can be over-simplistic, however, we refer to the abundance of observations for modern neural networks that are seen to be well-behaved with flat and wide minima [6,14]. Moreover, the assumption of well-behaved loss around tasks' local minima also forms the basis of a few other continual learning algorithms such as EWC [21] and OGD [8].
Formally, let w k be the model that we obtain at the end of the k-th task. We define the approximation of the sum of the first k empirical loss functions as where H k denotes the Hessian matrix ∇ 2 L k ( w k ) or its low rank approximation, and c k ∈ R is a constant that does not depend on w. We construct L k (w) at the end of task k, and when training on task k + 1, we minimize 1 k+1 ( L k (w) + L k+1 (w)). In the following, we name our algorithm SOLA, an acronym for Second-Order Loss Approximation.
As we can see, in the SOLA algorithm, after each task, if we choose to use the exact Hessian matrix, i.e., , H k = ∇ 2 L k ( w k ), it suffices to update A k and b k in memory, and thus the memory cost of the algorithm is O(d 2 ), which does not grow with the number of tasks. However, in practice, especially for overparameterized neural networks, the dimension d of the model is usually large, and thus the storage cost of memorizing the Hessian matrix can be high. Recent studies have shown that the Hessian matrices of the loss functions of a deep neural networks are usually approximately low rank [11]. If we choose H k as a rank-r approximation of ∇ 2 L k ( w k ), we need to keep accumulating the low rank approximations of the Hessian matrices in order to construct F k (w), and at the end of task k, the memory cost is O(krd), which in practice can be much smaller than that of using the exact Hessian matrices. We formally demonstrate our approach in Algorithm 1, and the methods that use exact Hessian matrices and low rank approximation of them are presented as options I and II, respectively. Moreover, we can use a recursive implementation for the low rank approximation, and the memory cost can be further reduced to O(rd), which does not grow with k. We present the details of the recursive implementation in Section 6.
Algorithm 1 Continual learning with second-order loss approximation (SOLA) 1: Input: initial weights w 0 , learning rate η, the number of tasks K, the rank of Hessian approximation r (for option II) 2: for k = 1, 2, . . . , K do 3: access training data for the k-th task (x k,i , y k,i ), i ∈ [n k ], w ← w k−1

4:
while termination condition not satisfied do 5: compute (stochastic) gradient of current loss ∇ L k (w) 6: compute gradient of loss function approximation

Theoretical analysis
In this section, we provide theoretical analysis of our algorithm. As we can see, the key idea in our algorithm is to approximate the loss functions of previous tasks using quadratic functions. This leads to the following theoretical question: By running gradient descent algorithm on an approximate loss function, can we still minimize the actual loss that we are interested in?
For the purpose of theoretical analysis, we make a few simplifications to our setup. Without loss of generality, we study the training process of the last task T K , and still use w k to denote the model parameters that we obtain at the end of the k-th task. We use the loss function approximation in (1), but for simplicity we ignore the finite-sample effect and replace the empirical loss function with the population loss function, i.e., we define where H k represents ∇ 2 L k ( w k ) or its low rank approximation. The reason for this simplification is that our focus is the optimization aspect of the problem, while the generalization aspect can be tackled by tools such as uniform convergence [30]. As discussed, during the training of the last task, we have access to the approximate loss function F (w) : , whereas the actual loss function that we care about is F (w) := 1 K K k=1 L k (w). We also focus on gradient descent instead of its stochastic counterpart. In particular, let w 0 := w K−1 be the initial model parameter for the last task. We run the following update for t = 1, 2, . . . , T : We use the following standard notions for differentiable function f : We make the assumptions that the loss functions are smooth and Hessian Lipschitz. We note that the Hessian Lipschitz assumption is standard in analysis of non-convex optimization [31,18].
We also assume that the error between the matrices H k and

Sufficient and worst-case necessary condition for one-step descent
We begin with analyzing a single step during training. Our goal is to understand by running a single step of gradient descent on F (w), whether we can minimize the actual loss function F (w). More specifically, we have the following result. Theorem 1. Under Assumptions 1 and 2, and suppose that in the t-th iteration, we observe and the learning rate satisfies η ≤ 2(1−1/c) µ , then we have We prove Theorem 1 in Appendix A. Here, we emphasize that this result does not assume any convexity of the loss functions. The theorem provides a sufficient condition (4), under which by running gradient descent on F , we can still minimize the true loss function F . Intuitively, this condition requires the gradient of F to be large enough, such that the magnitude of the gradient is larger than the error caused by the inexactness of the loss function. In Proposition 1 below, we will see that this condition is also necessary in the worst-case scenario, at least for the case where K = 2.
More specifically, we can construct cases in which (4) is violated and the gradients of F (w) and F (w) have opposite directions.
We prove Proposition 1 in Appendix B. In addition, we note that Theorem 1 also implies that as training going on and ∇ F (w t ) 2 decreasing, it is beneficial to decrease the learning rate η, since when c decreases, the upper bound on η that guarantees the decay of F (i.e., 2(1 − 1/c)/µ) also decreases. We notice that the importance of learning rate decay for continual learning has been observed in some empirical study recently [29].

Convergence analysis
Although the condition in (4) provides us with insights on the dynamics of the training algorithm, it is usually hard to check this condition in every step, since we may not have good estimates of δ and ρ. A practical implementation is to choose a constant learning rate along with an appropriate number of training steps. In this section, we provide bounds on the convergence behavior of our algorithm with a constant learning rate and T iterations, both for non-convex and convex loss functions. These results imply that early stopping can be helpful, and provide a theoretical treatment of the very intuitive fact that the more iterations one optimizes for the current task the more forgetting can happen for the previous ones. We begin with a convergence analysis for non-convex loss functions in Theorem 2, in which we use the common choice of learning rate 1/µ for gradient descent on smooth functions [4].
, and F * := min w∈W F (w). Then, under Assumptions 1 and 2, after running T iterations of the gradient descent update (3) with learning rate η = 1/µ, we have We prove Theorem 2 in Appendix C. Unlike standard optimization analysis, the average norm of the gradients does not always decrease as T increases, when δ = 0 or ρ = 0. Intuitively, as we move far from the points where we conduct Taylor expansion, the gradient of F becomes more and more inaccurate, and thus we need to stop early. In Section 7, we provide experimental evidence.
When the loss functions are convex, we can prove a better guarantee which does not have the O( √ T ) and O(T ) terms as in Theorem 2. More specifically, we have the following assumption and theorem. Assumption 3. L k (w) is convex and H k 0, ∀k ∈ [K].
Theorem 3 (convex). Suppose that Assumptions 1, 2, and 3 hold, and define F * = min w∈W F (w), w * ∈ arg min w∈W F (w), w * ∈ arg min w∈W F (w), and D := w 0 − w * 2 . After running T iterations of the gradient descent update (3) with learning rate η = 1/µ, we have where α = 2µ D 2 , and β = 1 We prove Theorem 3 in Appendix D. As we can see, if δ = 0 or ρ = 0, we still cannot guarantee the convergence to the true minimum of F , due to the inexactness of F . On the other hand, if the loss functions are quadratic and we save the full Hessian matrices, i.e., δ = ρ = 0, as we have full information about previous loss functions, we can recover the standard O(1/T ) convergence rate for gradient descent on convex and smooth functions.

Connection to EWC
The elastic weight consolidation (EWC) algorithm [21] for continual learning is proposed based on the Bayesian idea of estimating the posterior distribution of the model parameters. Interestingly, we notice that our algorithm has a connection with EWC, although their basic ideas are quite different. More specifically, we show that in some cases, the regularization technique that the EWC algorithm uses can be considered as a diagonal approximation of the Hessian matrix of the loss function. Suppose that in the k-th task, the data points are samples from a probabilistic model with the likelihood function being p k (x, y|w * k ), and we use negative log-likelihood as the loss function, i.e., k (w; x, y) = − log p k (x, y|w). Suppose that at the end of this task, we obtain the ground truth model parameter w * k . Then we know that ∇L k (w * k ) = 0, and that the Fisher information of the i-th coordinate of w * k is I(w * k,i ) = ∂ 2 i L k (w * k ). The EWC algorithm constructs a regularization term as a proxy of the loss function of the k-th task, and uses it in the following tasks. As we can see, in this case, the quadratic regularization in EWC is a diagonal approximation of the quadratic term in our loss function approximation approach.

A recursive implementation
As we have seen, one drawback of the SOLA algorithm with low rank approximation in Section 4 is that the memory cost grows with the number of tasks. In this section, we present a more practical and memory efficient implementation of SOLA with low rank approximation. Recall that L k (w) is the empirical loss function for the k-th task, k ∈ [K]. We then define the loss function approximation L k (w) in a recursive way. We begin with L 0 (w) ≡ 0, and for every k = 1, 2, . . . , K, we define where Q k− 1 2 ∈ R d×d is a rank-r approximation of the Hessian matrix ∇ 2 L k− 1 2 ( w k ). This means that at the end of task k, we compute the second-order Taylor expansion of the approximate loss function L k−1 (w) + L k (w) at w k , with the Hessian matrix being replaced by the low rank approximation Q k− 1 2 . Thus, we only need to store ∇ L k− 1 2 ( w k ) and Q k− 1 2 , and the memory cost is O(rd), which does not grow with k. We formally present this approach in Algorithm 2. In our experiments in Section 7, we use the recursive implementation for SOLA with low rank approximation.
Algorithm 2 Recursive implementation of SOLA with low rank approximation 1: Input: initial weights w 0 , learning rate η, the number of tasks K, rank r 2: for k = 1, 2, . . . , K do 3: access training data for the k-th task (x k,i , y k,i ), i ∈ [n k ], w ← w k−1

Experiments
We implement the experiments with TensorFlow [1]. When computing the exact or the low rank approximation of the Hessian matrix, we treat each tensor in the model independently; in other words, we compute the block diagonal approximation of the Hessian matrix. This technique has the benefit that the Hessian computation is independent of the model architecture and has been used in recent studies on second-order optimization [15]. We use the recursive implementation for SOLA with low rank approximation. In the following, we denote SOLA with exact Hessian matrix and low rank approximation by SOLA-exact and SOLA-prox, respectively. As for the calculation of the low rank matrix, we make use of Hessian-vector product and provide details in Appendix E.
Datasets. We use multiple standard continual learning benchmarks created based on MNIST [24] and CIFAR-10 [22] datasets, i.e., Permuted MNIST [13], Rotated MNIST [26], Split MNIST [47], and Split CIFAR (similar to a dataset in [5]). In Permuted MNIST, for each task, we choose a random permutation of the pixels of MNIST images, and reorder all the images according to the permutation. We use 5-task Permuted MNIST in the experiments. In Rotated MNIST, for each task, we rotate the MNIST images by a particular angle. In our experiments, we choose a 5-task Rotated MNIST, with the rotation angles being 0, 10, 20, 30, and 40 degrees. For Split MNIST, we Split the 10 labels Table 1: Average test accuracy (%) ± std. P-MNIST, R-MNIST, and S-MNIST represent Permuted, Rotated, and Split MNIST datasets, respectively. Boldface numbers correspond to best result among algorithms that do not store raw data points, i.e., excluding multi-task and A-GEM. For MLP [10,10], in OGD, we store 200 gradient samples for each task, and in SOLA-prox, we use r = 500; for other models, in OGD, we store 280 gradient samples for each task, and in SOLA-prox, we use r = 700.  of the MNIST dataset to disjoint subsets, and for each task, we use the MNIST data whose labels belong to a particular subset. In this paper, we use a 5-task Split MNIST, and the subsets of labels are Architecture. We use both multilayer perceptron (MLP) and convolutional neural network (CNN).
In most cases, we use MLP with two hidden layers, sometimes denoted by MLP[x, y], with x and y being the number of hidden units. We may use CNN-x to denote a CNN model with x convolutional layers, and provide details of the model in Appendix F. For Split MNIST and Split CIFAR, we use MLP and CNN models with a multi-head structure similar to what has been used in [5,8]. In the multi-head model, instead of having 10 logits in the output layer, we use separate heads for different tasks, and each head corresponds to the classes of the associated task. During training, for each task, we only optimize the cross-entropy loss over the logits and labels of the corresponding output head.
Baselines. We compare SOLA algorithm with several baselines: the vanilla algorithm which runs SGD over all the tasks without storing any side information; the multi-task algorithm which assumes access to all the training data of previous tasks; the repetition based A-GEM algorithm [5], which stores a subset of data samples from the previous tasks and forms constrained optimization algorithms when training on new tasks; the regularization based EWC algorithm [21] discussed in Section 5.3; and the orthogonal gradient descent (OGD) algorithm [8] that stores the gradients in previous tasks and forms a constrained optimization algorithm. 3 Among them, our algorithm, along with the vanilla, EWC, and OGD algorithms do not explicitly store the raw data samples. Following prior works [5,8,21], we choose a learning rate of 10 −3 and a batch size of 10. For all the results that we report, we present the average result over 10 independent runs, as well as the standard deviation (as the shaded areas in the figures).
Results. We provide a comprehensive comparison among SOLA and the baseline algorithms with a variety of combinations of datasets and models. Tables 1 and 2 present the results for MNIST-based datasets and Split CIFAR, respectively. We make a few notes before discussing the results. First, the multi-task algorithm uses all the data of previous tasks, which serves as an upper bound for the performance of continual learning algorithms. Second, since the A-GEM algorithm stores a subset of data samples from previous tasks, it is not completely fair to compare A-GEM with algorithms that do not store raw data. However, here we still report the results for A-GEM for reference, and in A-GEM we store 200 data points for each task. Third, since the performance of the algorithms depends on the number of epochs that we train for each task, we treat this quantity as a tuning parameter, and for each algorithm, we report the result corresponding to the best epoch choice for its performance. In particular, for MNIST-based datasets, we choose epoch from {1, 5, 10, 20}, and for Split CIFAR, we choose from {1, 5, 10, 20, 40}. Due to memory constraints, we only implement SOLA-exact on small models such as MLP [10,10] and CNN-2. We conclude from the results as follows: • If it is allowed to store raw data, repetition based algorithm such as A-GEM should be the choice.
This remarks the importance of the information contained in the raw data samples. In some cases we observe that SOLA outperforms A-GEM, e.g., on MLP [10,10]. However, we expect that the performance of A-GEM can be improved if more data are stored in memory.
• If it is not allowed to store raw data due to privacy concerns, then in many scenarios, SOLA outperforms the baseline algorithms. In particular, on MNIST-based datasets, SOLA-exact or SOLA-prox achieves the best performance in 6 out of 7 settings.
• On Split CIFAR, we observe mixed results. When the model is relatively small (CNN-2) and we can store the exact Hessian matrix, SOLA-exact achieves the best performance. On a relatively large CNN model (CNN-6), we observe that none of the continual learning algorithms (EWC, OGD, SOLA) significantly outperforms the vanilla algorithm. On a large MLP, we observe that OGD performs the best and the result for SOLA-prox becomes worse. We believe the reason is that since in this experiment we only use 200 eigenvectors to approximate a Hessian matrix with very high dimensions, the approximation error is so large that SOLA-prox cannot find a descent direction that is close to the true gradient. This remarks the importance of future study of SOLA on models with more complicated structure or higher dimensions.
Performance vs approximation. We study how the approximation of Hessian matrices affects the performance of SOLA-prox. In particular, we choose different values of the rank r in SOLA-prox and investigate its correlation with the final average test accuracy. Our theory implies that when the approximation of Hessian matrices is better, i.e., smaller δ, the final performance is better. Our experiments validate this point. Figure 1a and Figure 1b show that, as we increase r, i.e., using more eigenvectors to approximate the Hessian matrix, the average test accuracy over all tasks improves.
Early stopping. Our theoretical analysis in Section 5 implies that early stopping can be helpful for SOLA. Here, we discuss empirical evidence. As we can see from Figure 1c, on Permuted MNIST with MLP [10,10], the average test accuracy of SOLA-prox becomes worse if we train more than 5 epochs per task; similarly, from Figure 1a, we can also see that training each task for more epochs can hurt the performance of MLP[100, 100] on Permuted MNIST. However, this phenomenon is less severe on Rotated MNIST. In Figure 1d, for SOLA-prox with r = 100, we observe one case where the average test accuracy gradually decreases as we increase the number of epochs per task.
Moreover, we notice that we did not observe this phenomenon for SOLA-exact. Hence, we draw the conclusion that the importance of early stopping for SOLA depends on how different the tasks are and how well we approximate the Hessian matrix. In Permuted MNIST, the pixels are randomly shuffled when switching to new tasks, whereas in Rotated MNIST we only rotate the images by 10 degrees; thus early stopping is more important for Permuted MNIST. On the other hand, if we store the exact value of the Hessian matrix (δ = 0 in Theorem 1), the approximation error of the gradients can be small, and thus we can train more epochs on new tasks. In addition, we note that it has been observed that early stopping is typically helpful for other continual learning algorithms [8].

Conclusions
We propose the SOLA algorithm based on the idea of loss function approximation. We establish theoretical guarantees, make connections to the EWC algorithm, and present experimental results showing that in many scenarios, our algorithm outperforms several baseline algorithms, especially among the ones that do not explicitly store the raw data samples. Future directions include studying SOLA on broader classes of neural network architectures and parameter spaces with higher dimensions.

A Proof of Theorem 1
We first provide a bound for the difference between the gradients of ∇ F (w) and ∇F (w).
We prove Lemma 1 in Appendix A.1. Since the loss functions for all the tasks L k (w) are µ-smooth, we know that F (w) is also µ-smooth. Then we have Therefore, as long as Then we can complete the proof by combining (7) with Lemma 1.

A.1 Proof of Lemma 1
By the definition of F (w), for some ξ k ∈ [0, 1], k ∈ [K − 1], we have where the second equality is due to Lagrange's mean value theorem. Then, according to Assumptions 1 and 2, we have Then, according to triangle inequality, we obtain

B Proof of Proposition 1
We first note that it suffices to construct F (w) and F (w), as one can always choose L 2 (w) ≡ 0 and then the construction of F (w) and F (w) is equivalent to that of L 1 (w) and L 1 (w). Let w 1 = 0, One can easily check that F ( w 1 ) = F ( w 1 ), F ( w 1 ) = F ( w 1 ), and |F ( w 1 ) − F ( w 1 )| = δ 2 . In addition, since the second derivative of F (w) is always bounded in [0, 1], we know that F (w) is smooth. Since F (w) ≡ ρ, we know that F (w) is ρ-Hessian Lipschitz. Therefore, F (w) and F (w) satisfy all of our assumptions.

C Proof of Theorem 2
Similar to Appendix A, we define ∆(w) = ∇ F (w) − ∇F (w). According to Assumptions 1 and 2, we know that both F (w) and F are µ-smooth. By the smoothness of F (w) and using the fact that η = 1/µ, we get By averaging (9) over t = 1, . . . , T , we get By taking square root on both sizes, and using Cauchy-Schwarz inequality as well as the fact that We then proceed to bound ∆(w t−1 ) 2 2 . According to Lemma 1, we have does not depend on the iteration count t. By Cauchy-Schwarz inequality, we get Then we bound w t−1 − w 0 2 2 . By trianlge inequality and Cauchy-Schwarz inequality, we obtain and therefore On the other hand, since F (w) is also µ-smooth, we have for every t ≥ 1, where in the equality we use the fact that η = 1/µ. This implies that By combining (12) and (13), we obtain and combining (11) and (14), we get By averaging (15) over t = 1, . . . , T and plugging the result in (10), we obtain which completes the proof.

D Proof of Theorem 3
Let F * = min w∈W F (w). Since we run gradient descent with learning rate η = 1/µ on the convex and µ-smooth function F , according to standard results in convex optimization [4], we have We provide the following lemma that bounds the difference between F (w) and F (w). Lemma 2. For any w ∈ W, we have We prove Lemma 2 in Appendix D.1. Here, we proceed to analyze F (w T ) − F * . We have To bound F (w T ) − F (w T ), we use the fact that for any convex and smooth functions, when we run gradient descent with learning rate 1/µ, the iterates only move closer to the minimum of the function, i.e., w T − w * 2 ≤ D. Then, by triangle inequality, for any k ∈ [K − 1] w T − w k 2 ≤ w T − w * 2 + w * − w k 2 ≤ D + w 0 − w k 2 . Using the fact that for any two positive numbers a and b, (a + b) 2 ≤ 2(a 2 + b 2 ) and (a + b) 3 ≤ 4(a 3 + b 3 ), we obtain w T − w k By combining (18) and (19) with Lemma 2, we obtain We can use a similar argument to bound F * − F * . Note that F * − F * ≤ F (w * ) − F (w * ). Therefore, according to Lemma 2, we have Then we can complete the proof by combining (17), (20), and (21).

D.1 Proof of Lemma 2
By definition of F (w) and F (w), we have for some ξ k ∈ [0, 1]. Then according to (8), we get We note that for any fixed vector v and function F (w), ∇(v ∇F (w)) = ∇ 2 F (w)v. Therefore, it is easy to compute Hessian-vector products in TensorFlow using its automatic differentiation function. Then, we are able to compute the rank-r approximations of Hessian matrices via power method. We present details in Algorithm 3. In our experiments, we use 5 power iterations for each eigenvector.

F Experiment details
In this section, we provide details for the neural network architecture of the CNN models that we use in this paper. We denote a convolutional layer with x output channels, a kernel of size y, and a stride of z by conv(x, y, z). We denote a fully connected layer with x output units as full(x). We denote a max pooling layer with kernel size x and stride y by pool(x, y). We denote ReLU activation by ReLU.