# A Comparison of Handcrafted and Deep Neural Network Feature Extraction for
Classifying Optical Coherence Tomography (OCT) Images

Kuntoro Adi Nugroho  Department of Computer Engineering
Universitas Diponegoro
Semarang, Indonesia
Email<EMAIL_ADDRESS>
###### Abstract

Optical Coherence Tomography allows ophthalmologist to obtain cross-section
imaging of eye retina. Assisted with digital image analysis methods, effective
disease detection could be performed. Various methods exist to extract feature
from OCT images. The proposed study aims to compare the effectiveness of
handcrafted and deep neural network features. The evaluated dataset consist of
32339 instances distributed in four classes, namely CNV, DME, DRUSEN, and
NORMAL. The methods are Histogram of Oriented Gradient (HOG), Local Binary
Pattern (LBP), DenseNet-169, and ResNet50. As a result, the deep neural
network based methods outperformed the handcrafted feature with 88% and 89%
accuracy for DenseNet and ResNet compared to 50 % and 42 % for HOG and LBP
respectively. The deep neural network based methods also demonstrated better
result on the under represented class.

Keywords:  Optical Coherence Tomography, HOG, LBP, DenseNet, ResNet

##  1 Introduction

As a breakthrough of medical imaging technology, with most application in
ophthalmology, optical coherence tomography (OCT) allows in situ morphology
assessment without invasive procedure (surgical or removal of tissue sample).
The OCT imaging technique is based on light, in contrast to ultrasound B-mode
imaging that utilizes sound, thus allows more accurate resolution of $10\mu m$
[1]. This method is superior to fluorescein angiography in imaging all layers
of the retinal vasculature as reported by Spaide et al [2].

Digital image processing and analysis tasks in optical coherence tomography
include automatic image segmentation and detection of disease. Both tasks are
related given that segmentation can be viewed as a recognition task in an
$mxn$ pixels image. Various method are available to extract feature from OCT
images. The features can be handcrafted or obtained from automatic feature
learning such as deep neural network.

Some example of studies on OCT image segmentation based handcrafted features
are the study by Vermeer et al. and Lang et al [3] [4]. Vermeer et al.
proposed two types of features. The first is pixel value and the pixel above
and below it. The second is haar-like feature based on A-line scan [3]. Two
general categories of features are introduced by Lang et al. namely spatial
aware and context aware features which constitutes 27 features in total [4].
Other study by Gadde et al. demonstrated foveal avascular zone (FAZ) produce
segmentation using local fractal analysis [5].

(a) CNV

(b) DME

(c) DRUSEN

(d) NORMAL

Figure 1: Dataset samples [6]

Handcrafted features work relatively well in specific domain with small number
of data. With continuous improvements, feature learning based on deep neural
network is more favorable due to possibility of learning from larger data and
to learn more complex pattern. This is achieved by constructing the network
with deeper layer. Learning from small number of sample is also possible with
transfer learning.

OCT segmentation with automatic feature extraction are demonstrated in the
study by He et al. and Alonso-Caneiro et al [7] [8]. U-Net based Segmentation
Network (S-Net) followed by Regression Network (R-Net) were proposed to obtain
topology guaranteed segmentation [7]. In contrast to the S-Net and R-Net
framework, Convolutional Neural Network (CNN) was trained to compute
probability map as a sign of the presence of boundary at certain pixels [8].
Both studies demonstrated the superiority of automatic feature extraction
compared to the non-automatic counterpart.

Coherent with the presented OCT segmentation studies, automatically and non-
automatically extracted features are also evaluated in OCT image
classification studies.

Atlas based shape analysis with Linear Discriminant Analysis (LDA) was studied
by Lee et al. The method was evaluated for two experiments, the first was
discerning glaucoma versus normal eye and the second was normal versus suspect
eye (normal eye in pair with unilateral glaucoma). The proposed method relies
on fshape framework to obtain atlas estimation, a template object from which
each observation based on. With reliance on initial template to compute the
final mean template, the method requires a good initial template to start
with. The study was also limited in terms of the number of evaluated data [9].
Other reported study in OCT is Histogram of Oriented Gradient (HOG) feature
extractor and Support Vector Machine for classifying normal, diabetic macular
edema (DME), and dry age-related macular degeneration (AMD). Evaluated on
spectral domain-OCT data from 45 patients, the correctly classified fraction
of volumes are $86.67\%,100\%,$ and $100\%$ respectively [10].

Traditional computer vision application in object recognition usually employs
manually extracted feature. With deep neural network feature extraction,
limitation on the use of good initial template and relatively small number of
data sample as required in the previously discussed studies could be
mitigated.

Convolutional neural network has been studied at least before 1995 with
application in hand written word recognition [11]. A notable success of this
method is in ImageNet object recognition [12].

The success has been followed up in medical image analysis including
application in fundus autofluorescence (FAF) and OCT images. In FAF, for
example, deep convolutional neural network (DCNN) performed $96\%$ and $91\%$
in recognizing geographic atrophy (GA) versus normal and GA vs other diseases
respectively. Reported study in OCT image includes the proposed VGG-16 with
Xavier algorithm weight initialization which successfully identify age-related
macular degeneration with ROC area under curve of $92.78\%$ and accuracy of
$87.63\%$ [13]. Other study is the identification and quantification of
intraretinal cystoid fluid (IRC) and subretinal fluid (SRF) that resulted in a
high correlation between the method and the ground truth [14].

In spite of the recent development in deep neural network based methods,
comparison and evaluation of deep neural network and non-automatic
(handcrafted) feature extraction in OCT has not been thoroughly reported.
Therefore, the proposed study aims to evaluate image feature extractions
methods for classifying OCT images. Two types of methods are considered. The
first method belongs to the non-automatic feature extraction, while the second
belongs to deep convolutional neural network family.

##  2 Research Method

###  2.1 Dataset

OCT images from the study by Kermany et al. are used in this study [6]. The
OCT images belong to four class, namely Choroidal Neovascularization (CNV),
Diabetic Macular Edema (DME), Drusen (presents in eary Age-Related Macular
Degeneration), and Normal. Some of the samples are shown in Fig 1. The data
distribution is shown in Table 1.

Table 1: Dataset Instance Distribution | train | test | val |
---|---|---|---|---
CNV | 37205 | 242 | 8 | 37455
DME | 11348 | 242 | 8 | 11598
DRUSEN | 8616 | 242 | 8 | 8866
NORMAL | 26315 | 242 | 8 | 26565
| 83484 | 968 | 32 | 84484

###  2.2 Methodology

The proposed methodology consists of the following steps

  1. 1.

Data preprocessing

  2. 2.

Feature extraction

  3. 3.

Classifier training

###  2.3 Data Preprocessing

  1. 1.

Image rescaling and padding

The original images are varied in size. In order to extract the feature, all
of the images were adjusted into uniform size of 224 x 224 pixels. The longer
dimension (between height or width) will be rescaled into 224, and the shorter
will be adjusted with maintaining the image aspect ratio. After rescaling, the
image was padded with zero pixels located in the center. The procedure is
illustrated in Fig 2.

(a)

(b)

Figure 2: Rescaling and padding process: A single channel 469 x 768 image (a)
was rescaled into 144 x 224 and then padded with black pixels to form 224 x
224 pixels image (b)

  2. 2.

Image Resampling

There is a large gap between the number of training samples (83484) and
validation samples (32). Given the relatively large number of instances and
the imbalance between training and validation data as shown in Table 1, the
training and validation data were resampled and re-distributed. Approximately
$\pm 25\%$ of the training images were randomly sampled to form the new
training set. Another $\pm 12.5\%$ different partition of random samples were
taken from training set for the new validation set. The test set was not
altered. Therefore, only 38.278 % of the data were retained. The new
distribution is show in Table 2.

Table 2: Dataset Instance Distribution After Resampling | train | test | val |
---|---|---|---|---
CNV | 9261 | 242 | 4743 | 14246
DME | 2803 | 242 | 1441 | 4486
DRUSEN | 2106 | 242 | 1096 | 3444
NORMAL | 6571 | 242 | 3350 | 10163
| 20741 | 968 | 10630 | 32339

###  2.4 Feature Extraction

The training, validation, and testing images were converted into feature
vector using the following methods

  1. 1.

Histogram of Oriented Gradient (HOG) [15]

  2. 2.

Local Binary Pattern (LBP) [16] [17]

  3. 3.

Pretrained Convolutional Neural Network: Residual Network (ResNet50) [18]

  4. 4.

Pretrained Convolutional Neural Network: Densely Connected Network
(DenseNet-169) [19]

(a) Image sample

(b) HOG feature

(c) LBP feature

(d) DenseNet-169 feature

(e) ResNet50 feature

Figure 3: Feature extraction example (image: normal-336530-5.png)

HOG and LBP belong to non-automatic feature extraction method, while the other
two convolutional neural networks (ResNet50 and DenseNet-169) belong to the
automatic method. The parameter settings for HOG and LBP are summarized in
Table 3.

Both of the convolutional neural network feature extractors are based on Keras
implementation [20]. The network parameter and size is described in Table 4.
The weights are obtained by pretraining on Imagenet data.

Table 3: HOG and LBP Parameter Settings | Parameter | Value
---|---|---
HOG | Number of orientations | 8
Pixels per cell | 16 x 16
Cells per block | 2 x 2
LBP | Number of Points | 16
Radius | 2
Histogram block size | 16 x 16
| Method | Uniform
Table 4: Network Size | Number of Parameters | Size on Disk
---|---|---
DenseNet-169 | 14,307,880 | 57 MB
ResNet50 | 25,636,712 | 99 MB

The resulting feature dimension for each method is summarized in Table 5.

Table 5: Feature Size HOG | LBP | DenseNet | ResNet
---|---|---|---
5408 | 1960 | 1664 | 2048

###  2.5 Classification

The extracted features were then utilized to train a classifier. A classifier
was trained for each feature extraction type. The evaluated classifier is a
perceptron neural network with no hidden layer (logistic regression). Although
multilayer perceptron is a possible option to learn more complex (non-linear)
pattern, linear classifier is sufficient for benchmarking.

The classifier was implemented using one-vs-all strategy. Therefore, four
perceptron units were trained for each classifier as there are four classes
(CNV, DME, DRUSEN, NORMAL). The classifier uses softmax activation function.

The training parameters were summarized in Table 6.

Table 6: Training Parameters Parameter | Value
---|---
Number of epochs | 100
Learning rate | $10^{-4}$
Optimizer | Adam [21]
Loss function | Categorical crossentropy

##  3 Result and Discussion

###  3.1 Result

####  3.1.1 Training and Validation

Fig 4 shows per-epoch training and validation accuracy for each method. The
validation accuracy of the four methods against each other is depicted in Fig
5.

(a) HOG

(b) LBP

(c) DenseNet-169

(d) ResNet50

Figure 4: Training and Validation Accuracy in 100 Epochs. The gray straight
line and the black dashed line denote training and validation accuracy
respectively. Figure 5: Validation accuracy for each method against each other

####  3.1.2 Evaluation in Test Data

Evaluated on test data, the classifier accuracy is shown in Table 11. The
precision, recall, and F1-score for HOG, LBP, DenseNet, and ResNet are
summarized in Table 7, 8, 9, and 10 respectively.

