# Randomly Weighted Neuromodulation in Neural Networks Facilitates Learning of
Manifolds Common Across Tasks

Jinyung Hong1 Theodore P. Pavlic1,2
1School of Computing and Augmented Intelligence
2School of Life Sciences
Arizona State University
Tempe, AZ 85281
{ jhong53, tpavlic<EMAIL_ADDRESS>

###### Abstract

Geometric Sensitive Hashing functions, a family of Local Sensitive Hashing
functions, are neural network models that learn class-specific manifold
geometry in supervised learning. However, given a set of supervised learning
tasks, understanding the manifold geometries that can represent each task and
the kinds of relationships between the tasks based on them has received little
attention. We explore a formalization of this question by considering a
generative process where each task is associated with a high-dimensional
manifold, which can be done in brain-like models with neuromodulatory systems.
Following this formulation, we define _Task-specific Geometric Sensitive
Hashing (T-GSH)_ and show that a randomly weighted neural network with a
neuromodulation system can realize this function.111All source codes and
figures are available at https://github.com/PavlicLab/NeurIPS2023-UniReps-
Hong-TGSH-Randomly_Weighted_Neuromodulation_in_Neural_Networks.git.

##  1 Introduction

Deep Neural Networks (DNNs) have shown outstanding performance in various
fields due to their ability to learn to transform complex objects into useful
representations that can be easily separated in the embedding space. However,
due to their “black box” behavior, there is still little known about exactly
which information-rich features are captured by a DNN representation and what
information-poor variations are eliminated. To address the question, Dikkala
et al. [10] considered the manifold geometry of a generative process where
each class is associated with a manifold and parameter that shifts class
examples along that manifold. Using this perspective, they could show that
neural representations in supervised learning could be viewed as a kind of
Locality Sensitive Hash functions that they coined Geometric Sensitive Hashing
(GSH) [10].

In this work, we extend this geometric interpretation from one classification
task to a series of potentially related classification tasks. For example, in
a video clip of digits rotating counterclockwise, each possible digit within a
single frame is a class associated with its own manifold, as in GSH, and we
focus on the relationship between the manfolds induced by the classification
problem in one frame (i.e., one rotation angle) and the other frames. Such a
manifold-learning perspective allows us to compare representations of tasks
that are similar in principle so that we can identify representational kernels
that are shared among each individual classification task. In essence, for
this example digit-rotation task, we consider the possibility that rotated
digit images may in fact map to different points on the same induced digit
manifold that is then combined with a rotation manifold to represent the image
and its orientation.

The geometry of low- and high-level perceptual spaces has been studied in
cognitive science, and changes in those behavioral and cognitive states have
been associated with neuromodulatory systems in the brain. Neuromodulation has
been studied extensively in human and insect brains, where neuromodulatory
signals act as a kind of _switchboard_ that can remap a neural representation
in different ways so as to include adaptive control of behavior based on
internal state and environmental context [18, 16]. In response to the
importance of neuromodulatory systems in biological brains, neuromodulation-
inspired neural networks have been proposed in various domains, such as
reinforcement learning [4], modular networks [7], and adaptive behaviors [22].
For example, promising applications of neuromodulated neural networks to the
problem of to Continual Learning (CL) [7, 14, 17] can prevent catastrophic
forgetting by modulating synapsis that remain persistently stable.

In this work, we connect neuromodulation-inspired neural networks for
continual learning and geometric manifold learning in supervised learning. We
introduce the definition of _task-specific manifolds_ and demonstrate that
DNNs with neuromodulators can learn these manifolds.

(a) Data generating process.

(b) A DNN as a Task-specific Geometry Sensitive Hash (T-GSH) function.

Figure 1: (a) Each task $t$ consists of a set of points of classes on a simple
manifold. Furthermore, each point in the task can be characterized by three
latent parameters:
$\boldsymbol{\gamma}^{t},\boldsymbol{\delta},\boldsymbol{\theta}^{t}$.
Notably, we consider a task-specific manifold broader than the class-level one
in [10]. (b) Any set of points on the same task manifold map to
(approximately) the same representation (i.e., the penultimate layer feature
map), while a set of points from different task manifolds go to far away
representations.

##  2 Task-Specific Geometric Sensitive Hashing

###  2.1 Problem Definition and Proposed Randomly Weighted Neuromodulation
Method

Following [14], we consider a supervised continual-learning setup where
$\mathcal{T}$ tasks arrive to a learner in sequential order. Let
$\mathcal{D}_{t}=\\{\mathbf{x}_{i,t},y_{i,t}\\}_{t}^{n_{t}}$ be the dataset of
task $t$, composed of $n_{t}$ pairs of raw instances and corresponding labels.
For simplicity, we adopt a simple task-incremental continual learning setup
where $\mathcal{C}$ is the number of classes for every task. Next, we extend
the manifolds from Dikkala et al. [10] to represent _task-specific_ manifolds.
We assume that each input $\mathbf{x}_{i,t}\in\mathbb{R}^{d}$ is drawn from a
mixture distribution over $\mathcal{T}$ manifolds
$M_{1},\dots,M_{\mathcal{T}}$ sharing similar topologies (Fig. 1(a)), and each
label $y_{i,t}\in[\mathcal{C}]$ corresponds to the manifold of
$\mathbf{x}_{i,t}$. Moreover, each point $\mathbf{x}_{t}$ on manifold $M_{t}$
is determined by the latent vectors,
$\boldsymbol{\gamma}^{t}\in\mathbb{R}^{s}$,
$\boldsymbol{\delta}\in\mathbb{R}^{p}$, and
$\boldsymbol{\theta}^{t}\in\mathbb{R}^{k}$, where:

  * •

$\boldsymbol{\gamma}^{t}$ is the manifold identifier for the task $t$, and
there is a one-to-one correspondence between $\boldsymbol{\gamma}^{t}$ and
$M_{t}$.

  * •

$\boldsymbol{\delta}$ is a set of manifold identifiers for classes. This is an
explicitly shared manifold, representing canonical representation objects and
is shared across all task manifolds $\boldsymbol{\gamma}$.

  * •

$\boldsymbol{\theta}^{t}$ is the transformation in the task $t$. So, if we fix
$\boldsymbol{\gamma}^{t}$, the manifold $M_{t}$ can be generated by sampling
different values of $\boldsymbol{\theta}^{t}$.

Intuitively, $\boldsymbol{\gamma}^{t}$ represents the centroid of a cluster
for the task $t$, $\boldsymbol{\delta}$ includes explicitly shared clusters
for labels that every $\boldsymbol{\gamma}$ contains, and
$\boldsymbol{\theta}^{t}$ represents a small perturbation around the centroid.
So, the manifold $M_{\boldsymbol{\gamma}^{t}}$ is comprised of all inputs
$\mathbf{x}_{t}$ of the form
$\\{\boldsymbol{\gamma}^{t}+\boldsymbol{\delta}+\boldsymbol{\theta}^{t}|\;\lVert\boldsymbol{\theta}^{t}\rVert\leq\epsilon\\}$.
In this setting, we employ a Locality Sensitive Hashing (LSH) to map each
input to a hash bucket. With appropriate configuration, we empirically
demonstrate that DNNs exhibit LSH-like behavior on the family of manifolds
with a shared geometry defined by a set of analytic functions. In particular,
extending the work of Dikkala et al. [10], we show that the penultimate
“representation layer” $r$ (Fig. 1(b)) of an appropriately trained network
will satisfy the following property:

######  Definition 1 (Task-specific Geometric Sensitive Hashing (T-GSH))

We say $r$ is a T-GSH function with respect to a set of task manifolds if:

  1. a)

For any two points on the same task manifold $x_{1,t},x_{2,t}\in M_{t}$, the
distance $\lVert r(x_{1,t})-r(x_{2,t})\rVert$ is small regardless of the
associated classes.

  2. b)

For any two points on two well-separated task manifolds $x_{1,t}\in
M_{t}\;\text{and}\;x_{2,t^{\prime}}\in M_{t^{\prime}},\lVert
r(x_{1,t})-r(x_{2,t^{\prime}})\rVert$ is large regardless of the associated
classes.

That is, DNNs whose representations satisfy the T-GSH properly capture the
shared task manifold geometry similar to how LSH functions capture spatial
locality. Finally, by employing the recoverability of the manifold [10], we
give empirical evidence to support the recoverability of
$\boldsymbol{\gamma}^{t}$ via simple transformations on top of representations
learned by DNNs which behave as T-GSH functions.

We consider task manifolds as subsets of points in $\mathbb{R}^{d}$. Every
manifold $M_{\boldsymbol{\gamma}^{t}}$ has an associated latent vector
$\boldsymbol{\gamma}^{t}\in\mathbb{R}^{s}$ (where $s\leq d$) which acts as an
identifier of $M_{\boldsymbol{\gamma}^{t}}$. The task manifold is then defined
to be the set of points
$\mathbf{x}_{t}=\mathbf{f}(\boldsymbol{\gamma}^{t},\boldsymbol{\delta},\boldsymbol{\theta}^{t})=(f_{1}(\boldsymbol{\gamma}^{t},\boldsymbol{\delta},\boldsymbol{\theta}^{t}),\dots,f_{d}(\boldsymbol{\gamma}^{t},\boldsymbol{\delta},\boldsymbol{\theta}^{t}))\;\text{for}\;\boldsymbol{\delta}\in\Delta\subseteq\mathbb{R}^{p},p<d,\;\text{and
for}\;\boldsymbol{\theta}^{t}\in\Theta\subseteq\mathbb{R}^{k},k<d$. Here, the
manifold generating function $\mathbf{f}=\\{f_{i}(\cdot,\cdot)\\}_{i=1}^{d}$
where the $f_{i}$ are all analytic functions. Without significant loss of
generality, we assume our inputs $\mathbf{x}$ and $\boldsymbol{\gamma}$ are
normalized and lie on $S^{d-1},$ and $S^{s-1}$, the $d$ and $s$-dimensional
unit spheres, respectively. Given the above generative process, we leverage
the Assumption 1, which allows us to use a special form of the DNN
configuration used by Dikkala et al. [10].

######  Assumption 1 (Modified Invertibility [10])

There is an analytic function
$g(\cdot):\mathbb{R}^{d}\rightarrow\mathbb{R}^{s}$ with bounded norm Taylor
expansion such that for every point
$\mathbf{x}_{t}=\mathbf{f}(\boldsymbol{\gamma}^{t},\boldsymbol{\delta},\boldsymbol{\theta}^{t})\;\text{on}\;M_{\boldsymbol{\gamma}^{t}},g(\mathbf{x}_{t})=\boldsymbol{\gamma}^{t}$.

The conventional GSH neural network [10] is defined as a single-hidden-layer
network, $\mathbf{y}=A\cdot B\cdot\sigma(C\mathbf{x})$, where the input
$\mathbf{x}\in\mathbb{R}^{d}$ passes through a wide, randomly initialized,
fully connected, non-trainable layer $C\in\mathbb{R}^{D\times d}$ followed by
a ReLU activation $\sigma(\cdot)$. Then, there are two trainable, fully
connected layers $A\in\mathbb{R}^{\mathcal{C}\times T}$,
$B\in\mathbb{R}^{T\times D}$ with no non-linearity between them. Following the
above configuration, we define a special-case T-GSH neural network as follows:

$\hat{\mathbf{y}}_{t}=\mathbf{R}\cdot B^{t}\cdot\sigma(C\mathbf{x})$ (1)

for each task $t$. Thus, not only does a T-GSH differ from a GSH by the
introduction of multiple task labels, but the learnable matrix $A$ from a
standard GSH is replaced with another randomly weighted matrix $\mathbf{R}$
that acts as the explicitly shared manifold $\boldsymbol{\delta}$ in our
definition (i.e., it is shared across all tasks). Intuitively, a randomly
weighted matrix is nearly orthogonal, and so learning task-specific $B^{t}$
must be enforced to learn the commonality in the task $t$. In the upcoming
section, we connect a T-GSH neural network with one of the examples of
neuromodulation-inspired DNNs for continual learning.

###  2.2 Revising Configurable Random Weight Networks

Configurable Random Weight Networks (CRWNs) [14] were proposed as a simple yet
effective form of artificial neuromodulatory systems for continual learning.
The proposed approach consists of two types of artificial neuromodulation,
_Global_ and _Local_ modulation, and shows that their neuromodulatory signals
learn how to adjust long-lasting/unchanging random synaptic weights for
specific tasks, enabling task-specific learning. Furthermore, the proposed
neuromodulation can be applied in a layer-wise fashion. To improve the
efficiency of training and the applicability to any network architecture, Hong
and Pavlic [14] proposed two kinds of model architectures, which we summarize
in the single expression in Eq. (2),

$\begin{split}\hat{\mathbf{y}}_{t}&=\alpha_{t}\cdot\mathbf{R}\cdot(v^{t}\odot\sigma(C\mathbf{x}))\\\
&=\mathbf{R}\cdot((\alpha_{t}\cdot v^{t})\odot\sigma(C\mathbf{x}))\end{split}$
(2)

where:

  * •

$\mathbf{R}\in\mathbb{R}^{\mathcal{C}\times D}$ and $C\in\mathbb{R}^{D\times
d}$ are non-trainable, randomly weighted matrices.

  * •

$\alpha_{t}\in\mathbb{R}$ is a learnable constant acting as _global_
neuromodulation.

  * •

$v^{t}\in\mathbb{R}^{D}$ is a learnable vector mimicking _local_
neuromodulation.

In the following, we view Eq. (2) through the lens of T-GSH in Eq. (1).

#### CRWN as a T-GSH function.

CRWNs were first proposed as an efficient CL method by highlighting the
approach’s benefits concerning computational efficiency. However, we highlight
here that CRWNs are a T-GSH function achieved by the proposed artificial
neuromodulatory systems and thus can be used to shed light on the manifold
learning facilitated by neuromodulatory processes. Below, we compare Eq. (1)
and (2).

$\hat{\mathbf{y}}_{t}=\mathord{\overbrace{\mathbf{R}\cdot
B^{t}\cdot\sigma(C\mathbf{x})}^{\text{{Eq.~{}\eqref{eq:nn_tgsh}}}}}=\mathord{\overbrace{\mathbf{R}\cdot\mathord{(\underbrace{(\alpha_{t}\cdot
v^{t})}_{B^{t}}}\odot\sigma(C\mathbf{x}))}^{\text{{Eq.~{}\eqref{eq:crwn}}}}}$

