# Emergence of hierarchical modes from deep learning

Chan Li1 Haiping Huang1,2<EMAIL_ADDRESS>1PMI Lab, School of
Physics, Sun Yat-sen University, Guangzhou 510275, People’s Republic of China
2Guangdong Provincial Key Laboratory of Magnetoelectric Physics and Devices,
Sun Yat-sen University, Guangzhou 510275, People’s Republic of China

###### Abstract

Large-scale deep neural networks consume expensive training costs, but the
training results in less-interpretable weight matrices constructing the
networks. Here, we propose a mode decomposition learning that can interpret
the weight matrices as a hierarchy of latent modes. These modes are akin to
patterns in physics studies of memory networks, but the least number of modes
increases only logarithmically with the network width, and becomes even a
constant when the width further grows. The mode decomposition learning not
only saves a significant large amount of training costs, but also explains the
network performance with the leading modes, displaying a striking piecewise
power-law behavior. The modes specify a progressively compact latent space
across the network hierarchy, making a more disentangled subspaces compared to
standard training. Our mode decomposition learning is also studied in an
analytic on-line learning setting, which reveals multi-stage of learning
dynamics with a continuous specialization of hidden nodes. Therefore, the
proposed mode decomposition learning points to a cheap and interpretable route
towards the magical deep learning.

Introduction.— Deep neural networks are dominant tools with a broad range of
applications in not only image and language processing, but also scientific
researches Goodfellow et al. (2016); Carleo et al. (2019). These networks are
parameterized by a huge amount of trainable weight matrices, thereby consuming
an expensive training cost. However, these weight matrices are hard to
interpret, and thus mechanisms underlying the macroscopic performance of the
networks remain a big mystery in theoretical studies of neural networks Huang
(2022); Roberts et al. (2022).

To save the computational cost, previous studies of deep networks applied
singular value decomposition to the weight matrices Jaderberg et al. (2014);
Yang et al. (2020); Giambagli et al. (2021); Chicchi et al. (2021). This
decomposition requires the orthogonality condition for the singular vectors
and positive singular values. The training also involves a carefully-designed
structure for the trainable decomposition scheme Giambagli et al. (2021);
Chicchi et al. (2021). These constraints and designs make the training process
complicated, and thus a concise physics interpretation is still lacking. In
addition, previous studies of recurrent memory networks showed that the
network weight can be decomposed into separate random orthogonal patterns with
corresponding importance scores Jiang et al. (2021); Zhou et al. (2021).
Inspired by these studies, we conjecture that the learning in deep networks is
shaped by a hierarchy of latent modes, which are not necessarily orthogonal,
and the weight matrix can be expressed by these modes.

The mode decomposition learning (MDL) leads to a progressively compact latent
mode space across the network hierarchy, and meanwhile the subspaces
corresponding to different types of input are strongly disentangled,
facilitating discrimination. The least number of latent modes achieving the
comparable performance with the costly standard methods grows only
logarithmically with the network width and even could be a constant, thereby
reducing significantly the training cost. The mode spectrum exhibits an
intriguing piecewise power-law behavior. In particular, these properties do
not depend on details of the training setting. Therefore, our proposed MDL
calls for a rethinking of conventional weight-based deep learning through the
lens of cheap and interpretable mode-based learning.

Model.— To show the effectiveness of the MDL scheme, we train a deep network
to implement a classification task of handwritten digits mni . The deep
network has $L$ layers ($L-2$ hidden layers) with $N_{l}$ neurons in the
$l$-th layer. The weight value of the connection from the neuron $i$ at the
upstream layer $l$ to the neuron $j$ at the downstream layer $l+1$ is
specified by $w_{ij}^{l}$. The activation of the neuron $j$ at the downstream
layer $h_{j}^{l+1}=f(z_{j}^{l+1})=\max(0,z_{j}^{l+1})$, where the pre-
activation $z_{j}^{l+1}=\sum_{i}w_{ij}^{l}h_{i}^{l}$. For the output layer,
the softmax function $h_{k}=e^{z_{k}}/\sum_{i}e^{z_{i}}$ is chosen to specify
the probability over all classes of the input images. The cross entropy
$\mathcal{C}=-\sum_{i}\hat{h}_{i}\ln h_{i}$ is used as the cost function for
the supervised learning, and $\hat{h}_{i}$ is the target label (one-hot form).
After training (the cross entropy is repeatedly averaged over mini-batches of
training examples), we evaluate the generalization performance of the network
on an unseen test dataset.

Single weight values are not interpretable. According to our hypothesis,
latent patterns would emerge from training in each layer. We call these
patterns hierarchical modes for deep learning. Therefore, the relationship
between the modes and weight values is expressed by the following mode
decomposition,

$\mathbf{w}^{l}=\bm{\hat{\xi}}^{l}\bm{\Sigma}^{l}(\bm{{\xi}}^{l+1})^{\operatorname{T}},$
(1)

where there are $p^{l}$ upstream modes
$\bm{\hat{\xi}}^{l}\in\mathbb{R}^{N_{l}\times p^{l}}$, and the same number of
downstream modes $\bm{\xi}^{l+1}\in\mathbb{R}^{N_{l+1}\times p^{l}}$. The
importance of each pair of adjacent modes is specified by the diagonal of the
importance matrix $\bm{\Sigma}^{l}\in\mathbb{R}^{p^{l}\times p^{l}}$. These
modes may not be orthogonal with each other, and the importance score can take
a real value. This setting allows for more degrees of freedom for learning
features of input-output mappings. We will detail their geometric and physical
interpretations below.

Figure 1:  A simple illustration of the mode decomposition learning. (a) A
deep neural network of three layers, including one hidden layer with three
hidden nodes, for a classification task of non-linearly separable data. The
weight matrix
$w_{ij}^{l}=\sum_{\alpha=1}^{p}\hat{\xi}_{i,\alpha}^{l}\Sigma_{\alpha}^{l}{\xi}_{j,\alpha}^{l}$,
where $p=3$. The distribution of input data is modeled as a Gaussian mixture
(see the main text) from which samples are assigned to labels $t=\pm 1$ based
on the corresponding mixture component. The training performance is measured
by the mean-squared-error loss function
$\ell_{\operatorname{MSE}}(y,t)=\|y-t\|^{2}/2$. (b) The representation of
hidden neurons $\mathbf{h}$ plotted in the 3D space, displaying the geometric
separation. (c) The successive mappings from input sample $\mathbf{x}$ (grey)
to $(\bm{\hat{\xi}}^{1})^{\operatorname{T}}\mathbf{x}$ (dark red), followed by
$\bm{\Sigma}^{1}(\bm{\hat{\xi}}^{1})^{\operatorname{T}}\mathbf{x}$ (green),
and finally
$\bm{\xi}^{2}\bm{\Sigma}^{1}(\bm{\hat{\xi}}^{1})^{\operatorname{T}}\mathbf{x}$
(blue).

A geometric interpretation of Eq. (1) in a simple learning task is shown in
Fig. 1. We use a three-layer network with three hidden neurons. The input data
is sampled from a four-component Gaussian mixture Fischer et al. (2022),

$\mathbb{P}(\mathbf{x},t)=P(t)\sum_{\pm}P_{\pm}\mathcal{N}\left(\mathbf{x}|\mu_{x}^{t,\pm},\Sigma_{x}^{t,\pm}\right),$
(2)

where $\mathcal{N}\left(\mathbf{x}|\mu_{x}^{t,\pm},\Sigma_{x}^{t,\pm}\right)$
denotes a Gaussian distribution with mean $\mu_{x}^{t,\pm}$ and covariances
$\Sigma_{x}^{t,\pm}$, and $P(t)=P_{\pm}=\frac{1}{2}$. For the label $t=+1$,
$\mu_{x}^{t=+1,\pm}=\pm(0.5,0.5)^{\operatorname{T}}$, while for $t=-1$,
$\mu_{x}^{t=-1,\pm}=\pm(-0.5,0.5)^{\operatorname{T}}$. Covariances are
isotropic throughout with $\Sigma_{x}^{t,\pm}=0.05\mathds{1}$. The input
samples $\mathbf{x}\in\mathbb{R}^{2}$ are first projected to the input pattern
space spanned by $(\hat{\bm{\xi}}^{1})^{\operatorname{T}}_{i}$ ($i=1,2,3$).
Then all three directions of this projection get expanded or contracted via
$\bm{\Sigma}^{1}(\bm{\hat{\xi}}^{1})^{\operatorname{T}}\mathbf{x}$. Finally
the geometrically modified representation is re-mapped to the downstream
representation space of a higher dimensionality, as
$\bm{{\xi}}^{2}\bm{\Sigma}^{1}(\bm{\hat{\xi}}^{1})^{\operatorname{T}}\mathbf{x}$
[Fig. 1 (c)]. The non-linearity of the transfer function is then applied to
the last linear transformation, leading to the geometric separation [Fig. 1
(b)]. We conclude that the MDL provides rich angles to look at the geometric
transformation of the input information along the hierarchy of deep networks.

Rather than the conventional weight values in standard backpropagation (BP)
algorithms Goodfellow et al. (2016), the trainable parameters are latent
patterns in the MDL. The training is implemented by stochastic gradient
descent in the mode space
$\bm{\theta}^{l}=(\bm{\hat{\xi}}^{l},\bm{\Sigma}^{l},\bm{{\xi}}^{l+1})$ SM ,

$\displaystyle\Delta\xi^{l+1}_{j\alpha}$
$\displaystyle\equiv-\eta\frac{\partial\mathcal{L}}{\partial\xi^{l+1}_{j\alpha}}=-\eta\mathcal{K}_{j}^{l+1}\Sigma_{\alpha}^{l}\sum_{i}\hat{\xi}_{i\alpha}^{l}h_{i}^{l},$
(3) $\displaystyle\Delta\Sigma_{\alpha}^{l}$
$\displaystyle\equiv-\eta\frac{\partial\mathcal{L}}{\partial\Sigma^{l}_{\alpha}}=-\eta\sum_{j}\mathcal{K}_{j}^{l+1}\xi_{j\alpha}^{l+1}\sum_{i}\hat{\xi}_{i\alpha}^{l}h_{i}^{l},$
$\displaystyle\Delta\hat{\xi}^{l}_{i\alpha}$
$\displaystyle\equiv-\eta\frac{\partial\mathcal{L}}{\partial\hat{\xi}_{i\alpha}^{l}}=-\eta\Sigma_{\alpha}^{l}h_{i}^{l}\sum_{j}\mathcal{K}_{j}^{l+1}\xi_{j\alpha}^{l+1},$

where $\mathcal{L}$ denotes the cost function (e.g., cross-entropy or mean-
squared error) over a mini-batch of training data, $\eta$ denotes the learning
rate, and $\mathcal{K}_{j}^{l+1}\equiv\partial\mathcal{L}/\partial
z_{j}^{l+1}$ denotes the error term, which could back-propagate from the top
layer where $\mathcal{K}_{j}^{L}=-\hat{h}_{j}^{L}\left(1-h_{j}^{L}\right)$ for
$\mathcal{L}=\mathcal{C}$ (cross entropy). Based on the chain rule, the error
backpropagation equation can be derived as
$\mathcal{K}_{i}^{l}=\sum_{j}\mathcal{K}_{j}^{l+1}\sum_{\alpha}\xi_{i\alpha}^{l+1}\Sigma_{\alpha}^{l}\hat{\xi}_{j\alpha}^{l}f^{\prime}(z_{i}^{l})$
SM . To ensure the pre-activation is independent of the upstream-layer width,
we take the initialization scheme that
$[\bm{{\xi}}^{l+1}\bm{{\Sigma}}^{l}(\hat{\bm{\xi}}^{l})^{\operatorname{T}}]_{ij}\sim\mathcal{O}(\frac{1}{\sqrt{N_{l}}})$
Jiang et al. (2021). To avoid the ambiguity of choosing patterns (e.g., scaled
by a factor), we impose an identical regularization with strength $10^{-4}$
for all trainable parameters. However, our result does not change
qualitatively with the specific values of regularization SM .