Table 7: Classifier with HOG Feature | Precision | Recall | F1-score
---|---|---|---
CNV | 0.40 | 0.96 | 0.57
DME | 0.60 | 0.04 | 0.07
DRUSEN | 0.59 | 0.04 | 0.08
NORMAL | 0.66 | 0.96 | 0.78
average | 0.56 | 0.50 | 0.37
Table 8: Classifier with LBP Feature | Precision | Recall | F1-score
---|---|---|---
CNV | 0.36 | 0.95 | 0.52
DME | 0.00 | 0.00 | 0.00
DRUSEN | 0.00 | 0.00 | 0.00
NORMAL | 0.55 | 0.74 | 0.63
average | 0.23 | 0.42 | 0.29
Table 9: Classifier with DenseNet-169 Feature | Precision | Recall | F1-score
---|---|---|---
CNV | 0.72 | 0.99 | 0.84
DME | 0.98 | 0.79 | 0.87
DRUSEN | 0.98 | 0.76 | 0.86
NORMAL | 0.94 | 0.98 | 0.96
average | 0.90 | 0.88 | 0.88
Table 10: Classifier with ResNet50 Feature | Precision | Recall | F1-score
---|---|---|---
CNV | 0.75 | 1.00 | 0.86
DME | 1.0 | 0.80 | 0.89
DRUSEN | 0.95 | 0.79 | 0.86
NORMAL | 0.94 | 0.98 | 0.96
average | 0.91 | 0.89 | 0.89
Figure 6: Recall Across The Four Methods Table 11: Test Accuracy HOG | LBP | DenseNet | ResNet
---|---|---|---
0.5010330 | 0.4235537 | 0.880165 | 0.8925619

###  3.2 Discussion

Fig 4 shows that the training and validation accuracy converges by the end of
the training process. The progress is almost static at the final epochs,
although LBP still shows relatively higher gradient compared to other methods.

HOG and LBP show almost no gap between training and validation accuracy, while
a larger gap is shown by the CNN based methods, especially ResNet50. The gap
is related to the bias-variance characteristic of the methods. The larger the
gap, the more the model variance. Against each other, the CNN based methods
consistently demonstrated higher accuracy than HOG and LBP, as shown in Fig 5.
As shown by the lower error (bias on training and validation data), the CNN
based method produced better feature compared to HOG and LBP.

The performance on the test data is consistent with the performance on the
training data. Overall, ResNet and DenseNet features outperformed HOG and LBP
in terms of accuracy, as shown in Table 11. The rank for each class recall is
almost consistent across the methods, as visualized in Fig 6. The recall on
CNV is the highest for all methods (0.96, 0.95, 0.99, 1.0 for each HOG, LBP,
DenseNet, and ResNet), followed by the recall on NORMAL. The result is very
likely due to the distribution of the data. CNV and NORMAL class are the
classes with the larger number of instances compared to the rest, as shown in
Table 1. Although DME and DRUSEN slightly under represented, DenseNet and
ResNet still performed much higher than other methods, greater than 75 % in
terms of recall.

In terms of classification performance, ResNet50 slightly outperformed
DenseNet-169. However, considering the feature and network parameters,
DenseNet-169 performed almost equally with ResNet50 with much less parameter
and network size ($\pm$ 14 M compared to $\pm$ 25 M, see Table 4).

##  4 Conclusion

In conclusion, the classifiers trained on features extracted from deep neural
network demonstrated the best performance. The result indicates that deep
neural network based methods generates better feature compared to HOG and LBP
in OCT image classification.

## References

  * [1] James G Fujimoto, Mark E Brezinski, Guillermo J Tearney, Stephen A Boppart, Brett Bouma, Michael R Hee, James F Southern, and Eric A Swanson.  Optical biopsy and imaging using optical coherence tomography.  1(9):970–972.
  * [2] Richard F. Spaide, James M. Klancnik, and Michael J. Cooney.  Retinal vascular layers imaged by fluorescein angiography and optical coherence tomography angiography.  133(1):45.
  * [3] K. A. Vermeer, J. van der Schoot, H. G. Lemij, and J. F. de Boer.  Automated segmentation by pixel classification of retinal layers in ophthalmic OCT images.  2(6):1743.
  * [4] Andrew Lang, Aaron Carass, Matthew Hauser, Elias S. Sotirchos, Peter A. Calabresi, Howard S. Ying, and Jerry L. Prince.  Retinal layer segmentation of macular OCT images using boundary classification.  4(7):1133.
  * [5] Santosh G. K. Gadde, Neha Anegondi, Devanshi Bhanushali, Lavanya Chidambara, Naresh Kumar Yadav, Aruj Khurana, and Abhijit Sinha Roy.  Quantification of vessel density in retinal optical coherence tomography angiography images using local fractal dimension.  57(1):246.
  * [6] Daniel S Kermany, Michael Goldbaum, Wenjia Cai, Carolina CS Valentim, Huiying Liang, Sally L Baxter, Alex McKeown, Ge Yang, Xiaokang Wu, Fangbing Yan, and others.  Identifying medical diagnoses and treatable diseases by image-based deep learning.  172(5):1122–1131.
  * [7] Yufan He, Aaron Carass, Bruno M. Jedynak, Sharon D. Solomon, Shiv Saidha, Peter A. Calabresi, and Jerry L. Prince.  Topology guaranteed segmentation of the human retina from OCT using convolutional neural networks.
  * [8] David Alonso-Caneiro, Scott A Read, Jared Hamwood, Stephen J Vincent, and Michael J Collins.  Automatic segmentation of retinal and choroidal thickness in OCT images using convolutional neural networks.  page 1.
  * [9] Sieun Lee, Nicolas Charon, Benjamin Charlier, Karteek Popuri, Evgeniy Lebed, Marinko V. Sarunic, Alain Trouvé, and Mirza Faisal Beg.  Atlas-based shape analysis and classification of retinal optical coherence tomography images using the functional shape (fshape) framework.  35:570–581.
  * [10] Pratul P. Srinivasan, Leo A. Kim, Priyatham S. Mettu, Scott W. Cousins, Grant M. Comer, Joseph A. Izatt, and Sina Farsiu.  Fully automated detection of diabetic macular edema and dry age-related macular degeneration from optical coherence tomography images.  5(10):3568.
  * [11] Yoshua Bengio, Yann LeCun, and Donnie Henderson.  Globally trained handwritten word recognizer using spatial representation, convolutional neural networks, and hidden markov models.  page 8.
  * [12] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton.  ImageNet classification with deep convolutional neural networks.  60(6):84–90.
  * [13] Cecilia S. Lee, Doug M. Baughman, and Aaron Y. Lee.  Deep learning is effective for classifying normal versus age-related macular degeneration OCT images.  1(4):322–327.
  * [14] Thomas Schlegl, Sebastian M. Waldstein, Hrvoje Bogunovic, Franz Endstraßer, Amir Sadeghipour, Ana-Maria Philip, Dominika Podkowinski, Bianca S. Gerendas, Georg Langs, and Ursula Schmidt-Erfurth.  Fully automated detection and quantification of macular fluid in OCT using deep learning.  125(4):549–558.
  * [15] N. Dalal and B. Triggs.  Histograms of oriented gradients for human detection.  In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 1, pages 886–893. IEEE.
  * [16] T. Ojala, M. Pietikainen, and D. Harwood.  Performance evaluation of texture measures with classification based on kullback discrimination of distributions.  In Proceedings of 12th International Conference on Pattern Recognition, volume 1, pages 582–585. IEEE Comput. Soc. Press.
  * [17] T. Ojala, M. Pietikainen, and T. Maenpaa.  Multiresolution gray-scale and rotation invariant texture classification with local binary patterns.  24(7):971–987.
  * [18] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778. IEEE.
  * [19] Gao Huang, Zhuang Liu, Laurens van der Maaten, and Kilian Q. Weinberger.  Densely connected convolutional networks.
  * [20] François Chollet and others.  Keras.
  * [21] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.


# Analysis of Catastrophic Forgetting for Random Orthogonal Transformation
Tasks in the Overparameterized Regime

Daniel Goldfarb
Khoury College of Computer Sciences
Northeastern University
Boston, MA 02115
<EMAIL_ADDRESS>
&Paul Hand
Dept. of Mathematics and Khoury College of Computer Sciences
Northeastern University
Boston, MA 02115
<EMAIL_ADDRESS>

###### Abstract

Overparameterization is known to permit strong generalization performance in
neural networks. In this work, we provide an initial theoretical analysis of
its effect on catastrophic forgetting in a continual learning setup. We show
experimentally that in permuted MNIST image classification tasks, the
generalization performance of multilayer perceptrons trained by vanilla
stochastic gradient descent can be improved by overparameterization, and the
extent of the performance increase achieved by overparameterization is
comparable to that of state-of-the-art continual learning algorithms. We
provide a theoretical explanation of this effect by studying a qualitatively
similar two-task linear regression problem, where each task is related by a
random orthogonal transformation. We show that when a model is trained on the
two tasks in sequence without any additional regularization, the risk gain on
the first task is small if the model is sufficiently overparameterized.

##  1 Introduction

Continual learning is the ability of a model to learn continuously from a
stream of data, building on what was previously learned and retaining
previously learned skills without the need for retraining. A major obstacle
for neural networks to learn continually is the catastrophic forgetting
problem: the abrupt drop in performance on previous tasks upon learning new
ones. Modern neural networks are typically trained to greedily minimize a loss
objective on a training set, and without any regularization, the model’s
performance on a previously trained task may degrade. Techniques for
mitigating catastrophic forgetting fall under three main groups: generative
replay, parameter isolation, and regularization methods [7]. Generally, the
goal of regularization methods is to determine important parameters from
previous tasks and constrain them so that they do not get modified too much
while training subsequent tasks. Two common regularization methods are
Synaptic Intelligence (SI) [22] and Elastic Weight Consolidation (EWC) [13].
See Appendix A for a detailed description of them.

It is well-known that strong generalization performance for neural networks is
typically obtained in the overparameterized regime, where the number of
learnable parameters is greater than the number of training examples. Work on
overparameterized machine learning has led to research on the so-called double
descent phenomenon, where test error improves as model complexity increases
beyond the level needed to fit the training data, outperforming all
underparameterized versions of the model [2]. One of the first observations of
this behavior in modern neural networks was in extremely wide ResNet18 models
that generalize better than their underparameterized counterparts on CIFAR-10
despite fitting to label noise [17]. This model-wise double descent phenomenon
has been demonstrated analytically in a variety of machine learning models
[10, 3, 1], including some as simple as linear regression. Such linear models
will be the basis of theoretical analysis in the present paper.

While investigating the relationship of overparameterization with catastrophic
forgetting, we made an interesting experimental observation. We considered 10
permuted-MNIST tasks, where each task has training data given by random
permutations of the original MNIST images [15]. We trained 2-layer multilayer
perceptrons (MLPs) with a variety of layer widths given by $[400w,400w]$,
where $w=1,3,5,7,9$, using vanilla stochastic gradient descent and the
continual learning algorithms, SI and EWC. We compared the test accuracy on
all seen tasks. As expected, and as shown in Figure 1, average accuracy with
SGD drops significantly after learning multiple tasks, and that drop is
mitigated by using SI or EWC. Interestingly, we observe that a significant
fraction of the accuracy gain achieved by SI or EWC can be obtained with
vanilla SGD by simply overparameterizing the model. This can be seen by
comparing the $w=1$ and $w=9$ curves with SGD to the curves with SI and EWC.
See Appendix B for more details on the experiments.

Figure 1: Results of permuted MNIST experiment. Red curves denote performance
of SI, yellow curves denote performance of EWC, blue curves denote performance
of Vanilla SGD. Bolder saturation of lines corresponds to larger width
parameters (more overparameterization). Specific hyperparameters are reported
in Appendix B. Curves for $w=7,9$ for EWC are omitted due to computational
constraints on Fisher matrix estimates.