Thus, $\alpha_{t}\cdot v^{t}$ from the CRWN in Eq. (2) plays the role of
$B^{t}$ in the T-GSH in Eq. (1), and so a CRWN is a simple realization of a
T-GSH. Following Dikkala et al. [10], it is important to use a proper
regularization term on the weights $B^{t}$ in the classification loss for
training a T-GSH function, and the learnable constant $\alpha_{t}$ in CRWNs
acts as a similar regularizer for $v^{t}$.

##  3 Experiments

(a) A comparison of cosine similarities of the points. (Intra similarity): the
cosine similarities of the points with the _different labels_ on the _same
task manifold_. (Inter similarity): the cosine similarities of the points with
the _same label_ on the _different task manifolds_.

(b) A confusion matrix of intra (same task manifold) VS inter (different task
manifolds) cosine similarity of task representations trained on
_RotationMNIST_. Cosine similarity between task context vectors before (left)
and after training (right). Notably, adjacent tasks are expected to be the
most similar to each other as they represent the smallest angular difference
in images.

Figure 2: Experiment results on RotationMNIST.

In this section, we demonstrate the CRWN is a T-GSH function using variants of
the MNIST [15] dataset, including _RotationMNIST_ , and our newly proposed
experimental setup, _ShiftMNIST_ and _AugmentMNIST_.

###  3.1 Experiments for showing that CRWN is a T-GSH function.

We first validate the CRWN on RotationMNIST to show it is a T-GSH function.
RotationMNIST is one of the popular CL benchmark datasets. Following Hong and
Pavlic [14], we define the dataset so that each of the 36 tasks corresponds to
images counterclockwise rotated by a multiple of 10 degrees.

#### CRWNs as a GSH and T-GSH function.

Using experimental results of CRWNs for continual learning, we can empirically
demonstrate that CRWNs are T-GSH functions. In particular, Hong and Pavlic
[14, Table 1] showed that CRWNs achieved about $95\%$ test accuracy average
over all 36 RotationMNIST tasks (FlyNet: $94.9\%$ and NeuroModNet: $95.5\%$,
respectively). This indicates that CRWNs perform well as a GSH function that
successfully classifies all labels for every task.

Figure 2(a) depicts a comparison of cosine similarities of latent CRWN task
representations. _Intra similarity_ indicates the cosine similarities of the
points with different labels on the same task manifold (i.e., same
RotationMNIST rotation angle). In contrast, _Inter similarity_ means the
cosine similarity of the points with the same label, but on the different task
manifolds (i.e., different RotationMNIST rotational angles). This shows CRWN
is a T-GSH function in Fig. 1(b), showing that the similarity of the points
having different labels, but on the same task manifold must be larger than one
of the points with the same label, but on the different task manifolds.

#### Satisfaction of Assumption 1.

Figure 3: Empirical evidence of approximate data generating process with our
defined invertibility assumption on RotationMNIST. (Top Row) Reconstructed
samples of digit “3”. (Bottom Row) Reconstructed samples of digit “1”. (Left
Column) Reconstructed samples of digits from the task manifold $M^{T1}$, which
is $0$°rotation. (Middle Column) the samples of digits from the task manifold
$M^{T5}$, which is $40$ °counterclockwise rotation. (Right Column) the samples
of digits from the task manifold $M^{T10}$, which is $90$ °counterclockwise
rotation.

We can also show that the inverse function of CRWNs can mimic the data-
generating process, satisfying our Assumption 1. Figure 3 depicts the
reconstructed image samples from CRWNs from the learned task manifolds.

Because of the ReLU function, we cannot obtain the complete inverse function
of CRWNs. Thus, we omitted ReLU and computed the pseudo-inverse of the weights
of the learned CRWNs. Because we fixed the manifolds for classes
$\boldsymbol{\delta}$ using the non-trainable matrix $\mathbf{R}$, each row in
$\mathbf{R}$ can be seen as the canonical representations of objects. As shown
in Fig. 3, the canonical representation for the digit “3” is the 4th row, and
for the digit “1” is the 2nd row. Thus, let $\mathbf{f}$ be the inverse CRWNs
omitting the ReLU function. We can reconstruct a digit sample with class $i$
as follows:

$\bar{\mathbf{x}}_{i,t}=\mathbf{f}(\textbf{s})\;\text{where}\;\mathbf{s}\sim\mathcal{N}(\mu_{i,t},\sigma_{i,t}),\mu_{i,t}=R_{i}^{\top}\odot(\alpha_{t}\cdot
v_{t}),\sigma_{i,t}=\mathbf{1}_{D}\cdot 1/D$ (3)

$R_{i}$ is $i$-th row of the matrix $\mathbf{R}$ and $\mathbf{1}_{D}$ is the
vector containing ones with the size $D$. In Fig. 3, despite approximate
reconstruction, the shape of digits and the rotation can be shown, indicating
the CRWN satisfies Assumption 1.

###  3.2 Experiments For Task Similarities

In this section, due to the architectural characteristic of a T-GSH function,
we demonstrate unique features of the CRWN that allow for computing the
relationship among tasks. In particular, we introduce the scalar–vector
product $c_{t}\triangleq\alpha_{t}\cdot v^{t}$ so that the CRWN expression
becomes:

$\hat{\textbf{y}}_{t}=\mathbf{R}\cdot((\alpha_{t}\cdot
v^{t})\odot\sigma(C\mathbf{x}))=\mathbf{R}(c_{t}\odot\sigma(C\mathbf{x})).$

We refer to $c_{t}$ as a _context vector._ The following experiments on
RotationMNIST, ShiftMNIST, and AugmentMNIST demonstrate the semantic
interpretation of the context vector $c_{t}$ among tasks.

#### RotationMNIST.

Figure 2(b) shows the pairwise cosine similarity222We compute
$1-\texttt{scipy.spatial.distance.cosine}$. between context-vector pairs
$(c_{j},c_{k})$ for every $(j,k)$ pair of RotationMNIST tasks, shown both
before training and after training. In the task ordering, adjacent tasks are
expected to be the most similar to each other as they represent the smallest
angular difference in images. In other words, networks trained for task $i$
should have minimal degradation in performance for test data from tasks $i-1$
and $i+1$. Consistent with this functional expectation, the mechanistic
pattern Fig. 2(b) shows greatest similarity between adjacent context vectors
after training, and so trained context vectors exist in a representational
space that captures fundamental similarities among tasks.

#### ShiftMNIST.

(a) ShiftMNIST-1px tasks

(b) ShiftMNIST-1px results

Figure 4: ShiftMNIST-1px. (4(a)) The difference between the datasets of the
two adjacent tasks, $i$ and $i+1$, is that the images of task $i+1$ are the
images of task $i$ shifted a single pixel to the right. Thus, the total number
of tasks to learn is 28. (4(b)) Cosine similarity between the context vectors
for each task before (left) and after training (right).

(a) ShiftMNIST-2px tasks

(b) ShiftMNIST-2px results

Figure 5: ShiftMNIST-2px. (5(a)) The difference between the datasets of the
two adjacent tasks, $i$ and $i+1$, is that the images of task $i+1$ are the
images of task $i$ shifted two pixels to the right. Thus, the total number of
tasks to learn is 14. (5(b)) Cosine similarity between the context vectors for
each task before (left) and after training (right).

To test whether geometric adjacency generally implies similarity of trained
vectors, we introduce a new sequence of tasks, _ShiftMNIST_ , as shown in Figs
4(a) and 5(a). Whereas RotationMNIST rotates each of the MNIST characters,
ShiftMNIST translates them laterally by a number of pixels (with periodic
boundaries). Adjacent ShiftMNIST-1px tasks in Fig. 4(a) are separated by 1
pixel, and so an image in Task 3 is the same as an image in Task 2 shifted to
the right by 1 pixel. Similarly, adjacent ShiftMNIST-2px tasks in Fig. 5(a)
are separated by 2 pixels. Because MNIST images are each 28 pixels wide, there
are 28 distinct ShiftMNIST-1px tasks and 14 distinct ShiftMNIST-2px tasks. We
conduct two experimental setups separately using different random seeds and
check whether meaningful relationships with each other can be derived. The
cosine similarities of context vectors before and after training for
ShiftMNIST-1px and ShiftMNIST-2px are shown in Figs 4(b) and 5(b),
respectively. The diagonal pattern in both cases matches that of Fig. 2(b),
which confirms that the similarity vectors for ShiftMNIST have encoded
geometric relationships among tasks in a similar way as with RotationMNIST.
Furthermore, Fig. 5(b) appears as a more coarsely subsampled version of Fig.
4(b), which demonstrates that cosine similarity differences are repeatable and
a function of the tasks themselves and not an artifact of the architecture of
the network. Thus, training these compact context vectors extracts meaningful
semantic information about individual tasks.

#### AugmentMNIST.

Here, we propose AugmentMNIST, which employs a sequence of 8 off-the-shelf,
commonly used data-augmentation tasks as shown in Fig. 6 (e.g., Random Erasing
[24] that is provided by PyTorch). We list the details of the setting in Table
1. After training on each of the 8 tasks, we use hierarchical agglomerative
clustering to sort the tasks so that adjacent tasks tend to have highest
similarity.

Table 1: PyTorch routines used for AugmentMNIST. Here, TT is shorthand for torchvision.transforms. AugmentMNIST Task: | PyTorch method and configuration
---|---
Horizontal Flip (T2): | TT.RandomHorizontalFlip(p=1.0)
Vertical Flip (T3): | TT.RandomVerticalFlip(p=1.0)
Gaussian Blur (T4): | TT.GaussianBlur(kernel_size=5, sigma=2)
Perspective (T5): | TT.RandomPerspective(distortion_scale=0.5, p=1.0)
Random Erasing (T6): | TT.RandomErasing(p=1.0)
Invert (T7): | TT.RandomInvert(p=1.0)
RandomResizedCrop (T8): | TT.RandomResizedCrop(size=28)

Fig. 7 shows the resulting cosine similarities as well as a dendrogram of the
hierarchical clustering by cosine similarity.

Figure 6: AugmentMNIST task. The list of variants of augmentation methods
applied to produce the task. See the implementation details in Table 1. Figure
7: AugmentMNIST results. Sorted cosine similarity between the context vectors
for each task after training (left) and task-similarity clustering (right).

The results show that context-vector representations of Random Erasing (T6),
Gaussian Blur (T4), and Perspective (T5) are very close to that of the
original image (T1); in other words, the four tasks that maintain the shape,
color, and orientation of the image are grouped together. Furthermore, the
inverted-image task (T7) is an outgroup as it has negative cosine similarity
with all other tasks, which is consistent with it being the only task
operating on a white background. Furthermore, the cosine similarities suggest
a closer relationship between T1 and T7 than between T2 and T7, indicating
that the horizontal flip combined with the color inversion has an additive
effect in terms of context vector dissimilarity. Thus, the geometry of the
trained context-vector space encodes semantic information about comparisons
between tasks themselves.

##  4 Related Work

There have been various works on studying classification problems as manifold
learning [20, 3]. Furthermore, manifold learning has been adapted in various
literature, such as feature learning with different modalities [19],
explainability [11], speech recognition [21] and neuroscience [5]. We leverage
the concept of Locality Sensitive Hashing, and this property has been
extensively studied in insect-brain-like architecture [8], and several works
have been proposed in novelty detection [9] and relational learning [13, 12].
The connection between DNNs and Hash Functions was explored before [23].
Another related work is the benefits of wide non-linear layers [6] and over-
parameterized networks [2, 1]. We empirically demonstrate that T-GSH holds for
certain architectures under the manifold data assumption given a series of
classification tasks.

##  5 Discussion

We studied the problem of a sequence of multiple supervised classification
tasks as a task manifold-learning problem under a specific generative process
wherein the manifolds share geometry. We demonstrated empirically that
properly trained DNNs with a neuromodulation-inspired architecture satisfy the
T-GSH property by recovering each task’s semantically meaningful latent
representation $\gamma$. Our work provides new geometric intuitions about the
functioning of randomly weighted neuromodulation systems found throughout
nature and allows us to better understand how randomness helps exploit finite
representation spaces. We plan to extend our work to provide theoretical
evidence in various continual-learning setups, such as class-incremental and
domain-incremental learning.

## Acknowledgements

This work was supported in part by NSF award 2223839 and USACE ERDC award
W912HZ-21-2-0040.