We remark that for each hidden layer, there exist two types of pattern
($\bm{\xi}^{l}\neq\hat{\bm{\xi}}^{l}$). Equation (3) is used to learn these
patterns. We call this case 1L2P. If we assume
$\bm{\xi}^{l}=\hat{\bm{\xi}}^{l}$, the training can be further simplified as
in SM , and we call this case 1L1P. The nature of this mode-based-computation
can be understood as an expanded linear-nonlinear layered computation, as
$f(z_{j}^{l+1})=f(\sum_{\alpha}c_{\alpha j}\kappa_{\alpha})$ where the linear
field $\kappa_{\alpha}=\sum_{i}\hat{\xi}_{i\alpha}^{l}h_{i}^{l}$ and the
equivalent weight $c_{\alpha j}=\xi_{j\alpha}^{l+1}\Sigma_{\alpha}^{l}$.
Therefore, the number of modes acts as the linear-layer width. We leave a
systematic exploration of this linear-nonlinear structure by statistical
mechanics in forthcoming works.

On-line learning dynamics in a shallow network.— The MDL can be analytically
understood in an on-line learning setting, where we consider one-hidden-layer
architecture. The on-line learning can be considered as a special case of the
above mini-batch learning (i.e., the batch size is set to one, and the sample
is visited by the learning only once). The training dataset consists of $n$
pairs $\left\\{\mathbf{x}^{\nu},y^{\nu}\right\\}_{\nu=1}^{n}$. Each training
example is independently sampled from a probability distribution
$\mathbb{P}(\mathbf{x},y)=\mathbb{P}(y|\mathbf{x})\mathbb{P}(\mathbf{x})$,
where $\mathbb{P}(\mathbf{x})$ is a standard Gaussian distribution, and the
scalar label ${y}^{\nu}$ is generated by the neural network of $k$ hidden
neurons, (i.e., teacher, indicated by the symbol $*$ below). Given an input
$\mathbf{x}^{\nu}\in\mathbb{R}^{d}$, the corresponding label is created by

$y^{\nu}=\frac{1}{k}\sum_{r=1}^{k}\sigma\left(\frac{[\bm{{\xi}}^{*}\bm{{\Sigma}^{*}}(\bm{\hat{\xi}^{*}})^{\operatorname{T}}]_{r}\mathbf{x}^{\nu}}{\sqrt{d}}\right)=\frac{1}{k}\sum_{r=1}^{k}\sigma\left(\lambda_{r}^{*\nu}\right),$
(4)

where
$[\bm{\xi}^{*}\bm{\Sigma}^{*}(\hat{\bm{\xi}}^{*})^{\operatorname{T}}]_{r}$
denotes the $r$-th row of the matrix
$\bm{\xi}^{*}\bm{\Sigma}^{*}(\hat{\bm{\xi}}^{*})^{\operatorname{T}}$, and
$\lambda_{r}^{*\nu}=[\bm{{\xi}}^{*}\bm{{\Sigma}^{*}}(\bm{\hat{\xi}^{*}})^{\operatorname{T}}]_{r}\mathbf{x}^{\nu}/{\sqrt{d}}$
represents the $r$-th element of the teacher local field vector
$\bm{\lambda}^{*\nu}\in\mathbb{R}^{k}$. The teacher network is quenched as
$[\bm{{\xi}}^{*}\bm{{\Sigma}^{*}}(\bm{\hat{\xi}^{*}})^{\operatorname{T}}]_{ij}\sim\mathcal{O}(1)$.
Here, we focus on the non-linear transfer function
$\sigma(x)=\operatorname{erf}(x/\sqrt{2})$. In addition, we train the other
shallow network called the student network, by minimizing the loss function
$\mathcal{L}(y,\hat{f}(\mathbf{x},\Theta))$ over the training data (labels are
given by the teacher network), where $\Theta$ denotes the trainable
parameters. The student’s prediction for a fresh sample $\mathbf{x}$ is given
by

$\hat{f}(\mathbf{x},\bm{\hat{\xi}},\bm{{\Sigma}},\bm{{\xi}})=\frac{1}{m}\sum_{r=1}^{m}\sigma\left(\frac{[\bm{{\xi}}\bm{{\Sigma}}(\bm{\hat{\xi}})^{\operatorname{T}}]_{r}\mathbf{x}}{\sqrt{d}}\right)=\frac{1}{m}\sum_{r=1}^{m}\sigma\left(\lambda_{r}\right),$
(5)

where $\lambda_{r}$ denotes the $r$-th component of the student local field
$\bm{\lambda}=\bm{{\xi}}\bm{{\Sigma}}(\bm{\hat{\xi}})^{\top}\mathbf{x}$, and
the student has $m$ hidden neurons. The student is supplied with data samples
in sequence (one sample each time step). We next use $\nu$ to indicate the
time step as well.

The mean-squared-error can be evaluated as

$\ell_{\mathrm{MSE}}(\bm{\Omega})=\frac{1}{2}\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\bm{\Omega}\right)}\left[\left(\hat{f}(\bm{\lambda})-f\left(\bm{\lambda}^{*}\right)\right)^{2}\right],$
(6)

where $f(\cdot)$ indicates the teacher’s output, and we have replaced the
expectation $\mathbb{E}_{\mathbf{x},y\sim\mathbb{P}(\mathbf{x},y)}[\cdot]$ by
$\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\bm{\Omega}\right)}[\cdot]$, because of the central-limit theorem and the
i.i.d. setting we consider Biehl and Schwarze (1995); Saad and Solla (1995);
Goldt et al. (2019). The covariance of the local field
$\bm{\Omega}^{\nu}\in\mathbb{R}^{(k+m)\times(k+m)}$ can be specified as
follows,

$\bm{\Omega}^{\nu}\equiv\left[\begin{array}[]{cc}\mathbf{Q}^{\nu}&\mathbf{M}^{\nu}\\\
(\mathbf{M}^{\nu})^{\operatorname{T}}&\mathbf{P}\end{array}\right],$ (7)

where
$\mathbf{Q}^{\nu}\equiv\mathbb{E}_{\mathbf{x},y\sim\mathbb{P}(\mathbf{x},y)}\left[\bm{\lambda}^{\nu}(\bm{\lambda}^{\nu})^{\operatorname{T}}\right]$,
$\mathbf{M}^{\nu}\equiv\mathbb{E}_{\mathbf{x},y\sim\mathbb{P}(\mathbf{x},y)}\left[\bm{\lambda}^{\nu}(\bm{\lambda}^{*\nu})^{\operatorname{T}}\right]$,
and
$\mathbf{P}^{\nu}\equiv\mathbb{E}_{\mathbf{x},y\sim\mathbb{P}(\mathbf{x},y)}\left[\bm{\lambda}^{*\nu}(\bm{\lambda}^{*\nu})^{\operatorname{T}}\right]$.
By definition, $\mathbf{P}$ is fixed, while $\mathbf{Q}^{\nu}$ and
$\mathbf{M}^{\nu}$ evolve according to the gradient updates, following a set
of deterministic ordinary differential equations (ODEs) as the input dimension
$d\to\infty$ SM . These matrices are exactly the order parameters in physics.
For simplicity, we consider $\bm{\xi}=\bm{\xi}^{*}$ and
$\bm{\Sigma}=\bm{\Sigma}^{*}$, i.e., only the upstream patterns are learned.

Figure 2: Test performance and mode hierarchy of MDL in deep neural networks.
(a) Training trajectories of a four-layer network, indicated by
$784$-$100$-$100$-$10$, where each number indicates the corresponding layer
width. The number of modes $p^{l}=p$ for layer $l$, where $l=1,\ldots,L$.
$p=20$, or $p=30$. Networks are trained on the full MNIST dataset ($6\times
10^{4}$ images) and tested on an unseen dataset containing $10^{4}$ images.
The fluctuation is computed over five independent runs. (b) Testing accuracy
versus $p$ (the number of modes is the same for all layers). The same
architecture as (a) is used. The error bar characterizes the fluctuation
across five independently trained networks, and each marker denotes the
average result. The least number of modes is indicated by the dash-dot line.
(c) The performance changes with the network width. The inset shows the least
number of modes versus the layer width $N$ (in the logarithmic scale). The
network architecture is given by $784$-$N$-$N$-$10$. The dash-dot line in the
inset separates the piecewise logarithmic increase ($\propto\ln N$) regions.
The result is obtained from five independent runs. (d) The averaged Euclidean
distance (dispersion) from the pattern-cloud center
$(\frac{1}{p}\sum_{\alpha}\bm{\xi}^{l}_{\alpha})$ as a function of layer
index. The network architecture is specified by $784$-$100$-$100$-$100$-$10$
($p=30$). (e-f) Subspace overlap (principal angle) versus layer. The overlap
is averaged with five independent runs, and seven-layer networks with hidden-
layer width $100$ are trained ($p=30$).

Results.— MDL can reach a similar test accuracy with that of BP performed in
the weight space, when $p$ is sufficiently large [Fig. 2 (a)]. The
computational cost of the BP scales with $N_{l}^{2}$. In contrast, MDL works
in the mode space, requiring a training cost of only the order of $pN_{l}$.
Note that $p$ is much smaller than $N_{l}$ (or
$\lim_{N_{l}\to\infty}p^{l}/N_{l}=0$), and our MDL does not need any
additional training constraints (compared to other matrix factorization
algorithms SM ). Remarkably, when $p=30$, the performance of MDL already
matches that of BP [Fig. 2 (b)], but only utilizes $40\%$ of the full sets of
parameters that are consumed by the BP. In fact, each hidden layer can have
two different types of latent pattern (1L2P) due to the mode decomposition.
But if we assume that $\bm{{\xi}}^{l}=\bm{\hat{\xi}}^{l}$, i.e., each layer
share a single type of pattern (1L1P), we can further reduce the computational
cost by an amount of $\sum_{l}p^{l}N_{l}$, without sacrificing the test
accuracy [Fig. 2 (b)]. Varying the network width, we reveal a logarithmic
increase of the least number of modes [Fig. 2 (c)], which is a novel property
of deep learning in the mode space, in stark contrast to a linear number of
memory patterns in previous studies Jiang et al. (2021). When the network
width further grows, the least number can even become a constant. We argue
that this manifests three separated phases of poor-good-saturated performance
with increasing layer width (see Fig. S9 in SM ).

Figure 3: The robustness properties of well-trained four-layer MDL models with
the architecture $784$-$100$-$100$-$10$. The case of 1L1P is considered with
$p=70$ in the hidden layers. (a) Effects of removing modes through two
protocols: removing modes with weak measure $\tau$ first (solid line) and
removing modes randomly (dashed line). The fluctuation is computed over ten
independent runs. (b) The rescaled $\ell_{2}$ norms $\gamma\|\xi\|_{2}$,
$\gamma\|\hat{\xi}\|_{2}$ and the absolute values of $\Sigma$ versus their
rank (in descending order) in the hidden layers, where
$\gamma=\sum_{\alpha}|{\Sigma_{\alpha}}|/\sum_{\alpha}(\|{\xi}_{\alpha}\|_{2}+\|{\hat{\xi}}_{\alpha}\|_{2})$.
The inset shows a log-log plot of the $\tau$ measure, displaying a piecewise
power-law behavior. The error bar is computed over five independent runs. The
marked percentage indicates the generalization accuracy after removing the
corresponding side of modes.