The goal of the present paper is to analytically illustrate the effect
overparameterization can have on catastrophic forgetting. As with the
illustrations of double descent in [10, 3, 1], we choose to study the effect
with a linear regression problem for simplicity and mathematical convenience.
While the experiments above study the case of tasks related by random
permutations, our analysis will instead study the qualitatively similar case
where two tasks are related by a random orthogonal transformation, which
results in simpler mathematical analysis, as we discuss in Section 3. The
relation given by the random orthogonal transformation gives two data feature
spaces corresponding to each task that are approximately but not exactly
orthogonal. We construct two tasks, A and B. Let task A be defined by data
matrix $X_{A}\in\mathbb{R}^{n\times p}$, with rows being $p$ noisy random
projections of some low-dimensional latent features, and responses
$y\in\mathbb{R}^{n}$ that are noiseless and linear in the latent features. Let
$O$ be a random $p\times p$ orthogonal matrix. Then task B is defined by data
$X_{B}=X_{A}O^{\top}$ and the same responses $y$. Learning tasks A and B
involve estimating a $\beta\in\mathbb{R}^{p}$ for the predictor
$f_{\beta}:x\mapsto x^{\top}\beta$ to fit the data $(X_{A},y)$ and
$(X_{B},y)$, respectively.

We analyze the increase in statistical risk on task A between an estimator
trained on task $A$ by minimizing square loss, with initialization at zero,
and one sequentially trained on task $A$ and then task $B$ with no explicit
regularization. Let $R(f_{\beta})$ be the risk on task A of an estimator $f$
with parameters $\beta$. Let $\hat{\beta}_{A}$ be the parameters of the model
that is trained on task A. Let $\hat{\beta}_{BA}$ be the parameters of a model
that is initialized at $\hat{\beta}_{A}$ and then trained on task B. Our main
result is that if there are more training examples than the intrinsic (latent)
dimensionality of the data and if there is not too much noise in the observed
features, then

$\displaystyle
R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})\lesssim\sqrt{\frac{n}{p}}$ (1)

with high probability. The result asserts that under our linear model, the
extent of catastrophic forgetting is arbitrarily small if the
overparameterization ratio, $p/n$, is sufficiently large. We thus see an
analytical illustration that catastrophic forgetting can be ameliorated by
overparametization in the case of a suitable linear model. The full theorem is
stated in Section 2.4 and its proof is provided in Appendix E.

The contributions of this paper are:

  * •

We empirically observe that overparameterization can account for a majority of
the performance drop due to catastrophic forgetting in a permuted image task
using a multi-layer perceptron.

  * •

We provide a linear regression problem that exhibits a corresponding effect
for overparameterization and continual learning.

  * •

We establish a non-asymptotic bound on the performance drop of this linear
model in an orthogonal transformation task setting using results from random
matrix theory. This result provides a formal illustration that continual
learning can in some cases be ameliorated by overparameterization.

##  2 Analysis of Catastrophic Forgetting in a Linear Model

In this section, we present a latent space model for linear regression that we
will analyze in order to illustrate that overparameterization can ameliorate
catastrophic forgetting. Our single task model is the latent space model of
[10] without label noise. Then, we present the analogy between this linear
model and neural networks. Next, we empirically demonstrate that under this
model, overparameterization ameliorates catastrophic forgetting. Finally, we
present a theorem that establishes that observation with high probability.

###  2.1 Latent Space Models for Two Linear Regression tasks

Let $\mathcal{Z}=\mathbb{R}^{d}$, which we call the latent feature space.
Consider data for regression generated by a noiseless linear response to
standard Gaussian latent features. That is, for some
$\theta\in\mathbb{R}^{d}$, let an example be given by

$\displaystyle z$ $\displaystyle\sim\mathcal{N}(0,I_{d}),$ (2) $\displaystyle
y$ $\displaystyle=z^{\top}\theta.$ (3)

Let $\mathcal{X}=\mathbb{R}^{p}$, which we call the observed feature space. We
consider the case where, for each example, we have access only to $p$ observed
features, given by noisy random projections of the latent features:

$\displaystyle x$ $\displaystyle=Wz+u$ (4)

where $W\in\mathbb{R}^{p\times d}$ and $u\sim\mathcal{N}(0,I_{p})$. We could
take $W$ to have i.i.d. $\mathcal{N}(0,\gamma)$ entries, but for mathematical
convenience, we will instead study the idealization in which $W$ has columns
that form a scaled orthonormal basis of a random $d$-dimensional subspace of
$\mathbb{R}^{p}$. Namely, $W^{\top}W=p\gamma I_{d}$. For large $p$, this
idealization is approximately satisfied under the above Gaussian model for
$W$.

We consider two tasks, $A$ and $B$, both with $n$ examples. Task $A$ has data
$(X_{A},y)\in\mathbb{R}^{n\times p}\times\mathbb{R}^{n}$ where each of the $n$
rows of $X_{A}$ and entries of $y$ are sampled independently by (2) - (4). Let
$O$ be a random $p\times p$ orthogonal matrix. Task B has data $(X_{B},y)$
where $X_{B}=X_{A}O^{\top}$.

We study estimators that are linear in the observed features $x$:

$\displaystyle f_{\hat{\beta}}:x\mapsto x^{\top}\hat{\beta},$ (5)

and we will sometimes refer to the parameters $\hat{\beta}$ as the estimator.
We estimate the parameters of this model by gradient descent with a square
loss. We are interested in the case of $d<n<p$. As $n<p$, the solution to this
problem depends on initialization and solves the following optimization
problem:

$\displaystyle\arg\min_{\hat{\beta}}\frac{1}{2}\|\hat{\beta}-\beta_{0}\|^{2}\text{
s.t. }y=X\hat{\beta},$ (6)

where $\beta_{0}$ is the initialization, and $X$ is either $X_{A}$ or $X_{B}$,
depending on the task being solved. To study the sequential training of tasks
A and B, we define the following estimators:

  * •

$\hat{\beta}_{A}$ is the solution to task $A$ when initialized at $0$,

  * •

$\hat{\beta}_{B}$ is the solution to task $B$ when initialized at $0$,

  * •

$\hat{\beta}_{BA}$ is the solution to task $B$ when initialized at
$\hat{\beta}_{A}$.

These parameters are found by solving the following optimization problems:

$\displaystyle\hat{\beta}_{A}$
$\displaystyle=\arg\min_{\hat{\beta}}\frac{1}{2}\|\hat{\beta}\|^{2}\text{ s.t.
}y=X_{A}\hat{\beta},$ (7) $\displaystyle\hat{\beta}_{B}$
$\displaystyle=\arg\min_{\hat{\beta}}\frac{1}{2}\|\hat{\beta}\|^{2}\text{ s.t.
}y=X_{B}\hat{\beta},$ (8) $\displaystyle\hat{\beta}_{BA}$
$\displaystyle=\arg\min_{\hat{\beta}}\frac{1}{2}\|\hat{\beta}-\hat{\beta}_{A}\|^{2}\text{
s.t. }y=X_{B}\hat{\beta}.$ (9)

The optimization problem in (6) has the following closed form solution when
$X$ has rank $n$:

$\displaystyle\hat{\beta}$
$\displaystyle=\beta_{0}+X^{\top}(XX^{\top})^{-1}y-X^{\top}(XX^{\top})^{-1}X\beta_{0}$
(10)
$\displaystyle=\beta_{0}+X^{\top}(XX^{\top})^{-1}y-\mathcal{P}_{X^{\top}}\beta_{0},$
(11)

where $\mathcal{P}_{X^{\top}}$ is the orthogonal projector onto the range of
$X^{\top}$. As $n<p$, $X_{A}$ and $X_{B}$ have rank $n$ with probability $1$,
and this gives the following closed forms for
$\hat{\beta}_{A},\hat{\beta}_{B},\hat{\beta}_{BA}$:

$\displaystyle\hat{\beta}_{A}$
$\displaystyle=X_{A}^{\top}(X_{A}X_{A}^{\top})^{-1}y,$ (12)
$\displaystyle\hat{\beta}_{B}$
$\displaystyle=X_{B}^{\top}(X_{B}X_{B}^{\top})^{-1}y,$ (13)
$\displaystyle\hat{\beta}_{BA}$
$\displaystyle=\hat{\beta}_{A}+\hat{\beta}_{B}-\mathcal{P}_{X_{B}^{\top}}\hat{\beta}_{A}.$
(14)

We evaluate these estimators on task A. The risk on task A of an estimator $f$
with parameters $\hat{\beta}$ is given by

$\displaystyle
R(f_{\hat{\beta}})=\sigma^{2}+(\hat{\beta}-\beta)^{\top}\Sigma(\hat{\beta}-\beta)$
(15)

where

$\displaystyle\Sigma$ $\displaystyle=WW^{\top}+I_{p}$ (16)
$\displaystyle\beta$ $\displaystyle=(I+WW^{\top})^{-1}W\theta$ (17)
$\displaystyle\sigma^{2}$
$\displaystyle=\theta^{\top}(W^{\top}W+I_{d})^{-1}\theta$ (18)

See Appendix D for the derivation of (15)–(18). It follows from showing that
the latent space model described above is equivalent to an anisotropic
regression model where $X_{A}$ has i.i.d. rows
${X_{A}}_{i}\sim\mathcal{N}(0,\Sigma)$ and labels $y=X_{A}\beta+\epsilon$
where $\epsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$.

Figure 2: The solid black line depicts the span of $W$. The true parameters
corresponding to tasks A and B are given by $\beta\in W$ and $O\beta$. The
gray lines depict the set of solutions to $X_{A}\beta=y$ and $X_{B}\beta=y$.
The estimators $\hat{\beta}_{A},\hat{\beta}_{B},\hat{\beta}_{BA}$ are given by
orthogonal projections of an initialization on the respective consistent
solutions. The red and blue ellipses depict lines of constant risk for tasks A
and B, respectively.

We aim to bound $R(f_{\hat{\beta}_{BA}})$ relative to
$R(f_{\hat{\beta}_{A}})$. Figure 2 illustrates the estimators
$\hat{\beta}_{A},\hat{\beta}_{B},\hat{\beta}_{BA}$ and curves of constant
risk. As depicted, if $p$ is large enough, $\hat{\beta}_{BA}$ has low risk on
Task A and Task B simultaneously.

###  2.2 Analogy of Linear Model to Neural Networks

The linear model we study is intended to be a mathematically tractable
idealization of a neural network, and it is meant to analytically illustrate
that overparameteriation can ameliorate catastrophic forgetting. The analogy
of this linear model and neural network training on image data is as follows:

Natural images in a neural network’s training distribution can be
(approximately) modeled as being on a nonlinear manifold and having a low-
dimensional latent representation. Instead of observing the latent
representation of an image, the neural network only sees a high-dimensional
representation either directly in pixel space or perhaps in a representation
computed from pixel space. Either of these representations contain noise in
the features used for prediction. Responses can be approximated by a neural
network.

In our linear model, the low-dimensional representation of an input image is
in a $d$-dimensional latent feature linear space. The responses are linear in
the latent features. We assume the response is noiseless for the sake of
simplicity, though our results could be extended to the noisy case. In our
linear model, predictions are made off of a $p$-dimensional model given by
noisy random projections of the latent features. We constrain $W$ to have
orthonormal columns which is a mathematical idealization of Gaussian
measurements. We study two tasks with the same responses like in the
permutation task setup, but for mathematical convenience we study tasks that
are related by a random orthogonal transformation instead.

###  2.3 Numerical Experiment

Before we establish our theoretical result about the system described in
Section 2.1, we provide empirical evidence that the latent space linear
regression model above exhibits the phenomenon that overparameterization can
ameliorate catastrophic forgetting. Specifically, we provide empirical
evidence that $R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$ decreases with
$p$.