## References

  * Allen-Zhu et al. [2019] Z. Allen-Zhu, Y. Li, and Y. Liang.  Learning and generalization in overparameterized neural networks, going beyond two layers.  _Advances in neural information processing systems_ , 32, 2019.
  * Arora et al. [2019] S. Arora, S. Du, W. Hu, Z. Li, and R. Wang.  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  In _International Conference on Machine Learning_ , pages 322–332. PMLR, 2019.
  * Belkin et al. [2006] M. Belkin, P. Niyogi, and V. Sindhwani.  Manifold regularization: A geometric framework for learning from labeled and unlabeled examples.  _Journal of machine learning research_ , 7(11), 2006.
  * Botvinick et al. [2020] M. Botvinick, J. X. Wang, W. Dabney, K. J. Miller, and Z. Kurth-Nelson.  Deep reinforcement learning and its neuroscientific implications.  _Neuron_ , 107(4):603–616, 2020.
  * Busch et al. [2023] E. L. Busch, J. Huang, A. Benz, T. Wallenstein, G. Lajoie, G. Wolf, S. Krishnaswamy, and N. B. Turk-Browne.  Multi-view manifold learning of human brain-state trajectories.  _Nature Computational Science_ , 3(3):240–253, 2023.
  * Daniely et al. [2016] A. Daniely, R. Frostig, and Y. Singer.  Toward deeper understanding of neural networks: The power of initialization and a dual view on expressivity.  _Advances in neural information processing systems_ , 29, 2016.
  * Daram et al. [2020] A. Daram, A. Yanguas-Gil, and D. Kudithipudi.  Exploring neuromodulation for dynamic learning.  _Front. Neurosci._ , page 928, 2020.
  * Dasgupta et al. [2017] S. Dasgupta, C. F. Stevens, and S. Navlakha.  A neural algorithm for a fundamental computing problem.  _Science_ , 358(6364):793–796, 2017.
  * Dasgupta et al. [2018] S. Dasgupta, T. C. Sheehan, C. F. Stevens, and S. Navlakha.  A neural data structure for novelty detection.  _PNAS USA_ , 115(51):13093–13098, 2018.
  * Dikkala et al. [2021] N. Dikkala, G. Kaplun, and R. Panigrahy.  For manifold learning, deep neural networks can be locality sensitive hash functions.  _arXiv preprint arXiv:2103.06875_ , 2021.
  * Han et al. [2022] H. Han, W. Li, J. Wang, G. Qin, and X. Qin.  Enhance explainability of manifold learning.  _Neurocomputing_ , 500:877–895, 2022.
  * Hong and Pavlic [2022a] J. Hong and T. Pavlic.  Representing prior knowledge using randomly, weighted feature networks for visual relationship detection.  In _The First International Workshop on Combining Learning and Reasoning: Programming Languages, Formalisms, and Representations (CLeaR) at the AAAI Conference on Artificial Intelligence_ , 2022a.
  * Hong and Pavlic [2021] J. Hong and T. P. Pavlic.  An insect-inspired randomly, weighted neural network with random fourier features for neuro-symbolic relational learning.  In _Proceedings of the 15th International Workshop on Neuro-Symbolic Learning and Reasoning (NeSy 20/21)_ , October 25–27, 2021.
  * Hong and Pavlic [2022b] J. Hong and T. P. Pavlic.  Learning to modulate random weights: Neuromodulation-inspired neural networks for efficient continual learning.  _arXiv preprint arXiv:2204.04297_ , 2022b.
  * LeCun [1998] Y. LeCun.  The MNIST database of handwritten digits.  _http://yann.lecun.com/exdb/mnist/_ , 1998.
  * Mei et al. [2022] J. Mei, E. Muller, and S. Ramaswamy.  Informing deep neural networks by multiscale principles of neuromodulatory systems.  _Trends Neurosci._ , 2022.
  * Miconi et al. [2020] T. Miconi, A. Rawal, J. Clune, and K. O. Stanley.  Backpropamine: training self-modifying neural networks with differentiable neuromodulated plasticity.  _arXiv:2002.10585_ , 2020.
  * Modi et al. [2020] M. N. Modi, Y. Shuai, and G. C. Turner.  The _drosophila_ mushroom body: from architecture to algorithm in a learning circuit.  _Annu. Rev. Neurosci._ , 43:465–484, 2020.
  * Nguyen et al. [2022] N. D. Nguyen, J. Huang, and D. Wang.  A deep manifold-regularized learning model for improving phenotype prediction from multi-modal data.  _Nature computational science_ , 2(1):38–46, 2022.
  * Tenenbaum et al. [2000] J. B. Tenenbaum, V. d. Silva, and J. C. Langford.  A global geometric framework for nonlinear dimensionality reduction.  _science_ , 290(5500):2319–2323, 2000.
  * Tomar and Rose [2013] V. S. Tomar and R. C. Rose.  Efficient manifold learning for speech recognition using locality sensitive hashing.  In _2013 IEEE international conference on acoustics, speech and signal processing_ , pages 6995–6999. IEEE, 2013.
  * Vecoven et al. [2020] N. Vecoven, D. Ernst, A. Wehenkel, and G. Drion.  Introducing neuromodulation in deep neural networks to learn adaptive behaviours.  _PloS ONE_ , 15(1):e0227922, 2020.
  * Wang et al. [2017] J. Wang, T. Zhang, N. Sebe, H. T. Shen, et al.  A survey on learning to hash.  _IEEE transactions on pattern analysis and machine intelligence_ , 40(4):769–790, 2017.
  * Zhong et al. [2020] Z. Zhong, L. Zheng, G. Kang, S. Li, and Y. Yang.  Random erasing data augmentation.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(07):13001–13008, 2020.  doi: 10.1609/aaai.v34i07.7000.


* [32] L. B. Anderson, J. Gray, A. Lukas, and E. Palti, “Heterotic Line Bundle Standard Models”, JHEP 06 (2012)113, arXiv:1202.1757 [hep-th].
  * [33] L. B. Anderson, A. Constantin, J. Gray, A. Lukas, and E. Palti, “A Comprehensive Scan for Heterotic SU(5) GUT models”, JHEP 01 (2014)047, arXiv:1307.4787 [hep-th].
  * [34] S. Groot Nibbelink, O. Loukas, F. Ruehle, and P. K. S. Vaudrevange, “Infinite number of MSSMs from heterotic line bundles?”, Phys. Rev. D 92 4, (2015)046002, arXiv:1506.00879 [hep-th].
  * [35] S. Groot Nibbelink, O. Loukas, and F. Ruehle, “(MS)SM-like models on smooth Calabi-Yau manifolds from all three heterotic string theories”, Fortsch. Phys. 63 (2015)609–632, arXiv:1507.07559 [hep-th].
  * [36] V. Braun, Y.-H. He, B. A. Ovrut, and T. Pantev, “Vector bundle extensions, sheaf cohomology, and the heterotic standard model”, Adv. Theor. Math. Phys. 10 4, (2006)525–589, arXiv:hep-th/0505041.
  * [37] V. Braun, Y.-H. He, and B. A. Ovrut, “Stability of the minimal heterotic standard model bundle”, JHEP 06 (2006)032, arXiv:hep-th/0602073.
  * [38] Z. Marshall, B. A. Ovrut, A. Purves, and S. Spinner, “Spontaneous $R$-Parity Breaking, Stop LSP Decays and the Neutrino Mass Hierarchy”, Phys. Lett. B 732 (2014)325–329, arXiv:1401.7989 [hep-ph].
  * [39] Z. Marshall, B. A. Ovrut, A. Purves, and S. Spinner, “LSP Squark Decays at the LHC and the Neutrino Mass Hierarchy”, Phys. Rev. D 90 1, (2014)015034, arXiv:1402.5434 [hep-ph].
  * [40] S. Dumitru, B. A. Ovrut, and A. Purves, “The $R$-parity Violating Decays of Charginos and Neutralinos in the B-L MSSM”, JHEP 02 (2019)124, arXiv:1810.11035 [hep-ph].
  * [41] S. Dumitru, B. A. Ovrut, and A. Purves, “$R$-parity Violating Decays of Wino Chargino and Wino Neutralino LSPs and NLSPs at the LHC”, JHEP 06 (2019)100, arXiv:1811.05581 [hep-ph].
  * [42] S. Dumitru, C. Herwig, and B. A. Ovrut, “$R$-parity Violating Decays of Bino Neutralino LSPs at the LHC”, JHEP 12 (2019)042, arXiv:1906.03174 [hep-ph].
  * [43] R. Deen, B. A. Ovrut, and A. Purves, “The minimal SUSY B $-$ L model: simultaneous Wilson lines and string thresholds”, JHEP 07 (2016)043, arXiv:1604.08588 [hep-ph].
  * [44] M. Ambroso and B. Ovrut, “The B-L/Electroweak Hierarchy in Heterotic String and M-Theory”, JHEP 10 (2009)011, arXiv:0904.4509 [hep-th].
  * [45] B. A. Ovrut, A. Purves, and S. Spinner, “Wilson Lines and a Canonical Basis of SU(4) Heterotic Standard Models”, JHEP 11 (2012)026, arXiv:1203.1325 [hep-th].
  * [46] B. A. Ovrut, A. Purves, and S. Spinner, “A statistical analysis of the minimal SUSY B–L theory”, Mod. Phys. Lett. A 30 18, (2015)1550085, arXiv:1412.6103 [hep-ph].
  * [47] B. A. Ovrut, A. Purves, and S. Spinner, “The minimal SUSY B–L model: from the unification scale to the LHC”, JHEP 06 (2015)182, arXiv:1503.01473 [hep-ph].
  * [48] V. Barger, P. Fileviez Perez, and S. Spinner, “Minimal gauged U(1)(B-L) model with spontaneous R-parity violation”, Phys. Rev. Lett. 102 (2009)181802, arXiv:0812.3661 [hep-ph].
  * [49] P. Fileviez Perez and S. Spinner, “Spontaneous R-Parity Breaking in SUSY Models”, Phys. Rev. D 80 (2009)015004, arXiv:0904.2213 [hep-ph].
  * [50] P. Fileviez Perez and S. Spinner, “The Minimal Theory for R-parity Violation at the LHC”, JHEP 04 (2012)118, arXiv:1201.5923 [hep-ph].
  * [51] G. Honecker and M. Trapletti, “Merging Heterotic Orbifolds and K3 Compactifications with Line Bundles”, JHEP 01 (2007)051, arXiv:hep-th/0612030.
  * [52] S. Groot Nibbelink, J. Held, F. Ruehle, M. Trapletti, and P. K. S. Vaudrevange, “Heterotic Z(6-II) MSSM Orbifolds in Blowup”, JHEP 03 (2009)005, arXiv:0901.3059 [hep-th].
  * [53] M. Blaszczyk, S. Groot Nibbelink, F. Ruehle, M. Trapletti, and P. K. S. Vaudrevange, “Heterotic MSSM on a Resolved Orbifold”, JHEP 09 (2010)065, arXiv:1007.0203 [hep-th].
  * [54] M. Blaszczyk, S. Groot Nibbelink, O. Loukas, and S. Ramos-Sanchez, “Non-supersymmetric heterotic model building”, JHEP 10 (2014)119, arXiv:1407.6362 [hep-th].
  * [55] M. Blaszczyk, S. Groot Nibbelink, O. Loukas, and F. Ruehle, “Calabi-Yau compactifications of non-supersymmetric heterotic string theory”, JHEP 10 (2015)166, arXiv:1507.06147 [hep-th].
  * [56] A. P. Braun, C. R. Brodie, and A. Lukas, “Heterotic Line Bundle Models on Elliptically Fibered Calabi-Yau Three-folds”, JHEP 04 (2018)087, arXiv:1706.07688 [hep-th].
  * [57] R. Deen, Y.-H. He, S.-J. Lee, and A. Lukas, “Machine Learning String Standard Models”, arXiv:2003.13339 [hep-th].
  * [58] M. Larfors and R. Schneider, “Explore and Exploit with Heterotic Line Bundle Models”, Fortsch. Phys. 68 5, (2020)2000034, arXiv:2003.04817 [hep-th].
  * [59] M. Larfors, D. Passaro, and R. Schneider, “Heterotic Line Bundle Models on Generalized Complete Intersection Calabi Yau Manifolds”, JHEP 05 (2021)105, arXiv:2010.09763 [hep-th].
  * [60] H. Otsuka and K. Takemoto, “Deep learning and k-means clustering in heterotic string vacua with line bundles”, JHEP 05 (2020)047, arXiv:2003.11880 [hep-th].
  * [61] V. Braun, Y.-H. He, and B. A. Ovrut, “Supersymmetric Hidden Sectors for Heterotic Standard Models”, JHEP 09 (2013)008, arXiv:1301.6767 [hep-th].
  * [62] A. Ashmore, S. Dumitru, and B. A. Ovrut, “Line Bundle Hidden Sectors for Strongly Coupled Heterotic Standard Models”, arXiv:2003.05455 [hep-th].
  * [63] S. Groot Nibbelin and F. Ruehle, “Line bundle embeddings for heterotic theories”, JHEP 04 (2016)186, arXiv:1601.00676 [hep-th].
  * [64] V. Braun, B. A. Ovrut, T. Pantev, and R. Reinbacher, “Elliptic Calabi-Yau threefolds with Z(3) x Z(3) Wilson lines”, JHEP 12 (2004)062, arXiv:hep-th/0410055.
  * [65] C. Schoen, “On fiber products of rational elliptic surfaces with section”, Mathematische Zeitschrift 197 2, (1988)177–200. http://eudml.org/doc/183725.
  * [66] A. Ashmore, S. Dumitru, and B. A. Ovrut, “Explicit Soft Supersymmetry Breaking in the Heterotic M-Theory $B-L$ MSSM”, arXiv:2012.11029 [hep-th].
  * [67] B. A. Ovrut, “Vacuum Constraints for Realistic Strongly Coupled Heterotic M-Theories”, Symmetry 10 12, (2018)723, arXiv:1811.08892 [hep-th].
  * [68] T. Banks and M. Dine, “Couplings and scales in strongly coupled heterotic string theory”, Nucl. Phys. B 479 (1996)173–196, arXiv:hep-th/9605136.
  * [69] R. Feger and T. W. Kephart, “LieART—A Mathematica application for Lie algebras and representation theory”, Comput. Phys. Commun. 192 (2015)166–195, arXiv:1206.6379 [math-ph].
  * [70] R. Feger, T. W. Kephart, and R. J. Saskowski, “LieART 2.0 – A Mathematica application for Lie Algebras and Representation Theory”, Comput. Phys. Commun. 257 (2020)107490, arXiv:1912.10969 [hep-th].
  * [71] S. Groot Nibbelink, M. Trapletti, and M. Walter, “Resolutions of C**n/Z(n) Orbifolds, their U(1) Bundles, and Applications to String Model Building”, JHEP 03 (2007)035, arXiv:hep-th/0701227.
  * [72] S. K. Donaldson, “Anti self-dual yang-mills connections over complex algebraic surfaces and stable vector bundles”, Proceedings of the London Mathematical Society s3-50 1, (1985)1–26. https://londmathsoc.onlinelibrary.wiley.com/doi/abs/10.1112/plms/s3-50.1.1.
  * [73] K. Uhlenbeck and S. T. Yau, “On the existence of hermitian-yang-mills connections in stable vector bundles”, Communications on Pure and Applied Mathematics 39 S1, (1986)S257–S293. https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.3160390714.
  * [74] L. B. Anderson, J. Gray, and B. A. Ovrut, “Transitions in the Web of Heterotic Vacua”, Fortsch. Phys. 59 (2011)327–371, arXiv:1012.3179 [hep-th].
  * [75] R. Donagi, A. Lukas, B. A. Ovrut, and D. Waldram, “Holomorphic vector bundles and nonperturbative vacua in M theory”, JHEP 06 (1999)034, arXiv:hep-th/9901009.
  * [76] M. Brandle, Aspects of branes in (heterotic) M-theory.  PhD thesis, Humboldt U., Berlin, 7, 2003.
  * [77] M. Dine, N. Seiberg, X. G. Wen, and E. Witten, “Nonperturbative Effects on the String World Sheet”, Nucl. Phys. B 278 (1986)769–789.
  * [78] M. Dine, N. Seiberg, and E. Witten, “Fayet-Iliopoulos Terms in String Theory”, Nucl. Phys. B 289 (1987)589–598.
  * [79] A. Lukas and K. S. Stelle, “Heterotic anomaly cancellation in five-dimensions”, JHEP 01 (2000)010, arXiv:hep-th/9911156.
  * [80] R. Blumenhagen, G. Honecker, and T. Weigand, “Loop-corrected compactifications of the heterotic string with line bundles”, JHEP 06 (2005)020, arXiv:hep-th/0504232.
  * [81] T. Weigand, “Compactifications of the heterotic string with unitary bundles”, Fortsch. Phys. 54 (2006)963–1077.
  * [82] L. B. Anderson, J. Gray, A. Lukas, and B. Ovrut, “Stability Walls in Heterotic Theories”, JHEP 09 (2009)026, arXiv:0905.1748 [hep-th].