To see how the latent patterns are transformed in geometry along the network
hierarchy, we first calculate the center of the pattern space. Then the
Euclidean distance from this center to each pattern is analyzed. We find that
the pattern space becomes progressively compact when going to deep layers
[Fig. 2 (d)]. To further characterize the geometric details, we define the
subspace spanned by the principal eigenvectors of the layer neural responses
to one type of inputs. Then the subspace overlap is calculated as the cosine
of the principal angle between two subspaces corresponding to two types of
inputs Bjoerck and Golub (1973); SM . We find that the hidden-layer
representation becomes more disentangled with layer in comparison with BP
[Fig. 2 (e,f)]. MDL shows great computational benefits of representation
disentanglement, thereby facilitating discrimination. A slight increase of the
overlap is observed for deeper layers, which is caused by the saturation of
the test performance (see more analyses in SM ).

Compared to other matrix factorization methods, MDL has no additional
constraints for the modes and importance scores, therefore being flexible for
feature extraction. We find that the interlayer patterns are less orthogonal
than the intralayer ones. The geometric transformation carried out by these
latent pattern matrices is not strictly a rotation for which the $\ell_{2}$
norm is preserved. This flexibility may be the key to make our method better
than other matrix factorization methods in both training cost and learning
performance (see details in SM ).

Figure 4: Mean-squared error dynamics in terms of $t=\frac{\nu}{d}$, where
$\nu$ denotes the on-line sample index, and $d$ is the input dimension. The
teacher and student networks share the same number of hidden neurons
($m=k=8$). Markers represent results of the simulation, while the solid lines
denote the theoretical predictions from solving the mean-field ODEs. The
number of modes $p^{*}=p=\alpha\ln d$ ($\alpha$ denotes the mode load here).
(a) Fixed $\alpha=1$. (b) Fixed $d=100$. The color deepens as $\alpha$ or $d$
increases. The insets display the evolving $\mathbf{M}$ matrix for $d=30$ and
$\alpha=1.0$, respectively.

We next ask whether some modes are more important than the others. Therefore,
we rank the modes according to the measure
$\tau_{\alpha}=\gamma\|\xi_{\alpha}\|_{2}+\gamma\|\hat{\xi}_{\alpha}\|_{2}+|\Sigma_{\alpha}|$,
where
$\gamma=\sum_{\alpha}|{\Sigma_{\alpha}}|/\sum_{\alpha}(\|{\xi}_{\alpha}\|_{2}+\|{\hat{\xi}}_{\alpha}\|_{2})$
to make comparable the magnitudes of the pattern and importance
($\bm{\Sigma}$) score. Removing modes with weak values of $\tau$ first yields
much higher accuracy than the random removal protocol [Fig. 3 (a)], suggesting
the existence of leading modes. Moreover, deeper layers are more robust.
Figure 3 (b) shows the measure as a function of rank in descending order,
which can be approximately captured by piecewise power-law behavior (a
transition point at the rank $10$). Ranking with only the importance scores
yields similar behavior SM . A small exponent is observed for the leading
measures, while the remaining measures bear a large exponent, thereby
revealing the coding hierarchy of latent modes in the deep networks. This
intriguing behavior does not change with the regularization strength or the
hidden-layer width SM .

Finally, the on-line mean-squared error dynamics of our model can be predicted
perfectly in a teacher-student setting. The number of modes strongly affects
the shape of the learning dynamics, and a large mode load can make the
plateaus disappear (Fig. 4). Moreover, during learning, the alignment between
receptive fields of the student’s hidden nodes and the teacher’s ones
continuously emerge, which is called the specialization transition Schwarze
(1993); Goldt et al. (2019).

Conclusion.— In this Letter, we propose a mode decomposition learning that
works in the mode space rather than the conventional weight space. This
learning scheme has three-fold technical and conceptual advances. First, the
learning can achieve the comparable performance with standard methods, with a
significant reduction of training costs. We also find that the least number of
modes grows only logarithmically with the network width and becomes even
independent of larger width, which is in stark contrast to a linear number of
patterns in recurrent memory networks. Second, the learning leads to
progressively compact pattern spaces, which promotes highly disentangled
hierarchical representations. The upstream pattern maps the activity into a
low-dimensional space, and then the resulting embedding is further expanded or
contracted. After that, the modified embedding is re-mapped into the high-
dimensional activity space. This sequence of geometric transformation can be
understood as a linear-nonlinear hidden structure. Third, all modes are not
equally important to the generalization ability of the network, showing an
intriguing piecewise power-law behavior. Finally, the mode learning dynamics
can be predicted by the mean-field ODEs, revealing the mode specialization
transition. Therefore, the MDL inspires a rethinking of conventional deep
learning, offering a faster, more interpretable training framework. Future
works along this direction will be inspired. For example, the impact of other
structured dataset, mode dynamics in over-parameterized or recurrent networks,
and the origin of adversarial vulnerability of deep networks in terms of
geometry of the mode space.

##  I Acknowledgments

This research was supported by the National Natural Science Foundation of
China for Grant number 12122515, and Guangdong Provincial Key Laboratory of
Magnetoelectric Physics and Devices (No. 2022B1212010008), and Guangdong Basic
and Applied Basic Research Foundation (Grant No. 2023B1515040023).

## Supplemental Material

##  Appendix A Derivation of learning equations

In this section, we show how to derive the updating equations for the mode
parameters
$\bm{\theta}^{l}=(\bm{\hat{\xi}}^{l},\bm{\Sigma}^{l},\bm{{\xi}}^{l+1})$ where
the superscript $l$ indicates the layer index in the range from $1$ to $L$.
The loss function is the cross entropy $\mathcal{C}=-\sum_{i}\hat{h}_{i}\ln
h_{i}$ averaged over all training examples (divided into mini-batches in
stochastic gradient descent), where $\hat{h}_{i}$ is defined as the target
label (one-hot representation as common in machine learning). After training
the network on the training dataset with size $T$, we evaluate the
generalization performance of the network on the unseen dataset with size $V$.

In our framework of mode decomposition learning, the weight is decomposed into
the form as follows,

$\mathbf{w}^{l}=\bm{\hat{\xi}}^{l}\bm{\Sigma}^{l}(\bm{{\xi}}^{l+1})^{\operatorname{T}}.$
(S1)

The mode parameters are updated according to gradient descent of the loss
function,

$\Delta\bm{\theta}_{ij}^{l}=-\eta\mathcal{K}_{j}^{l+1}\frac{\partial
z_{j}^{l+1}}{\partial\bm{\theta}_{ij}^{l}},$ (S2)

where $\eta$ denotes the learning rate, and the error propagation term
$\mathcal{K}_{j}^{l+1}\equiv\partial\mathcal{C}/\partial z_{j}^{l+1}$. On the
top layer, $\mathcal{K}_{j}^{l+1}$ can be computed with the result
$\mathcal{K}_{j}^{L}=-\hat{h}_{j}^{L}\left(1-h_{j}^{L}\right)$. For lower
layers, the term $\mathcal{K}_{i}^{l}$ can be iteratively computed using the
chain rule. More precisely,

$\displaystyle\mathcal{K}_{i}^{l}$ $\displaystyle=\partial\mathcal{C}/\partial
z_{i}^{l}=\sum_{j}\frac{\partial\mathcal{C}}{\partial
z_{j}^{l+1}}\frac{\partial z_{j}^{l+1}}{\partial z_{i}^{l}}$ (S3)
$\displaystyle=\sum_{j}\mathcal{K}_{j}^{l+1}\sum_{\alpha}\xi^{l+1}_{i\alpha}\Sigma^{l}_{\alpha}\hat{\xi}_{j\alpha}^{l}f^{\prime}(z_{i}^{l}).$

The explicit expressions of gradient steps for the three sets of mode
parameters are given as follows,

$\displaystyle\Delta\xi^{l+1}_{j\alpha}$
$\displaystyle=-\eta\frac{\partial\mathcal{C}}{\partial\xi^{l+1}_{j\alpha}}=-\eta\mathcal{K}_{j}^{l+1}\sum_{i}\Sigma_{\alpha}^{l}\hat{\xi}_{i\alpha}^{l}h_{i}^{l},$
(S4) $\displaystyle\Delta\Sigma_{\alpha}^{l}$
$\displaystyle=-\eta\frac{\partial\mathcal{C}}{\partial\Sigma^{l}_{\alpha}}=-\eta\sum_{j}\mathcal{K}_{j}^{l+1}\sum_{i}\xi_{j\alpha}^{l+1}\hat{\xi}_{i\alpha}^{l}h_{i}^{l},$
$\displaystyle\Delta\hat{\xi}^{l}_{i\alpha}$
$\displaystyle=-\eta\frac{\partial\mathcal{C}}{\partial\hat{\xi}_{i\alpha}^{l}}=-\eta\sum_{j}\mathcal{K}_{j}^{l+1}\xi_{j\alpha}^{l+1}\Sigma_{\alpha}^{l}h_{i}^{l}.$

The above learning equations apply to the case of 1L2P case.

Next, we consider the 1L1P case ($\bm{{\xi}}^{l}=\bm{\hat{\xi}}^{l}$). Apart
from the single input pattern for the first layer $\bm{\hat{\xi}}^{1}$ and the
single output pattern for the last layer $\bm{\xi}^{L}$, two types of pattern
in each hidden layer $[\bm{{\xi}}^{l},\bm{\hat{\xi}}^{l}]$ take the same form,
and we denote $\bm{{\xi}}^{l}=\bm{\hat{\xi}}^{l}=\bm{\Xi}^{l}$. The expression
of $\mathcal{K}_{i}^{l}$ remains unchanged, and we can then update
$[\bm{\hat{\xi}}^{1},\bm{\Sigma}^{l},\bm{{\xi}}^{L}]$ according to Eq. (S4).
Next, we give the gradient descent equation for $\bm{\Xi}^{l}$ where
$l=2,...,L-1$ as follows

$\displaystyle\Delta\Xi^{l}_{j\alpha}$
$\displaystyle=-\eta\frac{\partial\mathcal{C}}{\partial\Xi^{l}_{j\alpha}}=-\eta\mathcal{K}_{j}^{l}\sum_{i}\Sigma_{\alpha}^{l-1}\hat{\xi}_{i\alpha}^{l-1}h_{i}^{l-1}-\eta\sum_{i}\mathcal{K}_{i}^{l+1}\xi_{i\alpha}^{l+1}\Sigma_{\alpha}^{l}h_{j}^{l},$
(S5)

where two terms contribute to the gradient—the first one comes from the
contribution of $\bm{\xi}^{l}$, while the second one originates from the fact
that the same pattern can act as $\bm{\hat{\xi}}^{l}$.

To ensure the weighted sum in the pre-activation is independent of the
upstream layer width and the number of modes $p^{l}$, we choose the
initialization scheme such that
$[\bm{{\xi}}^{l+1}\bm{{\Sigma}}^{l}(\bm{\hat{\xi}}^{l})^{\operatorname{T}}]_{ij}\sim\mathcal{O}(\frac{1}{\sqrt{N_{l}}})$.
This scaling is inspired by studies of Hopfield models Jiang et al. (2021). In
practice, we independently and identically sample the initial elements
$\xi^{l+1}_{i\alpha},\Sigma^{l}_{\alpha},\hat{\xi}^{l}_{j\alpha}$ from the
standard Gaussian distribution, and then the weight values are multiplied by a
factor of $1/\sqrt{N_{l}\ln N_{l}}$. Note that the number of modes are assumed
to be proportional to $\ln N_{l}$. But if the number is a constant denoted by
$P^{l}$, then the factor could be $1/\sqrt{P^{l}N_{l}}$.