Let $d=20$, $n=100,\gamma=1$, $\beta_{0}=\vec{0}$, and
$\theta\sim\mathcal{N}(0,I_{d})$. We plot $R(f_{\hat{\beta}_{BA}})$,
$R(f_{\hat{\beta}_{A}})$, $R(f_{\beta_{0}})$ as a function of $p\in(n,2000)$
averaged over 100 samplings of $W,X_{A},O,u$.

Figure 3: Risk as a function of model complexity $p$ for model (2)–(14). The
dashed line depics the null risk, corresponding to the zero estimator. Note
the log scale on the vertical axis. Result of simulated numerical experiment
for random orthogonal transformation tasks. Dotted black line denotes risk of
null estimator, blue line denotes risk of estimator trained on task A, orange
line denotes risk of estimator trained on task A then task B.

Figure 3 shows the results of the experiment. We first note that both
$\hat{\beta}_{A}$ and $\hat{\beta}_{BA}$ outperform the null risk, given by
$\beta=0$. The null risk, $R(f_{\beta_{0}})$, defines a baseline that any
reasonable model must beat. We also observe that $R(f_{\hat{\beta}_{A}})$ and
$R(f_{\hat{\beta}_{BA}})$ are decreasing with $p$ in the overparameterized
regime, and that the difference between these risks appears to decrease for
increasing $p$. Note the log-scale of the vertical-axis. This provides
evidence that catastrophic forgetting is alleviated in the overparameterized
regime in our two-task learning setup.

###  2.4 Main Result

Our main result is an upper bound on the performance drop, defined as
$R(f_{\hat{\beta}_{A}})-R(f_{\hat{\beta}_{BA}})$, for the two-task latent
space linear regression model described above and inspired by the double
descent literature. As described in Section 2.1, we consider (2)–(14), where
$W$ satisfies the following assumption.

######  Assumption 2.1.

All non-zero singular values of $W$ are equal. Namely, $W^{\top}W=p\gamma
I_{d}$.

We begin with a proposition that defines the unlearned baseline for the
problem.

######  Proposition 2.2.

Fix $\theta\in\mathbb{R}^{d}$. Let $W\in\mathbb{R}^{p\times d}$ satisfy
Assumption 2.1. Then

$\displaystyle R(f_{0})=\|\theta\|^{2}.$ (19)

This risk calculation agrees with the numerical experiment in Section 2.3
where $\|\theta\|^{2}\approx d=20$. The result is formally stated and proven
in Lemma E.3. For our main result, we prove that if the number of examples
exceeds the problem’s latent dimensionality, if the number of parameters is
sufficiently large relative to the number of examples and relative to the
noise level of the observable features, then with high probability, the
performance drop is small.

######  Theorem 2.3.

Fix $\theta\in\mathbb{R}^{d}$. Let tasks $A,B$ be given by (2)–(14). Let
$W\in\mathbb{R}^{p\times d}$ satisfy Assumption 2.1 and $n\geq d,p\geq
max(17n,1/\gamma)$. Then there exists constant $c>0$ such that with
probability at least $1-10e^{-cd}$, the following holds:

$\displaystyle
R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})\leq\left(66\sqrt{\frac{n}{p}}+\frac{12}{p\gamma}\right)\|\theta\|^{2}$
(20)

Theorem 2.3 provides an upper bound on the amount of risk gained on task A
after subsequential training on task B given by two terms. The first term
provides dependence on the overparameterization ratio $p/n$ and decreases as
overparameterization becomes more extreme. The second term is given by the
signal to noise ratio of the noisy features. This term dominates only when
$\gamma\ll 1/\sqrt{np}$. Based on the theorem, we observe that the
overparameterization needs only to be linear in $n$ to achieve a negligible
performance drop in unregularized sequential task training compared to the
baseline of $\|\theta\|^{2}$ in Proposition 2.2. This shows that catastrophic
forgetting is ameliorated in the overparameterized regime. This result is
formally stated in Lemma E.11. A formal proof and supporting lemmas are
supplied in Appendix E. We provide a proof sketch here to outline the
techniques used.

###  2.5 Proof Sketch of Theorem 2.3

For readability, we write $X_{A}$ as $A$ and $X_{B}$ as $B$. As shown in
Appendix D, the latent space model described above is equivalent to an
anisotropic regression model where $A$ has i.i.d. rows
$A_{i}\sim\mathcal{N}(0,\Sigma)$ and labels $y=A\beta+\epsilon$ where
$\epsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$.

Performance drop is given by

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle=(\hat{\beta}_{BA}-\beta)^{\top}\Sigma(\hat{\beta}_{BA}-\beta)-(\hat{\beta}_{A}-\beta)^{\top}\Sigma(\hat{\beta}_{A}-\beta).$

After substituting the closed form solutions for
$\hat{\beta}_{A},\hat{\beta}_{BA}$, distributing terms, and applying simple
Cauchy-Schwarz and triangle inequalities, we get the following bound:

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq
8\frac{p\gamma\sqrt{p\gamma}}{p\gamma+1}\|\theta\|\|\mathcal{P}_{W}\hat{\beta}_{B}\|+14\sqrt{p\gamma}\|\theta\|\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|+12\frac{p\gamma}{(p\gamma+1)^{2}}\|\theta\|^{2}$
(21) $\displaystyle=I+II+III$ (22)

Lemmas E.5, E.6 establish results for orthogonal transformations to help bound
$\|\mathcal{P}_{W}\hat{\beta}_{B}\|,\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|$
respectively. $\mathcal{P}_{W}$ is a projection onto a $d$-dimensional space
which scales the norm in $I$ by $d/p$. $\mathcal{P}_{B^{\top}}$ is a
projection onto an $n$-dimensional space which scales the norm in $II$ by
$n/p$. As $d\leq n$ by assumption, $II$ dominates $I$ in the final bound.
Using these results and simplifying gives the following bound with probability
at least $1-10e^{-cd}$ for constant $c>0$:

$\displaystyle I+II\leq 66\sqrt{\frac{n}{p}}\|\theta\|^{2}$ (23)

We directly obtain

$\displaystyle III\leq\frac{12}{p\gamma}\|\theta\|^{2}$ (24)

Combining these bounds completes the proof.

##  3 Discussion

Overparameterization is a necessity for continual learning so that there can
exist an infinity of potential optima for each task [13]. This makes it likely
that there exists an optimum for some task B that is close to the solutions of
some task A. We provide experimental evidence that overparameterization can
provide additional benefits in combatting catastrophic forgetting for neural
networks solving permutation tasks. We use a linear model with clear analogies
to neural networks in order to study this behavior theoretically. In our
analysis of the linear model in the overparameterized regime, non-asymptotic
matrix estimates and results for orthogonal transformations provide bounds on
the performance drop. Our main result shows that, under our model,
catastrophic forgetting is ameliorated for sufficiently large
overparameterization. For the linear setting we study, the behavior we observe
can be explained geometrically: overparameterization causes the random
orthogonal transformation tasks to live in approximately orthogonal subspaces,
so training on subsequent tasks does not interrupt performance on learned
tasks.

We view the present work as helping to establish initial results for continual
learning theory. Before the field can rigorously understand machine learning
algorithms in practice, the behavior of simple systems should be well
understood. In particular, the behavior of linear systems with only vanilla
SGD is the most natural initial result. Our work remarks that future theory
should establish that continual learning algorithms beat not only a moderately
parameterized baseline, but also the performance of extremely
overparameterized models.

First we address the concern for using permutation tasks as realistic
benchmarks for continual learning methods. Researchers believe that
permutation tasks only provide a best-case for real world scenarios [9]. Also,
on a number of image classification datasets, MLPs do not experience
forgetting when only two permutation tasks are being learned [18]. Our
experiments confirm this effect while also showing that overparameterization
mitigates the observable forgetting on 10 task permuted MNIST. Despite these
critiques, we use permutation tasks as a launching point for theory because
each task is of the same ‘difficulty’ and is amenable to mathematical
analysis.

Next we discuss our choice to study the problem with a linear model. Linear
regression is the simplest setting, for which we know, that exhibits double
descent. The consensus of several works that study double descent in linear
models is that the risk of a model is monotonically decreasing in the
overparameterized regime with respect to number of parameters only if the data
has low effective dimension and high ambient dimension compared to the number
of training samples [6, 1, 10]. In order to have a model that has
monotonically decreasing performance drop for a particular continual learning
problem, it is a necessity that it exhibits monotonically decreasing risk on a
single task. Additionally in recent work, connections have been made between
neural networks and linear models using the so-called neural tangent kernel
(NTK) phenomenon [12]. The parameterization of a neural network can be so
large that training only changes its parameters slightly from its
initialization, resulting in functions that can be accurately approximated
linearly. Hence it is reasonable that the analysis of linear models can
explain the behavior of neural networks.

We now remark at a technical level two choices in our analysis. The first is
why we studied the case of random orthogonal transformation tasks instead of
permutation tasks. The empirical performance between orthogonal and
permutation tasks is similar; they both create tasks that are equally
’difficult’ for an MLP to learn, which spares us from needing to quantify
problem difficulty. Appendix C provides evidence that permutation and
orthogonal transformation tasks have the same difficulty in the linear
setting. Also, the mathematical analysis is easier when studying orthogonal
transformation tasks. With random orthogonal transformations, any subspace
gets mapped to a random subspace, for which the values of coefficients are
typically well spread out. With random permutations, some subspaces (e.g.
those aligned with the standard basis elements) do not exhibit the same
spreading effect, making the technical analysis more involved. Secondly, we do
not present a bound on $R(f_{\hat{\beta}_{A}})$, though it is expected to
approach zero for large $p$, as suggested by Figure 3. Whether or not this
risk goes to 0 in $p$, the performance drop goes to 0 in $p$ while the null
estimator remains with constant risk. So the regression problem is being
solved arbitrarily well for sufficiently large $p$.

With the growing popularity of continual learning, much of recent work is
focused on developing new algorithms to mitigate catastrophic forgetting [13,
22, 20, 16]. Only a few papers study the problem theoretically [14, 4, 8, 11,
5]. [14] uses set theory to prove that, in general, continual learning
problems are NP-hard, explaining why generative replay methods perform so
well. [4] uses the NTK regime to prove generalisation guarantees for an
existing continual learning method. [8] uses an NTK overlap matrix to define a
notion of task similarity and show that catastrophic forgetting is more severe
when tasks have high similarity. [11] studies a family of continual learning
methods that uses approximations of the Hessian to determine parameter
importance, presenting scenarios where continual learning provably fails and
succeeds. [5] shows that a number of regularization techniques that seem to be
derived from differing philosophies actually all study a variation of the
Fisher information matrix. While prior work has presented generalization
bounds for existing continual learning techniques, our work illustrates the
relationship between overparameterization and generalization.

A natural next step is to study the regimes in which catastrophic forgetting
is most problematic. This includes the setting where tasks do not have a
nearly orthogonal relationship but also when data does not necessarily live on
a low-dimensional manifold. We are also interested in understanding how the
ideas of this paper generalize to other continual learning benchmarks and for
more general neural network architectures. Prior work found experimental
evidence that catastrophic forgetting is most severe not when tasks are very
dissimilar but when they only have an intermediate level of similarity [19].
Using orthogonality as a proxy for task similarity, this agrees with our work
that shows that nearly orthogonal tasks are less prone to catastrophic
forgetting. An interesting future work would be to formalize this notion of
task similarity for our model. Moving forward, one goal of theory in continual
learning is to be able to analytically compare algorithms. Our work provides a
foundation of understanding this behavior in a simple linear regression
setting. In order to push this work forward, either non-linear models need to
be studied or tasks that are related by something more complex than
permutations.