$\mathbf{M}=\mathbf{U}^{\star}\Delta_{\mathbf{V}}^{\top}+\Delta_{\mathbf{U}}(\mathbf{V}^{\star})^{\top}$
for convenience, which can be thought of as arbitrary rank-$r$ matrices. Then
let $S$ be an $\epsilon$-net for all rank-$r$ matrices with Frobenius norm 1;
by Lemma 13 we have that $\log|S|\leq O(\max(t,d)r\log(\frac{1}{\epsilon}))$.
We set $\epsilon=\frac{1}{n^{3}}$ so $\log(\frac{1}{\epsilon})=O(\log n)$. Now
for any matrix $\mathbf{M}\in S$ we have that
$\frac{1}{n}\langle\mathbf{A}_{i},\mathbf{M}\rangle$ is a sub-gaussian random
variable with variance parameter at most $t\sigma^{2}\frac{\delta^{2}}{n}$.
Thus, using a sub-gaussian tail bound along with a union bound over the net
shows that uniformly over the $\mathbf{M}\in S$,

$\displaystyle|\frac{1}{n}\sum_{i=1}^{n}\langle\mathbf{M},\mathbf{A}_{i}\rangle
n_{i})|$ $\displaystyle\leq O\left(\sigma\delta\sqrt{\frac{t\max\\{t,d\\}r\log
n}{n}}\right),$

with probability at least $1-\frac{1}{\mathrm{poly}(d)}$. We now show how to
lift to the set of all $\mathbf{M}$. Note that with probability at least
$1-e^{-\Omega(n)}$ that $\|{\mathbf{n}}\|=O(\sqrt{t}\sigma\sqrt{n})$ by a sub-
gaussian tail bound (see for example Lemma 30). Let $\mathbf{M}$ be an
arbitrary element, and $\mathbf{M}^{\prime}$ its closest element in the cover;
then we have that
$\mathbf{z}_{i}=\langle\mathbf{A}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle\leq\frac{\delta}{n^{2}}$
using the precondition on $\mathbf{A}_{i}$. Combining and using a union bound
then shows that,

$\displaystyle|{\frac{1}{n}\sum_{i=1}^{n}n_{i}\langle\mathbf{A}_{i},\mathbf{M}\rangle}|\leq|{\frac{1}{n}\sum_{i=1}^{n}n_{i}\langle\mathbf{A}_{i},\mathbf{M}^{\prime}\rangle}|+|{\frac{1}{n}\sum_{i=1}^{n}n_{i}\langle\mathbf{A}_{i},\mathbf{M}-\mathbf{M}^{\prime}\rangle}|\leq$
$\displaystyle O\left(\sigma\delta\sqrt{\frac{t\max\\{t,d\\}r\log
n}{n}}\right)+\frac{\sqrt{t}\sigma\delta}{\sqrt{n}}\leq
O\left(\sigma\delta\sqrt{\frac{t\max(t,d)r\log n}{n}}\right).$

Rescaling and recalling the definition of $\delta$ gives the result. ∎

###  C.4 Task Diversity for the Landscape Analysis

Here we collect several useful results for interpreting the results of the
landscape analysis. Throughout this section we use the notation
$\mathbf{U}\in\mathbb{R}^{t\times r}$ and $\mathbf{V}\in\mathbb{R}^{d\times
r}$.

The first result allow us to convert a guarantee on error in Frobenius norm to
a guarantee in angular distance, assuming an appropriate diversity condition
on $\mathbf{U}$.

######  Lemma 16.

Suppose $\mathbf{V}$ and $\hat{\mathbf{V}}$ are orthonormal projection
matrices, that is $\mathbf{V}^{\top}\mathbf{V}=\mathbf{I}_{r}$, and
$\hat{\mathbf{V}}^{\top}\hat{\mathbf{V}}=\mathbf{I}_{r}$. Then, for any
$\epsilon>0$, if
$\|{\hat{\mathbf{U}}\hat{\mathbf{V}}^{\top}-\mathbf{U}\mathbf{V}^{\top}}\|_{\text{F}}^{2}\leq\epsilon$
for some $\hat{\mathbf{U}}$ and $\mathbf{U}$, then:

$\mathrm{dist}^{2}(\mathbf{V},\hat{\mathbf{V}})\leq\frac{\epsilon}{\nu t},$

where $\nu=\sigma_{r}(\mathbf{U}^{\top}\mathbf{U})/t$.

Here the distance function is the sine function of the principal angle; i.e.

$\mathrm{dist}(\mathbf{V},\hat{\mathbf{V}}):=\|{\mathbf{V}^{\top}\hat{\mathbf{V}}_{\perp}}\|,$

and $\nu=\sigma_{r}(\mathbf{U}^{\top}\mathbf{U})/t$ represents an analog of
the task diversity matrix.

###### Proof.

Define the function
$f(\tilde{\mathbf{U}})=\|{\tilde{\mathbf{U}}\hat{\mathbf{V}}^{\top}-\mathbf{U}\mathbf{V}^{\top}}\|_{\text{F}}^{2}$.
The precondition of the theorem states that there exists $\hat{\mathbf{U}}$ so
that
$\|{\hat{\mathbf{U}}\hat{\mathbf{V}}^{\top}-\mathbf{U}\mathbf{V}^{\top}}\|_{\text{F}}^{2}\leq\epsilon$.
This clearly implies the following:

$\min_{\tilde{\mathbf{U}}}f(\tilde{\mathbf{U}})\leq\epsilon.$ (21)

Setting the gradient $\mathrm{d}f/\mathrm{d}\tilde{\mathbf{U}}=0$, we have the
minimizer $\tilde{\mathbf{U}}^{\star}$ satisfies:

$(\tilde{\mathbf{U}}^{\star}\hat{\mathbf{V}}^{\top}-\mathbf{U}\mathbf{V}^{\top})\hat{\mathbf{V}}=0,$

which gives:

$\tilde{\mathbf{U}}^{\star}=\mathbf{U}\mathbf{V}^{\top}\hat{\mathbf{V}}.$

Plugging this back to Eq. (21) gives:

$\|{\mathbf{U}\mathbf{V}^{\top}(\hat{\mathbf{V}}\hat{\mathbf{V}}^{\top}-\mathbf{I})}\|_{\text{F}}^{2}\leq\epsilon.$

Finally, we have:

$\displaystyle\|{\mathbf{U}\mathbf{V}^{\top}(\hat{\mathbf{V}}\hat{\mathbf{V}}^{\top}-\mathbf{I})}\|_{\text{F}}^{2}=$
$\displaystyle\|{\mathbf{U}\mathbf{V}^{\top}\hat{\mathbf{V}}_{\perp}\hat{\mathbf{V}}_{\perp}^{\top}}\|_{\text{F}}^{2}=\|{\mathbf{U}\mathbf{V}^{\top}\hat{\mathbf{V}}_{\perp}}\|_{\text{F}}^{2}=\mathrm{tr}(\mathbf{U}^{\top}\mathbf{U}\mathbf{V}^{\top}\hat{\mathbf{V}}_{\perp}\hat{\mathbf{V}}_{\perp}^{\top}\mathbf{V})$
$\displaystyle\geq$
$\displaystyle\sigma_{r}(\mathbf{U}^{\top}\mathbf{U})\|{\mathbf{V}^{\top}\hat{\mathbf{V}}_{\perp}}\|_{\text{F}}^{2}\geq\sigma_{r}(\mathbf{U}^{\top}\mathbf{U})\|{\mathbf{V}^{\top}\hat{\mathbf{V}}_{\perp}}\|^{2}.$

The second last inequality follow since for any p.s.d. matrices $\mathbf{A}$
and $\mathbf{B}$, we have
$\mathrm{tr}(\mathbf{A}\mathbf{B})\geq\sigma_{\min}(\mathbf{A})\mathrm{tr}(\mathbf{B})$.
This concludes the proof. ∎

For the following let
$\mathbf{A}=(\bm{\alpha}_{1},\cdots,\bm{\alpha}_{t})^{\top}\in\mathbb{R}^{t\times
r}$ and denote the SVD of
$\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$. Next we remark that
our assumptions on task diversity and normalization implicit in the matrix
$\mathbf{A}$ are sufficient to actually imply an incoherence condition on
$\mathbf{U}$ (which is used in the matrix sensing/completion style analysis).

######  Lemma 17.

If $\mu=\frac{1}{r\sigma_{r}(\mathbf{A}^{\top}\mathbf{A}/t)}$ and
$\max_{i\in[t]}\|{\bm{\alpha}_{i}}\|^{2}\leq C$, then we have:

$\max_{i\in[t]}\|{\mathbf{e}_{i}^{\top}\mathbf{U}}\|^{2}\leq\frac{C\mu r}{t}.$

###### Proof.

Since $\max_{i\in[t]}\|{\bm{\alpha}_{i}}\|^{2}\leq C$, we have, for any
$i\in[t]$

$C\geq\|{\bm{\alpha}_{i}}\|^{2}=\|{\mathbf{e}_{i}^{\top}\mathbf{A}}\|^{2}=\|{\mathbf{e}_{i}^{\top}\mathbf{U}\mathbf{\Sigma}}\|^{2}\geq\|{\mathbf{e}_{i}^{\top}\mathbf{U}}\|^{2}\sigma_{\min}^{2}(\mathbf{\Sigma})=\|{\mathbf{e}_{i}^{\top}\mathbf{U}}\|^{2}\sigma_{r}(\mathbf{A}^{\top}\mathbf{A})=(t/\mu
r)\|{\mathbf{e}_{i}^{\top}\mathbf{U}}\|^{2},$

which finishes the proof. ∎

Note in the context of Assumption 2 the incoherence parameter corresponds to
the parameter $\bar{\kappa}\leq\kappa$ since under our normalization
$\mathrm{tr}(\mathbf{A}^{\top}\mathbf{A}/t)=\Theta(1)$. Further to quickly
verify the incoherence ball contains the true parameters it is important to
recall the scale difference $\mathbf{M}^{\star}$ and
$\mathbf{A}^{\top}\mathbf{A}/t$ by a factor of $\sqrt{t}$.

##  Appendix D Proofs for Section 4

Assuming we have obtained an estimate of the column space or feature set
$\hat{\mathbf{B}}$ for the initial set of tasks, such that
$\|{\hat{\mathbf{B}}^{\top}\mathbf{B}}\|\leq\delta$, we now analyze the
performance of the plug-in estimator (which explicitly uses the estimate
$\hat{\mathbf{B}}$ in lieu of the unknown $\mathbf{B}$) on a new task. Recall
we define the estimator for the new tasks by a projected linear regression
estimator: $\hat{\bm{\alpha}}=\mathop{\rm
argmin}_{\bm{\alpha}}\|{\mathbf{y}-\mathbf{X}\hat{\mathbf{B}}\bm{\alpha}}\|^{2}\implies\hat{\bm{\alpha}}=(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{y}$.

Analyzing the performance of this estimator requires first showing that the
low-dimensional empirical covariance and empirical correlation concentrate in
$\tilde{O}(r)$ samples and performing an error decomposition to compute the
bias resulting from using $\hat{\mathbf{B}}$ in lieu of $\mathbf{B}$ as the
feature representation. We measure the performance the estimator with respect
to its estimation error with respect to the underlying parameter
$\mathbf{B}\bm{\alpha}_{0}$; in particular, we use
$\|{\hat{\mathbf{B}}\hat{\bm{\alpha}}-\mathbf{B}\bm{\alpha}_{0}}\|^{2}$. Note
that our analysis can accommodate covariates $\mathbf{x}_{i}$ generated from
non-isotropic non-Gaussian distributions. In fact the only condition we
require on the design is that the covariates are sub-gaussian random vectors
in the following sense.

######  Assumption 3.

Each covariate vector $\mathbf{x}_{i}$ is mean-zero, satisfies
$\mathbb{E}[\mathbf{x}\mathbf{x}^{\top}]=\mathbf{\Sigma}$ such that
$\sigma_{\max}(\mathbf{\Sigma})\leq C_{\max}$ and
$\sigma_{\min}(\mathbf{\Sigma})\geq C_{\min}>0$ and is $\mathbf{\Sigma}$-sub-
gaussian, in the sense that
$\mathbb{E}[\exp(\mathbf{v}^{\top}\mathbf{x}_{i})]\leq\exp\left(\frac{\|\mathbf{\Sigma}^{1/2}\mathbf{v}\|^{2}}{2}\right)$.
Moreover, the additive noise $\epsilon_{i}$ is i.i.d. sub-gaussian with
variance parameter $1$ and is independent of $\mathbf{x}_{i}$.

In the context of the previous assumption we also define the conditioning
number as $C_{\text{cond}}=C_{\max}/C_{\min}$. Note that Assumption 1
immediately implies Assumption 3.

Throughout this section we will let $\hat{\mathbf{B}}$ and
$\hat{\mathbf{B}}_{\perp}$ be orthonormal projection matrices spanning
orthogonal subspaces which are rank $r$ and rank $d-r$ respectively—so that
$\text{range}{(}\hat{\mathbf{B}})\oplus\text{range}{(}\hat{\mathbf{B}}_{\perp})=\mathbb{R}^{d}$.

###### Proof of Theorem 4.

To begin we use the definition of

$\displaystyle\hat{\bm{\alpha}}=(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{y}=(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}\bm{\alpha}_{0}+(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\bm{\epsilon}$

to decompose the error as,

$\displaystyle(\hat{\mathbf{B}}\hat{\bm{\alpha}}-\mathbf{B}\bm{\alpha}_{0})=\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}\bm{\alpha}_{0}-\mathbf{B}\bm{\alpha}_{0}+\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\bm{\epsilon}.$

Now squaring both sides of the equation gives, so

$\displaystyle\|{\hat{\mathbf{B}}\hat{\bm{\alpha}}-\mathbf{B}\bm{\alpha}_{0}}\|^{2}\leq
2(\|{\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}\bm{\alpha}_{0}-\mathbf{B}\bm{\alpha}_{0}}\|^{2}+\|{\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\bm{\epsilon}}\|^{2}).$

The first bias term can be bounded by Lemma 18, while the the variance term
can be bounded by Lemma 19 . Combining the results and using a union bound
gives the result. ∎

We now present the lemmas which allow us to bound the variance terms in the
aforementioned error decomposition. For the following two results we also
track the conditioning dependence with respect $C_{\min}$ and $C_{\max}$. We
first control the term arising from the projection of the additive noise onto
the empirical covariance matrix.

######  Lemma 18.

Let the sequence of $n$ i.i.d. covariates $\mathbf{x}_{i}$ and $n$ i.i.d.
additive noise variables $\epsilon_{i}$ satisfy Assumption 3. Then if
$n\gtrsim C_{\text{cond}}^{2}r\log n$,

$\displaystyle\|{\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\bm{\epsilon}}\|^{2}\leq
O\left(\frac{r\log n}{C_{\min}n}\right),$

with probability at least $1-O(n^{-100})$.

###### Proof.

Since
$\|{\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\bm{\epsilon}}\|^{2}\leq\|{(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\bm{\epsilon}}\|^{2}$,
it suffices to bound the latter term. Consider
$\bm{\epsilon}^{\top}\underbrace{\frac{1}{n}\frac{\mathbf{X}\hat{\mathbf{B}}}{\sqrt{n}}(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})^{-2}\frac{\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}}{\sqrt{n}}}_{\mathbf{A}}\bm{\epsilon}$.
So applying the Hanson-Wright inequality [Vershynin, 2018, Theorem 6.2.1]
(conditionally on $\mathbf{X}$) to conclude that
$\mathbb{P}[|{\bm{\epsilon}^{\top}\mathbf{A}\bm{\epsilon}-\mathbf{E}[\bm{\epsilon}^{\top}\mathbf{A}\bm{\epsilon}]}|\geq
t]\leq
2\exp(-c\min(\frac{t^{2}}{\|{\mathbf{A}}\|_{F}^{2}},\frac{t}{\|{\mathbf{A}}\|}))$.
Hence
$\bm{\epsilon}^{\top}\mathbf{A}\bm{\epsilon}\leq\mathbf{E}[\bm{\epsilon}^{\top}\mathbf{A}\bm{\epsilon}]+O(\|{\mathbf{A}}\|_{F}\sqrt{\log(2/\delta_{1})})+O(\|{\mathbf{A}}\|\log(2/\delta_{1}))$
with probability at least $1-\delta_{1}$.