Figure S1: The comparison among the SVD training, MDL (1L2P) and traditional
BP in learning performance. The network structure is specified by
$[784,100,100,10]$ in all cases. The full MDL indicates the MDL with the same
number of parameters as that of the SVD training, while the blue dot (pruned
SVD) indicates the pruning of the full SVD model $60\%$ (the modes with small
$|s_{i}|$ ranked in descending order) modes off each layer (except the output
layer) to make the consuming parameter amount comparable with that of MDL with
$p=30$. Figure S2: The comparison among low-rank decomposition (LRD), MDL
(1L2P) and traditional BP. Both decomposition methods use $p=30$. The network
structure is $[784,100,100,10]$ in all cases.

##  Appendix B Comparison to other matrix factorization methods

Here, we compared our MDL method to other matrix factorization methods in
learning performance. These other methods include singular value decomposition
(SVD), low rank decomposition (LRD) and spectral training Yang et al. (2020);
Chicchi et al. (2021).

First, the SVD learning scheme is implemented by decomposing the weight of
each layer as

$\bm{W}^{l}=\bm{U}^{l}\operatorname{diag}(\bm{s}^{l})(\bm{V}^{l})^{\top},$
(S6)

where the diagonal matrix contains $\min(N_{l},N_{l+1})$ non-zero elements in
the diagonal, and the elements of $\mathbf{s}^{l}$ is constrained to be
positive. The orthogonality is forced by two regularization terms as

$L(\bm{U},\bm{s},\bm{V})=L_{T}+\lambda_{o}\sum_{l=1}^{D}L_{o}\left(\bm{U}_{l},\bm{V}_{l}\right)+\lambda_{s}\sum_{l=1}^{D}L_{s}\left(\bm{s}_{l}\right),$
(S7)

where $L_{T}$ is the original training loss,
$L_{o}(\bm{U},\bm{V})=\frac{1}{r^{2}}\left(\left\|\bm{U}^{T}\bm{U}-\bm{I}\right\|_{F}^{2}+\left\|\bm{V}^{T}\bm{V}-\bm{I}\right\|_{F}^{2}\right)$,
and
$L_{s}(\bm{s})=\frac{\|\bm{s}\|_{1}}{\|\bm{s}\|_{2}}=\frac{\sum_{i}\left|s_{i}\right|}{\sqrt{\sum_{i}s_{i}^{2}}}$.
$r$ is the rank of $\bm{U}$ and $\bm{V}$, $\|\bullet\|_{F}$ denotes the
Frobenius norm of a matrix. The regularization term $L_{o}$ forces $\bm{U}$
and $\bm{V}$ to be orthogonal, while $L_{s}$ adjusts the sparsity level of
$\bm{s}$. The gradients for each set of parameters are derived below,

$\displaystyle\frac{\partial L_{o}}{\partial\bm{U}}$
$\displaystyle=\frac{4}{r^{2}}\left(\bm{U}^{\top}\bm{U}-\bm{I}\right)^{\top}\times\bm{U}^{\top},$
(S8) $\displaystyle\frac{\partial L_{o}}{\partial\bm{V}}$
$\displaystyle=\frac{4}{r^{2}}\left(\bm{V}^{\top}\bm{V}-\bm{I}\right)^{\top}\times\bm{V}^{\top},$
$\displaystyle\frac{\partial L_{s}}{\partial s_{i}}$
$\displaystyle=\frac{\operatorname{sign}(s_{i})\sqrt{\sum_{i}s_{i}^{2}}-\sum_{i}|s_{i}|(\sum_{i}s_{i}^{2})^{-\frac{1}{2}}s_{i}}{\sum_{i}s_{i}^{2}}.$

For comparison, we carried out the SVD learning, with $L_{o}=100$,
$L_{s}=0.0$, and $L_{o}=100$, $L_{s}=5.0$, as shown in Fig. S1. We remark that
the training cost is larger for SVD models, which can be calculated as
$\sum_{l}[N_{l}\times
N_{l+1}+\min{(N_{l},N_{l+1})}^{2}+\min{(N_{l},N_{l+1})}]$. Taking
$[784,100,100,10]$ as an example, the learning needs $109710$ parameters in
total. However, for our MDL with $p=30$ which already reaches the traditional
BP performance, the learning only needs $35910$ parameters (but traditional BP
needs $89400$ parameters). In simulations, we prune the full SVD model $60\%$
((the modes with small $|s_{i}|$ ranked in descending order)) modes off each
layer (except the output layer) to make the number of trainable parameters
comparable with that of MDL with $p=30$. We conclude that the MDL consumes
less parameters, yet produces rapid learning with even better performances.

Figure S3: The comparison among the spectral learning, MDL (1L2P) and
traditional BP. The network structure is $[784,100,100,10]$ in all cases.

Next, we fix $\bm{\Sigma}=\mathbb{I}$ in our MDL, and this reduced form is
called the low rank decomposition as follows,

$\mathbf{W}^{l}=\bm{\hat{\xi}}^{l}(\bm{\xi^{l+1}})^{\top}.$ (S9)

In the simulation, we set $p=30$. We can see in Fig. S2 that the performance
of the LRD is much worse than that of MDL and traditional BP.

For the recently proposed spectral learning Chicchi et al. (2021), a
carefully-designed transformation matrix $\mathbf{A}^{k}$ (an $N\times N$
matrix, $N$ is the total number of units in the network, and $k$ is a layer
index) is used with a spectral decomposition. The eigenvalues and the
associated basis are optimized. However, this training performs worse compared
to our MDL in the examples shown in Fig. S3.

##  Appendix C Ranking the modes according to the importance matrix

Here, we rank the modes according to the diagonal of the importance matrix,
rather than the $\tau$ measure. We found that these two ranking schemes lead
to qualitatively identical results. Removing the most important modes
(according to either the $\tau$ measure or the importance score) will
significantly impair the generalization ability of the network. Details are
illustrated in Fig. S4. The non-smooth behavior can be attributed to the
existence of mode-contribution gap, i.e., the most important modes ($<15\%$
for the $\tau$ measure; $<30\%$ for the $\Sigma$ measure) dominate the
generalization capability of the network, while other modes capture irrelevant
noise in the data.

Figure S4: Ranking modes. The network structure is $[784,100,100,10]$, and we
analyze the 1L1P case here with $p=70$. The marked percentage in the inset
indicates the generalization accuracy after removing the corresponding side of
modes in the hidden layer. The piecewise power law behavior is retained for
both types of ranking.

##  Appendix D The qualitative behavior of the MDL does not change with the
regularization strength or the hidden-layer width

Further, our MDL is in essence a matrix factorization. Therefore, the pattern
and importance matrices are not unique. However, in practice, we impose the
$\ell_{2}$ norm level for these patterns and importance scores. In fact, we
find the intriguing properties of the MDL in deep learning do not change with
the regularization strength of the $\ell_{2}$ norm (denoted as $\lambda$, see
Fig. S8). Figure S5 shows an example for the behavior of the optimal number of
modes versus hidden-layer width, while Fig. S6 shows that the piecewise power
law behavior of the $\tau$ measure does not change with the regularization
strength. In addition, the piecewise power law behavior of the $\tau$ measure
does not change with the hidden-layer width as well (Fig. S7).

Figure S5: The piecewise increasing behavior of the least $p$ with the hidden-
layer width. The network has structure $[784,N,N,10]$, and we vary $N$ to get
the corresponding least $p$, which is defined as the least number of modes
that MDL needs to reach the performance of the traditional BP. 1L2P case is
considered. (a) and (b) are obtained under different regularization strengths.
(a) $\lambda=10^{-3}$. (b) $\lambda=10^{-4}$.  Figure S6: The piecewise power
law behavior of the $\tau$ measure does not change with the regularization
strength $\lambda$ in the 1L1P case. The network structure is specified by
$[784,100,100,10]$. (a) $\lambda=0.01$. (b) $\lambda=0.001$. (c)
$\lambda=0.0001$. Figure S7: The piecewise power law behavior of the $\tau$
measure does not change with the hidden-layer width in the 1L1P case. The
network structure is specified by $[784,N,N,10]$. (a) $N=100$. (b) $N=150$.
(c) $N=200$. (d) $N=300$.  Figure S8: The $\ell_{2}$ norm of parameters
$[\bm{\xi},\bm{\hat{\xi}},\bm{\Sigma}]$ under three regularization strengths—
$\lambda=10^{-2},10^{-3},10^{-4}$. The network structure is specified by
$[784,100,100,10]$ and $p=70$ for all layers, where the 1L1P case is
considered (the 1L2P case yields qualitatively the same results). (a, b, c)
are plotted for the first layer, the second layer, and the output layer,
respectively.

##  Appendix E Subspace overlap of layered response to pairs of stimuli

In this section, we provide details of estimating the average degree of
correlation between neural responses $\mathbf{h}^{l}$ to pairs of different
input stimuli (e.g., one stimulus contains images of the same class). The
covariance of neural response in each layer to the stimulus can be
diagonalized to specify a low-dimensional subspace. The subspace is spanned by
the first $K$ principal components. The subspace overlap can then be evaluated
via the cosine of the principal angle between these two subspaces
corresponding to two different stimuli. In practice, for neural responses in
each layer to the stimulus $s_{1}$ (e.g., many images of digit $0$), we first
identify the first $K$ principal components of the covariance of
$\mathbf{h}_{1}^{\ell}$, which explains over $80\%$ of the total variance, and
then reorganize the eigenvectors to an $N_{\ell}\times K$ matrix, namely
$\mathbf{Q}^{\ell}(s_{1})$. We repeat this procedure for another stimulus
$s_{2}$, and get another matrix $\mathbf{Q}^{\ell}(s_{2})$. Therefore, the
columns of $\mathbf{Q}^{\ell}(s_{1})$ and $\mathbf{Q}^{\ell}(s_{2})$ span two
subspaces corresponding to the neural responses to $s_{1}$ and $s_{2}$
respectively. The cosine of the principal angle between these two subspaces is
calculated as follows Bjoerck and Golub (1973)

$\cos\theta_{p}\left(s_{1},s_{2}\right)=\sigma_{\max}\left({\mathbf{Q}^{\ell}(s_{1})}^{\rm
T}\mathbf{Q}^{\ell}(s_{2})\right),$ (S10)

where $\sigma_{\max}(\mathbf{B})$ denotes the largest singular value of the
matrix $\mathbf{B}$. In simulations, we consider the classification task of
the MNIST dataset, where ten classes of digits are fed into a seven-layer
neural network. Specifically, we choose $K^{\ell}$ that can explain over
$80\%$ of the total variance for each stimulus and each layer, and therefore
the value of $K^{\ell}$ varies with layer and input stimulus.

In the main text, we observe a mild increase of the subspace overlap in deep
layers. Here, as shown in Fig. S9, we link this behavior to the saturation of
the test performance with increasing number of layers and network width. In
addition, the task we consider is relatively simple, and thus three hidden
layers (five layers in total) are sufficient to classify the digits with a
high accuracy. The subspace overlap under the MDL setting thus suggests a
consistent way to determine the optimal number of layers and the network width
in practical training.

Figure S9: The averaged subspace overlap versus layers. Different number of
layers with different hidden-layer widths are considered. The results are
averaged over five independent trainings. The inset shows the corresponding
test accuracy changing with the hidden-layer width. Figure S10: A simple
illustration of the teacher-student setup with i.i.d. standard Gaussian input.
The teacher network has $k$ hidden nodes and $p^{*}$ modes, while the student
network has $m$ hidden nodes and $p$ modes. The goal of the student network is
to predict the labels generated by the teacher network, minimizing the mean-
squared error. The weights to the output layer are set to one in the linear
readout for both teacher and student networks
($\mathbf{v}=\mathbf{I}_{m},\mathbf{v}^{*}=\mathbf{I}_{k}$, and
$\mathbf{I}_{x}$ is an all-one vector of length $x$).