## References

  * [1] Peter L Bartlett, Philip M Long, Gábor Lugosi, and Alexander Tsigler.  Benign overfitting in linear regression.  Proceedings of the National Academy of Sciences, 117(48):30063–30070, 2020.
  * [2] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
  * [3] Mikhail Belkin, Daniel Hsu, and Ji Xu.  Two models of double descent for weak features.  SIAM Journal on Mathematics of Data Science, 2(4):1167–1180, 2020\.
  * [4] Mehdi Abbana Bennani, Thang Doan, and Masashi Sugiyama.  Generalisation guarantees for continual learning with orthogonal gradient descent.  arXiv preprint arXiv:2006.11942, 2020.
  * [5] Frederik Benzing.  Unifying importance based regularisation methods for continual learning.  In International Conference on Artificial Intelligence and Statistics, pages 2372–2396. PMLR, 2022.
  * [6] Yehuda Dar, Vidya Muthukumar, and Richard G Baraniuk.  A farewell to the bias-variance tradeoff? an overview of the theory of overparameterized machine learning.  arXiv preprint arXiv:2109.02355, 2021.
  * [7] Matthias Delange, Rahaf Aljundi, Marc Masana, Sarah Parisot, Xu Jia, Ales Leonardis, Greg Slabaugh, and Tinne Tuytelaars.  A continual learning survey: Defying forgetting in classification tasks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [8] Thang Doan, Mehdi Abbana Bennani, Bogdan Mazoure, Guillaume Rabusseau, and Pierre Alquier.  A theoretical analysis of catastrophic forgetting through the ntk overlap matrix.  In International Conference on Artificial Intelligence and Statistics, pages 1072–1080. PMLR, 2021.
  * [9] Sebastian Farquhar and Yarin Gal.  Towards robust evaluations of continual learning.  arXiv preprint arXiv:1805.09733, 2018.
  * [10] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  arXiv preprint arXiv:1903.08560, 2019.
  * [11] Reinhard Heckel.  Provable continual learning via sketched jacobian approximations.  In International Conference on Artificial Intelligence and Statistics, pages 10448–10470. PMLR, 2022.
  * [12] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  arXiv preprint arXiv:1806.07572, 2018.
  * [13] James Kirkpatrick, Razvan Pascanu, Neil Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska-Barwinska, et al.  Overcoming catastrophic forgetting in neural networks.  Proceedings of the national academy of sciences, 114(13):3521–3526, 2017.
  * [14] Jeremias Knoblauch, Hisham Husain, and Tom Diethe.  Optimal continual learning has perfect memory and is np-hard.  In International Conference on Machine Learning, pages 5327–5337. PMLR, 2020.
  * [15] Yann LeCun.  The mnist database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.
  * [16] Zhizhong Li and Derek Hoiem.  Learning without forgetting.  IEEE transactions on pattern analysis and machine intelligence, 40(12):2935–2947, 2017.
  * [17] Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever.  Deep double descent: Where bigger models and more data hurt.  Journal of Statistical Mechanics: Theory and Experiment, 2021(12):124003, 2021.
  * [18] Benedikt Pfülb and Alexander Gepperth.  A comprehensive, application-oriented study of catastrophic forgetting in dnns.  arXiv preprint arXiv:1905.08101, 2019.
  * [19] Vinay V Ramasesh, Ethan Dyer, and Maithra Raghu.  Anatomy of catastrophic forgetting: Hidden representations and task semantics.  arXiv preprint arXiv:2007.07400, 2020.
  * [20] Hanul Shin, Jung Kwon Lee, Jaehong Kim, and Jiwon Kim.  Continual learning with deep generative replay.  Advances in neural information processing systems, 30, 2017.
  * [21] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027, 2010.
  * [22] Friedemann Zenke, Ben Poole, and Surya Ganguli.  Continual learning through synaptic intelligence.  In International Conference on Machine Learning, pages 3987–3995. PMLR, 2017.

##  Appendix A Description of Continual Learning Techniques

Synaptic Intelligence (SI) is a regularization technique that assigns to each
parameter of the network an estimate of importance for learned tasks [22].
This weight is determined in an online manner by tracking the amount that each
parameter contributed to the decrease in loss during training. The weight is
then used to penalize changes to the network parameters during subsequent
training in the form of a regularization term added to the loss function.

Elastic Weight Consolidation (EWC) is a regularization technique that
determines the importance of network weights using an estimation of the Fisher
Information Matrix [13]. Near a minimum of the loss function, the diagonals of
the Fisher matrix act as an estimate of the second order derivative of the
loss with respect to each parameter. The magnitude of this derivative is used
as a proxy for how sensitive the loss function is to fluctuation of the
parameter. Constraining parameters according to their corresponding Fisher
diagonal entries shows as an effective way of retaining the values of
important weights from previous tasks while training on new ones.

##  Appendix B MNIST Experiments

Table 1: Hyperparameters for the MNIST experiments Hyperparameter | SI | EWC | SGD
---|---|---|---
learning rate | 0.1 | 0.01 | 0.1
dropout | 0 | 0 | 0.5
batch size | 64 | 128 | 64
epochs / dataset | 5 | 20 | 5
c | 0.1 |  |
$\xi$ | 0.1 |  |
$\lambda$ |  | 150 |
fisher sample size |  | 1,000 |

Table 1 reports the hyperparameters used in the MNIST experiments. We adopted
the same hyperparameters for SI as in the original paper [22]. To our
surprise, EWC with default hyperparameters [13] did not compete with SI. A
basic grid search gave us a model that was more competitive. Blank entries
mean that the hyperparameter is not relevant for the particular method. Curves
for $w=7,9$ are omitted due to computational constraints in computing Fisher
matrix estimates.

##  Appendix C Permuted Numerical Experiment

Figure 4: Result of simulated numerical experiment for random permutation
tasks. Dotted black line denotes risk of null estimator, blue line denotes
risk of estimator trained on task A, orange line denotes risk of estimator
trained on task A then task B.

Figure 4 shows the result of the numerical experiment in Figure 3 but with a
random permutation matrix instead of a random orthogonal matrix. We observe
that the same behavior holds in this scenario.

##  Appendix D Equivalence of Models and Derivation of Risk

Recall in Section 2.1 where we defined the LSM model for linear regression. In
this section we show that LSM is equivalent to an anisotropic regression model
(ARM). We then use ARM to define the risk expression that we analyze
theoretically.

We begin by defining ARM. Define data matrix $X\in\mathbb{R}^{n\times p}$ and
responses $y=X\beta+\epsilon$ where
$\epsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$,
$\sigma^{2}=\theta^{\top}(W^{\top}W+I_{d})^{-1}\theta$, and
$\beta=(I+WW^{\top})^{-1}W\theta$ for some
$\theta\in\mathbb{R}^{d},W\in\mathbb{R}^{p\times d}$. Let rows $X_{i}$ be
independent random vectors in $\mathbb{R}^{p}$ with covariance
$\Sigma=WW^{\top}+I_{p}$. Then the model is defined by the distribution over
$(X,y)$.

Next we show that ARM is equivalent to LSM. First observe that for both
models, $(y_{i},x_{i}^{\top})\in\mathbb{R}^{p+1}$ are centered Gaussian
vectors. Thus to show that they induce the same distribution, it suffices to
show that they have the same covariance.

$\displaystyle
Cov((y_{i},x_{i}^{\top})^{\top})=\mathbb{E}(y_{i},x_{i}^{\top})^{\top}(y_{i},x_{i}^{\top})=\mathbb{E}\begin{bmatrix}y_{i}^{2}&y_{i}x_{i}^{\top}\\\
y_{i}x_{i}&x_{i}x_{i}^{\top}\end{bmatrix}$ (25)

We then compute the covariance matrices for each model.

Under LSM, we have:

$\displaystyle\mathbb{E}[y_{i}^{2}]$
$\displaystyle=\mathbb{E}(\theta^{\top}z_{i})(z_{i}^{\top}\theta)=\theta^{\top}I\theta$
(26) $\displaystyle\mathbb{E}[y_{i}x_{i}]$
$\displaystyle=\mathbb{E}(Wz_{i}+u_{i})(z_{i}^{\top}\theta)=W\theta$ (27)
$\displaystyle\mathbb{E}[x_{i}x_{i}^{\top}]$
$\displaystyle=\mathbb{E}(Wx_{i}+u_{i})(z_{i}^{\top}W^{\top}+u_{i})=WW^{\top}+I$
(28)

Plugging these quantities into (25) gives:

$\displaystyle
Cov((y_{i},x_{i}^{\top})^{\top})=\begin{bmatrix}\|\theta\|^{2}&(W\theta)^{\top}\\\
W\theta&I+WW^{\top}\end{bmatrix}$ (29)

Under ARM, we have:

$\displaystyle\mathbb{E}[y_{i}^{2}]$
$\displaystyle=\mathbb{E}(\beta^{\top}x_{i}+\epsilon_{i})(x_{i}^{\top}\beta+\epsilon_{i})=\beta^{\top}(I+WW^{\top})\beta+\sigma^{2}$
(30) $\displaystyle\mathbb{E}[y_{i}x_{i}]$
$\displaystyle=\mathbb{E}(x_{i}(x_{i}^{\top}\beta+\epsilon_{i}))=\mathbb{E}(x_{i}x_{i}^{\top}\beta+x_{i}\epsilon_{i})=(I+WW^{\top})\beta$
(31) $\displaystyle\mathbb{E}[x_{i}x_{i}^{\top}]$ $\displaystyle=I+WW^{\top}$
(32)

Plugging these quantities into (25) gives:

$\displaystyle
Cov((y_{i},x_{i}^{\top})^{\top})=\begin{bmatrix}\beta^{\top}(I+WW^{\top})\beta+\sigma^{2}&((I+WW^{\top})\beta)^{\top}\\\
(I+WW^{\top})\beta&I+WW^{\top}\end{bmatrix}$ (33)

We now show equivalence of the covariance matrices. Recall that for ARM,
$\beta=(I+WW^{\top})^{-1}W\theta$ and
$\sigma^{2}=\theta^{\top}(I+W^{\top}W)^{-1}\theta$. We first show equivalence
of the first row, first column entries of the covariance matrices:

$\displaystyle\beta^{\top}(I+WW^{\top})\beta+\sigma^{2}$
$\displaystyle=\theta^{\top}W^{\top}(I+WW^{\top})^{-1}W\theta+\theta^{\top}(I+W^{\top}W)^{-1}\theta$
(34)
$\displaystyle=\theta^{\top}(W^{\top}(I+WW^{\top})^{-1}W+(I+W^{\top}W)^{-1})\theta$
(35)

By Lemma E.12, $(I+WW^{\top})^{-1}W=W(I+W^{\top}W)^{-1}$. This gives

$\displaystyle\beta^{\top}(I+WW^{\top})\beta+\sigma^{2}$
$\displaystyle=\theta^{\top}(W^{\top}W(I+W^{\top}W)^{-1}+(I+W^{\top}W)^{-1})\theta$
(36) $\displaystyle=\theta^{\top}(W^{\top}W+I)(I+W^{\top}W)^{-1}\theta$ (37)
$\displaystyle=\theta^{\top}\theta=\|\theta\|^{2}$ (38)

Next we show equivalence of the second row, first column entries of the
covariance matrices:

$\displaystyle(I+WW^{\top})\beta$
$\displaystyle=(I+WW^{\top})(I+WW^{\top})^{-1}W\theta=W\theta$ (39)

The equivalence of the first row, second column entries also follows from this
equality. The equivalence of the second row, second column entries is trivial.