Now using cyclicity of the trace we have that
$\mathbb{E}[\bm{\epsilon}^{\top}\mathbf{A}\bm{\epsilon}]=\frac{1}{n}\mathrm{tr}[(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})^{-1}]$.
Similarly
$\|{\mathbf{A}}\|=\frac{1}{n}\|{(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})^{-1}}\|=\frac{1}{n}\|{(\mathbf{E}+\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}}\|$.
Applying Lemma 20 to the matrix $\mathbf{E}$ with $\delta=n^{-200}$ and
assuming $n\gtrsim C_{\text{cond}}^{2}r\log(1/\delta)\gtrsim
C_{\text{cond}}^{2}r\log n$ shows that
$\|{(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}\mathbf{E}}\|\leq\frac{1}{4}$.
Also note that on this event and this regime of sufficiently large $n$, this
concentration result shows that
$\sigma_{\min}(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})>C_{\min}/2$,
so the matrix is invertible. Hence an application of Lemma 25 shows that
$\|{\mathbf{A}}\|\leq\frac{1}{n}(\frac{1}{C_{\min}}\cdot(1+C_{\text{cond}}\sqrt{\frac{r\log
n}{n}}))\leq O(\frac{1}{C_{\min}n})$. Similarly since
$\frac{\mathbf{X}\hat{\mathbf{B}}}{\sqrt{n}}$ is rank $r$ and invertible on
this event, it follows
$\|{\mathbf{A}}\|_{\text{F}}\leq\sqrt{r}\|{\mathbf{A}}\|\leq
O(\frac{\sqrt{r}}{C_{\min}n})$ and that
$\frac{1}{n}\mathrm{tr}[(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})^{-1}]\leq\frac{r}{C_{\min}n}$.

Hence taking $\delta_{1}=n^{-200}$, and using the union bound, we conclude
that $\bm{\epsilon}^{\top}\mathbf{A}\bm{\epsilon}\leq\frac{1}{C_{\min}}\cdot
O(\frac{r}{n})+O(\frac{\sqrt{r\log n}}{n})+O(\frac{\log n}{n})\leq
O(\frac{r\log n}{C_{\min}n})$ with probability at least $1-O(n^{-100})$. ∎

We now control the error term which arises both from the variance in the
random design matrix $\mathbf{X}$ and the bias due to mismatch between
$\hat{\mathbf{B}}$ and $\mathbf{B}$.

######  Lemma 19.

Let the sequence of $n$ i.i.d. covariates $\mathbf{x}_{i}$ satisfy the design
assumptions in Assumption 3, and assume
$\sin(\hat{\mathbf{B}},\mathbf{B})\leq\delta\leq 1$. Then if $n\gtrsim
C_{\text{cond}}^{2}r\log n$,

$\displaystyle\|{\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}\bm{\alpha}_{0}-\mathbf{B}\bm{\alpha}_{0}}\|^{2}\leq
O(\|{\bm{\alpha}_{0}}\|^{2}C_{\text{cond}}^{2}\delta^{2}),$

with probability at least $1-O(n^{-100})$.

###### Proof.

To control this term we first insert a copy of the identity
$\mathbf{I}_{d}=\hat{\mathbf{B}}\hat{\mathbf{B}}^{\top}+\hat{\mathbf{B}}_{\perp}\hat{\mathbf{B}}_{\perp}^{\top}$
to allow the variance term in the design cancel appropriately in the span of
$\hat{\mathbf{B}}$; formally,

$\displaystyle\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}\bm{\alpha}_{0}-\mathbf{B}\bm{\alpha}_{0}=\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}(\hat{\mathbf{B}}\hat{\mathbf{B}}^{\top}+\hat{\mathbf{B}}_{\perp}\hat{\mathbf{B}}_{\perp}^{\top})\mathbf{B}\bm{\alpha}_{0}-\mathbf{B}\bm{\alpha}_{0}=$
$\displaystyle(\hat{\mathbf{B}}\hat{\mathbf{B}}^{\top}-\mathbf{I})\mathbf{B}\bm{\alpha}_{0}+\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}}_{\perp}\hat{\mathbf{B}}_{\perp}^{\top}\mathbf{B}\bm{\alpha}_{0}=$
$\displaystyle\hat{\mathbf{B}}_{\perp}\hat{\mathbf{B}}_{\perp}^{\top}\mathbf{B}\bm{\alpha}_{0}+\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}}_{\perp}\hat{\mathbf{B}}_{\perp}^{\top}\mathbf{B}\bm{\alpha}_{0}\implies$
$\displaystyle\|{\hat{\mathbf{B}}(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}\bm{\alpha}_{0}-\mathbf{B}\bm{\alpha}_{0}}\|^{2}\leq
2(\|{\bm{\alpha}_{0}}\|^{2}\delta^{2}+\|{(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}}_{\perp}}\|^{2}\delta^{2}\|{\bm{\alpha}_{0}}\|^{2}).$

We now turn to bounding the second error term,
$\|{(\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{B}}\|^{2}$.
Let
$\mathbf{E}_{1}=\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}}-\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}}$
and
$\mathbf{E}_{2}=\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\mathbf{B}-\hat{\mathbf{B}}\mathbf{\Sigma}\mathbf{B}$.
Applying Lemma 20 to the matrix $\mathbf{E}_{1}$ with $\delta=n^{-200}$ and
assuming $n\gtrsim C_{\text{cond}}^{2}r\log(1/\delta)\gtrsim
C_{\text{cond}}^{2}r\log n$ shows that
$\|{(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}\mathbf{E}_{1}}\|\leq\frac{1}{4}$
and $\|{\mathbf{E}_{1}}\|\leq O(C_{\max}\sqrt{\frac{r\log n}{n}})$ with
probability at least $1-O(n^{-100})$. A further application of Lemma 25 shows
that
$(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})^{-1}=(\mathbf{E}_{1}+\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}=(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}+\mathbf{F}$,
where
$\|{\mathbf{F}}\|\leq\frac{4}{3}\|{(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}}\|\|{\mathbf{E}_{1}(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}}\|$
on this event. Similarly, defining
$\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\mathbf{B}=\mathbf{E}_{2}+\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\mathbf{B}$
and applying Lemma 20 again but to the matrix $\mathbf{E}_{2}$ with
$\delta=n^{-200}$ and assuming $n\gtrsim
C_{\text{cond}}^{2}r\log(1/\delta)\gtrsim C_{\text{cond}}^{2}r\log n$,
guarantees that $\|{\mathbf{E}_{2}}\|\leq O(C_{\max}(\sqrt{\frac{r\log
n}{n}})$ with probability at least $1-O(n^{-100})$.

Hence on the intersection of these two events,

$\displaystyle\|{(\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\mathbf{B}}\|=\|{((\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}+\mathbf{F})(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\mathbf{B}+\mathbf{E}_{2})}\|\leq$
$\displaystyle\|{(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\mathbf{B}}\|+\|{(\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\hat{\mathbf{B}})^{-1}\mathbf{E}_{2}}\|+\|{\mathbf{F}\hat{\mathbf{B}}^{\top}\mathbf{\Sigma}\mathbf{B}}\|+\|{\mathbf{E}_{2}\mathbf{F}}\|\leq$
$\displaystyle C_{\text{cond}}+O(C_{\text{cond}}\sqrt{\frac{r\log
n}{n}})+O(C_{\text{cond}}^{2}\sqrt{\frac{r\log
n}{n}})+O(C_{\text{cond}}^{2}\frac{r\log n}{n})\leq$ $\displaystyle
C_{\text{cond}}+O(C_{\text{cond}}^{2}\sqrt{\frac{r\log
n}{n}})=O(C_{\text{cond}}),$

under the condition $n\gtrsim C_{\text{cond}}^{2}\frac{r\log n}{n}$. Taking a
union bound over the aforementioned events and combining terms gives the
result. ∎

Finally we present a concentration result for random matrices showing
concentration when the matrices are projected along two (potentially
different) subspaces.

######  Lemma 20.

Suppose a sequence of i.i.d. covariates $\mathbf{x}_{i}$ satisfy the design
assumptions in Assumption 3. Then, if $\mathbf{A}$ and $\mathbf{B}$ are both
rank $r$ orthonormal projection matrices,

$\displaystyle\|{(\mathbf{A}^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{n}\mathbf{B})-\mathbf{A}^{\top}\mathbf{\Sigma}\mathbf{B})}\|\leq
O(C_{\max}(\sqrt{\frac{r}{n}}+\frac{r}{n}+\sqrt{\frac{\log(1/\delta)}{n}}+\frac{\log(1/\delta)}{n})),$

with probability at least $1-\delta$.

###### Proof.

The result follows by a standard sub-exponential tail bound and covering
argument. First note that for any fixed
$\mathbf{u},\mathbf{v}\in\mathbb{S}^{r-1}$, we have that
$\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{x}_{i}$ and
$\mathbf{v}^{\top}\mathbf{B}^{\top}\mathbf{x}_{i}$ are both
$\text{sG}(\sqrt{C_{\max}})$. Hence for any fixed
$\mathbf{u},\mathbf{v}\in\mathbb{S}^{r-1}$,
$\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}\mathbf{B}\mathbf{v}-\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{\Sigma}\mathbf{B}\mathbf{v}$
is $\text{sE}(8C_{\max},8C_{\max})$.

Now, let $S$ denote a $\epsilon$-cover of $\mathbb{S}^{r-1}$ which has
cardinality at most $(\frac{3}{\epsilon})^{r}$ by a volume-covering argument.
Hence for $\epsilon=\frac{1}{5}$,

$\displaystyle\mathbb{P}[\sup_{\mathbf{u}\in S^{r-1},\mathbf{v}\in
S^{r-1}}\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}\mathbf{B}\mathbf{v}-\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{\Sigma}\mathbf{B}\mathbf{v}\geq
t]\leq 225^{r}\exp(-cn\min(t^{2}/C_{\max}^{2},t/C_{\max}))),$

using a union bound over the covers and a sub-exponential tail bound. Taking
$t=C\cdot
C_{\max}(\sqrt{\frac{r}{n}}+\frac{r}{n}+\sqrt{\frac{\log(1/\delta)}{n}}+\frac{\log(1/\delta)}{n})$
for sufficiently large $C$, shows that
$225^{r}\exp(-cn\min(t^{2}/C_{\max}^{2},t/C_{\max})))\leq\delta$. Finally a
standard Lipschitz continuity argument yields

$\displaystyle\sup_{\mathbf{u}\in\mathbb{S}^{r-1},\mathbf{v}\in\mathbb{S}^{r-1}}\mathbf{u}^{\top}\frac{1}{n}\sum_{i=1}^{n}\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}\mathbf{B}\mathbf{v}-\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{\Sigma}\mathbf{B}\mathbf{v}$
$\displaystyle\leq\frac{1}{1-3\epsilon}\sup_{\mathbf{u}\in
S^{r-1},\mathbf{v}\in
S^{r-1}}\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}\mathbf{B}\mathbf{v}-\mathbf{u}^{\top}\mathbf{A}^{\top}\mathbf{\Sigma}\mathbf{B}\mathbf{v},$

which gives the result. ∎

##  Appendix E Proofs for Section 5

We begin by presenting the proof of the main statistical lower bound for
recovering the feature matrix $\mathbf{B}$ and relevant auxiliary results.
Following this we provide relevant background on Grassmann manifolds.

As mentioned in the main text our main tool is to use is a non-standard
variant of the Fano method, along with suitable bounds on the cardinality of
the packing number and the distributional covering number, to obtain minimax
lower bound on the difficulty of estimating $\mathbf{B}$. We instantiate the
$f$-divergence based lower bound below (which we instantiate with
$\chi^{2}$-divergence). We restate this result for convenience.

######  Lemma 21.

[Guntuboyina, 2011, Theorem 4.1] For any increasing function
$\ell:[0,\infty)\to[0,\infty)$,

$\displaystyle\inf_{\hat{\theta}}\sup_{\theta\in\Theta}\mathbb{P}_{\theta}[\ell(\rho(\hat{\theta},\theta))\geq\ell(\eta/2)]\geq\sup_{\eta>0,\epsilon>0}\left\\{1-\left(\frac{1}{N(\eta)}+\sqrt{\frac{(1+\epsilon^{2})M_{C}(\epsilon,\Theta)}{N(\eta)}}\right)\right\\}.$