##  Appendix F Mean-field predictions of on-line learning dynamics

In this section, we give a detailed derivation of the mean-field ordinary
differential equations for the on-line dynamics. A sketch of the toy model
setting is shown in Fig. S10. The label for each sample $\mathbf{x}^{\nu}$
(i.i.d. standard Gaussian variable) is generated by the teacher network,

$y^{\nu}=f(\bm{\lambda}^{*\nu})=\frac{1}{k}\sum_{r=1}^{k}\sigma\left(\frac{[\bm{{\xi}}^{*}\mathbf{{\Sigma}^{*}}(\hat{\bm{\xi}}^{*})^{\operatorname{T}}]_{r}\mathbf{x}^{\nu}}{\sqrt{d}}\right)=\frac{1}{k}\sum_{r=1}^{k}\sigma\left(\lambda_{r}^{*\nu}\right),$
(S11)

where $[\mathbf{A}]_{r}$ denotes the $r$-th row of the matrix $\mathbf{A}$,
and
$\lambda_{r}^{*\nu}=[\bm{{\xi}}^{*}\bm{{\Sigma}}^{*}(\hat{\bm{\xi}}^{*})^{\operatorname{T}}]_{r}\mathbf{x}^{\nu}/{\sqrt{d}}$
represents the $r$-th element of the teacher local field vector
$\bm{\lambda}^{*\nu}\in\mathbb{R}^{k}$. To ensure the local field is
independent of the input dimension, we choose the initialization scheme for
the teacher network such that
$[\bm{{\xi}}^{*}\bm{\Sigma}^{*}(\hat{\bm{\xi}}^{*})^{\operatorname{T}}]_{ij}\sim\mathcal{O}(1)$.
More precisely, we set the elements
$\xi^{*}_{ik},\Sigma^{*}_{k},\hat{\xi}^{*}_{jk}$ to be independent standard
Gaussian variables, and then multiply the weight values by a factor of
$\frac{1}{\sqrt{\ln d}}$ for logarithmic increasing number of modes. This
scaling ensures that the magnitude of the weight values is of the order one.
Different forms of transfer function $\sigma(\cdot)$ can be considered, but we
choose the error function for the simplicity of the following theoretical
analysis. The prediction of the label by the student network for a new sample
$\mathbf{x}$ is given by

$\hat{f}(\mathbf{x},\bm{\hat{\xi}},\mathbf{{\Sigma}},\bm{{\xi}})=\frac{1}{m}\sum_{r=1}^{m}\sigma\left(\frac{[\bm{\xi}\bm{{\Sigma}}(\hat{\bm{\xi}})^{\operatorname{T}}]_{r}\mathbf{x}}{\sqrt{d}}\right)=\frac{1}{m}\sum_{r=1}^{m}\sigma\left(\lambda_{r}\right),$
(S12)

where $\lambda_{r}$ denotes the $r$-th component of the student local field
$\bm{\lambda}=\bm{\xi}\bm{\Sigma}\hat{\bm{\xi}}^{\operatorname{T}}\mathbf{x}$.
The student network has $m$ hidden nodes and $p$ patterns. For simplicity, we
assume $m=k$, $p=p^{*}$, and only the pattern $\hat{\bm{\xi}}$ is learned.

Training the student network with the one-pass gradient descent (on-line
learning) directly minimizes the following mean-squared error (MSE):

$\displaystyle\ell_{MSE}(\bm{\lambda},\bm{\lambda}^{*})=\frac{1}{2}\langle\left(\hat{f}(\bm{\lambda})-f\left(\bm{\lambda}^{*}\right)\right)^{2}\rangle,$
(S13)

where $\langle\cdot\rangle$ indicates the average over $\\{\mathbf{x,y}\\}$
that can be replaced by the average over local fields. For the Gaussian data
$\mathbb{P}(\mathbf{x})=\mathcal{N}(\mathbf{x}|0,\mathds{1})$, the dynamics of
$\ell_{MSE}$ can be completely determined by the following order parameters:
$\mathbf{Q}^{\nu}\equiv\mathbb{E}_{\mathbf{x},\mathbf{y}\sim\mathbb{P}(\mathbf{x},\mathbf{y})}\left[\bm{\lambda}^{\nu}(\bm{\lambda}^{\nu})^{\operatorname{T}}\right]=\frac{1}{d}\bm{\xi}^{\nu}\bm{\Sigma}^{\nu}(\bm{\hat{\xi}}^{\nu})^{\operatorname{T}}\bm{\hat{\xi}}^{\nu}\bm{\Sigma}^{\nu}(\bm{\xi}^{\nu})^{\operatorname{T}}$,
$\mathbf{M}^{\nu}\equiv\mathbb{E}_{\mathbf{x},\mathbf{y}\sim\mathbb{P}(\mathbf{x},\mathbf{y})}\left[\bm{\lambda}^{\nu}(\bm{\lambda}^{*\nu})^{\operatorname{T}}\right]=\frac{1}{d}\bm{\xi}^{\nu}\bm{\Sigma}^{\nu}(\bm{\hat{\xi}}^{\nu})^{\operatorname{T}}\bm{\hat{\xi}}^{*\nu}\bm{\Sigma}^{*\nu}(\bm{\xi}^{*\nu})^{\operatorname{T}}$,
and
$\mathbf{P}^{\nu}\equiv\mathbb{E}_{\mathbf{x},\mathbf{y}\sim\mathbb{P}(\mathbf{x},\mathbf{y})}\left[\bm{\lambda}^{*\nu}(\bm{\lambda}^{*\nu})^{\operatorname{T}}\right]=\frac{1}{d}\bm{\xi}^{*\nu}\bm{\Sigma}^{*\nu}(\bm{\hat{\xi}}^{*\nu})^{\operatorname{T}}\bm{\hat{\xi}}^{*\nu}\bm{\Sigma}^{*\nu}(\bm{\xi}^{*\nu})^{\operatorname{T}}$.
The corresponding matrix elements are denoted as
$q_{jl}^{\nu}\equiv\left[\mathbf{Q}^{\nu}\right]_{jl}$,
$m_{jr}^{\nu}\equiv\left[\mathbf{M}^{\nu}\right]_{jr}$ and
$\rho_{rs}\equiv[\mathbf{P}]_{rs}$. Then we can define the local-field
covariance matrix $\mathbf{\Omega}^{\nu}\in\mathbb{R}^{(k+m)\times(k+m)}$ at
time step $\nu$ as follows,

$\mathbf{\Omega}^{\nu}\equiv\left[\begin{array}[]{cc}\mathbf{Q}^{\nu}&\mathbf{M}^{\nu}\\\
(\mathbf{M}^{\nu})^{\operatorname{T}}&\mathbf{P}\end{array}\right],$ (S14)

where $\mathbf{P}$ is fixed by definition (parameters of the teacher network
are quenched), the sample index $\nu$ is also the time step in the on-line
learning setting, and the evolution of other order parameters
$\mathbf{Q}^{\nu}$ and $\mathbf{M}^{\nu}$ is driven by the gradient flow of
the mode parameters
$\bm{\theta}^{l}=(\hat{\bm{\xi}}^{l},\bm{\Sigma}^{l},\bm{{\xi}}^{l+1})$. The
loss is completely determined by the evolving order parameters,

$\displaystyle\ell_{\mathrm{MSE}}(\bm{\Omega})=\frac{1}{2}\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\bm{\Omega}\right)}\left[\left(\hat{f}(\bm{\lambda})-f\left(\bm{\lambda}^{*}\right)\right)^{2}\right]$
(S15)
$\displaystyle={\ell}_{\mathrm{t}}(\mathbf{P})+{\ell}_{\mathrm{s}}(\mathbf{Q})+\ell_{\mathrm{st}}(\mathbf{P},\mathbf{Q},\mathbf{M}),\
$

where

$\displaystyle{\ell}_{\mathrm{t}}(\mathbf{P})\equiv\mathbb{E}_{\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda}^{*}\mid
0,\mathbf{P}\right)}\left[f\left(\bm{\lambda}^{*}\right)^{2}\right]=\frac{1}{k^{2}}\sum_{r,s=1}^{k}\mathbb{E}_{\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda}^{*}\mid
0,\mathbf{P}\right)}\left[\sigma\left({\lambda}_{r}^{*}\right)\sigma\left({\lambda}_{s}^{*}\right)\right],$
(S16)
$\displaystyle{\ell}_{\mathrm{s}}(\mathbf{Q})\equiv\mathbb{E}_{\bm{\lambda}\sim\mathcal{N}(\bm{\lambda}\mid
0,\mathbf{Q})}\left[\hat{f}(\bm{\lambda})^{2}\right]=\frac{1}{m^{2}}\sum_{j,l=1}^{m}\mathbb{E}_{\bm{\lambda}\sim\mathcal{N}(\bm{\lambda}\mid
0,\mathbf{Q})}\left[\sigma\left(\lambda_{j}\right)\sigma\left(\lambda_{l}\right)\right]\text{,
}$
$\displaystyle\ell_{\mathrm{st}}(\mathbf{P},\mathbf{Q},\mathbf{M})\equiv\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\mathbf{\Omega}\right)}\left[\hat{f}(\bm{\lambda})f\left(\bm{\lambda}^{*}\right)\right]=-\frac{2}{mk}\sum_{j=1}^{m}\sum_{r=1}^{k}\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\mathbf{\Omega}\right)}\left[\sigma\left(\lambda_{j}\right)\sigma\left(\lambda_{r}^{*}\right)\right].$

To proceed, we define the integral
$I_{2}=\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\mathbf{\Omega}\right)}\left[\sigma\left(\bm{\lambda}^{\alpha}\right)\sigma\left(\bm{\lambda}^{\beta}\right)\right]$,
which has an analytic form for $\sigma(x)=\operatorname{erf}(x/\sqrt{2})$ as
follows,

$\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\mathbf{\Omega}\right)}\left[\sigma\left(\bm{\lambda}^{\alpha}\right)\sigma\left(\bm{\lambda}^{\beta}\right)\right]=\frac{2}{\pi}\arcsin\left(\frac{\Omega_{12}^{\alpha\beta}}{\sqrt{\left(1+\Omega_{11}^{\alpha\beta}\right)\left(1+\Omega_{22}^{\alpha\beta}\right)}}\right),$
(S17)

where $\Omega_{ij}^{\alpha\beta}$ denotes the element of the overlap matrix
for $\bm{\lambda}^{\alpha}$ and $\bm{\lambda}^{\beta}$, in which $\alpha$ and
$\beta$ indicate the attributes of the network—teacher or student. Therefore,
the generalization error can be estimated as follows,

$\displaystyle\begin{aligned}
\ell_{MSE}(\mathbf{\Omega})&=\frac{1}{k^{2}}\sum_{r,s=1}^{k}\frac{1}{\pi}\arcsin\left(\frac{\rho_{rs}}{\sqrt{\left(1+\rho_{rr}\right)\left(1+\rho_{ss}\right)}}\right)+\frac{1}{m^{2}}\sum_{j,l=1}^{m}\frac{1}{\pi}\arcsin\left(\frac{q_{jl}}{\sqrt{\left(1+q_{jj}\right)\left(1+q_{ll}\right)}}\right)\\\
&-\frac{2}{mk}\sum_{j=1}^{m}\sum_{r=1}^{k}\frac{1}{\pi}\arcsin\left(\frac{m_{jr}}{\sqrt{\left(1+q_{jj}\right)\left(1+\rho_{rr}\right)}}\right).\end{aligned}$
(S18)