Finally we derive the expression for the risk of ARM. By definition, the risk
of an estimator $f$ with parameters $\hat{\beta}$ has the following form:

$\displaystyle R(f_{\hat{\beta}})$
$\displaystyle=\mathbb{E}_{x,y}\|\hat{\beta}^{\top}x-y\|^{2}$ (40)
$\displaystyle=\mathbb{E}_{x,\epsilon}\|\hat{\beta}^{\top}x-\beta^{\top}x-\epsilon\|^{2}$
(41)
$\displaystyle=\mathbb{E}_{x}\|\hat{\beta}^{\top}x-\beta^{\top}x\|^{2}+\mathbb{E}_{\epsilon}\|\epsilon\|^{2}$
(42)
$\displaystyle=\mathbb{E}_{x}\|(\hat{\beta}-\beta)^{\top}x\|^{2}+\mathbb{E}_{\epsilon}\|\epsilon\|^{2}$
(43)
$\displaystyle=\mathbb{E}_{x}(\hat{\beta}-\beta)^{\top}xx^{\top}(\hat{\beta}-\beta)+\mathbb{E}_{\epsilon}\|\epsilon\|^{2}$
(44)
$\displaystyle=(\hat{\beta}-\beta)^{\top}\Sigma(\hat{\beta}-\beta)+\sigma^{2}$
(45)

where the third equality holds from independence of $\epsilon$ and the sixth
equality holds by definition of covariance.

We choose to study ARM with a slightly different but equivalent expression for
$\beta$. Using Lemma E.12,
$\beta=(I+WW^{\top})^{-1}W\theta=W(W^{\top}W+I)^{-1}\theta$.

##  Appendix E Supporting Lemmas

We begin with an assumption, inspired by [10], that all non-zero singular
values of $W$ are equal.

######  Assumption E.1.

All non-zero singular values of $W$ are equal. Namely, $W^{\top}W=p\gamma
I_{d}$.

######  Lemma E.2.

Assume $W\in\mathbb{R}^{p\times d}$ satisfies Assumption E.1. Then

$\displaystyle WW^{\top}=p\gamma\mathcal{P}_{W}$ (46)

where $\mathcal{P}_{W}$ is the orthogonal projection onto the range of $W$.

###### Proof.

We have that

$\displaystyle WW^{\top}=WW^{\top}\frac{p\gamma}{p\gamma}=p\gamma
W\left(\frac{1}{p\gamma}I_{p}\right)W^{\top}$ (47)

By Assumption E.1, we have

$\displaystyle WW^{\top}=p\gamma W(W^{\top}W)^{-1}W^{\top}$ (48)

$W^{\top}W$ has full rank with probability 1, so $W(W^{\top}W)^{-1}W^{\top}$
is given explicitly by $\mathcal{P}_{W}$, which completes the proof. ∎

######  Lemma E.3.

Let $\Sigma=WW^{\top}+I_{p}$ where $W\in\mathbb{R}^{p\times d}$ satisfies
Assumption E.1. For some $\theta\in\mathbb{R}^{d}$, let
$\beta=W(W^{\top}W+I_{d})^{-1}\theta$. Then

$\displaystyle R(f_{\vec{0}})=\|\theta\|^{2}$ (49)

###### Proof.

We have that

$\displaystyle R(f_{\vec{0}})$
$\displaystyle=(\vec{0}-\beta)^{\top}\Sigma(\vec{0}-\beta)+\sigma^{2}$ (50)
$\displaystyle=\beta^{\top}\Sigma\beta+\sigma^{2}$ (51)

By Lemma E.2, $\Sigma=p\gamma\mathcal{P}_{W}+I_{p}$ where $\mathcal{P}_{W}$ is
the orthogonal projection onto the range of $W$, which gives

$\displaystyle R(f_{\vec{0}})$
$\displaystyle=\beta^{\top}(p\gamma\mathcal{P}_{W}+I_{p})\beta+\sigma^{2}$
(52)

Since $\beta\in range(W)$,

$\displaystyle R(f_{\vec{0}})$
$\displaystyle=(p\gamma+1)\|\beta\|^{2}+\sigma^{2}$ (53)
$\displaystyle=(p\gamma+1)\theta^{\top}(W^{\top}W+I_{d})^{-1}W^{\top}W(W^{\top}W+I_{d})^{-1}\theta+\theta^{\top}(W^{\top}W+I_{d})^{-1}\theta$
(54) $\displaystyle=(p\gamma+1)\theta^{\top}((p\gamma+1)I_{d})^{-1}p\gamma
I_{d}((p\gamma+1)I_{d})^{-1}\theta+\theta^{\top}((p\gamma+1)I_{d})^{-1}\theta$
(55)
$\displaystyle=\left(\frac{p\gamma}{p\gamma+1}+\frac{1}{p\gamma+1}\right)\|\theta\|^{2}=\|\theta\|^{2}$
(56)

∎

######  Lemma E.4.

Let $x\sim\mathcal{N}(0,I_{d})$, and $\epsilon\leq 1$, then

$\mathbb{P}\left(d(1-\epsilon)\leq\|x\|_{2}^{2}\leq d(1+\epsilon)\right)\geq
1-e^{-c\epsilon^{2}d}$

where $c>0$ is an absolute constant.

###### Proof.

This statement follows from Corollary 5.17 in [21], concerning concentration
of sub-exponential random variables. ∎

######  Lemma E.5.

Assume $W\in\mathbb{R}^{p\times d}$ satisfies Assumption E.1. Let $O$ be a
random $p\times p$ orthogonal matrix. Fix $v\in\mathbb{R}^{p\times p}$. Then,
with probability at least $1-2e^{-c_{1}d}$,

$\displaystyle\|\mathcal{P}_{W}Ov\|^{2}\leq\frac{2d}{p}\|v\|^{2},$ (57)

for some universal constant $c_{1}>0$.

###### Proof.

Let $x=Ov$, and note that $\|x\|=\|v\|$, $\|x\|>0$ with probability 1 and
$\frac{x}{\|x\|}\sim\text{Uniform}(\mathbb{S}^{p-1})$. Letting
$z\sim\mathcal{N}(0,I_{p})$, we have that

$\displaystyle\|\mathcal{P}_{W}O\mathcal{P}_{A^{\top}}v\|$
$\displaystyle\leq\Bigl{\|}\mathcal{P}_{W}\frac{x}{\|x\|}\Bigr{\|}\|v\|\stackrel{{\scriptstyle
d}}{{=}}\Bigl{\|}\mathcal{P}_{W}\frac{z}{\|z\|}\Bigr{\|}\|v\|$ (58)

where the symbol $\stackrel{{\scriptstyle d}}{{=}}$ means equality in
distribution. Applying Lemma E.4 twice, we get that for any $\epsilon<1$, with
probability at least $1-e^{-c\epsilon^{2}p}-e^{-c\epsilon^{2}d}$,

$\displaystyle\|\mathcal{P}_{W}\frac{z}{\|z\|}\|\|v\|\leq\frac{\|\mathcal{P}_{W}z\|}{\sqrt{p}\sqrt{1-\epsilon}}\|v\|\leq\frac{\sqrt{d}\sqrt{1+\epsilon}}{\sqrt{p}\sqrt{1-\epsilon}}\|v\|$
(59)

for some universal constant $c>0$. By choosing suitable $\epsilon$, we obtain
that for $c_{1}=c\epsilon^{2}$, with probability at least $1-2e^{-c_{1}d}$,
$\|\mathcal{P}_{W}Ov\|^{2}\leq\frac{2d}{p}\|v\|^{2}$. ∎

######  Lemma E.6.

Define $A\in\mathbb{R}^{n\times p}$ with rows $A_{i}$ as independent random
vectors in $\mathbb{R}^{p}$ with covariance $\Sigma=WW^{\top}+I_{p}$ where
$W\in\mathbb{R}^{p\times d}$ satisfies Assumption E.1. Let $O$ be a random
$p\times p$ orthogonal matrix. Fix $v\in range(A^{\top})$. Then with
probability at least $1-2e^{-c_{1}n}$

$\displaystyle\|\mathcal{P}_{OA^{\top}}v\|^{2}\leq\frac{2n}{p}\|v\|^{2}$ (60)

for some universal constant $c_{1}>0$.

###### Proof.

We have that

$\displaystyle\|\mathcal{P}_{OA^{\top}}\mathcal{P}_{A^{\top}}v\|$
$\displaystyle=\|O\mathcal{P}_{A^{\top}}O^{\top}\mathcal{P}_{A^{\top}}v\|=\|\mathcal{P}_{A^{\top}}O^{\top}\mathcal{P}_{A^{\top}}v\|$
(61)
$\displaystyle=\|\mathcal{P}_{A^{\top}}\frac{O^{\top}\mathcal{P}_{A^{\top}}v}{\|O^{\top}\mathcal{P}_{A^{\top}}v\|}\|\|O^{\top}\mathcal{P}_{A^{\top}}v\|$
(62)
$\displaystyle\leq\|\mathcal{P}_{A^{\top}}\frac{O^{\top}\mathcal{P}_{A^{\top}}v}{\|O^{\top}\mathcal{P}_{A^{\top}}v\|}\|\|v\|$
(63) $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\|\mathcal{P}_{A^{\top}}\frac{z}{\|z\|}\|\|v\|$ (64)

where $z\sim\mathcal{N}(0,I_{p})$ and the last equality follows from the
rotational invariance of $O$. Applying Lemma E.4 twice, we get that for any
$\epsilon<1$, with probability at least
$1-e^{-c\epsilon^{2}p}-e^{-c\epsilon^{2}n}$,

$\displaystyle\|\mathcal{P}_{A^{\top}}\frac{z}{\|z\|}\|\|v\|\leq\frac{\|\mathcal{P}_{A^{\top}}z\|}{\sqrt{p}\sqrt{1-\epsilon}}\|v\|\leq\frac{\sqrt{n}\sqrt{1+\epsilon}}{\sqrt{p}\sqrt{1-\epsilon}}\|v\|$
(65)

for some universal constant $c>0$. By choosing suitable $\epsilon$, we obtain
that for $c_{1}=c\epsilon^{2}$, with probability at least $1-2e^{-c_{1}n}$,
$\|\mathcal{P}_{OA^{\top}}v\|^{2}\leq\frac{2n}{p}\|v\|^{2}$. ∎

######  Lemma E.7.

Let $a\in\mathbb{R}^{p}$ be generated by $\mathcal{N}(0,\Sigma)$,
$\Sigma=WW^{\top}+I_{p}$ where $W\in\mathbb{R}^{p\times d}$ satisfies
Assumption E.1. Then

$\displaystyle\mathbb{E}\|a\|_{2}^{2}=dp\gamma+p$ (66)

###### Proof.

It holds that $\mathbb{E}\|a\|_{2}^{2}=\|\Sigma^{1/2}\|_{F}^{2}$ [21].

$\displaystyle\|\Sigma^{1/2}\|_{F}^{2}$ $\displaystyle=\|\Sigma\|_{*}$ (67)

where $\|\|_{*}$ denotes the nuclear norm. By Lemma E.2,
$\Sigma=p\gamma\mathcal{P}_{W}+I_{p}$, which has $d$ singular values of
$p\gamma+1$ and $p-d$ singular values of 1. So

$\displaystyle\|\Sigma\|_{*}=d(p\gamma+1)+p-d=dp\gamma+p$ (68)

∎

######  Lemma E.8.