In the context of the previous result $N(\eta)$ denotes a lower bound on the
$\eta$-packing number of the metric space $(\Theta,\rho)$. Moreover,
$M_{C}(\epsilon,\Theta)$ is a positive real number for which there exists a
set $G$ with cardinality $\leq M_{C}(\epsilon,S)$ and probability measures
$Q_{\alpha}$, $\alpha\in G$ such that $\sup_{\theta\in S}\min_{\alpha\in
G}\chi^{2}(\mathbb{P}_{\theta},Q_{\alpha})\leq\epsilon^{2}$, where $\chi^{2}$
denotes the chi-squared divergence. In words, $M_{C}(\epsilon,S)$ is an upper
bound on the $\epsilon$-covering on the space
$\\{\mathbb{P}_{\theta}:\theta\in S\\}$ when distances are measured by the
square root of the $\chi^{2}$-divergence.

We obtain the other term capturing the difficulty of estimating
$\hat{\mathbf{B}}$ with respect to the task diversity minimum eigenvalue by
constructing a pair of feature matrices which are hard to distinguish for a
particular ill-conditioned task matrix $\mathbf{A}$. Using these two results
we provide the proof of our main lower bound.

###### Proof of Theorem 5.

For our present purposes we simply take $\ell(\cdot)$ to be the identity
function in our application of Lemma 21 as we obtain the second term in the
lower bound. Then by the duality between packing and covering numbers we have
that $\log N\geq\log M$ at the same scale (see for example Wainwright [2019,
Lemma 5.5]), so once again by Proposition 9 we have that $\log N(\eta)\geq
r(d-r)\log(\frac{c_{1}}{\eta})$. Then applying Lemma 22 we have that
$M_{C}(\epsilon,\Theta^{\prime})\leq(\frac{c_{2}n}{\log(1+\epsilon^{2})})^{r(d-r)/2}$.
For convenience we set $k=r(d-r)$ in the following. We now choose the pair
$\eta$ and $\epsilon$ appropriately in Lemma 21. The lower bound writes as,

$\displaystyle
1-\left(\frac{1}{N(\eta)}+\sqrt{\frac{(1+\epsilon^{2})M_{C}(\epsilon,\Theta)}{N(\eta)}}\right)\geq
1-\left((\frac{\eta}{c_{2}})^{k}+(\frac{\eta}{c_{2}})^{k/2}\cdot(c_{1}n)^{k/4}\frac{(1+\epsilon^{2})^{1/2}}{\log(1+\epsilon^{2})^{k/4}}\right),$

with the implicit constraint that
$\epsilon^{\prime}=\sqrt{\frac{2}{n}\log(1+\epsilon^{2})}<1$. A simple
calculus argument shows that
$\epsilon\to\sqrt{1+\epsilon^{2}}/(\log(1+\epsilon^{2}))^{k/4}$ is minimized
when $1+\epsilon^{2}=e^{k/2}$ (subject to
$\sqrt{\frac{1}{2n}\log(1+\epsilon^{2})}<1)$. This constraint can always be
ensured by taking $n>\frac{k}{4}$. We then have that the lower bound becomes

$\displaystyle
1-\left((\frac{\eta}{c_{2}})^{k}+(2e\frac{c_{1}}{c_{2}^{2}}\eta^{2}n/k)^{k/4}\right).$
(22)

We now take $\eta=C\sqrt{k/n}$ so the bound simplifies to
$1-\left((\frac{\eta}{c_{2}})^{k}+(2e\frac{c_{1}}{c_{2}^{2}}\eta^{2}n/k)^{k/4}\right)=1-\left((\frac{C}{c_{2}}\sqrt{k/n})^{k}+(2e\frac{c_{1}}{c_{2}^{2}}C^{2})^{k/4}\right)$.
By choosing $C$ to be sufficiently small and taking $n>k$ we ensure that
$1-\left((\frac{C}{c_{2}}\sqrt{k/n})^{k}+(2e\frac{c_{1}}{c_{2}^{2}}C^{2})^{k/4}\right)\geq\frac{99}{100}$.
Finally, under the condition $r\leq\frac{d}{2}$ we have that
$k\geq\frac{dr}{2}$. Combining with Lemma 21 gives the result for the second
term.

The first term is lower bounded using the LeCam two-point method in an
indepedent fashion as a consequence of Lemma 24. A union bound over the events
on which the lower bounds hold give the result. Note that a single choice of
$\mathbf{A}$ matrix can in fact be used for both lower bounds by simply opting
for the choice of $\mathbf{A}$ used in Lemma 24.

∎

To implement the lower bound we require an upper bound on the covering number
in the space of distributions of $\mathbb{P}_{\mathbf{B}}$. Throughout we use
standard properties of the $\chi^{2}$-divergence which can be found in
Tsybakov [2008, Section 2.4].

######  Lemma 22.

Suppose $n$ data points, $(\mathbf{x}_{i},y_{i})$, are generated from the
model in (1) with i.i.d. covariates
$\mathbf{x}_{i}\sim\mathcal{N}(0,\mathbf{I}_{d})$ and independent i.i.d. noise
$\epsilon_{i}\sim\mathcal{N}(0,1)$. Further, assume the task parameters
satisfy Assumption 2 with each task normalized to
$\|{\bm{\alpha}_{i}}\|=\frac{1}{2}$. Then if $r\leq\frac{d}{2}$,

$\displaystyle
M_{C}(\epsilon,\Theta^{\prime})\leq\left(\frac{cn}{\log(1+\epsilon^{2})}\right)^{r(d-r)/2},$

whenever $\sqrt{\frac{1}{2n}\log(1+\epsilon^{2})}<1$.

###### Proof.

We first upper bound the $\chi^{2}$ divergence between two data distributions
for two distinct $\mathbf{B}^{i}$ and $\mathbf{B}^{j}$. Now the joint
distribution over the observations for each measure
$\mathbb{P}_{\mathbf{B}^{i}}$ can be written as
$\mathbb{P}_{\mathbf{B}^{i}}\equiv\Pi_{k=1}^{t}p(\mathbf{X}_{k})p(\mathbf{y}|\mathbf{X}_{k},\mathbf{B}^{i},\bm{\alpha}_{j})$,
where $p(\mathbf{X}_{k})$ corresponds to the density of the Gaussian design
matrix, and $p(\mathbf{y}|\mathbf{X}_{k},\mathbf{B}^{i},\bm{\alpha}_{k})$ the
Gaussian conditionals of the observations $\mathbf{y}$. So using standard
properties of the $\chi^{2}$-divergence we find that,

$\displaystyle\chi^{2}(\mathbb{P}_{\mathbf{B}^{i}},\mathbb{P}_{\mathbf{B}^{j}})=\Pi_{k=1}^{t}(1+\mathbb{E}_{\mathbf{X}_{k}}[\chi^{2}(p(\mathbf{y}|\mathbf{X}_{k},\mathbf{B}^{i}),p(\mathbf{y}|\mathbf{X}_{k},\mathbf{B}^{j}))])-1=$
$\displaystyle\Pi_{k=1}^{t}\mathbb{E}_{\mathbf{X}_{k}}\left[\exp\left(\|{\mathbf{X}_{k}(\mathbf{B}^{i}\bm{\alpha}_{k}-\mathbf{B}^{j}\bm{\alpha}_{k})}\|^{2}\right)\right]-1.$

Now note that
$\|{(\mathbf{B}^{i}\bm{\alpha}_{k}-\mathbf{B}^{j}\bm{\alpha}_{k})}\|^{2}\leq
2\|{\bm{\alpha}_{k}}\|^{2}\sigma_{1}(\mathbf{I}_{r}-(\mathbf{B}^{i})^{\top}\mathbf{B}^{j})$.
Recognizing
$\sigma_{r}((\mathbf{B}^{i})^{\top}\mathbf{B}^{j})=\cos\theta_{1}(\mathbf{B}^{i},\mathbf{B}^{j})$
and $\|{(\mathbf{B}_{\perp}^{i})^{\top}\mathbf{B}^{j}}\|=\sin\theta_{1}$,
where $\theta_{1}$ is largest principal angle between the subspaces, we have
that
$1-\sigma_{r}((\mathbf{B}^{i})^{\top}\mathbf{B}^{j})=1-\sqrt{1-\|{(\mathbf{B}_{\perp}^{i})^{\top}\mathbf{B}^{j}}\|^{2}}\leq\|{(\mathbf{B}_{\perp}^{i})^{\top}\mathbf{B}^{j}}\|^{2}\leq
1$, using the elementary inequality $1-\sqrt{1-x^{2}}\leq x^{2}$ for $0\leq
x\leq 1$. Thus,
$\|{(\mathbf{B}^{i}\bm{\alpha}_{t}-\mathbf{B}^{j}\bm{\alpha}_{t})}\|^{2}\leq\frac{1}{2}$.

Now we use the identity that for
$\mathbf{x}\sim\mathcal{N}(0,\mathbf{I}_{d})$, and
$\|{\mathbf{v}}\|\leq\frac{1}{2}$ that
$\mathbb{E}_{\mathbf{x}}\exp((\mathbf{v}^{\top}\mathbf{x})^{2})=\mathbb{E}_{\mathbf{x}}[\exp(\|{\mathbf{v}}\|^{2}((\frac{\mathbf{v}^{\top}}{\|{\mathbf{v}}\|}\mathbf{x})^{2}-1)]\exp(\|{\mathbf{v}}\|^{2})=\frac{\exp(-\|{\mathbf{v}}\|^{2}))}{\sqrt{1-2\|{\mathbf{v}}\|^{2}}}\cdot\exp(\|{\mathbf{v}}\|^{2})\leq\exp(2\|{\mathbf{v}}\|^{4})\cdot\exp(\|{\mathbf{v}}\|^{2})\leq\exp(2\|{\mathbf{v}}\|^{2})$.
Hence combining the above two facts we obtain,

$\displaystyle\chi^{2}(\mathbb{P}_{\mathbf{B}^{i}},\mathbb{P}_{\mathbf{B}^{j}})\leq\exp\left(2\sum_{k=1}^{t}n_{t}\|{\mathbf{B}^{i}\bm{\alpha}_{k}-\mathbf{B}^{j}\bm{\alpha}_{k}}\|^{2}\right)-1\leq\exp(2n\cdot\sin^{2}\theta(\mathbf{B}^{i},\mathbf{B}^{j}))-1$

applying Lemma 23 in the inequality with a rescaling. Hence to ensure that
$\chi^{2}(\mathbb{P}_{\mathbf{B}^{i},\mathbf{B}^{j}})\leq\epsilon^{2}$ we take
$\mathbf{B}^{i}$ to be the closest element in a $\epsilon^{\prime}$ cover,
$S$, of $\text{Gr}_{r,d}(\mathbb{R})$ to $\mathbf{B}^{j}$. Since further
$\chi^{2}(\mathbb{P}_{\mathbf{B}^{i}},\mathbb{P}_{\mathbf{B}^{j}})\leq\exp(2n(\epsilon^{\prime})^{2})-1$
this is satisfied by taking
$\epsilon^{\prime}=\sqrt{\frac{1}{2n}\log(1+\epsilon^{2})}$ (with the
constraint we have $\epsilon^{\prime}<1$). Using Proposition 9, we then obtain
that,

$\displaystyle
M_{C}(\epsilon,\Theta^{\prime})\leq\left(\frac{cn}{\log(1+\epsilon^{2})}\right)^{r(d-r)/2},$

for a universal constant $c$. ∎

######  Lemma 23.

Let the task parameters $\bm{\alpha}_{i}$ each satisfy
$\|{\bm{\alpha}_{i}}\|=1$ with parameter
$\nu=\sigma_{r}(\frac{\mathbf{A}^{\top}\mathbf{A}}{t})>0$, and let
$\mathbf{B}^{i}$ and $\mathbf{B}^{j}$ be distinct, rank-$r$ orthonormal
feature matrices. Then,

$\displaystyle\sum_{k=1}^{t}n_{t}\|{\mathbf{B}^{i}\bm{\alpha}_{k}-\mathbf{B}^{j}\bm{\alpha}_{k}}\|^{2}\leq
n\sin^{2}\theta(\mathbf{B}^{i},\mathbf{B}^{j}).$

where $n=t\cdot n_{t}$.

###### Proof.

We can simplify the expression as follows,

$\displaystyle\sum_{k=1}^{t}n_{t}\|{\mathbf{B}^{i}\bm{\alpha}_{k}-\mathbf{B}^{j}\bm{\alpha}_{k}}\|^{2}=n\cdot\frac{1}{t}\sum_{j=1}^{t}2\cdot(\bm{\alpha}_{k}^{\top}\bm{\alpha}_{k}-\bm{\alpha}_{k}^{\top}(\mathbf{B}^{i})^{\top}\mathbf{B}^{j}\bm{\alpha}_{k})=n\cdot\mathrm{tr}\left(\left(\mathbf{I}_{r}-(\mathbf{B}^{i})^{\top}\mathbf{B}^{j}\right)\frac{\mathbf{A}^{\top}\mathbf{A}}{t}\right)$

using the fact that $\mathbf{B}^{\top}\mathbf{B}=\mathbf{I}_{r}$ for an
orthonormal feature matrix, the normalization of $\|{\bm{\alpha}_{k}}\|=1$,
and the cyclic property of the trace. Now use the fact that
$\mathrm{tr}\left(\left(\mathbf{I}_{r}-(\mathbf{B}^{i})^{\top}\mathbf{B}^{j}\right)\frac{\mathbf{A}^{\top}\mathbf{A}}{t}\right))\leq\sigma_{\max}(\mathbf{I}_{r}-(\mathbf{B}^{i})^{\top}\mathbf{B}^{j}))\cdot\mathrm{tr}(\frac{\mathbf{A}^{\top}\mathbf{A}}{t})\leq\frac{1}{4}(1-\cos\theta_{1})=\frac{1}{4}(1-\sqrt{1-\sin^{2}\theta_{1}})\leq\frac{1}{4}\sin^{2}\theta_{1}$,
using the elementary inequality $1-\sqrt{1-x^{2}}\leq x^{2}$ for $0\leq x\leq
1$. Note that $\mathrm{tr}(\mathbf{A}^{\top}\mathbf{A}/t)\leq\frac{1}{4}$
follows from the normalization of the $\bm{\alpha}_{j}$. ∎

We now provide brief background on $\text{Gr}_{r,d}(\mathbb{R})$ and establish
several pieces of notation relevant to the discussion. We denote the Grassmann
manifold, which consists of the the set of $r$-dimensional subspaces within
the underlying $d$-dimensional space, as $\text{Gr}_{r,d}(\mathbb{R})$.
Another way to define it is as the homogeneous space of the orthogonal group
$O(d)$ in the sense that,

$\displaystyle\text{Gr}_{r,d}(\mathbb{R})\cong O(d)/(O(r)\times O(d-r)),$