We next consider the evolution of the order parameters, which involves only
the update of $\bm{\hat{\xi}}$ in our toy model setting. Therefore, we derive
the evolution of order parameters $q^{\nu}_{jl}$ and $m^{\nu}_{jr}$ based on
the gradient of $\hat{\bm{\xi}}$:
$\Delta\hat{\xi}_{j\alpha}=-\eta\frac{\hat{f}(\bm{\lambda})-{f}(\bm{\lambda}^{*})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\alpha}\Sigma_{\alpha}x_{j}^{\nu}$.
In the high-dimensional limit ($d\to\infty$), we use the self-averaging
property of the order parameters considering the disorder average over the
input data distribution Biehl and Schwarze (1995); Saad and Solla (1995);
Goldt et al. (2019). Then we have the following expressions,

$\displaystyle q_{jl}^{\nu+1}-q_{jl}^{\nu}$
$\displaystyle=\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}(\hat{\xi}_{n\alpha}+\Delta\hat{\xi}_{n\alpha})(\hat{\xi}_{n\beta}+\Delta\hat{\xi}_{n\beta})\Sigma_{\beta}\xi_{l\beta}\right]-q_{jl}^{\nu},$
(S19) $\displaystyle m_{jr}^{\nu+1}-m_{jr}^{\nu}$
$\displaystyle=\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}(\hat{\xi}_{n\alpha}+\Delta\hat{\xi}_{n\alpha})\hat{\xi}^{*}_{n\beta}\Sigma^{*}_{\beta}\xi^{*}_{r\beta}\right]-m_{jr}^{\nu},$

where the expectation is carried out with respect to the data distribution.

Inserting the update equation of $\hat{\bm{\xi}}$ into the equation of the
order parameter $q_{jl}$, we get

$\displaystyle q_{jl}^{\nu+1}-q_{jl}^{\nu}$
$\displaystyle=\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}(\hat{\xi}_{n\alpha}+\eta\frac{{f}(\bm{\lambda}^{*})-\hat{f}(\bm{\lambda})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\alpha}\Sigma_{\alpha}x_{n}^{\nu})(\hat{\xi}_{n\beta}+\eta\frac{{f}(\bm{\lambda}^{*})-\hat{f}(\bm{\lambda})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\beta}\Sigma_{\beta}x_{n}^{\nu})\Sigma_{\beta}\xi_{l\beta}\right]$
(S20) $\displaystyle-q_{jl}^{\nu},$
$\displaystyle=\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}\hat{\xi}_{n\alpha}\eta\frac{{f}(\bm{\lambda}^{*})-\hat{f}(\bm{\lambda})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\beta}\Sigma_{\beta}x_{n}^{\nu}\Sigma_{\beta}\xi_{l\beta}\right]$
$\displaystyle+\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}\eta\frac{{f}(\bm{\lambda}^{*})-\hat{f}(\bm{\lambda})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\alpha}\Sigma_{\alpha}x_{n}^{\nu}\hat{\xi}_{n\beta}\Sigma_{\beta}\xi_{l\beta}\right]$
$\displaystyle+\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}\eta\frac{{f}(\bm{\lambda}^{*})-\hat{f}(\bm{\lambda})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\alpha}\Sigma_{\alpha}x_{n}^{\nu}(\eta\frac{{f}(\bm{\lambda}^{*})-\hat{f}(\bm{\lambda})}{m\sqrt{d}}\sum_{i}\sigma^{\prime}(\lambda_{i})\xi_{i\beta}\Sigma_{\beta}x_{n}^{\nu})\Sigma_{\beta}\xi_{l\beta}\right],$

where we have applied the definition of
$q_{jl}^{\nu}=\frac{1}{d}\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}\hat{\xi}_{n\alpha}\hat{\xi}_{n\beta}\Sigma_{\beta}\xi_{l\beta}$
to derive the second equality. Considering the definition of
$\hat{f}(\bm{\lambda})$, ${f}(\bm{\lambda}^{*})$, $\bm{\lambda}$, and
$\bm{\lambda}^{*}$, we recast Eq. (S20) as follows,

$\displaystyle q_{jl}^{\nu+1}-q_{jl}^{\nu}$
$\displaystyle=\frac{\eta}{dkm}\mathbb{E}\left[\sum_{\beta,r,i}\lambda_{j}\sigma\left(\lambda_{r}^{*\nu}\right)\sigma^{\prime}(\lambda_{i})\xi_{i\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right]-\frac{\eta}{dm^{2}}\mathbb{E}\left[\sum_{\beta,\hat{r},i}\lambda_{j}\sigma\left(\lambda_{\hat{r}}^{\nu}\right)\sigma^{\prime}(\lambda_{i})\xi_{i\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right]$
(S21)
$\displaystyle+\frac{\eta}{dkm}\mathbb{E}\left[\sum_{\alpha,r,i}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}\lambda_{l}\sigma\left(\lambda_{r}^{*\nu}\right)\sigma^{\prime}(\lambda_{i})\right]-\frac{\eta}{dm^{2}}\mathbb{E}\left[\sum_{\alpha,\hat{r},i}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}\lambda_{l}\sigma\left(\lambda_{\hat{r}}^{\nu}\right)\sigma^{\prime}(\lambda_{i})\right]$
$\displaystyle+\frac{\eta^{2}}{dm^{4}}\mathbb{E}\left[\sum_{\alpha,\beta,i,a,r,\hat{r}}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}\sigma^{\prime}(\lambda_{i})\sigma^{\prime}(\lambda_{a})\sigma\left(\lambda_{r}^{*\nu}\right)\sigma\left(\lambda_{\hat{r}}^{*\nu}\right)\xi_{{a\beta}}\Sigma_{\beta}^{2}\xi_{l\beta}\right]$
$\displaystyle+\frac{\eta^{2}}{dm^{2}k^{2}}\mathbb{E}\left[\sum_{\alpha,\beta,i,a,r,\hat{r}}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}\sigma^{\prime}(\lambda_{i})\sigma^{\prime}(\lambda_{a})\sigma\left(\lambda_{r}^{\nu}\right)\sigma\left(\lambda_{\hat{r}}^{\nu}\right)\xi_{{a\beta}}\Sigma_{\beta}^{2}\xi_{l\beta}\right]$
$\displaystyle-2\frac{\eta^{2}}{dm^{3}k}\mathbb{E}\left[\sum_{\alpha,\beta,i,a,r,\hat{r}}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}\sigma^{\prime}(\lambda_{i})\sigma^{\prime}(\lambda_{a})\sigma\left(\lambda_{r}^{*\nu}\right)\sigma\left(\lambda_{\hat{r}}^{\nu}\right)\xi_{{a\beta}}\Sigma_{\beta}^{2}\xi_{l\beta},\right].$

To proceed, we have to estimate the integral defined by
$I_{3}(\alpha,\beta,\eta)=\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\mathbf{\Omega}\right)}\left[\sigma^{\prime}\left(\bm{\lambda}^{\alpha}\right)\bm{\lambda}^{\beta}\sigma\left(\bm{\lambda}^{\eta}\right)\right]=\frac{2}{\pi}\frac{\Omega_{23}^{\alpha\beta\eta}\left(1+\Omega_{11}^{\alpha\beta\eta}\right)-\Omega_{12}^{\alpha\beta\eta}\Omega_{13}^{\alpha\beta\eta}}{\left(1+\Omega_{11}^{\alpha\beta\eta}\right)\sqrt{\left(1+\Omega_{11}^{\alpha\beta\eta}\right)\left(1+\Omega_{33}^{\alpha\beta\eta}\right)-\left(\Omega_{13}^{\alpha\beta\eta}\right)^{2}}}$
for our transfer function $\sigma(x)=\operatorname{erf}(x/\sqrt{2})$, where
$\Omega_{ij}^{\alpha\beta\eta}$ denotes the element of the field-covariance
matrix $\mathbf{\Omega}^{\alpha\beta\eta}$. We also have to estimate the
second integral defined by
$I_{4}(i,j,k,l)=\mathbb{E}_{\bm{\lambda},\bm{\lambda}^{*}\sim\mathcal{N}\left(\bm{\lambda},\bm{\lambda}^{*}\mid
0,\mathbf{\Omega}\right)}\left[\sigma^{\prime}(\lambda_{i})\sigma^{\prime}(\lambda_{j})\sigma(\lambda_{k})\sigma(\lambda_{l})\right]$,
which has a closed form as

$I_{4}(i,j,k,l)=\frac{4}{\pi^{2}}\frac{1}{\sqrt{\bar{\Omega}_{0}^{ijkl}}}\arcsin\left(\frac{\bar{\Omega}_{1}^{ijkl}}{\sqrt{\bar{\Omega}_{2}^{ijkl}\bar{\Omega}_{3}^{ijkl}}}\right),$
(S22)

where

$\begin{array}[]{c}\bar{\Omega}_{0}^{ijkl}\equiv\left(1+\Omega_{11}^{ijkl}\right)\left(1+\Omega_{22}^{ijkl}\right)-\left(\Omega_{12}^{ijkl}\right)^{2},\\\
\bar{\Omega}_{1}^{ijkl}\equiv\bar{\Omega}_{0}^{ijkl}\Omega_{34}^{ijkl}-\Omega_{23}^{ijkl}\Omega_{24}^{ijkl}\left(1+\Omega_{11}^{ijkl}\right)-\Omega_{13}^{ijkl}\Omega_{14}^{ijkl}\left(1+\Omega_{22}^{ijkl}\right)\\\
+\Omega_{12}^{ijkl}\Omega_{13}^{ijkl}\Omega_{24}^{ijkl}+\Omega_{12}^{ijkl}\Omega_{14}^{ijkl}\Omega_{23}^{ijkl},\\\
\bar{\Omega}_{2}^{ijkl}\equiv\bar{\Omega}_{0}^{ijkl}\left(1+\Omega_{44}^{ijkl}\right)-\left(\Omega_{24}^{ijkl}\right)^{2}\left(1+\Omega_{11}^{ijkl}\right)-\left(\Omega_{13}^{ijkl}\right)^{2}\left(1+\Omega_{22}^{ijkl}\right)+2\Omega_{12}^{ijkl}\Omega_{13}^{ijkl}\Omega_{23}^{ijkl},\\\
\bar{\Omega}_{3}^{ijkl}\equiv\bar{\Omega}_{0}^{ijkl}\left(1+\Omega_{44}^{ijkl}\right)-\left(\Omega_{24}^{ijkl}\right)^{2}\left(1+\Omega_{11}^{ijkl}\right)-\left(\Omega_{14}^{ijkl}\right)^{2}\left(1+\Omega_{22}^{ijkl}\right)+2\Omega_{12}^{ijkl}\Omega_{14}^{ijkl}\Omega_{24}^{ijkl}.\end{array}$
(S23)

In an analogous way, we can derive the mean-field evolution of $m_{jr}^{\nu}$
as follows,

$\displaystyle m_{jr}^{\nu+1}-m_{jr}^{\nu}$
$\displaystyle=\frac{1}{d}\mathbb{E}\left[\sum_{n,\alpha,\beta}\xi_{j\alpha}\Sigma_{\alpha}(\hat{\xi}_{n\alpha}+\Delta\hat{\xi}_{n\alpha})\hat{\xi}^{*}_{n\beta}\Sigma^{*}_{\beta}\xi^{*}_{r\beta}\right]-m_{jr}^{\nu}$
(S24)
$\displaystyle=\frac{\eta}{kmd}\mathbb{E}\left[\sum_{\alpha,i,a}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}{\sigma(\lambda_{a}^{*})}\sigma^{\prime}(\lambda_{i})\lambda_{r}^{*}\right]-\frac{\eta}{m^{2}d}\mathbb{E}\left[\sum_{\alpha,i,a}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}{\sigma(\lambda_{a})}\sigma^{\prime}(\lambda_{i})\lambda_{r}^{*}\right],$