Define $A\in\mathbb{R}^{n\times p}$ with rows $A_{i}$ independent random
vectors in $\mathbb{R}^{p}$ with covariance
$\Sigma=p\gamma\mathcal{P}_{W}+I_{p}$ where $W\in\mathbb{R}^{p\times d}$
satisfies Assumption E.1. Then with probability at least $1-2e^{-n}$,

$\displaystyle\sigma_{min}(A)^{2}\geq(\sqrt{p-d}-2\sqrt{n})^{2}$ (69)

###### Proof.

WLOG let $range(W)=span(e_{1},...,e_{d})$. Then we can decompose $A$ into two
pieces: $A_{(1)}\in\mathbb{R}^{n\times d}$ with i.i.d.
$\mathcal{N}(0,p\gamma+1)$ entries and $A_{(2)}\in\mathbb{R}^{n\times p-d}$
with i.i.d. $\mathcal{N}(0,1)$ entries. This gives

$\displaystyle\sigma_{min}(A)^{2}$
$\displaystyle=\sigma_{min}(AA^{\top})=\sigma_{min}(A_{(1)}A_{(1)}^{\top}+A_{(2)}A_{(2)}^{\top})$
(70)
$\displaystyle\geq\sigma_{min}(A_{(2)}A_{(2)}^{\top})=\sigma_{min}(A_{(2)})^{2}=\sigma_{min}(A_{(2)}^{\top})^{2}$
(71)

By Theorem 5.39 in [21],
$\sigma_{min}(A_{(2)}^{\top})^{2}\geq(\sqrt{p-d}-2\sqrt{n})^{2}$ with
probability at least $1-2e^{-n}$. ∎

######  Lemma E.9.

Define $A\in\mathbb{R}^{n\times p}$ with rows $A_{i}$ independent random
vectors in $\mathbb{R}^{p}$ with covariance $\Sigma=WW^{\top}+I_{p}$ where
$W\in\mathbb{R}^{p\times d}$ satisfies Assumption E.1. Let
$\epsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$ and
$\sigma^{2}=\theta^{\top}(W^{\top}W+I_{d})^{-1}\theta$ for some
$\theta\in\mathbb{R}^{d}$. Then with probability at least $1-2e^{-n}$,

$\displaystyle\|A^{\dagger}\epsilon\|^{2}\leq\frac{n\|\theta\|^{2}}{p\gamma(\sqrt{p-d}-2\sqrt{n})^{2}},$
(72)

where $A^{\dagger}$ is the pseudoinverse of $A$.

###### Proof.

We have that

$\displaystyle\|A^{\dagger}\epsilon\|^{2}\leq\|A^{\dagger}\|^{2}\|\epsilon\|^{2}$
(73)

We have that
$\|\epsilon\|^{2}=n\sigma^{2}=n\theta^{\top}(W^{\top}W+I_{d})^{-1}\theta$.
Under Assumption E.1, this gives
$\|\epsilon\|^{2}=\frac{n\|\theta\|^{2}}{p\gamma+1}$,

$\displaystyle\|A^{\dagger}\epsilon\|^{2}\leq\|A^{\dagger}\|^{2}\frac{n\|\theta\|^{2}}{p\gamma+1}$
(74)

It holds that $\|A^{\dagger}\|=1/\sigma_{min}(A^{\top})$ where
$\sigma_{min}(A^{\top})$ is the smallest singular value of $A^{\top}$. By
Lemma E.8, $\sigma_{min}(A)\geq\sqrt{p-d}-2\sqrt{n}$ with probability at least
$1-2e^{-n}$. This gives

$\displaystyle\|A^{\dagger}\epsilon\|^{2}\leq\frac{1}{(\sqrt{p-d}-2\sqrt{n})^{2}}\cdot\frac{n\|\theta\|^{2}}{p\gamma+1}\leq\frac{n\|\theta\|^{2}}{p\gamma(\sqrt{p-d}-2\sqrt{n})^{2}}$
(75)

∎

######  Lemma E.10.

Suppose $W\in\mathbb{R}^{p\times d}$ satisfies Assumption E.1. Let
$\beta=W(W^{\top}W+I_{d})^{-1}\theta$ for some $\theta\in\mathbb{R}^{d}$. If
$p\geq 1/\gamma$ and $p\geq 16n+d$, then

$\displaystyle\frac{\sqrt{n}\|\theta\|}{\sqrt{p\gamma}(\sqrt{p-d}-2\sqrt{n})}\leq\|\beta\|$
(76)

###### Proof.

We have that
$\|\beta\|^{2}=\theta^{\top}(W^{\top}W+I_{d})^{-1}W^{\top}W(W^{\top}W+I_{d})\theta$.
Using Assumption E.1, this gives
$\|\beta\|=\frac{\sqrt{p\gamma}}{p\gamma+1}\|\theta\|$. Suppose $p\geq
1/\gamma$, then we have that
$\|\beta\|\geq\frac{1}{2\sqrt{p\gamma}}\|\theta\|$. When $p\geq 16n+d$,
$\frac{\sqrt{n}}{\sqrt{p-d}-2\sqrt{n}}\leq\frac{1}{2}$, which gives

$\displaystyle\frac{\sqrt{n}\|\theta\|}{\sqrt{p\gamma}(\sqrt{p-d}-2\sqrt{n})}\leq\frac{1}{2\sqrt{p\gamma}}\|\theta\|\leq\|\beta\|$
(77)

∎

######  Theorem E.11.

Define data matrix $A\in\mathbb{R}^{n\times p}$ and responses
$y=A\beta+\epsilon$ where $\epsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$,
$\sigma^{2}=\theta^{\top}(W^{\top}W+I_{d})^{-1}\theta$, and
$\beta=W(W^{\top}W+I_{d})^{-1}\theta$ for some $\theta\in\mathbb{R}^{d}$. Let
rows $A_{i}$ be independent random vectors in $\mathbb{R}^{p}$ with covariance
$\Sigma=WW^{\top}+I_{p}$ where $W\in\mathbb{R}^{p\times d}$ follows Assumption
E.1 and $n\geq d,p\geq max(17n,1/\gamma)$. Let $O$ be a random $p\times p$
orthogonal matrix and $B=AO^{\top}$. Let $\hat{\beta}_{A}$ be the parameters
of the minimum norm estimator on $A$, and $\hat{\beta}_{BA}$ be the parameters
of the estimator on $B$ using $\hat{\beta}_{A}$ as initialization. Let
$R(f_{\hat{\beta}})$ be the risk on task A of an estimator with parameters
$\hat{\beta}$. Then there exists constant $c>0$ such that with probability at
least $1-10e^{-cd}$, the following holds:

$\displaystyle
R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})\leq\left(\frac{66\sqrt{n}}{\sqrt{p}}+\frac{12}{p\gamma}\right)\|\theta\|^{2}$
(78)

###### Proof.

We have that

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle=(\hat{\beta}_{BA}-\beta)^{\top}\Sigma(\hat{\beta}_{BA}-\beta)-(\hat{\beta}_{A}-\beta)^{\top}\Sigma(\hat{\beta}_{A}-\beta)$
(79)
$\displaystyle=(\hat{\beta}_{A}+\hat{\beta}_{B}-\mathcal{P}_{B^{\top}}\hat{\beta}_{A}-\beta)^{\top}\Sigma(\hat{\beta}_{A}+\hat{\beta}_{B}-\mathcal{P}_{B^{\top}}\hat{\beta}_{A}-\beta)$
$\displaystyle\qquad-(\hat{\beta}_{A}-\beta)^{\top}\Sigma(\hat{\beta}_{A}-\beta)$
(80)

Distributing terms with $\hat{\beta}_{B}$ and
$\mathcal{P}_{B^{\top}}\hat{\beta}_{A}$ gives

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle=(\hat{\beta}_{A}-\beta)^{\top}\Sigma(\hat{\beta}_{A}-\beta)-(\hat{\beta}_{A}-\beta)^{\top}\Sigma(\hat{\beta}_{A}-\beta)+2\hat{\beta}_{A}^{\top}\Sigma\hat{\beta}_{B}-2\hat{\beta}_{B}^{\top}\Sigma\beta$
$\displaystyle\qquad+\hat{\beta}_{B}^{\top}\Sigma\hat{\beta}_{B}-2\hat{\beta}_{A}^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}-2\hat{\beta}_{B}^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}+2\beta^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}$
$\displaystyle\qquad+(\mathcal{P}_{B^{\top}}\hat{\beta}_{A})^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}$
(81)
$\displaystyle=2\hat{\beta}_{A}^{\top}\Sigma\hat{\beta}_{B}-2\hat{\beta}_{B}^{\top}\Sigma\beta+\hat{\beta}_{B}^{\top}\Sigma\hat{\beta}_{B}-2\hat{\beta}_{A}^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}-2\hat{\beta}_{B}^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}$
$\displaystyle\qquad+2\beta^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}+(\mathcal{P}_{B^{\top}}\hat{\beta}_{A})^{\top}\Sigma\mathcal{P}_{B^{\top}}\hat{\beta}_{A}$
(82)

By Lemma E.2, $\Sigma=p\gamma\mathcal{P}_{W}+I_{p}$ where $\mathcal{P}_{W}$ is
the orthogonal projection onto the range of $W$. This implies that
$\|\Sigma\|=p\gamma+1$, giving

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle=2\hat{\beta}_{A}^{\top}(p\gamma\mathcal{P}_{W}+I_{p})\hat{\beta}_{B}-2\hat{\beta}_{B}^{\top}(p\gamma\mathcal{P}_{W}+I_{p})\beta+\hat{\beta}_{B}^{\top}(p\gamma\mathcal{P}_{W}+I_{p})\hat{\beta}_{B}$
$\displaystyle\qquad+(\mathcal{P}_{B^{\top}}\hat{\beta}_{A})^{\top}(p\gamma\mathcal{P}_{W}+I_{p})(\mathcal{P}_{B^{\top}}\hat{\beta}_{A}-2\hat{\beta}_{A}-2\hat{\beta}_{B}+2\beta)$
(83) $\displaystyle\leq
2p\gamma\hat{\beta}_{A}^{\top}\mathcal{P}_{W}\hat{\beta}_{B}+2\hat{\beta}_{A}^{\top}\hat{\beta}_{B}-2p\gamma\hat{\beta}_{B}^{\top}\mathcal{P}_{W}\beta-2\hat{\beta}_{B}^{\top}\beta+p\gamma\hat{\beta}_{B}^{\top}\mathcal{P}_{W}\hat{\beta}_{B}+\hat{\beta}_{B}^{\top}\hat{\beta}_{B}$
$\displaystyle\qquad+(p\gamma+1)\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}-2\hat{\beta}_{A}-2\hat{\beta}_{B}+2\beta\|$
(84)

Applying Cauchy-Schwarz and triangle inequality gives the following bound:

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq
2p\gamma\|\hat{\beta}_{A}\|\|\mathcal{P}_{W}\hat{\beta}_{B}\|+2\|\hat{\beta}_{A}\|\|\hat{\beta}_{B}\|+2p\gamma\|\beta\|\|\mathcal{P}_{W}\hat{\beta}_{B}\|+2\|\hat{\beta}_{B}\|\|\beta\|$
$\displaystyle\qquad+p\gamma\|\hat{\beta}_{B}\|\|\mathcal{P}_{W}\hat{\beta}_{B}\|+\|\hat{\beta}_{B}\|^{2}$
$\displaystyle\qquad+(p\gamma+1)\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|(\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|+2\|\hat{\beta}_{A}\|+2\|\hat{\beta}_{B}\|+2\|\beta\|)$
(85)