which defines its geometric structure. The underlying measure on the manifold
$G_{r,n}(\mathbb{R})$ is the associated, normalized invariant (or Haar)
measure.

Note that each orthonormal feature matrix $\mathbf{B}$, is contained in an
equivalence class (under orthogonal rotation) of an element in
$\text{Gr}_{r,d}(\mathbb{R})$. To define distance in
$\text{Gr}_{r,d}(\mathbb{R})$ we define the notion of a principal angle
between two subspaces $p$ and $q$. If $\mathbf{C}$ is an orthonormal matrix
whose columns form an orthonormal basis of $p$ and $\mathbf{D}$ is an
orthonormal matrix whose columns form an orthonormal basis of $q$, then the
singular values of the decomposition of
$\mathbf{C}^{\top}\mathbf{D}=\mathbf{U}\mathbf{D}\mathbf{V}^{\top}$ defines
the principal angles as follows:

$\displaystyle\mathbf{D}=\text{diag}(\cos\theta_{1},\cos\theta_{2},\ldots,\cos\theta_{k}),$

where $0\leq\theta_{k}\leq\ldots\leq\theta_{1}\leq\frac{\pi}{2}$. As shorthand
we let $\bm{\theta}=(\theta_{1},\theta_{2},\ldots,\theta_{k})$, and let $\sin$
and $\cos$ act element-wise on its components. The subspace angle distance
which is induced by $\ell_{\infty}$ norms on the vector $\sin\theta$. We refer
the reader to Pajor [1998] for geometric background on coding and
packing/covering bounds in the context of Grassmann manifolds relevant to our
discussion here.

In the following we let $M(\text{Gr}_{r,d}(\mathbb{R}),\sin\theta_{1},\eta)$
denote the $\eta$-covering number of $\text{Gr}_{r,d}(\mathbb{R})$ in the
subspace angle distance.

######  Proposition 9.

[Pajor, 1998, Adapted from Proposition 8] For any integers $1\leq
r\leq\frac{d}{2}$ and every $\epsilon>0$, we have that,

$\displaystyle r(d-r)\log(\frac{c_{1}}{\eta})\leq\log
M(\text{Gr}_{r,d}(\mathbb{R}),\sin\theta_{1},\eta)\leq
r(d-r)\log(\frac{c_{2}}{\eta}),$

for universal constants $c_{1},c_{2}>0$.

###### Proof.

Define for a linear operator $T$, $\sigma_{q}(T)=(\sum_{i\geq
1}|{s_{i}(T)}|^{p})^{1/p}$ for all $1\leq q\leq\infty$ where $s_{i}(T)$
denotes its $i$th singular value. Note that Proposition 8 in Pajor [1998]
states the result in the distance metric $d(E,F)=\sigma_{q}(P_{E}-P_{F})$
where $P_{E}$ and $P_{F}$ denotes the projection operator onto the subspace
$E$ and $F$ respectively, and $\sigma_{q}(T)=(\sum_{i\geq
1}|{s_{i}(T)}|^{p})^{1/p}$. However, as the computation in Proposition 6 of
Pajor [1998] establishes, we have that
$\sigma_{q}(P_{E}-P_{F})=(2\sum_{i=1}^{r}(1-\cos^{2}\theta_{i})^{q/2})^{1/q}$;
taking $q\to\infty$ implies $\sigma_{q}(P_{E}-P_{F})=\sin\theta_{1}$, and
hence directly translating the result gives the statement of the proposition.
∎

Finally, we include the proof of the lower bound which captures the dependence
on the task diversity parameter. The proof uses the LeCam two-point method
between two problem instances which are difficult to distinguish for a
particular, ill-conditioned task matrix.

######  Lemma 24.

Under the conditions of Theorem 5, for $n\geq\frac{1}{8\nu}$,

$\displaystyle\inf_{\hat{\mathbf{B}}}\sup_{\mathbf{B}\in\text{Gr}_{r,d}(\mathbb{R})}\sin\theta(\hat{\mathbf{B}},\mathbf{B})\geq\Omega\left(\sqrt{\frac{1}{\nu}}\sqrt{\frac{1}{n}}\right)$

with probability at least $\frac{3}{10}$.

###### Proof.

First consider an ill-conditioned task matrix where the sequence of
$\bm{\alpha}_{i}=\frac{1}{2}\mathbf{e}_{i}$ for $i\in[r-1]$ but then
$\bm{\alpha}_{r}=\frac{1}{2}(\sqrt{1-b^{2}}\mathbf{e}_{r-1}+b\mathbf{e}_{r})$
for $0<b<1$ where $\mathbf{e}_{i}$ are the standard basis in $\mathbb{R}^{r}$.
Now, consider two task models for two different subspaces $\mathbf{B}_{1}$ and
$\mathbf{B}_{2}$ which are distinct in a single direction. Namely we take
$\mathbf{B}_{1}=[\mathbf{e}_{1},\ldots,\mathbf{e}_{r}]$ and
$\mathbf{B}_{2}=[\mathbf{e}_{1},\ldots,\mathbf{e}_{r-1},\sqrt{1-a^{2}}\mathbf{e}_{r}+a\mathbf{e}_{r+1}]$,
for $0<a<1$, where $\mathbf{e}_{i}$ refer to the standard basis in
$\mathbb{R}^{d}$. Here we have
$\cos\theta_{1}=\|{\mathbf{B}_{2}^{\top}\mathbf{B}_{1}}\|=\sqrt{1-a^{2}}\implies\sin\theta_{1}=a$,
where $\theta_{1}$ refers to the largest principle angle between the two
subspaces.

Data is generated from the two linear models as,

$\displaystyle
y_{i}=\mathbf{x}_{i}^{\top}\mathbf{B}_{1}\bm{\alpha}_{j}+\epsilon_{i}\
i=1,\ldots,n$ $\displaystyle
y_{i}=\mathbf{x}_{i}^{\top}\mathbf{B}_{2}\bm{\alpha}_{j}+\epsilon_{i}\
i=1,\ldots,n$

with $n$ total samples generated evenly from each of $j$ in $[t]$ tasks
($n_{t}$ from each task) inducing two measures $\mathbb{P}_{1}$ and
$\mathbb{P}_{2}$ over their respective data. The LeCam two-point method (see
Wainwright [2019, Ch. 15] for example) shows that,

$\displaystyle\inf_{\hat{\mathbf{B}}}\sup_{\mathbf{B}\in\text{Gr}_{r,d}(\mathbb{R})}\mathbb{P}_{\mathbf{B}}[\sin\theta(\hat{\mathbf{B}},\mathbf{B})\geq
a]\geq\frac{1}{2}(1-\|{\mathbb{P}_{1}-\mathbb{P}_{2}}\|_{\text{TV}})$ (23)

for the $\mathbf{B}_{1}$ and $\mathbf{B}_{2}$ of our choosing as above.

We can now upper bound the total variation distance (via the Pinsker
inequality) similar to as in Lemma 22,

$\displaystyle\|{\mathbb{P}_{1}-\mathbb{P}_{2}}\|_{\text{TV}}^{2}\leq\frac{1}{2}\text{KL}(\mathbb{P}_{1}|\mathbb{P}_{2})=\frac{1}{4}\sum_{j=1}^{t}\sum_{i=1}^{n_{t}}\|{\mathbf{B}_{1}\bm{\alpha}_{j}-\mathbf{B}_{2}\bm{\alpha}_{j}}\|_{2}^{2}=\frac{1}{2}n_{t}\sum_{j=1}^{t}(\|{\bm{\alpha}_{j}}\|^{2}-\bm{\alpha}_{j}^{\top}\mathbf{B}_{1}^{\top}\mathbf{B}_{2}\bm{\alpha}_{j})=\frac{n}{2}\cdot(\frac{1}{4}-\mathrm{tr}(\mathbf{B}_{1}^{\top}\mathbf{B}_{2}\mathbf{C}))$
(24)

using cyclicity of the trace. Straightforward calculations show that given the
$\mathbf{A}$ matrix,
$\mathbf{C}=\frac{1}{4r}\begin{bmatrix}\mathbf{I}_{r-2}&\textbf{0}\\\
\textbf{0}&\mathbf{M}_{1}\end{bmatrix}$, where
$\mathbf{M}_{1}=\begin{bmatrix}2-b^{2}&b\sqrt{1-b^{2}}\\\
b\sqrt{1-b^{2}}&b^{2}\\\ \end{bmatrix}$. Similarly
$\mathbf{B}_{1}^{\top}\mathbf{B}_{2}=\begin{bmatrix}\mathbf{I}_{r-2}&\textbf{0}\\\
\textbf{0}&\mathbf{M}_{2}\end{bmatrix}$ where
$\mathbf{M}_{2}=\begin{bmatrix}1&0\\\ 0&\sqrt{1-a^{2}}\\\ \end{bmatrix}$.
Computing the trace term,

$\displaystyle\frac{1}{4}-\mathrm{tr}((\mathbf{B}_{1}^{\top}\mathbf{B}_{2})\mathbf{C})=\frac{1}{4}-\frac{1}{4}\left(\frac{r-2}{r}+\frac{1}{r}\cdot\mathrm{tr}\Big{(}\begin{bmatrix}1&0\\\
0&\sqrt{1-a^{2}}\\\ \end{bmatrix}\begin{bmatrix}2-b^{2}&b\sqrt{1-b^{2}}\\\
b\sqrt{1-b^{2}}&b^{2}\\\ \end{bmatrix}\Big{)}\right)=$ (25)
$\displaystyle\frac{1}{4r}b^{2}(1-\sqrt{1-a^{2}}).$ (26)

Hence,
$\|{\mathbb{P}_{1}-\mathbb{P}_{2}}\|_{\text{TV}}^{2}\leq\frac{n}{2r}b^{2}(1-\sqrt{1-a^{2}})\leq\frac{n}{2r}b^{2}a^{2}$.
Combining with the LeCam two-point lemma shows that,

$\displaystyle\inf_{\hat{\mathbf{B}}}\sup_{\mathbf{B}\in\text{Gr}_{r,d}(\mathbb{R})}\mathbb{P}_{\mathbf{B}}[\sin\theta(\hat{\mathbf{B}},\mathbf{B})\geq
a]\geq\frac{1}{2}(1-\|{\mathbb{P}_{1}-\mathbb{P}_{2}}\|_{\text{TV}})\geq\frac{1}{2}(1-\sqrt{\frac{n}{2r}}ba).$
(27)

Taking $a=\frac{1}{2}\sqrt{\frac{r}{n}}\frac{1}{b}<1$ suffices to ensure the
lower bound with probability at least $\frac{3}{10}$. This induces the
constraint $\frac{1}{2}\sqrt{\frac{r}{n}}\frac{1}{b}<1\implies
n>\frac{r}{4b^{2}}$. As a last remark note that the $\mathbf{C}$ matrix has
maximum and minimum eigenvalues $\frac{1}{4r}(1+\sqrt{1-b^{2}})$ and
$\frac{1}{4r}(1-\sqrt{1-b^{2}})$. So
$\nu=\frac{1}{4r}(1-\sqrt{1-b^{2}})\implies\sqrt{2}\sqrt{r\nu}\leq b\leq
2\sqrt{2}\sqrt{r\nu}$ for $0<b<1$ using the inequality $x^{2}/2\leq
1-\sqrt{1-x^{2}}\leq x^{2}$. Hence it follows
$a\geq\frac{1}{8}\frac{1}{\sqrt{\nu}}\sqrt{\frac{1}{n}}$ as well. Similarly
the constraint can reduce too $n>\frac{1}{8\nu}$. ∎

##  Appendix F Auxiliary Results

Here we collect several auxiliary results. We begin by stating a simple matrix
perturbation result.

######  Lemma 25.

Let $\mathbf{A}$ be a positive-definite matrix and $\mathbf{E}$ another matrix
which satisfies $\|{\mathbf{E}\mathbf{A}^{-1}}\|\leq\frac{1}{4}$, then,

$\displaystyle(\mathbf{A}+\mathbf{E})^{-1}=\mathbf{A}^{-1}+\mathbf{F},$

where
$\|{\mathbf{F}}\|\leq\frac{4}{3}\|{\mathbf{A}^{-1}}\|\|{\mathbf{E}\mathbf{A}^{-1}}\|$.

###### Proof.

$\displaystyle(\mathbf{A}+\mathbf{E})^{-1}=\mathbf{A}^{-1}(\mathbf{I}+\mathbf{E}\mathbf{A}^{-1})^{-1}.$

Under the condition, $\|{\mathbf{E}\mathbf{A}^{-1}}\|\leq\frac{1}{4}$,
$\mathbf{I}+\mathbf{E}\mathbf{A}^{-1}$ is invertible and has a convergent
power series expansion so

$\displaystyle(\mathbf{A}+\mathbf{E})^{-1}=\mathbf{A}^{-1}(\mathbf{I}-\mathbf{E}\mathbf{A}^{-1}+(\mathbf{E}\mathbf{A}^{-1})^{2}+\ldots)=\mathbf{A}^{-1}+\mathbf{F},$

where
$\mathbf{F}=\mathbf{A}^{-1}\cdot(-\mathbf{E}\mathbf{A}^{-1}+(\mathbf{E}\mathbf{A}^{-1})^{2}+\ldots)$.
Moreover,

$\displaystyle\|{\mathbf{F}}\|\leq\|{\mathbf{A}^{-1}}\|(\|{\mathbf{E}\mathbf{A}^{-1}}\|+\|{\mathbf{E}\mathbf{A}^{-1}}\|^{2}+\ldots)\leq\frac{\|{\mathbf{A}^{-1}}\|\|{\mathbf{E}\mathbf{A}^{-1}}\|}{1-\|{\mathbf{E}\mathbf{A}^{-1}}\|}\leq\frac{4}{3}\|{\mathbf{A}^{-1}}\|\|{\mathbf{E}\mathbf{A}^{-1}}\|.$

∎

We now present several results related to the concentration of measure.

######  Lemma 26.

Let $x,y$ be mean-zero random variables that are both sub-gaussian with
parameters $\kappa_{1}$ and $\kappa_{2}$ respectively. Then
$z=xy-\mathbb{E}[xy]\sim\text{sE}(8\kappa_{1}\kappa_{2},8\kappa_{1}\kappa_{2})$.

The proof is a standard argument and omitted. Next we prove a matrix
concentration result for the individually rescaled covariance matrices of
i.i.d. random variables. The proof uses a standard covering argument.

######  Lemma 27.

Let $\mathbf{X}\in\mathbb{R}^{n\times d}$ be a random matrix with rows
$a_{i}\mathbf{x}_{i}$, where $\mathbf{x}_{i}$ are i.i.d. random vectors
satisfying the design conditions in Assumption 1. Then,