where the definition of $m_{jr}$ has bee used. If we define $\tau\equiv\nu/d$,
and take the thermodynamic limit of $d\to\infty$, the time step becomes
continuous, and we can thus write down the following ODEs,

$\displaystyle\frac{\mathrm{d}q_{jl}}{\mathrm{d}\tau}$
$\displaystyle=\frac{\eta}{km}\left[\sum_{\beta,r^{*},i}I_{3}(i,j,r^{*})\xi_{i\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right]-\frac{\eta}{m^{2}}\left[\sum_{\beta,\hat{r},i}I_{3}(i,j,\hat{r})\xi_{i\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right]$
(S25)
$\displaystyle+\frac{\eta}{km}\left[\sum_{\alpha,r^{*},i}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{3}(i,l,r^{*})\right]-\frac{\eta}{m^{2}}\mathbb{E}\left[\sum_{\alpha,\hat{r},i}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{3}(i,l,\hat{r})\right]$
$\displaystyle+\frac{\eta^{2}}{m^{4}}\left[\sum_{\alpha,\beta,i,a,r^{*},\hat{r}^{*}}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{4}(i,a,r^{*},\hat{r}^{*})\xi_{a\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right]$
$\displaystyle+\frac{\eta^{2}}{m^{2}k^{2}}\left[\sum_{\alpha,\beta,i,a,r,\hat{r}}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{4}(i,a,r,\hat{r})\xi_{a\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right]$
$\displaystyle-2\frac{\eta^{2}}{m^{3}k}\left[\sum_{\alpha,\beta,i,a,r^{*},\hat{r}}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{4}(i,a,r^{*},\hat{r})\xi_{a\beta}\Sigma_{\beta}^{2}\xi_{l\beta}\right],$
$\displaystyle\frac{\mathrm{d}m_{jr^{*}}^{\nu}}{\mathrm{d}\tau}$
$\displaystyle=\frac{\eta}{km}\left[\sum_{\alpha,i,a}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{3}(i,r^{*},a^{*})\right]-\frac{\eta}{m^{2}}\left[\sum_{\alpha,i,a}\xi_{j\alpha}\Sigma_{\alpha}^{2}\xi_{i\alpha}I_{3}(i,r^{*},a)\right],$

where the index in $I_{3}$ or $I_{4}$ with the symbol $*$ labels the teacher’s
local-field.

## References

  * Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ (MIT Press, Cambridge, MA, 2016).
  * Carleo et al. (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019).
  * Huang (2022) H. Huang, _Statistical Mechanics of Neural Networks_ (Springer, Singapore, 2022).
  * Roberts et al. (2022) D. A. Roberts, S. Yaida, and B. Hanin, _The Principles of Deep Learning Theory: An Effective Theory Approach to Understanding Neural Networks_ (Cambridge University Press, Cambridge, 2022).
  * Jaderberg et al. (2014) M. Jaderberg, A. Vedaldi, and A. Zisserman, arXiv:1405.3866 (2014).
  * Yang et al. (2020) H. Yang, M. Tang, W. Wen, F. Yan, D. Hu, A. Li, H. Li, and Y. Chen, 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW) (2020).
  * Giambagli et al. (2021) L. Giambagli, L. Buffoni, T. Carletti, W. Nocentini, and D. Fanelli, Nature Communications 12, 1330 (2021).
  * Chicchi et al. (2021) L. Chicchi, L. Giambagli, L. Buffoni, T. Carletti, M. Ciavarella, and D. Fanelli, Phys. Rev. E 104, 054312 (2021).
  * Jiang et al. (2021) Z. Jiang, J. Zhou, T. Hou, K. Y. M. Wong, and H. Huang, Phys. Rev. E 104, 064306 (2021).
  * Zhou et al. (2021) J. Zhou, Z. Jiang, T. Hou, Z. Chen, K. Y. M. Wong, and H. Huang, Phys. Rev. E 104, 064307 (2021).
  * (11) Y. LeCun, The MNIST database of handwritten digits, retrieved from http://yann.lecun.com/exdb/mnist.
  * Fischer et al. (2022) K. Fischer, A. Ren’e, C. Keup, M. Layer, D. Dahmen, and M. Helias, arXiv:2202.04925 (2022).
  * (13) See the supplemental material at http://… for technical and experimental details.
  * Biehl and Schwarze (1995) M. Biehl and H. Schwarze, Journal of Physics A: Mathematical and General 28, 643 (1995).
  * Saad and Solla (1995) D. Saad and S. A. Solla, Phys. Rev. Lett. 74, 4337 (1995).
  * Goldt et al. (2019) S. Goldt, M. S. Advani, A. M. Saxe, F. Krzakala, and L. Zdeborová, arXiv:1901.09085 (2019).
  * Bjoerck and Golub (1973) A. Bjoerck and G. H. Golub, Mathematics of Computation 27, 579 (1973).
  * Schwarze (1993) H. Schwarze, J. Phys. A: Math. Gen. 26, 5781 (1993).


# Biological Valuation Map of Flanders: A Sentinel-2 Imagery Analysis

###### Abstract

In recent years, machine learning has become crucial in remote sensing
analysis, particularly in the domain of Land-use/Land-cover (LULC). The
synergy of machine learning and satellite imagery analysis has demonstrated
significant productivity in this field, as evidenced by several studies [1, 2,
3, 4]. A notable challenge within this area is the semantic segmentation
mapping of land usage over extensive territories, where the accessibility of
accurate land-use data and the reliability of ground truth land-use labels
pose significant difficulties. For example, providing a detailed and accurate
pixel-wise labeled dataset of the Flanders region, a first-level
administrative division of Belgium, can be particularly insightful. Yet there
is a notable lack of regulated, formalized datasets and workflows for such
studies in many regions globally. This paper introduces a comprehensive
approach to addressing these gaps. We present a densely labeled ground truth
map of Flanders paired with Sentinel-2 satellite imagery. Our methodology
includes a formalized dataset division and sampling method, utilizing the
topographic map layout ’Kaartbladversnijdingen,’ and a detailed semantic
segmentation model training pipeline. Preliminary benchmarking results are
also provided to demonstrate the efficacy of our approach.

Index Terms—  Machine learning, Remote sensing, Satellite imagery, Semantic
segmentation, Densely-labeled dataset

##  1 Introduction and background

###  1.1 Remote sensing with deep learning

Deep learning Remote sensing and satellite imagery analysis are key in
geoscience and have expanded into urban planning, plantation coverage, and
demographic monitoring. Mapping and evaluating land use globally presents a
complex challenge, requiring interdisciplinary collaboration. Traditionally,
this field relied on visual interpretation, spectral analysis, and the
Normalized Difference Vegetation Index (NDVI). However, the advent of machine
learning has revolutionized these methods. With the growing scale and
complexity of satellite imagery data, machine learning significantly improves
data processing capabilities, enhancing applications in various related
fields.

Semantic segmentation, state-of-the-art models and challenges Semantic
segmentation, a pivotal area in computer vision, involves assigning specific
labels to each pixel in an image, a process that is crucial for a detailed
understanding of images [5]. This approach goes beyond object detection, which
merely identifies object boundaries, and image classification, which
categorizes an entire image under a single label. Semantic segmentation stands
out by offering a detailed, pixel-level classification.

One of the most influential developments in this field is the introduction of
the U-Net architecture by Ronneberger et al. (2015) [6]. Initially used for
biomedical image segmentation, U-Net has since become a cornerstone in various
segmentation tasks due to its efficiency and effectiveness [7, 8]. U-Net sets
a benchmark in the wider field of image segmentation, which includes remote
sensing research.

While deep learning methods in remote sensing can provide promising results,
it is widely acknowledged that developing reliable segmentation dataset with
sufficient sample quantity for such training experiments poses a significant
challenge due to its inherent complexity.

###  1.2 Sentinel-2 database

As part of the Copernicus Program, Sentinel-2A, operated by the European Space
Agency (ESA), offers open-source earth observation imagery [9]. This satellite
system provides high-resolution images in 13 spectral bands, covering visible,
near-infrared, and shortwave infrared spectra with spatial resolutions of 10,
20, and 60 meters. Sentinel-2A’s swath width of 290 km allows frequent
coverage of Earth’s land surfaces, crucial for monitoring environmental
changes. For data processing, we utilized the OpenEO backend framework, which
facilitates easy access and processing of Earth observation data [10]. OpenEO
provides a standardized interface to various Earth observation cloud backends,
enhancing the accessibility of Sentinel-2A data products.

In our study, we used Sentinel-2A TOC (top of canopy) reflectance data from
Terrascope, excluding the SCL (scene classification layer) for cloud masking.
We focused on 11 data channels out of 13 for training our neural network
model, initially based on the 10m resolution of visible light, with other
channels resampled to 10m to serve as complementary features. This approach
leverages the rich dataset of Sentinel-2A for advanced neural network training
in environmental studies.

##  2 Constructing the BVM dataset

BVM (biological valuation map) is a densely labeled ground truth map
describing the land use of the Belgian Flemish region [11,
12].111https://www.vlaanderen.be/inbo/en-gb/biological-valuation-map/ The
original form of the dataset is different from what is ready for training use.
Multiple steps were taken to construct such a semantic segmentation dataset.

###  2.1 Polygon and label rasterization

Coordinate reference systems The BVM data was provided by the Flemish Research
Institute for Nature and Forestry (INBO) in the format of a geographic
shapefile containing labeled polygons covering the entire land area of
Flanders. The coordinate reference system of these polygons is under the code
EPSG:31370, the Sentinel-2 satellite imagery is provided in EPSG:32631.
Projecting both systems to the most commonly used EPSG:4326 system solves the
compatibility problem of the OpenEO backend.

Polygon rasterization In the BVM data map, all land areas are assigned labels
in the form of polygons. To be used as ground truth data in our machine
learning neural models, these polygons need to be rasterized so that each
pixel on the BVM map is assigned to a class label individually as shown in
Figure 1. There are in total 14 classes names listed in Table 1:

Habitat Type | # pixels | %
---|---|---
Coastal dune habitats | 183,904 | 0.14
Cultivated land | 49,007,068 | 37.95
Grasslands | 33,033,072 | 25.58
Heathland | 1,112,397 | 0.86
Inland marshes | 337,316 | 0.26
Marine habitats | 329,994 | 0.26
Pioneer vegetation | 940,283 | 0.73
Small Landscape Features - not Specified | 81,046 | 0.06
Small Non-woody Landscape Features | 192,329 | 0.15
Small Woody Landscape Features | 925,387 | 0.72
Unknown | 13,252 | 0.01
Urban areas | 39,130,401 | 30.30
Water bodies | 2,883,922 | 2.23
Woodland and shrub | 16,657,189 | 12.90
Table 1: Habitat Types, Values, and Percentages

One can notice that the data label is unevenly distributed, with classes such
as cultivated land and urban areas taking up the majority of the land area
while classes like inland marshes and coastal dune habitats take up much less.
We believe this class imbalance is often the case for most of the first-level
administrative divisions in other countries too. Either inventing new training
strategies or new metrics to address this problem in studies that adopt this
type of dataset can be an interesting and challenging task.

Fig. 1: Example of a rasterized BVM map (different colors represent different
pixel classes)

###  2.2 Pre-processing procedures

In Flanders, satellite imagery acquisition often suffers from frequent cloudy
weather, leading to heavily clouded Sentinel-2 data cubes. This significantly
hinders model training, as cloudy pixels do not contribute useful features. To
counter this, we utilized the SCL layer of Sentinel-2 data to mask out cloudy
pixels, focusing solely on cloud-free pixels for training.

Additionally, the rectangular patches from Sentinel-2 sometimes extend beyond
Flanders’ borders. Since the BVM map only includes land-use data within
Flanders, pixels outside this area are excluded, similar to cloudy pixels. To
prevent data loss from patches partially outside Flanders, we implemented a
strategy to mirror pixels near the Flemish border, creating a margin that
ensures all pixels in a patch are usable for training. This approach maximizes
data utility and maintains the integrity of our training dataset.

###  2.3 Data preparation method

In terms of preparation of the data, we introduce the topographic map sheet
cutouts, the ”kaartbladenversnijdingen (map sheet divisions)” (Figure 2), to
help us partition the BVM map further for training, validation, and test
dataset splits. As the following figure shows: for each number X indexed
block, we group X/1N, X/2N, X/1Z, X/2Z, X/5N, X/6N, X/5Z, X/6Z together as
training blocks; X/13N, X/4N, X/3Z, X/4Z as validation blocks and X/17N, X/8N,
X/7Z, X/8Z as test blocks. Each training or validation/test block covers a
land area of size 160 $km^{2}$.

Fig. 2: A demonstration of map sheet cutouts (Source: Digital Flanders Agency)

Using a kaartbladen grid to split up Flanders and choosing corresponding
blocks is a robust way to partition training, validation, and test datasets.
To ensure homogeneity in the class distribution of these three datasets, we
report the Chi-squared distance between them. Chi-squared distance is often
used to determine the similarity of two categorical distributions; it is
formulated as follows:
$D_{\chi^{2}}(P,Q)=\sum_{i}\frac{(P(i)-Q(i))^{2}}{P(i)+Q(i)}$, where P(i) and
Q(i) are the probabilities of the i-th element in distributions P and Q,
respectively.

The number of pixels for each class in three datasets and the total number of
pixels are taken. The Chi-squared distance can thus be calculated as 0.0038
between training and validation sets, 0.00392 between training and test sets,
and 0.0048 between validation and test sets. All three values are relatively
low, indicating that the class distributions of all three sets are very
similar to each other. Such close distances are desirable in our case as it
suggests that the validation and test sets are good representatives of the
training data. Details for class distribution are in Table 2.

Class | Train (%) | Val. (%) | Test (%)
---|---|---|---
Coastal dune habitats | 0.10 | 0.06 | 0.25
Cultivated land | 34.27 | 33.92 | 32.85
Grasslands | 23.07 | 22.92 | 22.14
Heathland | 0.54 | 1.08 | 0.95
Inland marshes | 0.24 | 0.23 | 0.22
Marine habitats | 0.26 | 0.04 | 0.34
Pioneer vegetation | 0.69 | 0.65 | 0.56
Small landscape features - not specified | 0.05 | 0.06 | 0.07
Small non-woody landscape features | 0.14 | 0.13 | 0.12
Small woody landscape features | 0.63 | 0.58 | 0.73
Unknown | 0.01 | 0.007 | 0.006
Urban areas | 26.27 | 27.05 | 28.56
Water bodies | 2.08 | 2.05 | 1.74
Woodland and shrub | 11.65 | 11.22 | 11.47
Table 2: Class distribution ratios in train, validation, and test Datasets

##  3 Segmentation training tests

###  3.1 Experiment specification

Model settings Our experiment is done on an Nvidia A100 40GB GPU with 256GB
RAM, the model we adopted is vanilla U-Net with VGG16 as the encoder backbone,
and the network is initiated with pre-trained weights on the Imagenet dataset
[13]. Depending on how many input channels are required in the input layer,
the first layer of the network is modified correspondingly, where features
from either 3 or 11 channels are merged together. The training is scheduled
for 20,000 steps; the height and width of each training image are 256x256; and
the batch size is set at 8.

Dataloader Our dataloader is designed for handling Sentinel-2 satellite
imagery data cubes in the Flanders region. It efficiently manages features
like cloud coverage and region-specific data handling. The dataloader is
initialized with user-defined parameters such as areas (kaartbladen), years,
months, patch size, and data split (train, validation, test). Cloud masking
using the SCL layer from Sentinel-2 data cubes has been applied. Region-
specific processing is a key aspect, as the dataloader filters out data
patches beyond Flanders’ boundaries, maintaining alignment with the BVM map’s
geographical scope. It also employs a mirroring function for patches partially
extending beyond Flanders’ border.

Given the dataset size and convolution nature in the contracting path, where a
target region is influenced by adjacent areas, our dataloader avoids simple
checkerboard cutting of the BVM map to preserve relational information between
nearby regions. It performs random cropping of satellite images and masks,
ensuring spatial consistency in the cropped patches. This involves validating
input types, applying padding, and determining cropping coordinates. This
random cropping is vital for creating diverse training samples while
maintaining the map’s overall class distribution. The details of our project,
including the dataloader and instructions on how to use it, can be found at
https://github.com/limingshi1994/GEOInformed/.

###  3.2 Results discussion

Our experiments were conducted in three modes: 3-channel RGB, 3-channel RGNIR,
and 11-channel, using different band combinations. The 3-channel RGB mode
utilizes bands B02, B03, B04, while 3-channel RGNIR uses B03, B04, B08 (near
infrared), and the 11-channel mode employs all available 11 bands. The network
employs an early fusion mechanism where the initial convolution block merges
features from all channels to produce feature embeddings in subsequent stages.

In terms of metrics, we monitored training and validation losses and overall
pixel accuracy (Figure 3). Overall accuracy is determined by the ratio of
correctly predicted pixels to all valid target pixels. Results indicate that
the basic U-Net model, with a simple encoder architecture, achieves an average
accuracy of 67% for RGB input and 69% for RGNIR input. The slight improvement
with RGNIR input suggests that spectral bands beyond visible RGB light provide
additional useful information for feature learning. Exploring more
combinations, backed by remote sensing and geoscience theories, could
potentially further enhance performance.

For the 11-channel input, we also achieved a 69% accuracy, similar to RGB.
This could be attributed to several factors. Firstly, our U-Net model is
relatively simple, with only 23 million trainable parameters, which could be
insufficient given our dataset’s size and the complexity of temporal and
spatial pixel features. Based off scaling laws proposed by Kaplan et al. for
language models [14], we hypothesize that our model might have reached its
performance limit simply inputting more features provided less improvement.
Secondly, the feature fusion in our model occurs early via simple downsampling
convolutions. While this reduces computational demands, it may not fully
exploit the rich features of the 11-channel input. Researchers with more
resources could develop more complex CNN architectures to better utilize the
dataset’s potential.

Fig. 3: Loss and accuracy of 3 modes

##  4 Conclusion

The results from our training experiments indicate that the BVM dataset, which
we specifically developed for deep neural network training, meets the
essential criteria for a comprehensive machine learning approach in the
context of remote sensing studies. Firstly, we have contributed a densely
labeled remote sensing semantic segmentation dataset, which covers an entire
first-level administrative division of Belgium. The number of samples and
credibility of labels makes it a valuable resource for future research.
Secondly, we introduced a systematic method for dividing the dataset into
training, validation, and test segments. This approach ensures a remarkable
consistency across all three categories, which is crucial for the reliability
and effectiveness of machine learning models. Thirdly, our initial experiments
with this dataset, employing a basic, unoptimized U-Net model, yielded
satisfying results, particularly in terms of overall accuracy. These findings
underscore the dataset’s potential even with relatively simple neural network
models.

Additionally, the dataset presents several challenges that are pivotal in
advancing model training towards enhanced performance. These challenges
include addressing class imbalance, managing the substantial computational
resources required, and the timing and methods of feature fusion in the
contracting path of the network.

## References

  * [1] Michele Volpi and Devis Tuia,  “Dense semantic labeling of subdecimeter resolution images with convolutional neural networks,”  IEEE Transactions on Geoscience and Remote Sensing, 2017.
  * [2] Lei Ma, Yu Liu, Xueliang Zhang, Yuanxin Ye, Gaofei Lin, and Brian Johnson,  “Deep learning in remote sensing applications: A meta-analysis and review,”  Isprs Journal of Photogrammetry and Remote Sensing, 2019.
  * [3] Qiangqiang Yuan, Huanfeng Shen, Tongwen Li, Zhiwei Li, Shuwen Li, Yun Jiang, Hongzhang Xu, Weiwei Tan, Qianqian Yang, Jiwen Wang, Jianhao Gao, and Liangpei Zhang,  “Deep learning in environmental remote sensing: Achievements and challenges,”  Remote Sensing of Environment, 2020.
  * [4] Xiaoxiang Zhu, D. Tuia, Lichao Mou, Gui-Song Xia, Liang-pei Zhang, Feng Xu, and F. Fraundorfer,  “Deep learning in remote sensing: A comprehensive review and list of resources,”  IEEE Geoscience and Remote Sensing Magazine, 2017.
  * [5] Evan Shelhamer, Jonathan Long, and Trevor Darrell,  “Fully convolutional networks for semantic segmentation,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017.
  * [6] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  Lecture Notes in Computer Science, 2015.
  * [7] Zongwei Zhou, Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang,  “Unet++: A nested u-net architecture for medical image segmentation,”  Lecture Notes in Computer Science, 2018.
  * [8] Özgün Çiçek, Ahmed Abdulkadir, Soeren S. Lienkamp, Thomas Brox, and Olaf Ronneberger,  “3d u-net: Learning dense volumetric segmentation from sparse annotation,”  Medical Image Computing and Computer Assisted Intervention, 2016.
  * [9] M. Drusch, U. Del Bello, S. Carlier, O. Colin, V. Fernandez, F. Gascon, B. Hoersch, C. Isola, P. Laberinti, P. Martimort, A. Meygret, F. Spoto, O. Sy, F. Marchese, and P. Bargellini,  “Sentinel-2: Esa’s optical high-resolution mission for gmes operational services,”  Remote Sensing of Environment, 2012.
  * [10] Matthias Schramm, Edzer Pebesma, Milutin Milenković, Luca Foresta, Jeroen Dries, Alexander Jacob, Wolfgang Wagner, Matthias Mohr, Markus Neteler, Miha Kadunc, Tomasz Miksa, Pieter Kempeneers, Jan Verbesselt, Bernhard Gößwein, Claudio Navacchi, Stefaan Lippens, and Johannes Reiche,  “The openEO API: Harmonising the use of earth observation cloud services using virtual data cube functionalities,”  Remote Sensing, vol. 13, no. 6, 2021.
  * [11] Steven De Saeger, Pieter Dhaluin, Rémar Erens, Robin Guelinckx, Dirk Hennebel, Indra Jacobs, Myriam Kumpen, Frank Van Oost, Toon Spanhove, An Leyssen, Patrik Oosterlynck, Guy Van Dam, Martine Van Hove, and Carine Wils,  Biologische Waarderingskaart en Natura 2000 Habitatkaart, uitgave 2023,  Number 31 in Rapporten van het Instituut voor Natuur- en Bosonderzoek. Instituut voor Natuur- en Bosonderzoek, België, 2023.
  * [12] Steven De Saeger, Patrik Oosterlynck, and Desiré Paelinckx,  “The biological valuation map (BVM): a field-driven survey of land cover and vegetation in the Flemish region of Belgium,”  Documents Phytosociologiques, , no. 6, pp. 373–382, 2018.
  * [13] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei,  “Imagenet large scale visual recognition challenge,”  International Journal of Computer Vision, 2015.
  * [14] Jared Kaplan, Samuel McCandlish, Thomas Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei,  “Scaling laws for neural language models,”  arXiv: Learning, 2020.