By definition in Section 2,
$\hat{\beta}_{A}=\mathcal{P}_{A^{\top}}\beta+A^{\dagger}\epsilon$ and
$\hat{\beta}_{B}=O\mathcal{P}_{A^{\top}}\beta+OA^{\dagger}\epsilon$ and it
holds that $\|\mathcal{P}_{A^{\top}}\beta\|\leq\|\beta\|$. By Lemma E.9,
$\|A^{\top}(AA^{\top})^{-1}\epsilon\|\leq\frac{\sqrt{n}\|\theta\|}{\sqrt{p\gamma}(\sqrt{p-d}-2\sqrt{n})}$
with probability at least $1-2e^{-n}$ (call this Event E). So by Lemma E.10 if
$p\geq max(17n,1/\gamma)$, then $\|\hat{\beta}_{A}\|\leq 2\|\beta\|$ and
$\|\hat{\beta}_{B}\|\leq 2\|\beta\|$, which gives

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq
8p\gamma\|\beta\|\|\mathcal{P}_{W}\hat{\beta}_{B}\|+14(p\gamma+1)\|\beta\|\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|+12\|\beta\|^{2}$
(86)

We have that
$\|\beta\|^{2}=\theta^{\top}(W^{\top}W+I_{d})^{-1}W^{\top}W(W^{\top}W+I_{d})\theta$.
Using Assumption E.1, this gives
$\|\beta\|=\frac{\sqrt{p\gamma}}{p\gamma+1}\|\theta\|$,

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq
8\frac{p\gamma\sqrt{p\gamma}}{p\gamma+1}\|\theta\|\|\mathcal{P}_{W}\hat{\beta}_{B}\|+14\sqrt{p\gamma}\|\theta\|\|\mathcal{P}_{B^{\top}}\hat{\beta}_{A}\|+12\frac{p\gamma}{(p\gamma+1)^{2}}\|\theta\|^{2}$
(87) $\displaystyle=I+II+III$ (88)

We will bound each of these terms separately, starting with $I$:

Substituting $\hat{\beta}_{B}=B^{\top}(BB^{\top})^{-1}y$ into this expression
and distributing accordingly, we get that
$\mathcal{P}_{W}\hat{\beta}_{B}=\mathcal{P}_{W}O\mathcal{P}_{A^{\top}}\beta+\mathcal{P}_{W}OA^{\top}(AA^{\top})^{-1}\epsilon$,

$\displaystyle I$
$\displaystyle=8\frac{p\gamma\sqrt{p\gamma}}{p\gamma+1}\|\theta\|\|\mathcal{P}_{W}O\mathcal{P}_{A^{\top}}\beta+\mathcal{P}_{W}OA^{\top}(AA^{\top})^{-1}\epsilon\|$
(89) $\displaystyle\leq
8\frac{p\gamma\sqrt{p\gamma}}{p\gamma+1}\|\theta\|\|\mathcal{P}_{W}O\mathcal{P}_{A^{\top}}\beta\|+8\frac{p\gamma\sqrt{p\gamma}}{p\gamma+1}\|\theta\|\|\mathcal{P}_{W}OA^{\top}(AA^{\top})^{-1}\epsilon\|$
(90)

By Lemma E.5, there exists constant $c_{1}>0$ such that
$\|\mathcal{P}_{W}O\mathcal{P}_{A^{\top}}\beta\|\leq
1.5\sqrt{\frac{d}{p}}\|\beta\|$ and
$\|\mathcal{P}_{W}OA^{\top}(AA^{\top})^{-1}\epsilon\|\leq
1.5\sqrt{\frac{d}{p}}\|A^{\top}(AA^{\top})^{-1}\epsilon\|$ with probability at
least $1-2e^{-c_{1}d}$ each. By Lemma E.9,
$\|A^{\top}(AA^{\top})^{-1}\epsilon\|\leq\frac{\sqrt{n}\|\theta\|}{\sqrt{p\gamma}(\sqrt{p-d}-2\sqrt{n})}$
(failure probability already accounted for on Event E). This gives the
following bound with probability at least $1-4e^{-c_{1}d}$:

$\displaystyle I$ $\displaystyle\leq
12\frac{p\gamma\sqrt{d\gamma}}{p\gamma+1}\|\theta\|\|\beta\|+12\frac{\gamma\sqrt{ndp}}{(p\gamma+1)(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}$
(91)

Substituting $\|\beta\|=\frac{\sqrt{p\gamma}}{p\gamma+1}\|\theta\|$ gives

$\displaystyle I\leq
12\frac{p\gamma^{2}\sqrt{dp}}{(p\gamma+1)^{2}}\|\theta\|^{2}+12\frac{\gamma\sqrt{ndp}}{(p\gamma+1)(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}$
(92)

Using $p\gamma+1>p\gamma$ gives

$\displaystyle I$ $\displaystyle\leq
12\frac{\sqrt{d}}{\sqrt{p}}\|\theta\|^{2}+12\frac{\sqrt{nd}}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}$
(93)

Now we bound term II:

Substituting $\hat{\beta}_{A}=A^{\top}(AA^{\top})^{-1}y$ into this expression
and distributing accordingly, we get that
$\mathcal{P}_{B^{\top}}\hat{\beta}_{A}=\mathcal{P}_{B^{\top}}\mathcal{P}_{A^{\top}}\beta+\mathcal{P}_{B^{\top}}A^{\top}(AA^{\top})^{-1}\epsilon$.
This gives the following bound:

$\displaystyle II$
$\displaystyle=14\sqrt{p\gamma}\|\theta\|\|\mathcal{P}_{B^{\top}}\mathcal{P}_{A^{\top}}\beta+\mathcal{P}_{B^{\top}}A^{\top}(AA^{\top})^{-1}\epsilon\|$
(94) $\displaystyle\leq
14\sqrt{p\gamma}\|\theta\|\|\mathcal{P}_{B^{\top}}\mathcal{P}_{A^{\top}}\beta\|+14\sqrt{p\gamma}\|\theta\|\|\mathcal{P}_{B^{\top}}A^{\top}(AA^{\top})^{-1}\epsilon\|$
(95)

By Lemma E.6, there exists constant $c_{2}>0$ such that
$\|\mathcal{P}_{B^{\top}}\mathcal{P}_{A^{\top}}\beta\|\leq
1.5\sqrt{\frac{n}{p}}\|\beta\|$ and
$\|\mathcal{P}_{B^{\top}}A^{\top}(AA^{\top})^{-1}\epsilon\|\leq
1.5\sqrt{\frac{n}{p}}\|A^{\top}(AA^{\top})^{-1}\epsilon\|$ with probability at
least $1-2e^{-c_{2}n}$ each. By Lemma E.9,
$\|A^{\top}(AA^{\top})^{-1}\epsilon\|\leq\frac{\sqrt{n}\|\theta\|}{\sqrt{p\gamma}(\sqrt{p-d}-2\sqrt{n})}$
(failure probability already accounted for on Event E). This gives the
following bound with probability $1-4e^{-c_{2}n}$:

$\displaystyle II$ $\displaystyle\leq
21\sqrt{n\gamma}\|\theta\|\|\beta\|+21\frac{n}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}$
(96)

Substituting $\|\beta\|=\frac{\sqrt{p\gamma}}{p\gamma+1}\|\theta\|$ gives

$\displaystyle II$ $\displaystyle\leq
21\frac{\gamma\sqrt{np}}{p\gamma+1}\|\theta\|^{2}+21\frac{n}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}$
(97)

Using $p\gamma+1>p\gamma$ gives

$\displaystyle II$ $\displaystyle\leq
21\frac{\sqrt{n}}{\sqrt{p}}\|\theta\|^{2}+21\frac{n}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}$
(98)

Lastly we bound term $III$. Using $p\gamma+1>p\gamma$ gives the following
bound:

$\displaystyle III\leq\frac{12}{p\gamma}\|\theta\|^{2}$ (99)

Putting all three terms together gives the following bound with probability
$1-10e^{-cd}$ where $c=min(c_{1},c_{2})$:

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq
12\frac{\sqrt{d}}{\sqrt{p}}\|\theta\|^{2}+12\frac{\sqrt{nd}}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}+21\frac{\sqrt{n}}{\sqrt{p}}\|\theta\|^{2}$
(100)
$\displaystyle\qquad+21\frac{n}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}+\frac{12}{p\gamma}\|\theta\|^{2}$
(101)
$\displaystyle=\frac{12\sqrt{d}+21\sqrt{n}}{\sqrt{p}}\|\theta\|^{2}+\frac{12\sqrt{nd}+21n}{\sqrt{p}(\sqrt{p-d}-2\sqrt{n})}\|\theta\|^{2}+\frac{12}{p\gamma}\|\theta\|^{2}$
(102)

Using $d\leq n$ gives the following bound:

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq\left(\frac{33\sqrt{n}}{\sqrt{p}}+\frac{33n}{\sqrt{p}(\sqrt{p-n}-2\sqrt{n})}+\frac{12}{p\gamma}\right)\|\theta\|^{2}$
(103)

By assumption, $p\geq 17n$, so
$\frac{n}{\sqrt{p}(\sqrt{p-n}-2\sqrt{n})}\leq\frac{\sqrt{n}}{\sqrt{p}}$, which
gives the following bound:

$\displaystyle R(f_{\hat{\beta}_{BA}})-R(f_{\hat{\beta}_{A}})$
$\displaystyle\leq\left(\frac{66\sqrt{n}}{\sqrt{p}}+\frac{12}{p\gamma}\right)\|\theta\|^{2}$
(104)

∎

######  Lemma E.12.

Let $W\in\mathbb{R}^{p\times d}$. Then

$\displaystyle(I+WW^{\top})^{-1}W=W(W^{\top}W+I)^{-1}$ (105)

###### Proof.

Let $W=USV$ be the SVD of $W$ where $U\in\mathbb{R}^{p\times
d},S\in\mathbb{R}^{d\times d},V\in\mathbb{R}^{d\times d}$. Then we have

$\displaystyle(I+WW^{\top})^{-1}W$
$\displaystyle=(I+USVV^{\top}SU^{\top})^{-1}USV$ (106)
$\displaystyle=(I+US^{2}U^{\top})^{-1}USV$ (107)

Let $\tilde{U}\in\mathbb{R}^{p\times p}$ have the first $d$ columns be $U$ and
the last $p-d$ columns be the rest of the orthonormal basis. Then we have

$\displaystyle(I+WW^{\top})^{-1}W$
$\displaystyle=(\tilde{U}\tilde{U}^{\top}+US^{2}U^{\top})^{-1}USV$ (108)
$\displaystyle=(\tilde{U}\tilde{U}^{\top}+\tilde{U}\begin{bmatrix}S^{2}&0\\\
0&0\end{bmatrix}\tilde{U}^{\top})^{-1}USV$ (109)
$\displaystyle=\left(\tilde{U}\begin{bmatrix}I+S^{2}&0\\\
0&I\end{bmatrix}\tilde{U}^{\top}\right)^{-1}USV$ (110)
$\displaystyle=\tilde{U}\left(\begin{bmatrix}I+S^{2}&0\\\
0&I\end{bmatrix}\right)^{-1}\tilde{U}^{\top}USV$ (111)
$\displaystyle=\tilde{U}\left(\begin{bmatrix}I+S^{2}&0\\\
0&I\end{bmatrix}\right)^{-1}(I_{d},0)^{\top}SV$ (112)
$\displaystyle=\tilde{U}((I+S^{2})^{-1},0)^{\top}SV$ (113)
$\displaystyle=U(I+S^{2})^{-1}SV$ (114) $\displaystyle=US(I+S^{2})^{-1}V$
(115) $\displaystyle=USVV^{\top}(I+S^{2})^{-1}V$ (116)
$\displaystyle=W(I+W^{\top}W)^{-1}$ (117)

∎