$\displaystyle\|{\frac{1}{n}\mathbf{X}^{\top}\mathbf{X}-\frac{1}{n}\sum_{i=1}^{n}a_{i}^{2}\mathbf{\Sigma}}\|\leq\|{\mathbf{\Sigma}}\|K^{2}\max(\delta,\delta^{2})\quad\text{
for }\delta=C(\sqrt{d/n}+t/\sqrt{n}),$

with probability at least $1-2\exp(-t^{2})$. Here $C$ denotes a universal
constant and $K=\max_{i}|{a_{i}}|$.

###### Proof.

First note that we bring all the vectors to isotropic position by rotating so
that
$\|{\frac{1}{n}\mathbf{X}^{\top}\mathbf{X}-\frac{1}{n}\sum_{i=1}^{n}a_{i}^{2}\mathbf{\Sigma}}\|\leq\|{\mathbf{\Sigma}}\|\|{\frac{1}{n}\sum_{i=1}^{n}a_{i}^{2}(\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i})(\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i})^{\top}-\frac{1}{n}\sum_{i=1}^{n}a_{i}^{2}\mathbf{I}_{d}}\|$.
Now by definition for any fixed
$\mathbf{v}\in\mathbb{R}^{d}:\|{\mathbf{v}}\|=1$, each
$\mathbf{v}^{\top}\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}$ is $\text{sG}(1)$ and
hence
$a_{i}^{2}(\mathbf{v}^{\top}\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i})(\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i})^{\top}\mathbf{v}$
is $\text{sE}(8a_{i}^{2},8a_{i}^{2})$ by Lemma 26. For the latter quantity
Vershynin [2018, Theorem 4.6.1, Eq. 4.22] proves the result when $a_{i}=1$
using a standard covering argument along with a sub-exponential tail bound. A
close inspection of the proof of Vershynin [2018, Theorem 4.6.1, Eq. 4.22]
shows that the aforementioned analogous statement holds when the sequence of
random vectors is scaled by $a_{i}$. ∎

We now include two useful results on operator norm bounds of higher-order
matrices. The results only require the condition of $O(1)$-L4-L2
hypercontractivity (which is directly implied by Assumption 1—see for example
the sub-gaussian moment bounds in Wainwright [2019, Theorem 2.6]). Formally,
we say a random vector $\mathbf{x}$ is $L$-L4-L2 hypercontractive if
$\mathbb{E}[\langle\mathbf{v},\mathbf{x}\rangle^{4}]\leq
L^{2}(\mathbb{E}[\langle\mathbf{v},\mathbf{x}\rangle^{2}])^{2}$ for all unit
vectors $\mathbf{v}$. Also note that if $\mathbf{x}$ is hypercontractive this
immediately implies that $P_{\mathbf{V}}\mathbf{x}$ is also hypercontractive
with the same constant where $P_{\mathbf{V}}$ is an orthogonal projection
operator.

######  Lemma 28.

Let $\mathbf{x}$ be a mean-zero random vector from a distribution that is
$L$-L$4$-L$2$ hypercontractive with covariance $\mathbf{\Sigma}$ and let
$P_{\mathbf{V}}$ be an orthogonal projection operator onto a rank-$r$
subspace. Then

$\|{\mathbb{E}[\|{\mathbf{x}}\|^{2}\mathbf{x}\mathbf{x}^{\top}]}\|\leq
L\mathrm{tr}(\mathbf{\Sigma})\|{\mathbf{\Sigma}}\|;\quad\|{\mathbb{E}[\|{P_{\mathbf{V}}\mathbf{x}}\|^{2}\mathbf{x}\mathbf{x}^{\top}]}\|\leq
Lr\|{\mathbf{\Sigma}}\|^{2};\quad\|{\mathbb{E}[\|{\mathbf{x}}\|^{2}P_{\mathbf{V}}\mathbf{x}(P_{\mathbf{V}}\mathbf{x})^{\top}}\|\leq
L\mathrm{tr}(\mathbf{\Sigma})\|{\mathbf{\Sigma}}\|.$ (28)

###### Proof of Lemma 28.

We introduce a vector $\mathbf{v}$ with $\|{\mathbf{v}}\|\leq 1$. Then,

$\displaystyle\mathbb{E}[\langle\mathbf{v},\|{\mathbf{x}}\|^{2}\mathbf{x}\mathbf{x}^{\top}\mathbf{v}\rangle]=\mathbb{E}[\|{\mathbf{x}}\|^{2}\langle\mathbf{v},\mathbf{x}\rangle^{2}]\leq(\mathbb{E}[\|{\mathbf{x}}\|^{4}])^{1/2}(\mathbb{E}[\langle\mathbf{v},\mathbf{x}\rangle^{4}])^{1/2},$
(29)

by the Cauchy-Schwarz inequality. For the first term we have
$(\mathbb{E}[\|{\mathbf{x}}\|^{4}])^{1/2}\leq\sqrt{L}\mathrm{tr}(\mathbf{\Sigma})$
by Lemma 29. For the second term once again using $L$-L$4$-L$2$
hypercontractivity we have
$(\mathbb{E}[\langle\mathbf{v},\mathbf{x}\rangle^{4}])^{1/2}\leq\sqrt{L}\mathbb{E}[\langle\mathbf{v},\mathbf{x}\rangle]^{2}\leq\sqrt{L}\|{\mathbf{\Sigma}}\|$.
Maximizing over $\mathbf{v}$ gives the result. The remaining statements follow
using an identical calculation and appealing to Lemma 29. ∎

######  Lemma 29.

Let $\mathbf{x}$ be a mean-zero random vector from a distribution that is
$L$-L$4$-L$2$ hypercontractive with covariance $\mathbf{\Sigma}$ and let
$P_{\mathbf{V}}$ be an orthogonal projection operator onto a rank-$r$
subspace. Then

$\displaystyle\mathbb{E}[\|{\mathbf{x}}\|^{4}]\leq
L(\mathrm{tr}\mathbf{\Sigma})^{2};\quad\mathbb{E}[\|{P_{\mathbf{V}}\mathbf{x}}\|^{2}]\leq
r\|{\mathbf{\Sigma}}\|;\quad\mathbb{E}[\|{P_{\mathbf{V}}\mathbf{x}}\|^{2}\|{\mathbf{x}}\|^{2}]\leq
Lr\|{\mathbf{\Sigma}}\|(\mathrm{tr}\mathbf{\Sigma});\|{\mathbb{E}[\|{P_{\mathbf{V}}\mathbf{x}}\|^{4}]}\|\leq
L\|{\mathbf{\Sigma}}\|^{2}r^{2}.$

###### Proof of Lemma 29.

A short computation using the Cauchy-Schwarz inequality and L4-L2 equivalence
shows that,

$\displaystyle\mathbb{E}[\|{\mathbf{x}}\|^{4}]=\mathbb{E}[(\sum_{i=1}^{d}\langle\mathbf{x},\mathbf{e}_{i}\rangle^{2})^{2}]=\mathbb{E}[\sum_{a,b}\langle\mathbf{x},\mathbf{e}_{a}\rangle^{2}\langle\mathbf{x},\mathbf{e}_{b}\rangle^{2}]\leq\sum_{a,b}(\mathbb{E}[\langle\mathbf{x},\mathbf{e}_{a}\rangle^{4}]\mathbb{E}[\langle\mathbf{x},\mathbf{e}_{b}\rangle^{4}])^{1/2}\leq$
$\displaystyle
L\sum_{a,b}\mathbb{E}[\langle\mathbf{x},\mathbf{e}_{a}\rangle^{2}]\mathbb{E}[\langle\mathbf{x},\mathbf{e}_{b}\rangle^{2}]\leq
L(\mathrm{tr}(\mathbf{\Sigma}))^{2}.$

The second statement follows since
$\mathbb{E}[\|{P_{\mathbf{V}}\mathbf{x}}\|^{2}]=\mathbb{E}[\mathrm{tr}(P_{\mathbf{V}}\mathbf{x}\mathbf{x}^{\top}]=\mathrm{tr}(P_{\mathbf{V}}\mathbf{\Sigma})\leq
r\|{\mathbf{\Sigma}}\|$ where the last line follows by the von Neumann trace
inequality and the fact $P_{\mathbf{V}}$ is a projection operator. The final
statements follow by combining the previous calculations. ∎

######  Lemma 30.

Let $\mathbf{x}$ be a random vector in $d$ dimensions from a distribution
satisfying Assumption 3. Then, we have:

$\|{\mathbf{x}}\|\leq O\left(\sqrt{C_{\max}}(\sqrt{d}+\sqrt{\log
1/\delta})\right),$

with probability at least $1-\delta$.

###### Proof of Lemma 30.

Note that by rotating the vectors into isotropic position
$\mathbf{y}=\mathbf{\Sigma}^{-1/2}\mathbf{x}$ is $\mathbf{I}_{d}$-subgaussian
in the sense of Assumption 3. Since
$\|{\mathbf{x}}\|\leq\sqrt{C_{\max}}\|{\mathbf{\Sigma}^{-1/2}\mathbf{x}}\|$ it
suffices to bound the norm of $\|{\mathbf{y}}\|$. First note that
$\mathbb{E}[\|{\mathbf{y}}\|]\leq\sqrt{\mathbb{E}[\|{\mathbf{y}}\|^{2}]}=\sqrt{d}$.
Now, take an $1/2$-net over the unit sphere, $G$; by a standard covering
argument the number of elements in $G$ can be upper bounded by $6^{d}$
[Wainwright, 2019, Chapter 5]. By definition of
$\mathbf{I}_{d}$-subgaussianity, for any $\mathbf{v}\in G$, we have that,

$\mathbb{P}[|{\langle\mathbf{v},\mathbf{y}\rangle}|\geq t]\leq
2\exp\left(-\frac{t^{2}}{2}\right)$

By a standard continuity argument it follows that $\max_{\mathbf{v}\in
S^{d-1}}|{\langle\mathbf{v},\mathbf{y}\rangle}|\leq 3\max_{\mathbf{v}\in
G}|{\langle\mathbf{v},\mathbf{y}\rangle}|$. So by a union bound,
$\mathbb{P}[\|{\mathbf{y}}\|\geq t]\leq(12)^{d}\exp(-t^{2}/20)$. Therefore, by
taking $t=C(\sqrt{d}+\sqrt{\log 1/\delta})$ for large-enough constant $C$ we
can ensure that $(12)^{d}\exp(-t^{2}/20)\leq\delta$. Re-arranging gives the
conclusion. ∎

Finally, we prove a truncated version of the matrix Bernstein inequality we
can apply to matrices that are unbounded in spectral norm. This is our primary
technical tool used to show concentration of the higher-order moments used in
the algorithm to recover the feature matrix $\mathbf{B}$.

######  Lemma 31.

Consider a truncation level $R>0$. If $Z_{i}$ is a sequence of symmetric
independent random matrices and if
$Z_{i}^{\prime}=Z_{i}\mathbbm{1}{[}\|{Z_{i}}\|\leq R]$, then

$\displaystyle\mathbb{P}[\|{\frac{1}{n}\sum_{i=1}^{n}Z_{i}-\mathbb{E}[Z_{i}]}\|\geq
t]\leq\mathbb{P}[\|{\frac{1}{n}\sum_{i=1}^{n}Z^{\prime}_{i}-\mathbb{E}[Z_{i}^{\prime}]}\|\geq
t-\Delta]+n\mathbb{P}[\|{Z_{i}}\|\geq R],$

where $\Delta\geq\|{\mathbb{E}[Z_{i}]-\mathbb{E}[Z_{i}^{\prime}]}\|$. Further,
for $t\geq\Delta$, we have that,

$\displaystyle\mathbb{P}[\|{\frac{1}{n}\sum_{i=1}^{n}Z^{\prime}_{i}-\mathbb{E}[Z_{i}^{\prime}]}\|\geq
t-\Delta]\leq
2d\exp\left(\frac{n^{2}(t-\Delta)^{2}}{\sigma^{2}+2Rn(t-\Delta)/3}\right),$

where
$\sigma^{2}=\|{\sum_{i=1}^{n}\mathbb{E}[(Z_{i}^{\prime}-\mathbb{E}[Z_{i}^{\prime}])^{2}]}\|\leq\|{\sum_{i=1}^{n}\mathbb{E}[Z_{i}^{2}]}\|$.

###### Proof.

The first statement follows by splitting on the event $\\{\|{Z_{i}}\|\leq
R:\forall i\in[n]\\}$, along with a union bound, and an application of the
triangle inequality to the first term. The second is simply a restatement of
the matrix Bernstein inequality in Tropp [2012] along with the almost sure
upper bound
$\|{Z_{i}^{\prime}-\mathbb{E}[Z_{i}^{\prime}]}\|\leq\|{Z_{i}^{\prime}}\|+\|{\mathbb{E}[Z_{i}^{\prime}]}\|\leq
R+\mathbb{E}[\|{Z_{i}^{\prime}}\|]\leq 2R$. The final bound on the matrix
variance follows from the facts that for the p.s.d. matrix
$\sum_{i=1}^{n}\mathbb{E}[(Z_{i}^{\prime}-\mathbb{E}[Z_{i}^{\prime}])^{2}]\preceq\sum_{i=1}^{n}\mathbb{E}[(Z_{i}^{\prime})^{2}]$
and that $(Z_{i}^{\prime})^{2}\preceq Z_{i}^{2}$. ∎

##  Appendix G Experimental Details

In our experiments we did find that gradient descent was able to decrease the
loss in 4, but the algorithm was slow to converge. In practice, we found using
the L-BFGS algorithm required no tuning and optimized the loss in 4 to high-
precision in far fewer iterations [Liu and Nocedal, 1989]. Hence we used this
first-order method throughout our experiments as our optimization routine. Our
implementation is in Python, and we leveraged the autograd package to compute
derivatives of the objective in 4, and the package Ray to parallelize our
experiments [Maclaurin et al., 2015, Moritz et al., 2018]. Each experiment is
averaged over 30 repetitions with error bars representing $\pm 1$ standard
deviation over the repetitions. All the experiments herein were run on
computer with 48 cores and 256 GB of RAM.

Note that after optimizing 4 directly using a first-order method in the
variable ($\mathbf{U},\mathbf{V}$), we can extract an estimate
$\hat{\mathbf{B}}$ of $\mathbf{B}$ by computing the column space of
$\mathbf{V}$ (for example using the SVD of $\mathbf{V}$ or applying the Gram-
Schmidt algorithm).