On the Dynamics Under the Unhinged Loss and Beyond

Xiong<EMAIL_ADDRESS>
Xianming LiuThe corresponding author<EMAIL_ADDRESS>
Hanzhang Wang<EMAIL_ADDRESS>
Deming Zhai<EMAIL_ADDRESS>
Junjun Jiang<EMAIL_ADDRESS>
School of Computer Science and Technology

Harbin Institute of Technology

Harbin, 150001, China
Xiangyang Ji<EMAIL_ADDRESS>
Department of Automation

Tsinghua University

Beijing, 100084, China

Samy Bengio

Recent works have studied implicit biases in deep learning, especially the behavior of last-layer features and classifier weights. However, they usually need to simplify the intermediate dynamics under gradient flow or gradient descent due to the intractability of loss functions and model architectures. In this paper, we introduce the unhinged loss, a concise loss function, that offers more mathematical opportunities to analyze the closed-form dynamics while requiring as few simplifications or assumptions as possible. The unhinged loss allows for considering more practical techniques, such as time-vary learning rates and feature normalization. Based on the layer-peeled model that views last-layer features as free optimization variables, we conduct a thorough analysis in the unconstrained, regularized, and spherical constrained cases, as well as the case where the neural tangent kernel remains invariant. To bridge the performance of the unhinged loss to that of Cross-Entropy (CE), we investigate the scenario of fixing classifier weights with a specific structure, (e.g., a simplex equiangular tight frame). Our analysis shows that these dynamics converge exponentially fast to a solution depending on the initialization of features and classifier weights. These theoretical results not only offer valuable insights, including explicit feature regularization and rescaled learning rates for enhancing practical training with the unhinged loss, but also extend their applicability to other loss functions. Finally, we empirically demonstrate these theoretical results and insights through extensive experiments.

implicit bias, neural collapse, gradient flow, gradient descent

§ INTRODUCTION

Deep learning with neural networks has achieved great success in a variety of tasks [LeCun et al., 2015],
which, however, is not entirely understood in the interpolation and generalization of the learned models [Zhang et al., 2017, Neyshabur et al., 2017, Nakkiran et al., 2019, Bubeck and Sellke, 2021, Mei and Montanari, 2021, Zhou et al., 2023, Langer et al., 2021, Qi et al., 2021]. Many modules, including loss functions [Lin et al., 2017, Hui and Belkin, 2021] and optimization algorithms [Auer et al., 2002, Duchi et al., 2011, Zeiler, 2012, Kingma and Ba, 2015], play a crucial role in the training of deep neural networks, but lack convincing explanations due to the complexity of multilayered architectures. Recent works are devoted to simplifying modeling to better understand the behavior of DNNs,
[Papyan et al., 2020, Mixon et al., 2022, Fang et al., 2021, Tian et al., 2021, Han et al., 2022] and then to gain insights for new algorithms, theoretical, and experimental investigations.

To better understand the implicit regularization that improves the generalization of trained models in deep learning, many studies have investigated the implicit bias of gradient descent [Hardt et al., 2016, Sekhari et al., 2021], with an emphasis on the behavior of linear predictors (or called classifiers) over linearly separable data [Soudry et al., 2018, Gunasekar et al., 2018, Nacson et al., 2019, Ji and Telgarsky, 2019, Ji et al., 2020, Shamir, 2021]. In particular, Soudry et al., 2018 demonstrated that gradient descent iterates under exponentially-tailed loss minimization on separable data are biased toward $\ell_2$-maximum-margin solutions and that continuing to optimize can still lead to performance improvements, even if the validation loss increases. Further, Ji et al., 2020 showed that the gradient descent path and the algorithm-independent regularization path converge to the same direction for general losses. Shamir, 2021 formally proved that standard gradient methods never overfit on separable data. These works impressively expose the implicit regularization induced by optimization algorithms and help to understand the generalization of the learned models, but they mainly focus on the behavior of linear classifiers that is only the last layer of neural networks, while the classifier actually interacts strongly with the features produced by many nonlinear layers and parameterized layers. Thus, the relevant conclusions do not always apply to deep learning. For example, in [Soudry et al., 2018], the convergence rate of gradient descent is rather slow, wherein for almost all datasets, the distance to the maximum-margin solution decreases only as $O(1/\log t)$, and in some degenerate datasets, the rate further slows down to $O(\log\log t/\log t)$.
However, the training of DNNs typically takes only a few hundred epochs. In this paper, we show that exponential convergence is more realistic.

This paper is closely related to another research line which emerged after the empirical discovery of Neural Collapse by Papyan et al., 2020. This phenomenon precisely characterizes a pervasive inductive bias of both features and linear classifiers at the terminal phase of training, and has opened a rich area of exploring this area with simplified mathematical frameworks [Mixon et al., 2022, Fang et al., 2021, Lu and Steinerberger, 2022, Galanti et al., 2021, Fang et al., 2021, Zhu et al., 2021, Hui et al., 2022, Tirer and Bruna, 2022, Lu and Steinerberger, 2022, Kothapalli et al., 2022, Zhou et al., 2022]. Neural collapse provides a clear view of how the last-layer features and linear classifiers behave after interpolation and enables us to understand the benefit of achieving zero training error in terms of generalization and robustness [Poggio and Liao, 2020, Wang et al., 2021, Kornblith et al., 2021, Thrampoulidis et al., 2022, Gao et al., 2023], but the intermediate dynamics have remained challenging to analyze due to the intractability of cross entropy (CE). To alleviate this issue, some studies [Mixon et al., 2022, Han et al., 2022, Zhou et al., 2022, Tirer and Bruna, 2022, Kothapalli et al., 2022, Xu et al., 2023] have explored the more tractable mean squared error (MSE) loss, which performs comparably to those trained with CE [Demirkaya et al., 2020, Hui and Belkin, 2021]. However, these studies still need to make additional simplifications or assumptions for intermediate dynamics. For instance, Mixon et al., 2022 formulate the gradient flow of the unconstrained feature model as a nonlinear ordinary differential equation and then linearize the equation by claiming that nonlinear terms are negligible for models initialized near the origin. Han et al., 2022 assume that the gradient flow is along the central path which requires the linear classifier to stay MSE-optimal for features throughout the dynamics. Furthermore, Tirer et al., 2023 delve into a more practical model, demonstrating the exponential decay of within-class variability. They assume that features remain in close proximity to predefined features under MSE loss function, but note that the gradient flow is still essential along the “central path". Therefore, MSE is still not simple enough to derive exact dynamics in certain mathematical frameworks, making it difficult to grasp and bridge the gap between the modeling and practical optimization.

Reference     Contribution     [c]Extra Simplification
or Assumption     [c]Loss

Papyan et al., 2020     Empirical results and
theoretical formulation         The CE loss

Fang et al., 2021     Global optimum in
regularized cases         The CE loss

Zhu et al., 2021     Landscape analysis in
regularized cases         The CE loss

Ji et al., 2022     Landscape analysis in
unconstrained cases         The CE loss

Tirer and Bruna, 2022     Global optimum with
extended unconstrained
feature model         The MSE loss

Mixon et al., 2022     Intermediate dynamics in
unconstrained cases         The MSE loss

Han et al., 2022     Intermediate dynamics in
weight-regularized cases         The MSE loss

Tirer et al., 2023     Intermediate dynamics by
perturbation analysis in
regularized cases         The MSE loss

This paper     Intermediate dynamics and
convergence analysis in
unconstrained, regularized,
spherical constrained cases
with time-varying
learning rates         The unhinged loss

Comparison of recent analysis for investigating the behavior of last-layer features and prototypes. Compared with prior work, this paper considers the time-varying learning rate, which is often used in practice, and provides intermediate dynamics and convergence analysis in unconstrained, regularized and spherical constrained cases while requiring fewer simplifications or assumptions.

In this paper, our objective is to analyze the closed-form dynamics under gradient descent within the layer-peeled model [Fang et al., 2021] (also known as the unconstrained features model [Mixon et al., 2022]) with minimal simplifications and assumptions. To achieve this, we introduce the unhinged loss, which possesses a concise form and intuitively expresses the classification objective. To further contextualize the performance of the unhinged loss in relation to the CE loss, we investigate scenarios involving fixing classifier weights with specific structures. Compared to previous works (as depicted in Table 1), the unhinged loss provides more mathematical opportunities for gaining insights into deep learning with closed-form dynamics, while demanding fewer simplifications or assumptions for intermediate dynamics. This equips us for more practical considerations and more rational designs. Our primary contributions are outlined as follows:

* We introduce the unhinged loss for analyzing closed-form dynamics in deep learning with as few simplifications or assumptions as possible. The unhinged loss with some auxiliary techniques can achieve comparative performance to the CE loss.

* We derive exact dynamics of last-layer features and prototypes in unconstrained, regularized, prototype-anchored and spherical constrained cases as well as the NTK-invariant case. For spherical constrained cases that do not exhibit convexity, Lipschitzness, and $\beta$-smoothness, we also prove that gradient descent biases the normalized features towards a global minimizer.

* We provide the corresponding convergence analysis, which shows that the features and classifier weights converge to a solution depending on the initialization rather than induce the neural collapse solution that forms a simplex equiangular tight frame, suggesting that not all losses under gradient descent would lead to neural collapse (as verified in <Ref>).

* We prove that the rate of convergence is exponential as a function of $\zeta(t)=\int_0^t\eta(\tau)\mathrm{d}\tau$, where $\eta(\tau)$ denotes the learning rate over time. This exponential convergence rate highlights the impact of the interaction of features and classifier weights.

* Moreover, we provide some insights and extensive verification for improvements in practical training with the unhinged loss and other losses (cf. <Ref>).

§ THE UNHINGED LOSS

In this paper, we primarily investigate the behavior of last-layer features and classifier weights in DNNs for classification tasks. We conduct our study to datasets comprising inputs from $C$ different classes, each with $N$ examples. The last-layer features $\mathbf{h}_{i,c} = \mathbf{f}_{\Theta}(\mathbf{x}_{i,c}) \in \mathbb{R}^p$ [For the closed-form solution of neural collapse, it is worth noting that Zhou et al., 2022 assume $p \geq C-1$, while Han et al., 2022 assume $p > C$, as last-layer features typically have a higher dimension than the number of classes. This work will explicitly address the relationship between $p$ and $C$ in some scenarios, covering all potential choices of feature dimension in others.] are obtained from the $i$-th example $\mathbf{x}_{i,c}$ through a series of parameterized layers $\mathbf{f}_{\Theta}:\mathcal{X}\rightarrow\mathbb{R}^{p}$, commonly treated as free optimization variables [Mixon et al., 2022, Fang et al., 2021, Han et al., 2022, Ji et al., 2022]. The final layer of the network, known as the linear classifier, is equipped with a class prototype $\mathbf{w}_c \in \mathbb{R}^p$ and a bias $b_c \in \mathbb{R}$ for each class $c \in [C]$. It predicts a label using the rule $\argmax_{c'}(\langle\mathbf{w}_{c'}, \mathbf{h}\rangle + b_{c'})$ [Unless specified otherwise, we denote vectors using boldface italicized letters, and elements within the vector are denoted with italicized letters and subscripts.].

To better understand the dynamics of last-layer features and prototypes based on gradient flow or gradient descent, we consider a concise loss that offers more mathematical opportunities than the hard-to-analyze CE loss and the MSE loss [The CE loss is known to have non-analytical intermediate dynamics, which makes it difficult to investigate the precise nature of the dynamics. In addition, the intermediate dynamics under MSE can be derived into closed-form, which however is quadratic and cannot be easily analyzed. Therefore, some additional simplifications or assumptions on CE and MSE are required to simplify the intermediate dynamics.]. Specifically, we investigate a generalized form of the unhinged loss [Van Rooyen et al., 2015] as follows:

\begin{equation}
\label{the-unhinged-loss}
L_{\gamma}(\mW\vh+\vb,y)=-\vw_y^\top \vh-b_y+\gamma\sum_{j\neq y}(\vw_j^\top\vh+b_j),
\end{equation}

where $\gamma>0$ is the trade-off parameter and $y$ denotes the class label of the feature $\vh$. When $\gamma=\frac{1}{C-1}$, the loss $\frac{1}{C-1}\sum_{j\neq y}[\vw_y^\top\vh+b_y-(\vw_j^\top\vh+b_j)]$ can be regarded as the unhinged version of the hinge loss[The multi-class hinge loss is formulated as $L_{\text{hinge}}(\vs,y;m)=\sum_{i\neq y}\max\{0, s_i-s_y+m\}$, where $m$ is the margin term.]  that removes the maximum operator and the margin term. We also note that the sample margin $m(\vh,y)=\vw_y^\top\vh+b_y-\max_{j\neq y}(\vw_j^\top\vh+b_j)$ [Koltchinskii and Panchenko, 2002, Cao et al., 2019, Zhou et al., 2022] is defined to measure the discriminativeness for a sample, which satisfies $m(\vh,y)\le \frac{1}{C-1}\sum_{j\neq y}[\vw_y^\top\vh+b_y-(\vw_j^\top\vh+b_j)]$, i.e.,  $L_{\gamma}(\mW\vh+\vb,y)$ with $\gamma=\frac{1}{C-1}$ averaging the margins over all non-target classes is the lower bound of $-m(\vh,y)$. Here, we replace $\frac{1}{C-1}$ with an additional parameter $\gamma$ that balances positive and negative logits to draw general conclusions. More clarification about the unhinged loss can be found in Appendix <ref>.

Intuitively, the unhinged loss in <ref> promotes the learned feature $\vh$ to increase the logit of the target class while decreasing the logits of the other classes. If we follow up the layer-peeled model [Fang et al., 2021] to restrict the norms of both features and prototypes, the global minimizer of $\frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C L_{\gamma}(\mW\vh_{i,c},y_{i,c})$ (the bias term $\vb$ is omitted) will lead to Neural Collapse [Papyan et al., 2020, Han et al., 2022]:

For norm-bounded prototypes and features, i.e., $\|\vw_c\|_2\le E_1$ and $\|\vh_{i,c}\|_2\le E_2$, $\forall i\in[N], \forall c\in[C]$, the global minimizer of $\frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C L_{\gamma}(\mW\vh_{i,c},y_{i,c})$ implies neural collapse when $p\ge C-1$. More specifically, the global minimizer is uniquely obtained at $\frac{\vw_i^\top\vw_j}{\|\vw_i
\|_2\|\vw_j\|_2}=-\frac{1}{C-1}$, $\forall i\neq j$, $\frac{\vw_{y_{i,c}}^\top\vh_{i,c}}{\|\vw_{y_{i,c}}\|_2\|\vh_{i,c}\|_2}=1$,  $\|\vw_c\|_2=E_1$, and $\|\vh_{i,c}\|_2=E_2$, $\forall i\in[N],$ $\forall c\in[C]$.

This lemma shows that the Neural Collapse solution[The Neural Collapse solution exhibits three critical properties [Papyan et al., 2020]: (i) Within-class variability collapse, the within-class variation of features becomes negligible as these features collapse to their class means; (ii) Convergence to self-duality, the class means and classifier weights converge to each other, up to rescaling; (iii) Convergence to a simplex equiangular tight frame (ETF), the vectors of class means converge to having equal length, forming equal-sized angles between any given pair, and being the maximally pairwise-distanced configuration. For instance, prototypes $\mW$ that forms a simplex ETF satisfies $\mW^\top\mW=\frac{C}{C-1}\mI-\frac{1}{C-1}\1\1^\top$ when $C<p$]. is the only global optimal solution to minimize $\frac{1}{CN}\sum_{i,c}L_{\gamma}(\mW\vh_{k,c},y_{i,c})$ in the norm-bounded case.
However, there exists an undesired direction to minimize the unhinged loss in unconstrained cases, since the norm of features and prototypes tends to grow to infinity. For example, we can directly scale up $\mW$ and $\vb$ to obtain a smaller loss if $L_{\gamma}(\mW\vh+\vb, y)<0$, which will happen analogously to CE [Liu et al., 2016, Wang et al., 2017, Zhou et al., 2022]. In this paper, we will analytically characterize the direction in which features $\mH$ and prototypes $\mW$ diverge. Specifically, we show that the gradient flow or gradient descent with the unhinged loss will exhibit an implicit bias associated with the initialization of features and prototypes. We further investigate the prototype-anchored case, wherein class prototypes are anchored with some specific structures, and the unhinged loss can perform comparative to CE.

§ MAIN THEORETICAL RESULTS

In this section, we begin with the commonly used assumption that treats last-layer features as free optimization variables. We then conduct a comprehensive analysis of the closed-form dynamics of last-layer features and prototypes under the unhinged loss in various scenarios, including unconstrained, regularized, prototype-anchored, spherical constrained, and NTK-invariant cases. Additionally, we provide convergence analyses for each case, revealing a surprising result: all cases exhibit exponential convergence. All proofs can be found in Appendix <ref>.

§.§ The Unconstrained Case

We first consider the unconstrained case [Mixon et al., 2022, Ji et al., 2022] in which there is no constraint or regularization on features and prototypes, i.e., learning with the following objective

\begin{equation}
\label{unconstrain-optimization}
\gL= \frac{1}{CN}\sum_{i,c} \big[-\vw_{y_{i,c}}^\top\vh_{i,c}-b_{y_{i,c}}+\gamma\sum_{j\neq y_{i,c}}(\vw_j^\top \vh_{i,c}+b_j)\big],
\end{equation}

which can be reformulated as

\begin{equation}
\label{objective}
\gL =\Tr(\mY^\top\mW^\top\mH)+\tfrac{\gamma C-\gamma-1}{C}\1_C^\top\vb,
\end{equation}

where $\mH=[\vh_{1,1},\ldots,\vh_{1,C},\vh_{2,1},\ldots,\vh_{N,C}]\in\mathbb{R}^{p\times CN}$ is the matrix resulting from stacking together the feature vectors as columns, $\mY=\frac{1}{CN}(\gamma\1_C\1_{CN}^\top-(1+\gamma)(\mI_C\otimes \1^\top_{N}))$, $\mW=[\vw_1,\vw_2,\ldots,\vw_C]\in\mathbb{R}^{p\times C}$ is the matrix resulting from stacking class prototypes as columns, $\otimes$ denotes the Kronecker product, $\mI_C$ is the identity matrix,  $\1_C$, $\1_N$, and $\1_{CN}$ are the length-$C$, -$N$, and -$CN$ vectors of ones, respectively. For brevity, we represent the label set $\{y_{i,c}\}_{1\le i\le N, 1\le c\le C}$ as the columns of the matrix $\mI_C\otimes \1_N^\top$.

Remark We follow the unconstrained features modeling perspective [Mixon et al., 2022, Ji et al., 2022] or the layer-peeled model [Fang et al., 2021] that treats $\mH$ as a free optimization variable. Within this model, we analyze the continuous dynamics of features $\mH$, prototypes $\mW$ and biases $\vb$ with gradient flow where  time-of-training is denoted by the variable $t$[Intuitively, we interpret $t=0$ as the initial state, that is $\mH(0)=\mH_0$, $\mW(0)=\mW_0$, and $\vb(0)=\vb_0$.].

Taking the partial derivative with respect to $\mH$, $\mW$, and $\vb$, respectively, we have the following:

\begin{equation}
\begin{aligned}
&\nabla_{\vb}\gL=\tfrac{\gamma C-\gamma-1}{C}\1_C,
\end{aligned}
\end{equation}

and the corresponding learning dynamics with respect to the gradient flow is

\begin{equation}
\label{gradient-flow}
\begin{aligned}
&\mH'(t)=\eta_1(t) \mW(t)\mM,\\
&\mW'(t)=\eta_2(t) \mH(t)\mM^\top, \\
&\vb'(t)=-\eta_2(t)\tfrac{\gamma C-\gamma-1}{C}\1_C,
\end{aligned}
\end{equation}

where $\mM=-\mY=\tfrac{1}{CN}((1+\gamma)(\mI_C\otimes \1^\top_{N})-\gamma\1_C\1_{CN}^\top)$, $\eta_1(t)$ and $\eta_2(t)$ are the time-varying part[In this paper, we are interested in investigating the effect of changing the learning rate over time. Thus, $\eta(t)$ is not the complete rate, but rather the time varying part of the learning rate. Specifically, for the example gradient flow $x'_t=\frac{\mathrm{d} x_t}{\mathrm{d} t}=\eta(t) \nabla f(x_t)$ considered in this paper is a continuous-time approximation of gradient descent in the limit of an infinitesimally small time step $\Delta_t$. Specifically, the corresponding  of gradient descent in discrete time can be formulated as $x_{t+\Delta_t}=x_{t}-\Delta_t \cdot \eta(t)\nabla f(x_t)$, where $\Delta_t\cdot \eta(t)$ denotes the complete learning rate, and $\Delta_t$ is the time step.] of the learning rate of the features $\mH$ and weights $(\mW,\vb)$, respectively.
The reason for introducing different learning rates is that the representation $\mH$ is a result of the interaction between a number of nonlinear layers rather than a completely free variable like network parameters. This implies that even if we use the same learning rate to optimize all network parameters, the feature $\mH$ assumed to be a free optimization variable is almost impossible to be optimized at this learning rate[In this paper, we mainly assume that $\eta_1(t_1)\eta_2(t_2)=\eta_1(t_2)\eta_2(t_1)$ for any pair of values of $t$, $t_1$ and $t_2$. This condition will be satisfied if and only if $\eta_1(t)$ is a scaled version of $\eta_2(t)$, i.e., $\eta_1(t)=s\cdot \eta_2(t)$]. Moreover, we consider time-varying rates $\eta_1=\eta_1(t)$ and $\eta_2=\eta_2(t)$ that are usually adopted in practical implementations, such as the cosine annealing decay [Loshchilov and Hutter, 2017]. As can be seen, the dynamics of features $\mH$ and prototypes $\mW$ are independent of the bias term $\vb$, thus we can analyze the dynamics of $\mH$ and $\mW$ jointly, and analyze $\vb$ independently:

Consider the continual gradient flow in gradient-flow, let $\mZ(t)=(\mH(t),\mW(t))$, if $\eta_1(t_1)\eta_2(t_2)=\eta_1(t_2)\eta_2(t_1)$ for any $t_1,t_2\ge0$, we have the following closed-form dynamics

\begin{equation}
\label{closed-form-zt}
\begin{aligned}
\mZ(t)=&\Pi_1^+\mZ_0 \left(\alpha_1^+(t)\mC(t)+\beta_1^+(t)\mI_{C(N+1)}\right)\\
+&\Pi_1^-\mZ_0 \left(\alpha_1^-(t)\mC(t)+\beta_1^-(t)\mI_{C(N+1)}\right)\\
+&\Pi_2^-\mZ_0 \left(\alpha_2^-(t)\mC(t)+\beta_2^-(t)\mI_{C(N+1)}\right)+\Pi_3\mZ_0,
\end{aligned}
\end{equation}

\begin{equation}
\vb(t)=\vb_0 +\frac{(1+\gamma-\gamma C)\zeta_2(t)}{C} \1_C,
\end{equation}

where $\alpha_1^\eps$, $\alpha_2^\eps$, $\beta_1^\eps$ and $\beta_2^\eps$ for $\eps\in\{\pm\}$ are the scalars that only depend on $C$, $N$, $\gamma$, $\eta_1(0)$ and $\eta_2(0)$ (where the detailed forms of these scalars can be seen in the appendix), $\mZ_0=(\mH_0,\mW_0)$, $\mC(t)=\left(\begin{smallmatrix}
\zeta_1(t)\mI_{CN} & 0\\
0 & \zeta_2(t)\mI_C\\
\end{smallmatrix}\right)$, $\zeta_1(t)=\int_0^t\eta_1(\tau)\mathrm{d}\tau$, $\zeta_2(t)=\int_0^t\eta_2(\tau)\mathrm{d}\tau$, $\Pi_1^{\eps}$, $\Pi_2^{\eps}$ and $\Pi_3$ for $\eps\in\{\pm\}$ are orthogonal projection operators onto the following respective eigenspaces:

\begin{equation}
\begin{aligned}
& \gE_1^\eps:=\{(\mH,\mW):\mH=\eps\cdot\tfrac{1}{\sqrt{N}}(\mW\otimes \1_N^\top),\mW\1_C=0\},\\
& \gE_2^\eps:=\{(\mH,\mW):\mH=\eps\cdot \tfrac{1}{\sqrt{N}}\vh\1^\top_{CN},\mW=\vh\1_C ^\top,\vh\in\mathbb{R}^p\},\\
& \gE_3:=\{(\mH,\mW): \mH(\mI_C\otimes\1_N)=0,\mW=0\}.
\end{aligned}
\end{equation}

[Train Loss]

[Train Accuracy]

Verification of the behavior of gradient descent iterates in <Ref> with $\gamma\in\{0, 0.05, 0.005, 0.001, \frac{1}{C-1}\}$, where we set $p=512$, $C=100$, $N=10$, $\eta_1(t)=\eta_2(t)=0.1$ (i.e., $s=\frac{\eta_1(0)}{\eta_2(0)}=1$, thus $\overline{\mZ}=\Pi_1^+\mZ_0$ according to Corollary <ref>), and then randomly initialize $\mH_0$ and $\mW_0$.   The red box in the figures represents the zoomed-in view of the last 2,000 epochs. (a) The training loss. (b) The training accuracy with the prediction rule $\argmax_c \vw_c^\top\vh$. As expected, the features align to their corresponding prototypes when $\gamma<\frac{2}{C-2}$. (c) The distance between $\hat{\mZ}(t)=\mZ(t)/\|\mZ(t)\|_2$ and $\hat{\overline{\mZ}}=\overline{\mZ}/\|\overline{\mZ}\|_2$. As expected in <ref>, the convergence rate is exponential when $0<\gamma<\frac{2}{C-2}$, and will be fastest if $\gamma=\frac{1}{C-1}$. (d) The norm of $\mZ(t)$, which increases exponentially. As can be noticed, $\hat{\mZ}(t)$ does not converge to $\hat{\overline{\mZ}}$ but tend to be orthogonal to $\hat{\overline{\mZ}}$ when $\gamma=0.05 >\frac{2}{C-2}$, that is, $\lim_{t\rightarrow\infty}\|\hat{\mZ}(t)-\hat{\overline{\mZ}}\|_2=\sqrt{2}$. Moreover, the exponential curve of $\|\boldsymbol Z(t)\|_2$ will be $O(e^{3.95 \times \sqrt{10}\times 10^{-4}t})$ when $\gamma=0.05$, which is faster than those curves of the rate roughly around $O(e^{\sqrt{10}\times 10^{-4}t})$ when $\gamma=0,0.001,0.005$ and $\frac{1}{C-1}=\frac{1}{99}$. This key distinction elucidates why the curves overlap when $\gamma=0, 0.001, 0.005$, and $\frac{1}{C-1}$.

Remark Note that $\gE^\eps_1$, $\gE^\eps_2$ and $\gE_3$ are orthogonal to each other,  $\gE^+_1$ (or $\gE^-_1$) denotes the subspace where all features are in the same (or opposite) direction of their corresponding prototypes while the mean of prototypes is zero, $\gE^\eps_2$ denotes the subspace where all features and prototypes collapse respectively into two scaled versions of the same unit vector, and $\gE_3$ denotes the subspace where the mean of all features from the same class is zero with all prototypes being zero. For classification tasks, we expect the features align to their corresponding prototypes with a cosine similarity of 1, i.e., the solution in $\gE=\{(\mH,\mW):\mH=k \mW\otimes\1_N^\top, \mW\in\sR^{p\times C}, k\in\sR^+\}$ that implies two manifestations of Neural Collapse: within-class variability collapse and convergence to self-duality [Papyan et al., 2020]. In the following, we prove that the normalized dynamics $\frac{\mZ(t)}{\|\mZ(t)\|}$ under the unhinged loss will converge to a solution in $\gE$:

Under the conditions and notation of Theorem <ref>, let $s=\frac{\eta_1(0)}{\eta_2(0)}$. If $0<\gamma < \frac{2}{C-2}$ (where $C>2$) or $C=2$, and $\lim_{t\rightarrow\infty}\zeta_1(t)=\infty$, the gradient flow (as in <ref>) will behave as:

\begin{equation}
\label{dynamics-of-Z(t)}
e^{-\frac{(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}\mZ(t)=\overline{\mZ} +\bm{\Delta}(t),
\end{equation}

where $\overline{\mZ}=\left(\tfrac{1+\sqrt{s}}{2}\mH_1^++\tfrac{1-\sqrt{s}}{2}\mH_1^-, \tfrac{1+\sqrt{s}}{2\sqrt{s}}\mW_1^+-\tfrac{1-\sqrt{s}}{2\sqrt{s}}\mW_1^-\right)$,  $(\mH_1^+,\mW_1^+)=\Pi_1^+\mZ_0$, $(\mH_1^-,\mW_1^-)=\Pi_1^-\mZ_0$, and the residual term $\bm{\Delta}(t)$ decreases as $\|\bm{\Delta}(t)\|=O\left(e^{\frac{\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\cdot\max\{-\gamma C,(C-2)\gamma-2\}}\right)$, and so the normalized $\mZ(t)$ converges to $\frac{\overline{\mZ}}{\|\overline{\mZ}\|}$ in

\begin{equation}
\label{convergence-rate-unconstrained-case}
\left\|\frac{\mZ(t)}{\|\mZ(t)\|}-\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\right\|=O\left(e^{\frac{\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\cdot\max\{-\gamma C,(C-2)\gamma-2\}}\right),
\end{equation}

which indicates $\lim_{t\rightarrow\infty}\frac{\mZ(t)}{\|\mZ(t)\|}=\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\in\gE$.
Moreover, if $\gamma\neq \tfrac{1}{C-1}$, $\lim_{t\rightarrow \infty}\frac{\max_i b_i(t)}{\min_i b_i(t)}=1$.

The corollary above shows that even without any mandatory constraints, the gradient flow under the unhinged loss will converge to a solution $\frac{\overline{\mZ}}{\|\overline{\mZ}\|}$ that belongs to $\gE$ and can be determined by the initialization of $\mZ_0$ and the ratio $s$. It is worth noting, however, that this solution does not conform to the geometric structure of Neural Collapse. Neural collapse typically entails the formation of a simplex equiangular tight frame, as exemplified in the case of CE [Papyan et al., 2020], MSE [Han et al., 2022], margin-based losses [Zhou et al., 2022],  and other losses [Zhou et al., 2022]. This result suggests that not all losses will lead to Neural Collapse solutions.

In addition, the rate of convergence is exponential as a function of the cumulative learning rates, i.e., $O\left(e^{\frac{\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\cdot\max\{-\gamma C,(C-2)\gamma-2\}}\right)$, which indicates that the convergence of updating both features and prototypes by gradient descent is much faster than $O(1/\log t)$ that only updates prototypes (linear predictors) on linearly separable data [Soudry et al., 2018]. In a sense, this convergence rate may help explain why training deep neural networks usually takes only several hundred or thousand epochs. Moreover, if $\gamma= \tfrac{1}{C-1}$, we can obtain the fastest convergence of <ref>, that is, $\left\|\frac{\mZ(t)}{\|\mZ(t)\|}-\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\right\|=O\left(e^{-\frac{\sqrt{\zeta_1(t)\zeta_2(t)}}{(C-1)\sqrt{N}}}\right)$. As shown in <ref>, if we set $\eta_1(t)=\eta$ as a constant learning rate, then $\zeta_1(t)=\eta t\rightarrow\infty$ as $t\rightarrow \infty$, and the gradient flow in <ref> exhibits an exponential convergence rate of $\frac{\mZ(t)}{\|\mZ(t)\|}$ to $\frac{\overline{\mZ}}{\|\overline{\mZ}\|}$. However, this will lead to an exponential increase in the norm of features and prototypes, with the rate $e^{\frac{(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}$. Such growth is almost unbearable for the practical training of DNNs. Therefore, in what follows, we consider to limit the excessive growth of these norms.

§.§ The Regularized Case

In this subsection, we focus on a regularized optimization problem that introduces $\ell_2$ regularization (also known as “weight decay”) on features, prototypes, and biases. This regularization term helps to prevent the norms of these variables from growing too large during training. The optimization problem is defined as follows:

\begin{equation}
\label{regularized-objective}
\min_{\mW,\mH,\vb}\Tr(\mY^\top \mW^\top\mH)+\tfrac{\gamma C-\gamma-1}{C}\1_C^\top\vb + \tfrac{\lambda_1}{2}\|\mH\|_F^2+\tfrac{\lambda_2}{2}\|\mW\|_F^2+\tfrac{\lambda_2}{2}\|\vb\|_2^2,
\end{equation}

where $\|\cdot\|_F$ denotes the Frobenius norm, $\lambda_1$ and $\lambda_2$ are trade-off parameters[To draw general conclusions, we consider a different regularization strength for each component.] that control the regularization strength.

Taking the partial derivative with respect to $\mH$, $\mW$, and $\vb$, we have

\begin{equation}
\notag
\begin{aligned}
&\nabla_{\vb}\gL=\tfrac{\gamma C-\gamma-1}{C}\1_C+\lambda_2\vb,
\end{aligned}
\end{equation}

and the corresponding learning dynamics following the gradient flow with time-vary learning rates, $\ell_2$ regularization can be formulated as

\begin{equation}
\label{regularized-gradient-flow}
\begin{aligned}
&\mH'(t)=\eta_1(t) \mW(t)\mM-\lambda_1\eta_1(t)\mH(t),\\
&\mW'(t)=\eta_2(t)\mH(t)\mM^\top-\lambda_2\eta_2(t)\mW(t), \\
&\vb'(t)=-\eta_2(t)\tfrac{\gamma C-\gamma-1}{C}\1_C-\lambda_2\eta_2(t)\vb(t).
\end{aligned}
\end{equation}

The dynamics of this regularized gradient flow can be proved as follows

Consider the continual gradient flow in regularized-gradient-flow, let $\mZ(t)=(\mH(t),\mW(t))$. If $\eta_1(t_1)\eta_2(t_2)=\eta_1(t_2)\eta_2(t_1)$ for any $t_1,t_2\ge0$, we have the following closed-form dynamics:

\begin{equation}
\begin{aligned}
\mZ(t)=&\Pi_1^+\mZ_0\begin{pmatrix}a_1^+(t)\mI_{CN} & 0\\0 & b_1^+(t)\mI_{C}\end{pmatrix}+\Pi_1^-\mZ_0\begin{pmatrix}a_1^-(t)\mI_{CN} & 0\\0 & b_1^-(t)\mI_{C}\end{pmatrix}\\
&+\Pi_2^+\mZ_0\begin{pmatrix}a_2^+(t)\mI_{CN} & 0\\0 & b_2^+(t)\mI_{C}\end{pmatrix}+\Pi_2^-\mZ_0\begin{pmatrix}a_2^-(t)\mI_{CN} & 0\\0 & b_2^-(t)\mI_{C}\end{pmatrix}\\
&+\Pi_3\mZ_0\begin{pmatrix}a_3(t)\mI_{CN} & 0\\0 & b_3(t)\mI_{C}\end{pmatrix},
\end{aligned}
\end{equation}

\begin{equation}
\vb(t)=\phi(t)\left(\vb_0+\tfrac{1+\gamma-\gamma C}{C}\psi(t)\1_C\right),
\end{equation}

where $\Pi_1^+\mZ_0$, $\Pi_1^-\mZ_0$, $\Pi_1^+\mZ_0$, $\Pi_1^-\mZ_0$, and $\Pi_3\mZ_0$ follow the definition in Theorem <ref>, $a_1^\eps$, $a_2^\eps$, $b_1^\eps$, $b_2^\eps$, $a_3$, and $b_3$ for $\eps\in\{\pm\}$ are the scalars that depend only on $C$, $N$, $\gamma$, $\lambda_1$, $\lambda_2$, $\eta_1$, and $\eta_2$ (where the detailed forms can be seen in <ref>), $\phi(t)=\exp(-\lambda_2\int_0^t\eta_2(\tau)\mathrm{d}\tau)$, and $\psi(t)=\int_0^t\zeta_2(\tau)\exp(\lambda_2\int_0^\tau\eta_2(s)\mathrm{d}s)\mathrm{d}\tau$.

The convergence under the regularized case can also be derived as:

Under the conditions and notation of Theorem <ref>, let $s=\frac{\eta_1(0)}{\eta_2(0)}$. If $0<\gamma<\frac{2}{C-2}$ (where $C>2$) or $C=2$, $\lambda_1=\lambda_2=\lambda$, and $\lim_{t\rightarrow\infty}\zeta_1(t)=\infty$, then there exist constants $\pi_h^+,\pi_h^-,\pi_w^+$, $\pi_w^-$, and $\omega$ only depending on $\lambda$, $\gamma$, $s$, $C$, and $N$, such that the gradient flow (as in <ref>) behaves as:

\begin{equation}
\label{convergence-of-Z-l2}
\left\|\frac{\mZ(t)}{\|\mZ(t)\|}-\frac{\mZ_\pi}{\|\mZ_\pi\|}\right\|=O\left(e^{-\omega\zeta_2(t)}\right),
\end{equation}

where $(\mH_1^+,\mW_1^+)=\Pi_1^+\mZ_0$, $(\mH_1^-,\mW_1^-)=\Pi_1^-\mZ_0$, and $\mZ_\pi=(\pi_h^+\mH_1^++\pi_h^-\mH_1^-,\pi_w^+\mW_1^++\pi_w^-\mW_1^-)$.

Furthermore, we have the following convergence results for $\mZ(t)$:

* If $\lambda > \frac{1+\gamma}{C\sqrt{N}}$, then $\lim_{t\rightarrow\infty}\|\mZ(t)\|=0$;

* If $\lambda =\frac{1+\gamma}{C\sqrt{N}}$, then $\lim_{t\rightarrow\infty}\mZ(t)=\left(\mH_1^+ + \frac{1-s}{1+s}\mH_1^-, \mW_1^+ - \frac{1-s}{1+s}\mW_1^-\right)$;

* If $\lambda < \frac{1+\gamma}{C\sqrt{N}}$, then $\lim_{t\rightarrow\infty}\|\mZ(t)\|=\infty$.

As can be seen, the features and classifier weights under the unhinged loss converge to the solution depending on the initialization with the form $\mZ_\pi(t)=(\pi_h^+ \mH_1^+ +\pi_h^- \mH_1^-, \pi_w^+ \mW_1^+ +\pi_w^- \mW_1^-)$, where $\mH_1^+=\frac{1}{\sqrt{N}}(\mW^+\otimes \1_N^\top)$ and $\mH_1^-=-\frac{1}{\sqrt{N}}(\mW^-\otimes \1_N^\top)$. Thus, the solution implies the property of within-class variability collapse. To full encompass all properties of Neural Collapse (when the class number and the feature dimensionality $p$ satisfies $C\le p-1$), it is essential for $\pi_w^+ \mW_1^+ +\pi_w^- \mW_1^-$ to form a Simplex ETF, and for $\frac{\pi_h^+}{\pi_w^+}=-\frac{\pi_h^-}{\pi_w^-}$ to hold, inducing the property of convergence to self-duality (i.e., the class means and linear classifiers converge to each other, up to rescaling). For example, when $\lambda=\frac{1+\gamma}{C\sqrt{N}}$, Corollary demonstrates that the solution satisfies the property of convergence to self-duality.

Moreover, the results in Corollary <ref> suggest that adding an appropriate weight decay on both features and prototypes can avoid impractical effects, since the norm of $\mZ(t)$ shrinking to 0 or diverging toward infinity will significantly affect the training of DNNs. Several recent works [Zhu et al., 2021, Zhou et al., 2022] described that the features are implicitly penalized, but this implicit penalization may be fragile when using the unhinged loss (as depicted in <ref>). Consequently, we emphasize the importance of adding explicit regularization to features, rather than relying solely on the implicit penalization attached by other components like batch normalization and weight decay. Explicit feature regularization also plays a role in mitigating minority collapse [Fang et al., 2021]. Firstly, minority collapse leads to features of the minority classes approaching $0$, since the minimization of the objective pays too much emphasis on enlarging the feature norms of the majority classes. In this context, explicit feature regularization can effectively restrain the excessive growth of feature norms of the majority classes. Secondly, explicit feature regularization in some sense intuitively reduces the energy of features in the optimization program [Fang et al., 2021], thereby shrinking the feasible domain of features. This further mitigates the imbalance of feature norms between majority and minority classes. In a nutshell, we need to limit the growth of feature norms.

§.§ The Prototype-anchored Case

We further consider the prototype-anchored case in which the class prototypes $\mW$ are fixed[That is, $\mW$ is not updated, which can be done by simply setting $W.\text{require\_grad} = False$ in PyTorch] into some specific structures (e.g., a simplex ETF) during the training process and the features $\mH$ are with $\ell_2$ regularization. Then, the dynamics of $\mH$ in <ref> will be formulated as the first-order non-homogeneous linear difference equation:

\begin{equation}
\label{eq:pal}
\mH'(t)=\eta(t)\mW\mM-\lambda \eta(t)\mH(t),
\end{equation}

and the solution to the linear difference equation can be easily derived as follows

Consider the continual gradient flow (<ref>) in which the prototypes $\mW$ are fixed, we have the closed-form dynamics:

\begin{equation}
\mH(t)=e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\mH(0)+\frac{1-e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}}{\lambda} \mW\mM,
\end{equation}

which further indicates that $\left\|\mH(t)-\frac{1}{\lambda}\mW\mM\right\|=O\left(e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\right)$ when $\lim_{t\rightarrow \infty}e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}=0$.

When the time-varying learning rate satisfies that $\lim_{t\rightarrow \infty}e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}=0$, then $\mH(t)$ converges to $\frac{1}{\lambda}\mW\mM$. That is, the unhinged loss in the prototype-anchored case coincides with a feature alignment task such that each feature $\vh_{i,c}$ in class $c$ aligns to $\frac{1}{\lambda CN}(\vw_c -\gamma\sum_{j\neq c}\vw_j)$. It's worth noting that when $\lambda=0$, the behavior of $\mH(t)$ can be described as $\mH(t)=\mH(0)+\int_0^t\eta(\tau) \mathrm{d}\tau \mW\mM$, which is predominantly influenced by the term $\mW\mM$ when the integral $\int_0^t\eta(\tau) d\tau$ becomes significantly large.

As depicted in <ref>, prototype-anchored learning demonstrates its effectiveness in mitigating training instability by transforming the classification problem into a feature alignment task. In this context, the unhinged loss with PAL yields results that are not only comparable but, in some cases, even superior to those achieved with the Cross-Entropy (CE) loss. This underscores the practicality of the unhinged loss as a valuable loss function for training classification models.

§.§ The Spherical Constrained Case

We consider another constrained case in which features are restricted on the $p$-sphere $\sS^{p-1}=\{\vx:\|\vx\|_2=1,\vx\in\sR^p\}$ by explicitly performing $\ell_2$ normalization[$\ell_2$ normalization can also prevent arithmetic overflow or underflow occurring in the training of DNNs.], and we fix the prototypes[The relevant studies are still few and often require some strict assumptions since the learning dynamics is very complicated when $\vw$ participates the optimization process with both feature and prototypes normalization. In this paper, we are going to try a more concise theoretical analysis with fixed prototypes.] to satisfy $\mW\1_C=0$[This aims to simplify <ref> as the objective of feature alignment, that is, $L_{\gamma}(\mW\hat{\vh},y)=-\vw_y^\top\hat{\vh}+\gamma\sum_{j\neq y}\vw_j^\top\hat{\vh}=\frac{(1+\gamma)\|\vw\|_2}{2}(\|\hat{\vh}-\hat{\vw}_y\|_2^2-2)$, and the global minimum will be obtained at $\hat{\vh}=\hat{\vw}_y$.], then the optimization problem in <ref> can be reformulated as

\begin{equation}
\label{spherical-objective}
\min_{\mH}-\frac{1+\gamma}{CN}\Tr((\mI_C\otimes \1_N)\mW^\top\hat{\mH}),
\end{equation}

where $\hat{\mH}=(\hat{\vh}_{1,1},...,\hat{\vh}_{c,N})$ and $\hat{\vh}=\frac{\vh}{\|\vh\|_2}$ denotes the $\ell_2$-normalized vector.

Take the partial derivative with respect to $\mH$, then the discrete dynamical system based on gradient descent is formulated as

\begin{equation}
\label{normalized-discrete-dynamics}
\mH(t+1)=\mH(t)+\tfrac{(1+\gamma)\eta(t)}{CN}\left(\tfrac{\partial \hat{\mH}}{\partial\mH}\big|_{\mH=\mH(t)}\right)^\top\mW(\mI_C\otimes \1_N^\top),
\end{equation}

where $(\frac{\partial \hat{\mH}}{\partial \mH})^\top$ is a vector-wise operator, and for any vector $\vh_{i,c}$ in $\mH'\in\sR^{p\times CN}$, we have $(\frac{\partial \hat{\mH}}{\partial \mH})^\top \vh_{i,c}=\left(\frac{1}{\|\vh_{i,c}\|_2}(\mI_p-\hat{\vh}_{i,c}\hat{\vh}_{i,c}^\top)\vh_{i,c}'\right)$.

Despite the fact that the optimization objective in <ref> does not show convexity, Lipschitzness, and $\beta$-smoothness on $\mH$ due to the $\ell_2$ normalization operator, we prove that the normalized features which obey the gradient descent iterates in <ref> can still converge to their corresponding normalized prototypes, i.e., achieve the global minimum of <ref>:

Considering the discrete dynamics in <ref>, if  $\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)>-1$ for any $i\in[N]$ and $c\in[C]$, the learning rate $\eta(t)$ satisfies that $\frac{\eta(t)}{\|\vh_{i,c}(t)\|_2}$ is non-increasing,  $\frac{\eta(0)(1+\gamma)}{CN\|\vh_{i,c}(0)\|_2}\le \frac{1}{\|\vw_c\|_2}$,  $\lim_{t\rightarrow\infty}\frac{\eta(t+1)}{\eta(t)}=1$, and there exists a constant $\varepsilon>0$, s.t., $\eta(t)>\varepsilon$, then we have

\begin{equation}
\lim_{t\rightarrow\infty}\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|=0,
\end{equation}

and further if $\lim_{t\rightarrow\infty}\|\mH(t)\|<\infty$, then there exists a constant $\mu>0$, such that the error above shows exponential decrease:

\begin{equation}
\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|= O(e^{-\mu t}).
\end{equation}

Moreover, if $\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)=-1$, then $\vh_{i,c}(t)=\vh_{i,c}(0)$.

Remark As shown in the theorem above, we prove that gradient descent exerts a bias that steers the normalized features toward the global minimizer $\hat{\mW}(\mI_C\otimes\1_N^\top)$, achieving exponential convergence under certain favorable conditions. This global solution shows two properties of Neural Collapse: within-class variability collapse and convergence to self-duality. If $\hat{W}$ also forms a simplex ETF, the global solution will exhibit all properties of Neural Collapse. Additionally, we establish that if the inner product between $\hat{\vw}_c$ and $\hat{\vh}_{i,c}(0)$ (i.e., the cosine similarity between $\vw_c$ and $\vh_{i,c}(0)$) is $-1$, then $\vh_{i,c}(t)=\vh_{i,c}(0)$, as the gradient induced by $\ell_2$ normalization will be 0. More details can be found in the proof of Theorem  <ref>. Therefore, we just analyze the case where the inner product $\hat{\vw}_c$ and $\hat{\vh}_{i,c}(0)$ is strictly greater than $-1$.

§.§ The NTK-Invariant Case

We further perform some analysis in the case where the neural tangent kernel $\nabla_{\bm{\Theta}}\mH^\top\nabla_{\bm{\Theta}}\mH$ [Jacot et al., 2018, Chizat et al., 2019, Yang and Hu, 2020, Yang and Hu, 2021] is assumed to be invariant during training. We prove that the unhinged loss has more potential to derive exact dynamics. Specifically, consider the last-layer feature $\vh_{i,c}=\vf_{\bm{\Theta}}(\vx_{i,c})$ extracted from the example $\vx_{i,c}$. When DNNs are trained using gradient descent to minimize the composition $\mathcal{L}\circ \boldsymbol{W}\circ \boldsymbol f_{\boldsymbol \Theta}$, we have

\begin{equation}
\begin{aligned}
\boldsymbol{\Theta}(t+1)-\boldsymbol{\Theta}(t)&=-\eta \nabla_{\boldsymbol{\Theta}}\mathcal{L}\big|_{\bm{\Theta}=\boldsymbol{\Theta}(t)},\\
\boldsymbol{W}(t+1)-\boldsymbol{W}(t)&=-\eta\nabla_{\boldsymbol{W}}\mathcal{L}\big|_{\bm{W}=\boldsymbol{W}(t)}.
% \\&\boldsymbol{b}(t+1)-\boldsymbol{b}(t)=-\eta\nabla_{\boldsymbol{b}}\mathcal{L}.
\end{aligned}
\end{equation}

According to the first-order Taylor expansion, we obtain

\begin{equation}
\boldsymbol H(t+1)-\boldsymbol H(t)\approx \nabla_{\boldsymbol \Theta} \boldsymbol H^\top [\boldsymbol{\Theta}(t+1)-\boldsymbol{\Theta}(t)].
\end{equation}

For the unhinged loss, we can obtain the following closed-form dynamics

Assume that the neural tangent kernel $\nabla_{\boldsymbol \Theta}\boldsymbol H ^\top\nabla_{\boldsymbol \Theta}\boldsymbol H$ remains invariant during iterations. Let $\vz(t)$ denote the row-first vectorization of $\begin{pmatrix}\boldsymbol{H}(t) & 0\\0 & \boldsymbol{W}(t)\end{pmatrix}$, $\boldsymbol B=\begin{pmatrix}0 & \boldsymbol M^\top\\\boldsymbol M & 0\end{pmatrix}$, and $\boldsymbol A=\begin{pmatrix}0 & \nabla_{\boldsymbol \Theta}\boldsymbol H ^\top\nabla_{\boldsymbol \Theta}\boldsymbol H\\\mathbf I_C & 0\end{pmatrix}$. Considering the eigendecomposition $\boldsymbol A=\boldsymbol U_{\boldsymbol A}\boldsymbol \Lambda_{\boldsymbol A} \boldsymbol U_{\boldsymbol A}^{-1}$ and $\boldsymbol B=\boldsymbol U_{\boldsymbol B}\boldsymbol \Lambda_{\boldsymbol B} \boldsymbol U_{\boldsymbol B}^{-1}$, we have

\begin{equation}
\mC\boldsymbol z(t)=\exp[(\boldsymbol \Lambda_{\boldsymbol A}\otimes \boldsymbol \Lambda_{\boldsymbol B})t]\mC\boldsymbol z(0),
\end{equation}

where $\mC=\boldsymbol U_{\boldsymbol B}^{-1}\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I$.

Though Theorem <ref> does not imply properties of Neural Collapse, it reveals a specific dynamics under the unhinged loss, which exhibits an exponential convergence rate, contingent upon the assumption that the neural tangent kernel remains constant. To demonstrate the mathematical simplicity of the unhinged loss, we consider the dynamics of the MSE loss $\mathcal L=\frac{1}{CN}\Vert\boldsymbol W^\top \boldsymbol H+\boldsymbol b \mathbf 1_{CN}^\top -\mathbf I\otimes \mathbf 1_N^\top \Vert_F^2$ for comparison, and similarly derive that

\begin{equation}
\begin{aligned}
\boldsymbol H'(t)&=\nabla_{\boldsymbol \Theta} \boldsymbol H^\top \nabla_{\boldsymbol \Theta} \boldsymbol H \boldsymbol W(t)\boldsymbol C(t),\quad
\boldsymbol W'(t) &= \boldsymbol H(t) \boldsymbol C(t)^\top,\quad\boldsymbol{b}'(t)&=\boldsymbol C(t) \mathbf 1_{CN},
\end{aligned}
\end{equation}

where $\boldsymbol C(t)=\frac{2\eta }{CN}(\boldsymbol W(t)^\top \boldsymbol H(t) + \boldsymbol b(t)\mathbf 1_{CN}^\top - \mathbf I_{C}\otimes \mathbf 1_N^\top)$. As can be seen, the above gradient flow is nonlinear and then intractable to analyze the exact dynamics though $\nabla_{\boldsymbol \Theta} \boldsymbol H^\top \nabla_{\boldsymbol \Theta} \boldsymbol H$ stays constant.

Through the above analysis, it is clear to see that the unhinged loss offers more mathematical opportunities than the MSE loss to analyze the closed-form dynamics under the assumption of an invariant neural tangent kernel.

§ INSIGHTS AND EXPERIMENTS

In this section, we provide some insights into better training DNNs according to the conclusions in <ref>. We then corroborate our theoretical results and insights with extensive experiments. More specifically, in Section <ref>, we propose to use prototype-anchored learning (PAL) as a means of resolving the instability issues that arise during training with the unhinged loss. In Section <ref>, we conduct experiments to highlight the benefits of explicit feature regularization on imbalanced learning and out-of-distribution (OOD) detection. In these experiments, we employ the unhinged loss with PAL and the CE loss as training objective. In Section <ref>, we propose the rescaling learning rates (RLR) with feature norms for the spherical case to address the problem of slow convergence resulting from feature normalization, which may also have implications for improving other methods of performing feature normalization. More details and results can be found in Appendix <ref>.

§.§ The Unhinged Loss with Prototype-Anchored Learning

Since directly using the unhinged loss will lead to volatile effects, which is mainly reflected in the rapid increase of feature norms and the imbalance between class prototypes when training DNNs with the stochastic gradient method, as shown in <ref>. Inspired by recent works [Zhou et al., 2022, Kasarla et al., 2022, Yang et al., 2022] that use the Neural Collapse structure as an inductive bias (also called prototyping-anchored learning, PAL), we fix prototypes $\mW$ as a simplex ETF during training, i.e., $\mW^\top\mW=\frac{C}{C-1}\mI-\frac{1}{C-1}\1\1^\top$[These prototypes $\mW$ can be obtained using one of two methods: (i) By minimizing the objective $\sum_{i=1}^C\log \frac{\exp\left(s\hat{\vm}_i^\top\hat{\vh}_i\right)}{\sum_{j=1}^C\exp\left(s\hat{\vm}_j^\top\hat{\vh}_i\right)}$ and setting $\mW=\hat{\mM}$, where $s$ is a scale factor [Zhou et al., 2022]; or (ii) By deriving them from the standard simplex ETF $\mM=\sqrt{\frac{C}{C-1}}(\mI_C-\frac{1}{C}\1_C\1_C^\top)$, i.e., $\mW=s\mU\mM\in\sR^{p\times C}$, where $\mU\in\sR^{p\times C}$ ($p\ge C$) is a partial orthogonal matrix [Papyan et al., 2020].].

Experimental Results We conduct experiments on widely-used classification datasets including CIFAR-10, CIFAR-100, and ImageNet-100. To mitigate training instability under the unhinged loss, we employ the Prototype-Anchored Learning (PAL) and Feature-Normalized PAL (FNPAL). As depicted in <ref>, the results obtained by the unhinged loss with PAL and FNPAL variants demonstrate comparable or even better performance in comparison to CE. This substantiates the feasibility of utilizing the unhinged loss as a practical training objective for standard classification tasks.

Validation accuracy of different loss functions on CIFAR-10, CIFAR-100, and ImageNet, where $*$ denotes training with explicit feature regularization. PAL and FNPAL denote the model trained with prototype-anchored learning (PAL) and feature-normalized and prototype-anchored learning (FNPAL) [Zhou et al., 2022]. The curve in the red box represents the zoomed-in curve of the last 50 epochs. As can be seen, DNNs trained with the unhinged loss can achieve comparative or even better performance than those of CE.

Validation accuracies on long-tailed CIFAR-10/-100 with CE and different explicit feature regularization under Cross-Entropy (CE) and the unhinged loss (UL) with PAL. The imbalance ratio $\rho=\frac{\max_i n_i}{\min_i n_i}$ is the ratio between sample sizes of the most frequent and least frequent classes, and $\rho=1$ denotes the original CIFAR-10/-100. $\lambda=0$ denotes the model training with CE. All values are percentages. Bold numbers indicate the results that are better than the baseline vanilla CE or UL. The best results are underlined.

Dataset     5c|Long-tailed CIFAR-10     5cLong-tailed CIFAR-100

Imbalance Ratio      100     50     20     10     1     100     50     20     10     1

Vanilla CE     67.81     72.93     83.97     88.37     95.28     33.37     39.40     42.96     56.38     75.42

CE ($\lambda=5e-6$)     67.84     72.85     83.17     89.06     95.27     36.00     41.92     50.75     60.13     76.48

CE ($\lambda=1e-5$)     67.74     76.14     84.17     89.19     95.23      36.61     42.36     49.21     58.91     77.34

CE ($\lambda=5e-5$)     69.74     77.29     84.92     88.64     95.39     34.88     42.74     54.72     60.84     76.19

Vanilla UL     72.60     77.20     85.46     89.29     95.23     39.20     45.81     54.66     59.92     75.37

UL ($\lambda=1e-7$)     71.85     78.80     85.76     89.68     95.33     41.39     46.13     55.25     61.15     75.18

UL ($\lambda=5e-7$)     73.47     79.54     86.28     89.36     95.14     41.85     47.58     54.00     61.31     75.51

UL ($\lambda=1e-6$)     73.71     79.40     86.12     89.43     95.50     40.85     47.07     55.71     61.47     75.84

UL ($\lambda=5e-6$)     71.13     78.90     85.33     89.20     95.10     41.87     45.83     55.92     60.87     77.20

§.§ Explicit Feature Regularization

In this paper, we directly consider explicit feature regularization to avoid excessive growth of feature norms, i.e., adding the regularization term $\lambda \sum_{\vx\in\gB}\|\vf_{\boldsymbol{\Theta}}(\vx)\|_2^2$ in the objective. Explicit regularization on features can significantly remedy over-confidence and even improve generalization. Moreover, adding explicit feature regularization can speed up the convergence of $\hat{\mH}(t)$ to $\hat{\mW\mM}$ according to Theorem <ref>, as verified in <ref>.

Experimental Results To validate the role of explicit feature regularization, we conduct experiments on two tasks: (i) Long-tailed recognition on benchmarks CIFAR-10-LT and CIFAR-100-LT with artificially created long-tailed settings; (ii) Out-of-distribution (OOD) detection between SVHN and CIFAR-10/-100. For long-tailed classification, we follow the controllable class imbalance strategy in [Cao et al., 2019] by reducing the number of training examples per class and keeping the validation set unchanged. As shown in Table <ref>, explicit feature regularization effectively improves the performance on long-tailed classification in most cases, even for normal classification. For ODD detection, we train ResNet-18 and ResNet-34 on in-distribution datasets CIFAR-10 and CIFAR-100, respectively, and then use SVHN as the OOD dataset to evaluate the performance.

[FPR95: 43.04]

[FPR95: 27.87]

[FPR95: 89.84]

[FPR95: 81.41]

[FPR95: 52.10]

[FPR95: 24.94]

[FPR95: 95.54]

[FPR95: 87.83]

Distribution of energy scores (a-d) [Liu et al., 2020]  and feature norms (e-h) from classification models trained without (a & c & e & g) or with (b & d & f & h) explicit feature regularization (EFR) ($\lambda=1e-5$). (a & b & e & f) and (c & d & g & h) are from ResNet-18 [He et al., 2016] trained on CIFAR-10 and from ResNet-34 trained on CIFAR-100, respectively. As can be seen, EFR can improve the performance of OOD detection by alleviating the over-confidence of OOD samples and making the energy scores of ID samples more concentrated. More intuitively, comparing (f) to (e) and (h) to (g), EFR effectively limits the growth of feature norms and significantly improves the distinction between ID samples and OOD samples in feature norm.

OOD detection performance using softmax-based [Hendrycks and Gimpel, 2016], energy-based [Liu et al., 2020], and feature norm-based approaches while model training with explicit feature regularization (EFL) (where $\lambda=1e-5$ for CE and $\lambda=5e-7$ for the unhinged loss (UL)). We use ResNet-18 and ResNet-34 to train on the in-distribution datasets CIFAR-10 and CIFAR-100, respectively. We then use SVHN [Netzer et al., 2011] as the OOD dataset to evaluate the performance of OOD detection. All values are percentages. $\uparrow$ indicates large values are better, and $\downarrow$ indicates smaller values are better. The best results are underlined.

Dataset $\gD_{\text{in}}^{\text{test}}$     Method     FPR95 $\downarrow$      AUROC $\uparrow$     AUPR $\uparrow$

5*CIFAR-10     4cSoftmax-based / Energy-based / Feature Norm-based

CE     52.09 / 43.04 / 52.10     91.67 / 91.94 / 89.54     84.11 / 82.80 / 77.06

CE + EFL     37.39 / 27.87 / 24.94     93.90 / 94.60 / 94.17     85.48 / 85.34 / 83.15

UL     53.29  / 53.33 / 52.34     87.72 / 87.73 / 89.45     74.86 / 74.83 / 79.05

UL+EFL     27.27  / 27.24 / 26.62     94.27 / 94.46 / 95.12     86.51 / 86.60 / 88.46

5*CIFAR-100     4cSoftmax-based / Energy-based / Feature Norm-based

CE        87.75 / 89.84 / 95.54     71.01 / 71.94 / 59.54     55.42 / 56.69 / 43.21

CE+EFL     81.48 / 81.41 / 87.83     77.02 / 78.03 / 73.91     62.92 / 63.66 / 58.81

UL     80.43  / 76.86 / 76.23     66.97 / 74.83 / 76.19     37.83 / 45.07 / 47.51

UL+EFL     75.74  / 74.64 / 73.59     75.21 / 75.36 / 76.80     45.16 / 45.27 / 48.34

§.§ Rescaled Learning Rate for the Spherical Case

In this subsection, we propose a strategy of rescaling the learning rate for the spherical case in which the features are performed $\ell_2$ normalization.

As shown in Theorem <ref> and its proof, the convergence behavior in the spherical constrained case is characterized by the term $\frac{\eta(t)}{\|\vh_{i,c}(t)\|_2}$, when $\|\vh_{i,c}\|_2$ is monotonically increasing and thus leads to a slower convergence rate. This result underscores the influence of the feature norm on the gradient dynamics in spherical case.  Also, we observe that the training plots in the spherical case exhibit a slower rate in comparison to the CE loss in Figure 2. To bridge the performance of the unhinged loss to that of CE, we propose a straightforward approach—rescaling the learning rates based on feature norms, that is, scaling up the learning rate $\eta(t)$ with the feature norm $\|\vh_{i,c}\|_2$ for each example[We can also implement this strategy by rescaling the loss, but the multiplicative feature norm stops gradients.]. Though we change the learning rates, we still guarantee the convergence analysis in Theorem <ref> if $\eta(t)$ is non-increasing and satisfies $\frac{\eta(0)(1+\gamma)}{CN}\le \frac{1}{\|\vw_c\|_2}$.

Experimental Results To confirm the efficacy of the rescaled learning rate in accelerating convergence for the spherical case, we experiment with the unhinged loss (FNPAL). As illustrated in Figure <ref>, the utilization of the unhinged loss (FNPAL) with the rescaled learning rate (RLR) clearly demonstrates accelerated convergence in comparison to its counterpart without RLR.

§ CONCLUSION

In this paper, we introduced the unhinged loss as a surrogate to analyze the behavior of last-layer features and prototypes. Due to the conciseness of the unhinged loss, we derived exact dynamics under gradient descent in various scenarios, including unconstrained, regularized, and spherical constrained cases, as well as the case with an invariant neural tangent kernel. Furthermore, we demonstrated that these dynamics converge exponentially to a specific solution depending on the initialization. Inspired by these results, we proposed additional insights for improvements, such as using prototype-anchored learning with the unhinged loss to bridge its performance to that of cross-entropy loss, employing explicit feature regularization to address over-confidence, and implementing a rescaled learning rate to accelerate convergence in the spherical case. Finally, we validated these theoretical results and insights through extensive experiments, covering numerical analysis, visual classification, imbalanced learning, and out-of-distribution detection. We anticipate that the unhinged loss will serve as a valuable tool for the community to gain a deeper understanding of the behavior of deep neural networks, extending beyond the scope of this paper.

§ ACKNOWLEDGEMENTS

This work was supported in part by National Natural Science Foundation of China under Grants 92270116 and 62071155, and in part by the Fundamental Research Funds for the Central Universities (Grant No.HIT.OCEF.2023043 and HIT.DZJJ.2023075). Xiong Zhou thanks the support from MEGVII Technology and Qiyuan Lab.

Appendix for “On the Dynamics under the Unhinged Loss and Beyond”

§ CLARIFICATION FOR THE UNHINGED LOSS

In this section, we provide some clarification for the unhinged loss as follows.

* About the contribution. As stated in the body of this paper, the unhinged loss serves as an alternative object of study, which has the advantage of being linear in the outputs of the model. Thus, it can provide more opportunities to analyze more concrete behaviors and explore the shortcomings of the existing modeling. The main contribution of this work is to glimpse into deep learning through a concise form under the unhinged loss, which enjoys the great advantage to simplify theoretical research very much. Moreover, it is not just a tool for theoretical analysis. To make the unhinged loss practical, we further provide several insights theoretically and empirically, including using prototype-anchored learning with the unhinged loss (Sec. <ref>), applying explicit regularization on features (Sec. <ref>), rescaling the learning rates for the spherical constrained case (Sec. <ref>), and other potential insights (Sec. <ref>).

* About the criterion of correct classification. The value of the unhinged loss in <ref> cannot be directly used as the criterion to judge whether a sample is classified correctly since the unhinged loss can be arbitrary negative in unconstrained cases, since the gradient of the unhinged loss if $\frac{\partial L_{\gamma}}{\partial \vh}=-(\vw_y-\gamma\sum_{i\neq y}\vw_i)$ cannot converge to zero unless $\vw_y-\gamma\sum_{i\neq y}\vw_i=0$ for all $y$. This is also the main reason why the unhinged loss cannot be used directly for practical training. To overcome this drawback, we turn to limit the rapid growth of norms. For instance, in the spherical constrained case, the features are constrained to the unit sphere, and class prototypes are fixed as a simplex equiangular tight frame (ETF). Moreover, under this case, the unhinged loss would be lower bounded by $-(1+\gamma)$, and thus can be used in practical training (as shown in <ref>). We can also derive that $L_{\gamma}<-(1+\gamma)\sqrt{\frac{C-2}{2(C-1}}\Rightarrow \vw_y^\top \vh>\max_{i\neq y}\vw_i^\top\vh$ when these prototypes are fixed as a simplex ETF, so the sample can be correctly classified if the value of the unhinged loss is less than $(1+\gamma)\sqrt{\frac{C-2}{2(C-1}}$.

* About other losses.

1) For the CE and multi-binary CE loss, we have the following inequality.

\begin{aligned}
\le&\min\bigg\{\log (1+\exp(-\vw_y^\top \vh-b_y))+\gamma \sum_{j\neq y}\log (1+\exp(\vw_j^\top \vh+b_j)),\\
&(\gamma(C-1)-1)(\vw_y^\top+b_y) \vh-\gamma(C-1)\log \frac{\exp(\vw_y^\top \vh+b_y)}{\sum_{i=1}^C\exp(\vw_i^\top \vh+b_i)}\bigg\}.
\end{aligned}

2) For the Wasserstein loss in Arjovsky et al., 2017, we can find that the empirical form of the unhinged loss $\frac{1}{N}\sum_{i=1}^N(\vw_y^\top\vh+b_y-\gamma\sum_{i\neq y}(\vw_i^\top\vh+b_i))$ resembles the Wasserstein loss $\E_{x\sim P_r}[f_w(x)]-\E_{x\sim P_\theta}[f_w(x)]$ if we regard positive logits $\vw_y^\top\vh+b_y$ as scores from $P_r$ and negative logits $\vw_i^\top\vh+b_i$ ($i\neq y$) as scores from $P_{\theta}$. The main difference is that in the unhinged loss the two logits $\vw_y^\top\vh+b_y$ and $\vw_i^\top\vh+b_i$ ($i\neq y$) are computed from $\vh$ (i.e., the same input), rather than being sampled separably from two distributions as done in Wasserstein loss. Moreover, the unhinged loss does not constrain the values of these scores, in contrast the Wasserstein loss requires $f_w$ to be Lipschitz continuous.

§ PROOFS FOR LEMMAS, THEOREMS, PROPOSITIONS AND COROLLARIES

§.§ Proof of Lemma <ref>

Lemma <ref> (Neural Collapse of The Unhinged Loss).

For norm-bounded prototypes and features, i.e., $\|\vw_c\|_2\le E_1$ and $\|\vh_{i,c}\|_2\le E_2$, $\forall i\in[N], \forall c\in[C]$, the global minimizer of $\frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C L_{\gamma}(\mW\vh_{i,c},y_{i,c})$ implies neural collapse when $p\ge C-1$. More specifically, the global minimizer is uniquely obtained at $\frac{\vw_i^\top\vw_j}{\|\vw_i
\|_2\|\vw_j\|_2}=-\frac{1}{C-1}$, $\forall i\neq j$, $\frac{\vw_{y_{i,c}}^\top\vh_{i,c}}{\|\vw_{y_{i,c}}\|_2\|\vh_{i,c}\|_2}=1$,  $\|\vw_c\|_2=E_1$, and $\|\vh_{i,c}\|_2=E_2$, $\forall i\in[N],$ $\forall c\in[C]$.

The proof is based on lower bounding the objective $\frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C L_{\gamma}(\mW\vh_{i,c},y_{i,c})$ by a sequence of inequalities that holds if and only if the solution forms Neural Collapse [Papyan et al., 2020]. Let $\hat \vw=\frac{1}{C}\sum_{c=1}^C\vw_c$, according to the definition of $L_{SM}$, we have

\begin{equation}
\notag
\begin{aligned}
&\frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C L_{SM}(\mW\vh_{i,c},y_{i,c})\\
=& \frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C(-\vw_{y_{i,c}}^\top\vh_{i,c}+\gamma\sum_{j\neq y_{i,c}}\vw_j^\top\vh_{i,c})\\
\ge& \frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C(-E_1E_2+\gamma(C\hat{\vw}-\vw_{y_{i,c}})^\top\vh_{i,c})\\
\ge& -\frac{\gamma E_2}{CN}\sum_{i=1}^N\sum_{c=1}^C\|C\hat{\vw}-\vw_{y_{i,c}}\|_2-E_1 E_2\\
\ge& -\gamma E_2\sqrt{\frac{1}{CN}\sum_{i=1}^N\sum_{c=1}^C\|C\hat{\vw}-\vw_{y_{i,c}}\|_2^2}-E_1 E_2\\
=&-\gamma E_2\sqrt{\frac{1}{C}\sum_{c=1}^C \|\vw_{c}\|_2^2-C^2\|\hat{\vw}\|_2^2}-E_1E_2\\
\ge& -(1+\gamma) E_1 E_2
\end{aligned},
\end{equation}

where  the first and second inequalities are based on the facts that $\vw_{y_{i,c}}^\top\vh_{i,c}\le E_1E_2$ and $(C\hat{\vw}-\vw_{y_{i,c}})^\top\vh_{i,c}\ge -E_2\|C\hat{\vw}-\vw_{y_{i,c}}\|_2$, respectively. In the third equality, we used the Cauchy-Schwarz inequality, and the last inequality we use the facts that $\|\vw_c\|_2\le E_1$ and $\|\hat{\vw}\|_2\ge 0$.

According the above derivation, the equality holds if and only if $\forall i\in [N]$, $\forall c\in [C]$, $\vw_{y_{i,c}}^\top\vh_{i,c}= E_1E_2$, $(C\hat{\vw}-\vw_{y_{i,c}})^\top\vh_{i,c}= -E_2\|C\hat{\vw}-\vw_{y_{i,c}}\|_2$, $\|C\hat{\vw}-\vw_{c}\|_2=\|C\hat{\vw}-\vw_{C}\|_2$, $\|\vw_c\|_2= E_1$, and $\|\hat{\vw}\|_2= 0$. These equations can be simplified as $\frac{\vw_i^\top\vw_j}{\|\vw_i
|_2\|\vw_j\|_2}=-\frac{1}{C-1}$, $\forall i\neq j$, $\frac{\vw_{y_{i,c}}^\top\vh_{i,c}}{\|\vw_{y_{i,c}}\|_2\|\vh_{i,c}\|_2}=1$,  $\|\vw_c\|_2=E_1$, and $\|\vh_{i,c}\|_2=E_2$, $\forall i\in[N],$ $\forall c\in[C]$, which also implies neural collapse.

§.§ Proof of Theorem <ref>

In this section, we will provide the proof of Theorem <ref>. Our analysis will actually rely on the eigenvalues and eigenspaces of five subspaces $\gE_1^+$, $\gE_1^-$, $\gE_2^+$, $\gE_2^-$ and $\gE_3$ in Theorem <ref>. Their concrete projection operator can be found in <ref>. In the following, we show that these five subspaces are orthogonal:

The following five subspaces are orthogonal to each other and satisfy $\sR^{p\times C(N+1)}=\gE_1^+\oplus\gE_1^-\oplus\gE_2^+\oplus\gE_2^-\oplus \gE_3$:

\begin{equation}
\begin{aligned}
& \gE_1^\eps:=\{(\mH,\mW):\mH=\eps\cdot\tfrac{1}{\sqrt{N}}(\mW\otimes \1_N^\top),\mW\1_C=0, \mW\in\sR^{p\times C}\},\\
& \gE_2^\eps:=\{(\mH,\mW):\mH=\eps\cdot \tfrac{1}{\sqrt{N}}\vh\1^\top_{CN},\mW=\vh\1_C ^\top,\vh\in\mathbb{R}^p\},\\
& \gE_3:=\{(\mH,\mW): \mH(\mI_C\otimes\1_N)=0,\mW=0, \mH\in\sR^{p\times CN}\}.
\end{aligned}
\end{equation}

where $\eps\in\{\pm 1\}$, and $k\neq 0$.

For $(\mH_1,\mW_1)=(\frac{1}{\sqrt{N}}(\mW_1\otimes \1_N^\top), \mW_1)\in\gE_1^+$ and $(\mH_2,\mW_2)=(-\frac{1}{\sqrt{N}}(\mW_2\otimes \1_N^\top), \mW_2)\in\gE_1^-$, we have

\begin{equation}
\notag
\mH_1\mH_2^\top + \mW_1\mW_2^\top=-\frac{1}{N}(\mW_1\otimes \1_N^\top)(\mW_2^\top\otimes \1_N) + \mW_1\mW_2^\top=0.
\end{equation}

For $(\mH_1,\mW_1)=(\frac{1}{\sqrt{N}}\vh_1\1_{CN}^\top, \vh_1\1_C^\top)\in\gE_2^+$ and $(\mH_2,\mW_2)=(-\frac{1}{\sqrt{N}}\vh_2\1_{CN}^\top, \vh_2\1_C^\top)\in\gE_2^-$, we have

\begin{equation}
\notag
\mH_1\mH_2^\top + \mW_1\mW_2^\top=-\frac{1}{N}\vh_1\1_{CN}^\top \1_{CN}\vh_2^\top + \vh_1\1_{C}^\top \1_{C}\vh_2^\top=0.
\end{equation}

For $(\mH_1,\mW_1)=(\eps_1\frac{1}{\sqrt{N}}(\mW_1\otimes \1_N^\top), \mW_1)\in\gE_1^\eps$ and $(\mH_2,\mW_2)=(\eps_2\frac{1}{\sqrt{N}}\vh_2\1_{CN}^\top, \vh_2\1_C^\top)\in\gE_2^\eps$, since $\mW_1\1_C=0$ and $(\mW_1\otimes\1_N^\top)\1_{CN}=N\mW_1\1_C=0$, we have

\begin{equation}
\notag
\mH_1\mH_2^\top + \mW_1\mW_2^\top=\frac{\eps_1\eps_2}{N}(\mW_1\otimes\1_N^\top)\1_{CN}\vh_2^\top+\mW_1\1_C\vh_2^\top =0.
\end{equation}

For $(\mH_1,\mW_1)=(\eps\frac{1}{\sqrt{N}}(\mW_1\otimes \1_N^\top), \mW_1)\in\gE_1^\eps$ and $(\mH_2,\mW_2)=(\mH_2,0)\in\gE_3$, we have

\begin{equation}
\notag
\mH_1\mH_2^\top + \mW_1\mW_2^\top=\frac{\eps}{\sqrt{N}}(\mW_1\otimes \1_N^\top)\mH_2^\top=\frac{\eps}{\sqrt{N}}\mW_1(\mH_2(\mI_C\otimes\1_N))^\top = 0.
\end{equation}

For $(\mH_1,\mW_1)=(\eps\frac{1}{\sqrt{N}}\vh_1\1_{CN}^\top,\vh_1\1_C^\top)\in\gE_2^\eps$ and $(\mH_2,\mW_2)=(\mH_2,0)\in\gE_3$, since $\mH_2(\mI_C\otimes\1_N)=0$, we have

\begin{equation}
\notag
\mH_1\mH_2^\top + \mW_1\mW_2^\top=\frac{\eps}{\sqrt{N}}\vh_1\1_{CN}^\top\mH_2^\top=\frac{\eps}{\sqrt{N}}\vh_1(\mH_2(\mI_C\otimes\1_N)
\1_C^\top)^\top = 0
\end{equation}

To sum up, we prove that the five subspaces $\gE_1^+,\gE_1^-,\gE_2^+,\gE_2^-,\gE_3$ are orthogonal to each other. Moreover, we have

\mathop{\text{dim}} \gE_1^+=\mathop{\text{dim}} \gE_1^-=p(C-1), \quad \mathop{\text{dim}}\gE_2^+= \mathop{\text{dim}}\gE_2^-=p,\quad \mathop{\text{dim}}\gE_3=pC(N-1).

Since these dimensions sum to $pC(N+1)=\mathop{\text{dim}}(\mathbb{R}^{p\times CN}\oplus\mathbb{R}^{p\times C})$, then $\sR^{p\times C(N+1)}=\gE_1^+\oplus\gE_1^-\oplus\gE_2^+\oplus\gE_2^-\oplus \gE_3$.

Theorem <ref> (Dynamics of Features, Class Prototypes and Biases without Constraints).

Consider the continual gradient flow in gradient-flow, let $\mZ(t)=(\mH(t),\mW(t))$, if $\eta_1(t_1)\eta_2(t_2)=\eta_1(t_2)\eta_2(t_1)$ for any $t_1,t_2\ge0$, we have the following closed-form dynamics

\begin{equation}
\begin{aligned}
\mZ(t)=&\Pi_1^+\mZ_0 \left(\alpha_1^+(t)\mC(t)+\beta_1^+(t)\mI_{C(N+1)}\right)\\
+&\Pi_1^-\mZ_0 \left(\alpha_1^-(t)\mC(t)+\beta_1^-(t)\mI_{C(N+1)}\right)\\
+&\Pi_2^-\mZ_0 \left(\alpha_2^-(t)\mC(t)+\beta_2^-(t)\mI_{C(N+1)}\right)+\Pi_3\mZ_0,
\end{aligned}
\end{equation}

\begin{equation}
\vb(t)=\vb_0 +\frac{(1+\gamma-\gamma C)\zeta_2(t)}{C} \1_C,
\end{equation}

where $\alpha_1^\eps$, $\alpha_2^\eps$, $\beta_1^\eps$ and $\beta_2^\eps$ for $\eps\in\{\pm\}$ are the scalars that only depend on $C$, $N$, $\gamma$, $\eta_1$ and $\eta_2$ (where the detailed forms of these scalars can be seen in the appendix), $\mZ_0=(\mH_0,\mW_0)$, $\mC(t)=\left(\begin{smallmatrix}
\zeta_1(t)\mI_{CN} & 0\\
0 & \zeta_2(t)\mI_C\\
\end{smallmatrix}\right)$, $\zeta_1(t)=\int_0^t\eta_1(\tau)\mathrm{d}\tau$, $\zeta_2(t)=\int_0^t\eta_2(\tau)\mathrm{d}\tau$, $\Pi_1^{\eps}$, $\Pi_2^{\eps}$ and $\Pi_3$ for $\eps\in\{\pm\}$ are orthogonal projection operators onto the following respective eigenspaces:

\begin{equation}
\begin{aligned}
& \gE_1^\eps:=\{(\mH,\mW):\mH=\eps\cdot\tfrac{1}{\sqrt{N}}(\mW\otimes \1_N^\top),\mW\1_C=0\},\\
& \gE_2^\eps:=\{(\mH,\mW):\mH=\eps\cdot \tfrac{1}{\sqrt{N}}\vh\1^\top_{CN},\mW=\vh\1_C ^\top,\vh\in\mathbb{R}^p\},\\
& \gE_3:=\{(\mH,\mW): \mH(\mI_C\otimes\1_N)=0,\mW=0\}.
\end{aligned}
\end{equation}

Writing $\mZ(t)=(\mH(t),\mW(t))$, then the unsolved portion of the system is given by

\begin{equation}
\mZ'(t)=\mZ(t)\begin{pmatrix}0 & \mM^\top\\\mM & 0\\\end{pmatrix}\begin{pmatrix}
\eta_1(t)\mI_{CN} & 0\\
0 & \eta_2(t)\mI_C\\
\end{pmatrix},
\end{equation}

where $\mM=\frac{1}{CN}((1+\gamma)(\mI_C\otimes \1^\top_{N})-\gamma\1_C\1_{CN}^\top)$.

Let $\mA(t)=\begin{pmatrix}0 & \mM^\top\\\mM & 0\\\end{pmatrix} \begin{pmatrix}\eta_1(t)\mI_{CN} & 0\\0 & \eta_2(t)\mI_C\\\end{pmatrix}=\begin{pmatrix}0 & \eta_2(t)\mM^\top\\\eta_1(t)\mM & 0\\\end{pmatrix}$, then the equation above can be reformulated as the initial-value problem associated with the linear ordinary differential equation:

\begin{equation}
\mZ'(t)=\mZ(t)\mA(t),\quad \mZ(0)=\mZ_0.
\end{equation}

For any $t_1,t_2$, we have the matrix commutator of $\mA(t_1)$ and $\mA(t_2)$

\begin{aligned}
=&\begin{pmatrix}0 & \eta_2(t_1)\mM^\top\\\eta_1(t_1)\mM & 0\end{pmatrix}\begin{pmatrix}0 & \eta_2(t_2)\mM^\top\\\eta_1(t_2)\mM & 0\end{pmatrix}-\mA(t_2)\mA(t_1)\\
=&\begin{pmatrix}\eta_2(t_1)\eta_1(t_2)\mM^\top\mM & 0\\ 0 & \eta_1(t_1)\eta_2(t_2)\mM\mM^\top\end{pmatrix}-\mA(t_2)\mA(t_1)\\
=&\begin{pmatrix}(\eta_2(t_1)\eta_1(t_2)-\eta_2(t_2)\eta_1(t_1))\mM^\top\mM & 0\\ 0 & (\eta_1(t_1)\eta_2(t_2)-\eta_2(t_1)\eta_1(t_2))\mM\mM^\top\end{pmatrix}\\
\end{aligned}

where the last equality is based on the fact that $\eta_2(t_1)\eta_1(t_2)=\eta_2(t_2)\eta_1(t_1)$. Therefore, according to Magnus approach, we have

\begin{equation}
\mZ(t)=\mZ_0\exp\left(\int_0^t \mA(\tau)\mathrm{d}\tau\right)=\mZ_0\exp\begin{pmatrix}
0 & \int_0^t\eta_2(\tau) \mathrm{d}\tau\mM^\top\\
\int_0^t\eta_1(\tau) \mathrm{d}\tau\mM & 0\\
\end{pmatrix}.
\end{equation}

Let $\zeta_1(t)=\int_0^t\eta_1(\tau)\mathrm{d}\tau$, $\zeta_2(t)=\int_0^t\eta_2(\tau)\mathrm{d}\tau$, $\mB=\begin{pmatrix}
0 & \mM^\top\\
\mM & 0\\
\end{pmatrix}$, $\mC(t)=\begin{pmatrix}
\zeta_1(t)\mI_{CN} & 0\\
0 & \zeta_2(t)\mI_C\\
\end{pmatrix}$, and $\mL(t)=\mB\mC(t)$,  we have

\begin{equation}
\mZ(t)=\mZ_0\exp(\mL(t))=\mZ_0\sum_{k=0}^\infty \frac{(\mL(t))^k}{k!}.
\end{equation}

Moreover, we have

\begin{equation}
\begin{aligned}
0 & \zeta_2(t)\mM^\top\\
\zeta_1(t)\mM & 0\\
\end{pmatrix}\begin{pmatrix}
0 & \zeta_2(t)\mM^\top\\
\zeta_1(t)\mM & 0\\
\end{pmatrix}\\
\zeta_1(t)\zeta_2(t)\mM^\top\mM & 0\\
0 & \zeta_1(t)\zeta_2(t)\mM\mM^\top\\
\end{pmatrix}\\
\end{aligned}
\end{equation}

thus we obtain

\begin{equation}
\begin{aligned}
\mZ(t)=&\mZ_0\left(\sum_{k=0}^\infty \frac{(\mL(t))^{2k+1}}{(2k+1)!}+\sum_{k=0}^\infty \frac{(\mL(t))^{2k}}{(2k)!}\right)\\
=&\mZ_0\left(\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k+1}\mC(t)}{(2k+1)!}+\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k}}{(2k)!}\right)\\
=&\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k+1}\mC(t)}{(2k+1)!}+\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k}}{(2k)!}
\end{aligned}
\end{equation}

Looking at the above equation, we just need to analyze the eigenspaces and eigenvalues of $\mB$.

Considering the following five subspaces:

\begin{equation}
\begin{aligned}
& \gE_1^\eps:=\{(\mH,\mW):\mH=\eps\cdot\tfrac{1}{\sqrt{N}}(\mW\otimes \1_N^\top),\mW\1_C=0\},\\
& \gE_2^\eps:=\{(\mH,\mW):\mH=\eps\cdot \tfrac{1}{\sqrt{N}}\vh\1^\top_{CN},\mW=\vh\1_C ^\top,\vh\in\mathbb{R}^p\},\\
& \gE_3:=\{(\mH,\mW): \mH(\mI_C\otimes\1_N)=0,\mW=0\}.
\end{aligned}
\end{equation}

where $\eps\in\{\pm\}$. According to Lemma <ref>, these five subspaces are orthogonal to each other and satisfy $\sR^{p\times C(N+1)}=\gE_1^+\oplus\gE_1^-\oplus\gE_2^+\oplus\gE_2^-\oplus \gE_3$.

In the following, we will prove that $\gE_{1}^\eps$, $\gE_{2}^\eps$ and $\gE_3$ are five eigenspaces of $\mB$. More specifically, each nonzero member of each claimed eigenspace is an eigenvector, and the claimed eigenspaces have distinct eigenvalues.

Note that for $(\mH,\mW)\in \sR^{p\times CN}\oplus\sR^{p\times C}$, we have $(\mH,\mW)\mB=(\mW\mM^\top,\mH\mW)$.

For $(\mH,\mW)\in\gE_1^\eps$, we have $\mH=\frac{\eps}{\sqrt{N}}\mW\otimes \1_N^\top$ and  $\mW\1_C=0$, thus

\begin{equation}
\notag
\begin{aligned}
\mW\mM&=\tfrac{1}{CN}\mW((1+\gamma)(\mI_C\otimes \1^\top_{N})-\gamma\1_C\1_{CN}^\top)\\
&=\tfrac{(1+\gamma)}{CN}\mW\otimes \1_N^\top\\
&=\tfrac{\eps(1+\gamma)}{C\sqrt{N}} \mH,\\
\mH\mM^\top&=\tfrac{1}{CN}[\eps\cdot\tfrac{1}{\sqrt{N}}(\mW\otimes \1_N^\top)][(1+\gamma)(\mI_C\otimes\1^\top_{N})-\gamma\1_C\1_{CN}^\top]^\top\\
&=\tfrac{\eps}{CN\sqrt{N}}[(\mW\otimes \1_N^\top)((1+\gamma)(\mI_C\otimes \1_N)-\gamma(\1_C\otimes \1_N)\1_C^\top)]\\
&=\tfrac{\eps(1+\gamma)}{C\sqrt{N}} \mW,
\end{aligned}
\end{equation}

i.e., $(\mH,\mW)$ is an eigenvector of $\mB$ with eigenvalue $\tfrac{\eps(1+\gamma)}{C\sqrt{N}}$.

For $(\mH,\mW)\in\gE_2^\eps$, we have $\mH=\frac{\eps}{\sqrt{N}}\vh\1_{CN}^\top$ and $\mW=\vh\1_C^\top$, thus

\begin{equation}
\notag
\begin{aligned}
\mW\mM&=\tfrac{1}{CN}\vh\1_C^\top ((1+\gamma)(\mI_C\otimes\1^\top _{N})-\gamma\1_C\1_{CN}^\top )\\
&=\tfrac{1}{CN}\vh((1+\gamma)\1_C(\mI_C \otimes \1_N^\top )-\gamma C\1_{CN}^\top )\\
&=\tfrac{(1+\gamma-\gamma C)}{CN}\vh \1_{CN}^\top \\
&=\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}} \mH,\\
\mH\mM^\top &=\tfrac{1}{CN\sqrt{N}}(\eps\cdot \vh \1_{CN}^\top )((1+\gamma)(\mI_C\otimes\1^\top _{N})-\gamma\1_C\1_{CN}^\top)^\top \\
&=\tfrac{\eps}{CN}\vh((1+\gamma)\1_{CN}^\top(\mI_C\otimes \1_N)-\gamma\1_{CN\sqrt{N}}^\top\1_{CN}\1_C^\top)\\
&=\tfrac{\eps}{CN\sqrt{N}}\vh((1+\gamma)N\1_{C}^\top-\gamma CN\1_C^\top) \\
&=\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}}\vh\1_C ^\top \\
&=\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}}\mW,
\end{aligned}
\end{equation}

i.e., $(\mH,\mW)$ is an eigenvector of $\mB$ with eigenvalue $\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}}$.

For $(\mH,\mW)\in\gE_3$, we have $\mH(\mI_C\otimes\1_N)=0$ and $\mW=0$, thus

\begin{equation}
\notag
\begin{aligned}
\mW\mM&=\tfrac{1}{CN}\cdot0((1+\gamma)(\mI_C\otimes\1^\top _{N})-\gamma\1_C\1_{CN}^\top)=0,\\
\mH\mM^\top &=\tfrac{1}{CN}\mH((1+\gamma)(\mI_C\otimes\1^\top_{N})-\gamma\1_C\1_{CN}^\top )^\top \\
&=\tfrac{1}{CN}\mH((1+\gamma)(\mI_C\otimes\1_{N})-\gamma\1_{CN}\1_C^\top ) \\
&=-\tfrac{\gamma}{CN}\mH (\mI_C\otimes\1_{N})\1_C\1_C^\top \\
\end{aligned}
\end{equation}

i.e., $(\mH,\mW)$ is an eigenvector of $\mB$ with eigenvalue $0$.

Overall, letting $\Pi_i^\eps$ denote orthogonal projection onto $\gE_i^\eps$, we have the spectral decomposition

\begin{equation}
\mB=\tfrac{1}{C\sqrt{N}}\left[(1+\gamma)(\Pi_1^+-\Pi_1^-)+(1+\gamma-\gamma C)(\Pi_2^+-\Pi_2^-)\right].
\end{equation}

We then provide the concrete formulation of $\mZ(t)=\mZ_0\exp(\mL(t))$ by the orthogonal projection of $\mZ_0$ onto each eigenspace of $\mB$, i.e.,

\mZ_0=\Pi_1^+\mZ_0+\Pi_1^-\mZ_0+\Pi_2^+\mZ_0+\Pi_2^-\mZ_0+\Pi_3\mZ_0.

Decomposition along $\Pi_1^\eps \mZ_0$.
First, $\Pi_1^\eps \mZ_0 \mB=\frac{\eps(1+\gamma)}{C\sqrt{N}}\Pi_1^\eps\mZ_0$, so $ \Pi_1^\eps \mZ_0 \mB^k = \left(\frac{\eps(1+\gamma)}{C\sqrt{N}}\right)^k\Pi_1^\eps\mZ_0$ for $k\ge0$, then

\begin{equation}
\begin{aligned}
&\Pi_1^\eps\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k+1}\mC(t)}{(2k+1)!}\\
=&\Pi_1^\eps\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\left(\frac{\eps(1+\gamma)}{C\sqrt{N}}\right)^{2k+1}\mC(t)}{(2k+1)!}\\
=&\frac{\Pi_1^\eps\mZ_0 \mC(t)}{\sqrt{\zeta_1(t)\zeta_2(t)}}\sum_{k=0}^\infty \frac{\left(\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)^{2k+1}}{(2k+1)!}\\
=&\frac{\Pi_1^\eps\mZ_0 \mC(t)}{2\sqrt{\zeta_1(t)\zeta_2(t)}}\left(e^{\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}-e^{-\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}\right),
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
& \Pi_1^\eps\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k}}{(2k)!}\\
=&\Pi_1^\eps\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\left(\frac{\eps(1+\gamma)}{C\sqrt{N}}\right)^{2k}}{(2k)!}\\
\end{aligned}
\end{equation}

which is based on the facts that $\frac{e^x-e^{-x}}{2}=\sum_{k=0}^\infty \frac{x^{2k+1}}{(2k+1)!}$ and $\frac{e^x+e^{-x}}{2}=\sum_{k=0}^\infty \frac{x^{2k}}{(2k)!}$. Thus we have

\begin{equation}
\Pi_1^\eps\mZ_0 \exp(\mL(t))=\Pi_1^\eps\mZ_0(\alpha_1^\eps\mC(t)+\beta_1^\eps\mI_{C(N+1)}),
\end{equation}

\begin{equation}
\begin{aligned}
&\alpha_1^\eps(t)=\frac{\exp\left(\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)-\exp\left(-\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)}{2\sqrt{\zeta_1(t)\zeta_2(t)}},\\ &\beta_1^\eps(t)=\frac{\exp\left(\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)+\exp\left(-\frac{\eps(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)}{2}.
\end{aligned}
\end{equation}

Decomposition along $\Pi_2^\eps \mZ_0$. Similarly, for $\Pi_2^\eps\mZ_0 \mB=\frac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}}\Pi_2^\eps$, we have

\Pi_2^\eps\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k+1}\mC(t)}{(2k+1)!}=\frac{\Pi_2^\eps\mZ_0 \mC(t)}{2\sqrt{\zeta_1(t)\zeta_2(t)}}\left(e^{\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}-e^{-\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}\right),

\Pi_2^\eps\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k}}{(2k)!}=\frac{\Pi_2^\eps\mZ_0}{2}\left(e^{\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}+e^{-\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}\right).

Thus we have

\begin{equation}
\Pi_2^\eps\mZ_0 \exp(\mL(t))=\Pi_2^\eps\mZ_0(\alpha_2^\eps\mC(t)+\beta_2^\eps\mI_{C(N+1)}),
\end{equation}

\begin{equation}
\begin{aligned}
&\alpha_2^\eps(t)=\frac{\exp\left(\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)-\exp\left(-\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)}{2\sqrt{\zeta_1(t)\zeta_2(t)}},\\
&\beta_2^\eps(t)=\frac{\exp\left(\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)+\exp\left(-\frac{\eps(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\right)}{2}.
\end{aligned}
\end{equation}

Decomposition along $\Pi_3 \mZ_0$. Since each vector in $\gE_3$ is a eigenvector of $\mB$ with eigenvalue 0, then we have

\begin{equation}
\Pi_3\mZ_0\left(\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k+1}\mC(t)}{(2k+1)!}+\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k}}{(2k)!}\right)=\Pi_3\mZ_0
\end{equation}

Note that $\gE_1^+$, $\gE_1^-$, $\gE_2^+$, $\gE_2^-$ and $\gE_3$ are orthogonal subspace of $\sR^{p\times CN}\oplus\sR^{C\times p}$, thus

\begin{equation}
\begin{aligned}
\mZ(t)=&\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k+1}\mC(t)}{(2k+1)!}+\mZ_0\sum_{k=0}^\infty \frac{(\zeta_1(t)\zeta_2(t))^k\mB^{2k}}{(2k)!}\\
=&\Pi_1^+\mZ_0 \left(\alpha_1^+(t)\mC(t)+\beta_1^+(t)\mI_{C(N+1)}\right)+\Pi_1^-\mZ_0 \left(\alpha_1^-(t)\mC(t)+\beta_1^-(t)\mI_{C(N+1)}\right)+\Pi_3\mZ_0\\
&+\Pi_2^+\mZ_0\left(\alpha_2^+(t)\mC(t)+\beta_2^+(t)\mI_{C(N+1)}\right)+\Pi_2^-\mZ_0 \left(\alpha_2^-(t)\mC(t)+\beta_2^-(t)\mI_{C(N+1)}\right)
\end{aligned}
\end{equation}

Moreover, since $\vb'(t)=-\eta_2(t)\tfrac{\gamma C-\gamma-1}{C}\1_C$, we obtain

\begin{equation}
\vb(t)=\vb(0)+\int_0^t-\eta_2(\tau)\frac{\gamma C-\gamma-1}{C} \mathrm{d}\tau \1_C=\vb_0+\frac{(1+\gamma-\gamma C)\zeta_2(t)}{C} \1_C,
\end{equation}

with $\zeta_2(t)=\int_0^t\eta_2(\tau)\mathrm{d}\tau$.

§.§ Proof of Corollary <ref>

Corollary <ref>

Under the conditions and notation of Theorem <ref>, let $s=\frac{\eta_1(0)}{\eta_2(0)}$, if $0<\gamma < \frac{2}{C-2}$ (where $C>2$) or $C=2$, and $\lim_{t\rightarrow\infty}\zeta_1(t)=\infty$, then the gradient flow (as in <ref>) will behave as:

\begin{equation}
e^{-\frac{(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}}\mZ(t)=\overline{\mZ} +\bm{\Delta}(t),
\end{equation}

where $\overline{\mZ}=\left(\tfrac{1+\sqrt{s}}{2}\mH_1^++\tfrac{1-\sqrt{s}}{2}\mH_1^-, \tfrac{1+\sqrt{s}}{2\sqrt{s}}\mW_1^+-\tfrac{1-\sqrt{s}}{2\sqrt{s}}\mW_1^-\right)$,  $(\mH_1^+,\mW_1^+)=\Pi_1^+\mZ_0$, $(\mH_1^-,\mW_1^-)=\Pi_1^-\mZ_0$, and the residual term $\bm{\Delta}(t)$ decreases as $\|\bm{\Delta}(t)\|=O\left(e^{\frac{\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\cdot\max\{-\gamma C,(C-2)\gamma-2\}}\right)$, and so the normalized $\mZ(t)$ converges to $\frac{\overline{\mZ}}{\|\overline{\mZ}\|}$ in

\begin{equation}
\left\|\frac{\mZ(t)}{\|\mZ(t)\|}-\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\right\|=O\left(e^{\frac{\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}\cdot\max\{-\gamma C,(C-2)\gamma-2\}}\right),
\end{equation}

which further indicates $\lim_{t\rightarrow\infty}\frac{\mZ(t)}{\|\mZ(t)\|}\in\gE$.
Moreover, if $\gamma\neq \tfrac{1}{C-1}$, then $\lim_{t\rightarrow \infty}\frac{\max_i b_i(t)}{\min_i b_i(t)}=1$.

Let $(\mH_1^\eps,\mW_1^\eps)=\Pi_1^\eps\mZ_0$, $(\mH_2^\eps,\mW_2^\eps)=\Pi_2^\eps\mZ_0$ and $(\mH_3,\mW_3)=\Pi_3\mZ_0$ for $\eps\in\{\pm1\}$, according to Theorem <ref>, we have

\begin{equation}
\label{decomposition-H-W}
\begin{aligned}
&\mH(t)=\sum_{\substack{i\in\{1,2\}\\\eps\in\{\pm\}}} (\alpha_i^\eps(t)\zeta_1(t) + \beta_i^\eps(t))\mH_i^\eps + \mH_3,\\
&\mW(t)=\sum_{\substack{i\in\{1,2\}\\\eps\in\{\pm\}}} (\alpha_i^\eps(t)\zeta_2(t) + \beta_i^\eps(t))\mW_i^\eps + \mW_3.
\end{aligned}
\end{equation}

Since $\eta_1(t_1)\eta_2(t_2)=\eta_1(t_2)\eta_2(t_1)$, then $\eta_1(t)=\frac{\eta_1(0)}{\eta_2(0)}\eta_2(t)$. Let $s=\frac{\eta_1(0)}{\eta_2(0)}$, $p(t)=\frac{(1+\gamma)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}$, and $q(t)=\frac{(1+\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}$, we have $\zeta_1(t)=\int_0^t\eta_1(\tau)\mathrm{d}\tau=s\int_0^t\eta_2(\tau)\mathrm{d}\tau= s\zeta_2(t)$, and then

\begin{equation}
\notag
\begin{aligned}
&\alpha_1^\eps(t)\zeta_1(t) + \beta_1^\eps(t)
=\tfrac{1+\sqrt{s}}{2}e^{\eps p(t)}+\tfrac{1-\sqrt{s}}{2}e^{-\eps p(t)}=\tfrac{1+\eps\sqrt{s}}{2}e^{p(t)}+O(e^{-p(t)}),\\
&\alpha_1^\eps(t)\zeta_2(t) + \beta_1^\eps(t)
=\tfrac{1+\sqrt{s}}{2\sqrt{s}}e^{\eps p(t)} -\tfrac{1-\sqrt{s}}{2\sqrt{s}}e^{-\eps p(t)}=\tfrac{\eps+\sqrt{s}}{2\sqrt{s}}e^{p(t)}+O(e^{-p(t)}),\\
&\alpha_2^\eps(t)\zeta_1(t) + \beta_2^\eps(t)
=\tfrac{1+\sqrt{s}}{2}e^{\eps q(t)}+\tfrac{1-\sqrt{s}}{2}e^{-\eps q(t)}=\tfrac{1+\eps\sqrt{s}}{2}e^{q(t)}+O(e^{-q(t)}),\\
&\alpha_2^\eps(t)\zeta_2(t) + \beta_2^\eps(t)
=\tfrac{1+\sqrt{s}}{2\sqrt{s}}e^{\eps q(t)} -\tfrac{1-\sqrt{s}}{2\sqrt{s}}e^{-\eps q(t)}=\tfrac{\eps+\sqrt{s}}{2\sqrt{s}}e^{q(t)}+O(e^{-q(t)}).\\
\end{aligned}
\end{equation}

Since $0<\gamma<\frac{2}{C-2}$ ( where $C>2$) or $C=2$, then we have $p(t)-q(t)=\frac{\gamma C\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}>0$, $p(t)+q(t)=\frac{(2+2\gamma-\gamma C)\sqrt{\zeta_1(t)\zeta_2(t)}}{C\sqrt{N}}>0$, and substitute these results into <Ref> to obtain

\begin{equation}
\begin{aligned}
\end{aligned}
\end{equation}

where $\|\bm{\Delta}_1(t)\|=O(e^{\max\{q(t)-p(t),-q(t)-p(t)\}}$ and $\|\bm{\Delta}_2(t)\|=O(e^{\max\{q(t)-p(t),-q(t)-p(t)\}}$. Therefore, we have

\begin{equation}
e^{-p(t)}\mZ(t)=\left(\tfrac{1+\sqrt{s}}{2}\mH_1^++\tfrac{1-\sqrt{s}}{2}\mH_1^-, \tfrac{1+\sqrt{s}}{2\sqrt{s}}\mW_1^+-\tfrac{1-\sqrt{s}}{2\sqrt{s}}\mW_1^-\right) +\bm{\Delta}(t),
\end{equation}

where $\bm{\Delta}(t)=(\bm{\Delta}_1(t),\bm{\Delta}_2(t))$ and $\|\bm{\Delta}(t)\|\le \|\bm{\Delta}_1(t)\|+\|\bm{\Delta}_2(t)\|=O(e^{\max\{q(t)-p(t),-q(t)-p(t)\}}$.

Let $\overline{\mZ}=\left(\tfrac{1+\sqrt{s}}{2}\mH_1^++\tfrac{1-\sqrt{s}}{2}\mH_1^-, \tfrac{1+\sqrt{s}}{2\sqrt{s}}\mW_1^+-\tfrac{1-\sqrt{s}}{2\sqrt{s}}\mW_1^-\right)$, we have

\begin{equation}
\begin{aligned}
\left\|\frac{\mZ(t)}{\|\mZ(t)\|}-\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\right\|=\left\|\frac{\overline{\mZ}+\bm{\Delta}(t)}{\|\overline{\mZ}+\bm{\Delta}(t)\|}-\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\right\|\le\frac{2\|\overline{\mZ}\|\|\bm{\Delta}(t)\|}{\|\overline{\mZ}+\bm{\Delta}(t)\|\|\overline{\mZ}\|}= \frac{2\|\bm{\Delta}(t)\|}{\|\overline{\mZ}+\bm{\Delta}(t)\|},
\end{aligned}
\end{equation}

thus $\left\|\frac{\mZ(t)}{\|\mZ(t)\|}-\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\right\|=O(e^{\max\{q(t)-p(t),-q(t)-p(t)\}}$, and further $\lim_{t\rightarrow\infty}\frac{\mZ(t)}{\|\mZ(t)\|}=\frac{\overline{\mZ}}{\|\overline{\mZ}\|}$ when $\lim_{t\rightarrow \infty}\zeta_1(t)=\infty$.

According to the definition of $\gE_1$ and $\gE_2$, we have

\begin{equation}
\begin{aligned}
\overline{\mZ}=\left(\tfrac{\sqrt{s}}{\sqrt{N}}\left(\tfrac{1+\sqrt{s}}{2\sqrt{s}}\mW_1^+-\tfrac{1-\sqrt{s}}{2\sqrt{s}}\mW_1^-\right)\otimes \1_N^\top, \tfrac{1+\sqrt{s}}{2\sqrt{s}}\mW_1^+-\tfrac{1-\sqrt{s}}{2\sqrt{s}}\mW_1^- \right),
\end{aligned}
\end{equation}

thus we have $\overline{\mZ}\in\gE$, then $\lim_{t\rightarrow\infty}\frac{\mZ(t)}{\|\mZ(t)\|}=\frac{\overline{\mZ}}{\|\overline{\mZ}\|}\in\gE$.

Moreover, we have $\vb(t)=\vb_0 +\frac{(1+\gamma-\gamma C)\zeta_2(t)}{C} \1_C$, then $\forall i,j$,

\begin{equation}
\notag
\lim_{t\rightarrow \infty}\frac{b_i(t)}{b_j(t)}=\lim_{t\rightarrow \infty}\frac{b_{i}(0)+\frac{1+\gamma -\gamma C}{C}\zeta_2(t)}{b_{j}(0)+\frac{1+\gamma -\gamma C}{C}\zeta_2(t)}=1,
\end{equation}

thus $\lim_{t\rightarrow \infty}\frac{\max_i b_i(t)}{\max_i b_i(t)}$=1.

§.§ Proof of Theorem <ref>

Theorem <ref> (Dynamics of Features and Prototypes Under Weight Decay)

Consider the continual gradient flow in regularized-gradient-flow, let $\mZ(t)=(\mH(t),\mW(t))$. If $\eta_1(t_1)\eta_2(t_2)=\eta_1(t_2)\eta_2(t_1)$ for any $t_1,t_2\ge0$, we have the following closed-form dynamics:

\begin{equation}
\begin{aligned}
\mZ(t)=&\Pi_1^+\mZ_0\begin{pmatrix}a_1^+(t)\mI_{CN} & 0\\0 & b_1^+(t)\mI_{C}\end{pmatrix}+\Pi_1^-\mZ_0\begin{pmatrix}a_1^-(t)\mI_{CN} & 0\\0 & b_1^-(t)\mI_{C}\end{pmatrix}+\\
&\Pi_2^+\mZ_0\begin{pmatrix}a_2^+(t)\mI_{CN} & 0\\0 & b_2^+(t)\mI_{C}\end{pmatrix}+\Pi_2^-\mZ_0\begin{pmatrix}a_2^-(t)\mI_{CN} & 0\\0 & b_2^-(t)\mI_{C}\end{pmatrix}+\\
&\Pi_3\mZ_0\begin{pmatrix}a_3(t)\mI_{CN} & 0\\0 & b_3(t)\mI_{C}\end{pmatrix}
\end{aligned}
\end{equation}

\begin{equation}
\vb(t)=\phi(t)\left(\vb_0+\tfrac{1+\gamma -\gamma C}{C}\psi(t)\1_C\right),
\end{equation}

where $\Pi_1^+\mZ_0$, $\Pi_1^-\mZ_0$, $\Pi_1^+\mZ_0$, $\Pi_1^-\mZ_0$, and $\Pi_3\mZ_0$ follow the definition in Theorem <ref>, $a_1^\eps$, $a_2^\eps$, $b_1^\eps$, $b_2^\eps$, $a_3$, and $b_3$ for $\eps\in\{\pm\}$ are the scalars that depend only on $C$, $N$, $\gamma$, $\lambda_1$, $\lambda_2$, $\eta_1$, and $\eta_2$ (where the detailed forms can be seen in <ref>), $\phi(t)=\exp(-\lambda\int_0^t\eta_2(\tau)\mathrm{d}\tau)$, and $\psi(t)=\int_0^t\zeta_2(\tau)\exp(\lambda\int_0^\tau\eta_2(s)\mathrm{d}s)\mathrm{d}\tau$.

According to the gradient flow in <ref> and the notations in the proof of Theorem <ref>, we have:

\begin{equation}
\mZ'(t)=\mZ(t)\mA(t)-\mZ(t)\begin{pmatrix}
\lambda_1\eta_1(t)\mI_{CN} & 0\\
0 & \lambda_2\eta_2(t)\mI_{C}
\end{pmatrix},
\end{equation}

i.e., $\mZ'(t)=\mZ(t)\mA_\lambda(t)$, where $\mA_\lambda(t)=\mA(t)-\bm{\Lambda}(t)$ and $\bm{\Lambda}(t)=\left(\begin{smallmatrix}
\lambda_1\eta_1(t)\mI_{CN} & 0\\
0 & \lambda_2\eta_2(t)\mI_{C}
\end{smallmatrix}\right)$.

For any $t_1,t_2$, we have the matrix commutator of $\mA_\lambda(t_1)$ and $\mA_\lambda(t_2)$

\begin{equation}
\begin{aligned}
=& \bm{\Lambda}(t_2)\mA(t_1)-\bm{\Lambda}(t_1)\mA(t_2)+\mA(t_2)\bm{\Lambda}(t_1)-\mA(t_1)\bm{\Lambda}(t_2) \\
0 & \lambda_1[\eta_1(t_2)\eta_2(t_1)-\eta_1(t_1)\eta_2(t_2)]\mM^\top\\
\lambda_2[\eta_2(t_2)\eta_1(t_2)-\eta_2(t_1)\eta_1(t_2)]\mM & 0
\end{pmatrix}\\
\end{aligned}
\end{equation}

where the last equality is based on the fact that $\eta_2(t_1)\eta_1(t_2)=\eta_2(t_2)\eta_1(t_1)$. Therefore, according to Magnus approach, we have

\begin{equation}
\mZ(t)=\mZ_0\exp\left(\int_0^t \mA_\lambda(\tau)\mathrm{d}\tau\right)=\mZ_0\exp\begin{pmatrix}
-\lambda_1\zeta_1(t)\mI_{CN} & \zeta_2(t)\mM^\top\\
\zeta_1(t)\mM & -\lambda_2\zeta_2(t)\mI_{C}\\
\end{pmatrix},
\end{equation}

where $\zeta_1(t)=\int_0^t\eta(\tau)\mathrm{d}\tau$ and $\zeta_2(t)=\int_0^t\eta(\tau)\mathrm{d}\tau$.

Let $\mB_{\lambda} = \mB-\left(\begin{smallmatrix}\lambda_1\mI_{CN} & 0\\0 & \lambda_2\mI_{C}\end{smallmatrix}\right)$, we have

\begin{equation}
\mZ(t)=\mZ_0\exp(\mB_\lambda\mC(t))=\mZ_0\sum_{k=0}^\infty \frac{(\mB_\lambda\mC(t))^k}{k!}.
\end{equation}

We again consider the orthogonal decomposition of $\mZ_0$, i.e., $\mZ_0=(\Pi_1^++\Pi_1^-+\Pi_2^++\Pi_2^-+\Pi_3)\mZ_0$. As mentioned in the proof of Theorem <ref>, we have

\begin{equation}
\begin{aligned}
\Pi_1^\eps \mZ_0 \mB&=\tfrac{\eps(1+\gamma)}{C\sqrt{N}}\Pi_1^\eps\mZ_0,\\
\Pi_2^\eps\mZ_0 \mB&=\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}}\Pi_2^\eps\mZ_0, \text{ and }\\
\Pi_3\mZ_0\mB&=0.
\end{aligned}
\end{equation}

Therefore, for any $ \mD=(\mH,\mW)\in\{\Pi_1^\eps\mZ_0,\Pi_2^\eps\mZ_0,\Pi_3\mZ_0\}$ (where $\mH\in\sR^{p\times CN}$ and $\mW\in\sR^{p\times C}$) and the corresponding eigenvalue $\sigma\in\{\frac{\eps(1+\gamma)}{C\sqrt{N}}, \frac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}}, 0\}$, we have

\begin{equation}
\mD\mB=\sigma\mD,\ \mH\mM^\top = \sigma\mW,\text{ and } \mW\mM = \sigma\mH.
\end{equation}

In the following, we will prove that there exist two scalars $a(t)$ and $b(t)$, such that $\mD\exp(\mB_\lambda\mC(t))=\mD\begin{pmatrix}
a(t)\mI_{CN} & 0\\ 0 & b(t)\mI_C
\end{pmatrix}$ for any $\mD=(\mH,\mW)\in\{\Pi_1^\eps\mZ_0,\Pi_2^\eps\mZ_0,\Pi_3\mZ_0\}$.

First, we prove that $\mD(\mB_\lambda\mC(t))^k$ can be represented as $\mD(\mB_\lambda\mC(t))^k=(a_k(t)\mH,b_k(t)\mW)$ by induction, where $a_k(t),b_k(t)\in\sR$.

For $k=0$, we have $\mD(\mB_\lambda\mC(t))^0=(\mH,\mW)$, i.e., $a_0=b_0=1$. Assume that $\mD(\mB_\lambda\mC(t))^n=(a_n(t)\mH,b_n(t)\mW)$ for $k=n$. Then for $k=n+1$, we have

\begin{equation}
\begin{aligned}
=&(a_n(t)\mH,b_n(t)\mW)\begin{pmatrix}-\lambda_1\mI_{CN} & \mM^\top\\ \mM & -\lambda_2\mI_{C}\end{pmatrix} \begin{pmatrix}\zeta_1(t)\mI_{CN} & 0\\ 0 & \zeta_2(t)\mI_{C}\\\end{pmatrix}\\
=& (b_n(t)\mW\mM-\lambda_1 a_n(t)\mH, a_n(t)\mH\mM^\top-\lambda_2 b_n(t)\mW)\begin{pmatrix}\zeta_1(t)\mI_{CN} & 0\\ 0 & \zeta_2(t)\mI_{C}\\\end{pmatrix}\\
=& (\zeta_1(t)(\sigma b_n(t)-\lambda_1 a_n(t))\mH, \zeta_2(t)(\sigma a_n(t)-\lambda_2 b_n(t))\mW)
\end{aligned},
\end{equation}

thus $a_{n+1}(t)=\zeta_1(t)(\sigma b_n(t)-\lambda_1 a_n(t))$ and $b_{n+1}(t)=\zeta_2(t)(\sigma a_n(t)-\lambda_2 b_n(t))$.

To sum up, we have shown by induction that $\mD(\mB_\lambda\mC(t))^k$ can be represented as $\mD(\mB_\lambda\mC(t))^k=(a_k(t)\mH,b_k(t)\mW)=\mD\begin{pmatrix}a_k(t)\mI_{CN}&0\\0 & b_k(t)\mI_C\end{pmatrix}$, and $\begin{pmatrix}
a_k(t) \\ b_k(t)
\end{pmatrix}$ satisfies

\begin{equation}
\begin{pmatrix}a_{k}(t) \\ b_{k}(t)\end{pmatrix}=\begin{pmatrix}-\lambda_1\zeta_1(t) & \sigma\zeta_1(t) \\ \sigma\zeta_2(t) & -\lambda_2\zeta_2(t)\end{pmatrix}\begin{pmatrix}a_{k-1}(t) \\ b_{k-1}(t)\end{pmatrix} \quad \text{with}\quad \begin{pmatrix}a_{0} \\ b_{0}\end{pmatrix}=\begin{pmatrix}1\\1\end{pmatrix},
\end{equation}

i.e., $\begin{pmatrix}a_{k}(t) \\ b_{k}(t)\end{pmatrix} = \left(\mS(\sigma,\lambda_1,\lambda_2,\zeta_1(t),\zeta_2(t))\right)^k\begin{pmatrix}1\\1\end{pmatrix}$, where $\mS(\sigma,\lambda_1,\lambda_2,\zeta_1,\zeta_2)=\begin{pmatrix}-\lambda_1\zeta_1 & \sigma\zeta_1 \\ \sigma\zeta_2 & -\lambda_2\zeta_2\end{pmatrix}$.

Therefore, we have

\begin{equation}
\mD\exp(\mB_\lambda\mC(t))=\mD\sum_{k=0}^\infty \frac{(\mB_\lambda\mC(t))^k}{k!}=\mD \begin{pmatrix}
a(t)\mI_{CN} & 0\\
0 & b(t)\mI_{C}\\
\end{pmatrix},
\end{equation}

with $a(t)=\sum_{k=0}^\infty\frac{a_k(t)}{k!}$ and $b(t)=\sum_{k=0}^\infty\frac{b_k(t)}{k!}$, i.e.,

\begin{equation}
\begin{aligned}
\begin{pmatrix}a(t) \\ b(t)\end{pmatrix}&=\sum_{k=0}^\infty \frac{1}{k!}\begin{pmatrix}a_k(t) \\ b_k(t)\end{pmatrix}=\sum_{k=0}^\infty \frac{(\mS(\sigma,\lambda_1,\lambda_2,\zeta_1(t),\zeta_2(t)))^k}{k!}\begin{pmatrix}1\\1\end{pmatrix}\\
\end{aligned}
\end{equation}

Next, we are going to derive a concrete expression of $\exp(\mS(\sigma,\lambda_1,\lambda_2,\zeta_1(t),\zeta_2(t)))$. Let the determinant $|\mS(\sigma,\lambda_1,\lambda_2,\zeta_1(t),\zeta_2(t)-\theta\mI|=0$, we can derive that the eigenvalues of $\mS(\sigma,\lambda_1,\lambda_2,\zeta_1(t),\zeta_2(t))$ are

\begin{equation}
\begin{aligned}
&\theta_1= \tfrac{-(\lambda_1\zeta_1(t)+\lambda_2\zeta_2(t))-\sqrt{(\lambda_1\zeta_1(t)-\lambda_2\zeta_2(t))^2+4\sigma^2\zeta_1(t)\zeta_2(t)}}{2}<0\text{ and } \\
% =\tfrac{(\sqrt{(\lambda_1s-\lambda_2)^2+4s\sigma^2}-(s\lambda_1+\lambda_2))\zeta_2(t)}{2},
\end{aligned}
\end{equation}

and the corresponding eigenvectors are

\vv_1=\begin{pmatrix}1\\ \frac{\lambda_1\zeta_1(t)+\theta_1}{\sigma\zeta_1(t)}\end{pmatrix}, \text{ and } \vv_2=\begin{pmatrix}1\\ \frac{\lambda_1\zeta_1(t)+\theta_2}{\sigma\zeta_1(t)}\end{pmatrix}.

Let $\mP=(\vv_1,\vv_2)$, we have

\begin{equation}
\mS(\sigma,\lambda_1,\lambda_2,\zeta_1(t),\zeta_2(t)) = \mP\begin{pmatrix}
\theta_1 & 0\\
0 & \theta_2\\
\end{pmatrix}\mP^{-1},\  \mP^{-1}=\tfrac{\sigma\zeta_1(t)}{\theta_2-\theta_1}\begin{pmatrix}
\frac{\lambda_1\zeta_1(t)+\theta_2}{\sigma\zeta_1(t)} & -1\\
-\frac{\lambda_1\zeta_1(t)+\theta_1}{\sigma\zeta_1(t)} & 1
\end{pmatrix},
\end{equation}

\begin{equation}
\begin{aligned}
e^{\theta_1} & 0\\
0 & e^{\theta_2}\\
\end{pmatrix}\mP^{-1}\\
1 & 1\\
\frac{\lambda_1\zeta_1(t)+\theta_1}{\sigma\zeta_1(t)} & \frac{\lambda_1\zeta_1(t)+\theta_2}{\sigma\zeta_1(t)}\\
\end{pmatrix}
\begin{pmatrix}
e^{\theta_1} & 0\\
0 & e^{\theta_2}\\
\end{pmatrix}
\begin{pmatrix}
\frac{\lambda_1\eta_1(t)+\theta_2}{\sigma\zeta_1(t)} & -1\\
-\frac{\lambda_1\zeta_1(t)+\theta_1}{\sigma\zeta_1(t)} & 1
\end{pmatrix}
\\
(\lambda_1\zeta_1(t)+\theta_2)e^{\theta_1}-(\lambda_1\zeta_1(t)+\theta_1)e^{\theta_2} & \sigma \zeta_1(t)e^{\theta_2}-\sigma\zeta_1(t)e^{\theta_1}\\
(\lambda_1\zeta_1(t)+\theta_2)e^{\theta_2}-(\lambda_1\zeta_1(t)+\theta_1)e^{\theta_1} & \sigma\zeta_2(t)e^{\theta_2}-\sigma\theta_2(t)e^{\theta_1}
\end{pmatrix}
\end{aligned},
\end{equation}

\begin{equation}
\begin{aligned}
& a(t)=\frac{e^{\theta_2}}{\theta_2-\theta_1}\left[(\sigma\zeta_1(t)-\lambda_1\zeta_1(t)-\theta_1)-(\sigma\zeta_1(t)-\lambda_1\zeta_1(t)-\theta_2)e^{\theta_1-\theta_2}\right],\\
& b(t)=\frac{e^{\theta_2}}{\theta_2-\theta_1}\left[(\lambda_1\zeta_1(t)+\theta_2+\sigma\zeta_2(t))+(\lambda_1\zeta_1(t)+\theta_1-\sigma\zeta_2(t))e^{\theta_1-\theta_2}\right].\\
\end{aligned}
\end{equation}

Finally, since $\gE_1^+$, $\gE_1^-$, $\gE_2^+$, $\gE_2^-$ and $\gE_3$ are orthogonal subspace of $\sR^{p\times CN}\oplus\sR^{C\times p}$, we obtain that

\begin{equation}
\begin{aligned}
\mZ(t)=&(\Pi_1^+\mZ_0+\Pi_1^-\mZ_0+\Pi_2^+\mZ_0+\Pi_2^-\mZ_0+\Pi_3\mZ_0)\exp(\mB_\lambda \mC(t))\\
=&\Pi_1^+\mZ_0\begin{pmatrix}a_1^+(t)\mI_{CN} & 0\\0 & b_1^+(t)\mI_{C}\end{pmatrix}+\Pi_1^-\mZ_0\begin{pmatrix}a_1^-(t)\mI_{CN} & 0\\0 & b_1^-(t)\mI_{C}\end{pmatrix}+\\
&\Pi_2^+\mZ_0\begin{pmatrix}a_2^+(t)\mI_{CN} & 0\\0 & b_2^+(t)\mI_{C}\end{pmatrix}+\Pi_2^-\mZ_0\begin{pmatrix}a_2^-(t)\mI_{CN} & 0\\0 & b_2^-(t)\mI_{C}\end{pmatrix}+\\
&\Pi_3\mZ_0\begin{pmatrix}a_3(t)\mI_{CN} & 0\\0 & b_3(t)\mI_{C}\end{pmatrix},
\end{aligned}
\end{equation}

where $a_1^\eps(t)$, $b_1^\eps(t)$, $a_2^\eps(t)$, $b_2^\eps(t)$, $a_3(t)$ and $b_3(t)$ satisfy

\begin{equation}
\begin{aligned}
&\begin{pmatrix}a_1^\eps(t) \\ b_1^\eps(t)\end{pmatrix}=\exp\left(\mS\left(\tfrac{\eps(1+\gamma)}{C\sqrt{N}},\lambda_1, \lambda_2, \zeta_1(t),\zeta_2(t)\right)\right)\begin{pmatrix}1\\ 1\end{pmatrix},\\
&\begin{pmatrix}a_2^\eps(t) \\ b_2^\eps(t)\end{pmatrix}=\exp\left(\mS\left(\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}},\lambda_1, \lambda_2, \zeta_1(t),\zeta_2(t)\right)\right)\begin{pmatrix}1\\ 1\end{pmatrix},\\
& \begin{pmatrix}a_3(t) \\ b_3(t)\end{pmatrix}=\exp\left(\mS\left(0,\lambda_1, \lambda_2, \zeta_1(t),\zeta_2(t)\right)\right)\begin{pmatrix}1\\ 1\end{pmatrix}.\\
\end{aligned}
\end{equation}

Moreover, since $\vb'(t)=\eta_2(t)\frac{1+\gamma -\gamma C}{C}\1_C-\lambda_2\eta_2(t)\vb(t)$ is a first-order linear differential equation, then we have

\begin{equation}
\vb(t)=\phi(t)\left(\vb_0+\tfrac{1+\gamma-\gamma C}{C}\psi(t)\1_C\right),
\end{equation}

where $\phi(t)=\exp(-\lambda_2\int_0^t\eta_2(\tau)\mathrm{d}\tau)$ and $\psi(t)=\int_0^t\zeta_2(\tau)\exp(\lambda_2\int_0^\tau\eta_2(s)\mathrm{d}s)\mathrm{d}\tau$.

§.§ Proof of Corollary <ref>

Corollary <ref>.

Under the conditions and notation of Theorem <ref>, let $s=\frac{\eta_1(0)}{\eta_2(0)}$. If $0<\gamma<\frac{2}{C-2}$ (where $C>2$) or $C=2$, $\lambda_1=\lambda_2=\lambda$, and $\lim_{t\rightarrow\infty}\zeta_1(t)=\infty$, then there exist constants $\pi_h^+,\pi_h^-,\pi_w^+$, $\pi_w^-$, and $\omega$ only depending on $\lambda$, $\gamma$, $s$, $C$, and $N$, such that the gradient flow (as in <ref>) behaves as:

\begin{equation}
\left\|\frac{\mH(t)}{\|\mH(t)\|}-\frac{\pi_h^+\mH_1^++\pi_h^-\mH_1^-}{\|\pi_h^+\mH_1^++\pi_h^-\mH_1^-\|}\right\|+\left\|\frac{\mW(t)}{\|\mW(t)\|}-\frac{\pi_w^+\mW_1^++\pi_w^-\mW_1^-}{\|\pi_w^+\mH_1^++\pi_w^-\mH_1^-\|}\right\|=O(e^{-\omega\zeta_2(t)}),
\end{equation}

where $(\mH_1^+,\mW_1^+)=\Pi_1^+\mZ_0$, $(\mH_1^-,\mW_1^-)=\Pi_1^-\mZ_0$.
Furthermore, we have the following results:

* If $\lambda > \frac{1+\gamma}{C\sqrt{N}}$, then $\lim_{t\rightarrow\infty}\|\mZ(t)\|=0$;

* If $\lambda =\frac{1+\gamma}{C\sqrt{N}}$, then $\lim_{t\rightarrow \infty}\mH(t)=\mH_1^+ + \frac{1-s}{1+s}\mH_1^-, \lim_{t\rightarrow \infty}\mW(t)=\mW_1^+ - \frac{1-s}{1+s}\mW_1^-$;

* If $\lambda < \frac{1+\gamma}{C\sqrt{N}}$, then $\lim_{t\rightarrow\infty}\|\mZ(t)\|=\infty$.

Since $\zeta_1(t)=s\zeta_2(t)$, then the eigenvalues of $\mS(\sigma,\lambda,\zeta_1(t),\zeta_2(t))$ satisfy

\begin{equation}
\begin{aligned}
\theta_1 &=-\frac{\lambda(s+1)+\sqrt{\lambda^2(s-1)^2+4s\sigma^2}}{2}\zeta_2(t)\rightarrow -\infty \text{ as } t\rightarrow\infty, \text{ and } \\
\theta_2 &=\frac{(\sqrt{\lambda^2(s-1)^2+4s\sigma^2}-\lambda (s+1))}{2}\zeta_2(t).
\end{aligned}
\end{equation}

Let $\omega_1(\sigma,\lambda,s)=-\tfrac{\lambda(s+1)+\sqrt{\lambda^2(s-1)^2+4s\sigma^2}}{2}<0$ and $\omega_2(\sigma,\lambda,s)=\tfrac{(\sqrt{\lambda^2(s-1)^2+4s\sigma^2}-\lambda (s+1))}{2}$. For brevity, let $\omega_1$ and $\omega_2$ denote $\omega_1(\sigma,\lambda,s)$ and $\omega_2(\sigma,\lambda,s)$, respectively, and then we can reformulate $a(t)$ and $b(t)$ as

\begin{equation}
\begin{aligned}
&=\tfrac{s\sigma-s\lambda-\omega_1}{\omega_2-\omega_1}e^{\omega_2\zeta_2(t)} + O(e^{(\omega_1)\zeta_2(t)}),\\
&=\tfrac{s\lambda+\omega_2+\sigma}{\omega_2-\omega_1}e^{\omega_2\zeta_2(t)} + O(e^{\omega_1\zeta_2(t)}).
\end{aligned}
\end{equation}

Moreover, according to Theorem <ref>, we have

\begin{equation}
\begin{aligned}
i\in\{1,2\}}}a_i^\eps(t)\mH_i^\eps + a_3(t)\mH_3,\quad
\mW(t)=\sum_{\substack{\eps\in\{\pm\}\\
i\in\{1,2\}}}b_i^\eps(t)\mW_i^\eps + b_3(t)\mW_3,\\
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
&\begin{pmatrix}a_1^\eps(t) \\ b_1^\eps(t)\end{pmatrix}=\exp\left(\mS\left(\tfrac{\eps(1+\gamma)}{C\sqrt{N}},\lambda, \zeta_1(t),\zeta_2(t)\right)\right)\begin{pmatrix}1\\ 1\end{pmatrix},\\
&\begin{pmatrix}a_2^\eps(t) \\ b_2^\eps(t)\end{pmatrix}=\exp\left(\mS\left(\tfrac{\eps(1+\gamma-\gamma C)}{C\sqrt{N}},\lambda, \zeta_1(t),\zeta_2(t)\right)\right)\begin{pmatrix}1\\ 1\end{pmatrix},\\
& \begin{pmatrix}a_3(t) \\ b_3(t)\end{pmatrix}=\exp\left(\mS\left(0,\lambda, \zeta_1(t),\zeta_2(t)\right)\right)\begin{pmatrix}1\\ 1\end{pmatrix}.\\
\end{aligned}
\end{equation}

Since $0<\gamma<\frac{2}{C-2}$ (where $C>2$) or $C=2$, thus we have $\frac{1+\gamma }{C\sqrt{N}} > \frac{|1+\gamma -\gamma C|}{C\sqrt{N}}$, and then $\omega_2\left(\frac{1+\gamma}{C\sqrt{N}},\lambda,s\right)>\omega_2\left(\frac{|1+\gamma-\gamma C|}{C\sqrt{N}},\lambda,s\right)$. When $t\rightarrow \infty$, the dominant terms in $\mH(t)$ and $\mW(t)$ are the ones whose coefficient contains $\exp\left(\omega_2\left(\frac{1+\gamma}{C\sqrt{N}},\lambda,s\right)\zeta_2(t)\right)$, i.e., $a_1^+(t)$, $a_1^-(t)$, $b_1^+(t)$ and $b_1^-(t)$. Let $(\mH_1^\eps,\mW_1^\eps)=\Pi_1^\eps\mZ_0$, $(\mH_2^\eps,\mW_2^\eps)=\Pi_2^\eps\mZ_0$ and $(\mH_3,\mW_3)=\Pi_3\mZ_0$ for $\eps\in\{\pm1\}$, thus we have

\begin{equation}
\begin{aligned}
&e^{-\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)\zeta_2(t)}\mH(t)= \pi_h^+(\lambda, \gamma, s,C,N)\mH_1^+ + \pi_h^-(\lambda, \gamma, s,C,N)\mH_1^- + \bm{\Delta}_1,\\
&e^{-\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)\zeta_2(t)}\mW(t)= \pi_w^+(\lambda, \gamma, s,C,N)\mW_1^+ + \pi_w^-(\lambda, \gamma, s,C,N)\mW_1^- + \bm{\Delta}_2,
\end{aligned}
\end{equation}

where $\bm{\Delta}_1$ and $\bm{\Delta}_2$ decrease to zero as least as $O\left(e^{\frac{\left(\omega_2\left(\frac{|1+\gamma -\gamma C|}{C\sqrt{N}}, \lambda,s\right)-\omega_2\left(\frac{1+\gamma}{C\sqrt{N}}, \lambda,s\right)\right)\zeta_2(t)}{C\sqrt{N}}}\right)$, and

\begin{equation}
\begin{aligned}
&\pi_h^+(\lambda, \gamma, s,C,N)=\frac{s\frac{1+\gamma}{C\sqrt{N}}-s\lambda-\omega_1\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)}{\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)-\omega_1\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)},\\
&\pi_h^-(\lambda, \gamma, s,C,N) = \frac{-s\frac{1+\gamma}{C\sqrt{N}}-s\lambda-\omega_1\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)}{\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)-\omega_1\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)},\\
&\pi_w^+(\lambda, \gamma, s,C,N)=\frac{s\lambda+\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)+\tfrac{1+\gamma}{C\sqrt{N}}}{\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)-\omega_1\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)},\\
&\pi_w^-(\lambda, \gamma, s,C,N) = \frac{s\lambda+\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)-\tfrac{1+\gamma}{C\sqrt{N}}}{\omega_2\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)-\omega_1\left(\tfrac{1+\gamma}{C\sqrt{N}},\lambda,s\right)}.\\
\end{aligned}
\end{equation}

Therefore, we have

\begin{equation}
\begin{aligned}
&\lim_{t\rightarrow \infty}\frac{\mH(t)}{\|\mH(t)\|} = \frac{\pi_h^+(\lambda, \gamma, s,C,N)\mH_1^+ + \pi_h^-(\lambda, \gamma, s,C,N)\mH_1^-}{\|\pi_h^+(\lambda, \gamma, s,C,N)\mH_1^+ + \pi_h^-(\lambda, \gamma, s,C,N)\mH_1^-\|},\\
&\lim_{t\rightarrow \infty}\frac{\mW(t)}{\|\mW(t)\|} = \frac{\pi_w^+(\lambda, \gamma, s,C,N)\mW_1^+ + \pi_w^-(\lambda, \gamma, s,C,N)\mW_1^-}{\|\pi_w^+(\lambda, \gamma, s,C,N)\mW_1^+ + \pi_w^-(\lambda, \gamma, s,C,N)\mW_1^-\|},
\end{aligned}
\end{equation}

and the rate of convergence is $O\left(e^{\frac{\left(\omega_2\left(\frac{|1+\gamma -\gamma C|}{C\sqrt{N}}, \lambda,s\right)-\omega_2\left(\frac{1+\gamma}{C\sqrt{N}}, \lambda,s\right)\right)\zeta_2(t)}{C\sqrt{N}}}\right)$.

Moreover, we have the following conclusions:

* If $\lambda=\frac{1+\gamma}{C\sqrt{N}}$, we have $\omega_2=0$, $\omega_1=-\lambda(s+1)$, $\pi_h^+=1$, $\pi_h^-=\frac{1-s}{1+s}$, $\pi_w^+=1$ and $\pi_w^-=-\frac{1-s}{1+s}$. Since $\lim_{t\rightarrow \infty}\zeta_2(t)=\infty$, we have

\begin{equation}
\lim_{t\rightarrow \infty}\mH(t)=\mH_1^+ + \frac{1-s}{1+s}\mH_1^-, \lim_{t\rightarrow \infty}\mW(t)=\mW_1^+ - \frac{1-s}{1+s}\mW_1^-.
\end{equation}

* If $\lambda > \frac{1+\gamma}{C\sqrt{N}}$, we have $\omega_2 > 0$, and then $\lim_{t\rightarrow \infty}\|\mZ(t)\|=0$ since $\lim_{t\rightarrow \infty}\zeta_2(t)=\infty$.

* If $\lambda <\frac{1+\gamma}{C\sqrt{N}}$, we have $\omega_2 < 0$, and then $\lim_{t\rightarrow \infty}\|\mZ(t)\|=\infty$ since $\lim_{t\rightarrow \infty}\zeta_2(t)=\infty$.

So far the proof has been completed.

§.§ Proof of Theorem <ref>

For $\vh(t), \vw\in\sR^p$, $\eta(t)>0$, let $\hat{\vv}=\frac{\hat{\vv}}{\|\vv\|_2}$ denote the $\ell_2$-normalized vector of $\vv$, considering the discrete dynamical system $\vh(t+1)=\vh(t)+\frac{\eta(t)}{\|\vh(t)\|_2}\left(\mI_p-\hat{\vh}(t)\hat{\vh}^\top(t)\right)\vw$, if $\hat{\vw}^\top\hat{\vh}(0)>-1$, the learning rate $\eta(t)$ satisfies that $\lim_{t\rightarrow\infty}\frac{\eta(t+1)}{\eta(t)}=1$, $\frac{\eta(t)}{\|\vh(t)\|_2}$ is non-increasing with $\frac{\eta(0)}{\|\vh(0)\|_2}<\frac{1}{\|\vw\|_2}$, and there exists a constant $\varepsilon>0$, s.t., $\eta(t)>\varepsilon$, s.t., $\eta(t)>\varepsilon$, then we have

\begin{equation}
\lim_{t\rightarrow \infty}\left\|\frac{\vh(t)}{\|\vh(t)\|_2}-\frac{\vw}{\|\vw\|_2}\right\| = 0.
\end{equation}

For brevity, let $\alpha_t=\frac{\eta(t)\|\vw\|_2}{\|\vh(t)\|_2^2}$, $\xi_t=\frac{\eta(t+1)}{\eta(t)}$, and $\beta_t=\hat{\vw}^\top\hat{\vh}(t)$ denote the cosine similarity between $\hat{\vw}$ and $\hat{\vh}(t)$, then we can easily derive that $\alpha_t>0$ and $\beta_{t}>-1$ for all $t\ge 0$.

We will show that $\alpha_t$ is monotonically decreasing and $\beta_t$ is monotonically increasing. Note that $\vh(t)$ is orthogonal with $\left(\mI_p-\hat{\vh}(t)\hat{\vh}^\top(t)\right)\vw$, thus we have

\begin{equation}
\notag
\|\vh(t+1)\|_2^2=\|\vh(t)\|_2^2 + \frac{\eta^2(t)}{\|\vh(t)\|_2^2} \left\|\left(\mI_p-\hat{\vh}(t)\hat{\vh}^\top(t)\right)\vw\right\|_2^2\ge \|\vh(t)\|_2^2,
\end{equation}

which indicates $\|\vh(t)\|_2$ is monotonically increasing as a function of $t$, and then $\alpha_t$ is monotonically decreasing since $\frac{\eta(t)\|\vw\|_2}{\|\vh(t)\|_2}$ is non-increasing.

Moreover, we can rearrange the discrete dynamics and formulate $\vh(t+1)$ as a positive combination of $\vh(t)$ and $\vw$:

\begin{equation}
\label{discrete-dynamics-of-H}
\vh(t+1)=\left(1-\frac{\eta(t)\vw^\top\hat{\vh}(t)}{\|\vh(t)\|_2}\right)\vh(t)+\frac{\eta(t)}{\|\vh(t)\|_2}\vw,
\end{equation}

so $\hat{\vw}^\top\hat{\vh}(t+1)\ge\hat{\vw}^\top\hat{\vh}(t)$, i.e., $\beta_{t}$ is monotonically increasing, which is based on the facts that $1-\frac{\eta(t)\vw^\top\hat{\vh}(t)}{\|\vh(t)\|_2}\ge 1-\frac{\eta(t)\|\vw\|_2}{\|\vh(t)\|_2}\ge 1-\frac{\eta(0)\|\vw\|_2}{\|\vh(0)\|_2}>0$, $\frac{\eta(t)}{\|\vh(t)\|_2}>0$, and $\frac{\vy^\top(\vx+k\vy )}{\|\vx+k\vy\|_2}\ge \frac{\vy^\top\vx}{\|\vx\|_2}$ holds for all $k>0$ and $\vx,\vy\neq 0$.

We can formulate the discrete iterations of $\alpha_t$ and $\beta_t$ from the <ref> as follows:

\begin{equation}
\label{discrete_dynamical-system}
\begin{aligned}
\beta_{t+1}&=\frac{{\vw}^\top{\vh}(t+1)}{\|\vw\|_2\|\vh(t+1)\|_2}\\
&=\frac{\vw^\top\left(\vh(t)+\frac{\eta(t)}{\|\vh(t)\|_2}\left(\mI_p-\hat{\vh}(t)\hat{\vh}^\top(t)\right)\vw\right)}{\|\vw\|_2\sqrt{\|\vh(t)\|_2^2 + \frac{\eta^2(t)}{\|\vh(t)\|_2^2} \left\|\left(\mI_p-\hat{\vh}(t)\hat{\vh}^\top(t)\right)\vw\right\|_2^2}}\\
\alpha_{t+1}&= \frac{\eta(t+1)\|\vw\|_2}{\|\vh(t+1)\|_2^2}\\
&=\frac{\xi_t\eta(t)\|\vw\|_2}{\|\vh(t)\|_2^2 + \frac{\eta^2(t)}{\|\vh(t)\|_2^2} \left\|\left(\mI_p-\hat{\vh}(t)\hat{\vh}^\top(t)\right)\vw\right\|_2^2}\\
\end{aligned}
\end{equation}

with $\beta_0=\hat{\vw}^\top\hat{\vh}(0)>-1$, $\alpha_0=\frac{\eta(0)\|\vw\|_2}{\|\vh(0)\|_2^2}>0$, $\xi_t\le 1$ and $\lim_{t\rightarrow\infty}\xi_t=1$.

To prove $\lim_{t\rightarrow \infty}\left\|\frac{\vh(t)}{\|\vh(t)\|_2}-\frac{\vw}{\|\vw\|_2}\right\| = 0$, we just need prove $\lim_{t\rightarrow \infty} \beta_t=1$. Note that $\alpha_t$ is monotonic decreasing and lower bounded by $0$, then the sequence $(\alpha_t)$ is convergent. Similarly, the sequence $(\beta_t)$ is convergent. Let $a = \lim_{t\rightarrow} \alpha_t$ and $b=\lim_{t\rightarrow \beta_t}$, we obtain

\begin{equation}
\lim_{t\rightarrow\infty}\alpha_{t+1}=\lim_{t\rightarrow\infty} \frac{\xi_t\alpha_t}{1+\alpha_t^2(1-\beta_t^2)}\quad \Rightarrow\quad a=\frac{a}{1+\lim_{t\rightarrow\infty}\alpha_t^2(1-\beta_t^2)},
\end{equation}

thus $a=0$ or $\lim_{t\rightarrow \infty}\alpha_t^2(1-\beta_t^2)=0$, i.e., $a=0$ or $b=1$. Therefore, the limits of $\alpha_t$ and $\beta_t$ exist if and only if $\lim_{t\rightarrow \infty} \alpha_t=0$ or $\lim_{t\rightarrow \infty} \beta_t=1$. In the following, we will prove that the limit of $\beta_t$ must be equal to 1.

Firstly, we prove a simpler result when $\beta_0>0$:

For the discrete dynamics in <ref>, if $\beta_0\ge0$, then $\lim_{t\rightarrow \infty} \beta_t=1$.

As aforementioned, due to the existence of the limit of $\alpha_t$, we have  $\lim_{t\rightarrow\infty}\alpha_t=0$ or $\lim_{t\rightarrow\infty}\beta_t=1$. Thus, we just need to prove that $\lim_{t\rightarrow\infty}\beta_t=1$ as $\lim_{t\rightarrow\infty}\alpha_t=0$.

When $\lim_{t\rightarrow \infty}\alpha_t=0$, then there exists $\tau$, such that $\forall t>\tau$, $\alpha_t\le 1$.

According to the iterations in <ref>, we can derive that

\begin{equation}
\notag
\begin{aligned}
\frac{1-\beta_{t+1}^2}{\alpha_{t+1}}
&=\frac{1-\beta_t^2}{\xi_t\alpha_t}\cdot (1+\alpha_t^2\beta_t^2-2\alpha_t\beta_t)\\
&=\frac{1-\beta_t^2}{\alpha_t}\cdot \frac{(1-\alpha_t\beta_t)^2}{\xi_t},
\end{aligned}
\end{equation}

then $\forall t>\tau$,

\begin{equation}
\begin{aligned}
&=  \alpha_{t+1}\cdot \frac{\eta(0)(1-\beta_0^2)}{\alpha_0\eta(t+1)}\cdot \prod_{i=0}^\tau (1-\alpha_i\beta_i)^2 \cdot\prod_{i=\tau+1}^t (1-\alpha_i\beta_i)^2\\
&\le\alpha_{t+1}\cdot \frac{\eta(0)(1-\beta_0^2)}{\varepsilon\alpha_0}\cdot \prod_{i=0}^\tau (1-\alpha_i\beta_i)^2,
\end{aligned}
\end{equation}

where the inequality is based on the fact that $1-\alpha_i\beta_i\in(0,1]$ when $0\le\beta_0\le\beta_i\le1$, $\eta(t+1)\ge \varepsilon$, and $\alpha_i\le 1$ for $i>\tau$. Since $\frac{\eta(0)(1-\beta_0^2)}{\varepsilon\alpha_0}\cdot \prod_{i=0}^\tau (1-\alpha_i\beta_i)^2$ is a constant, we obtain

\begin{equation}
\notag
\lim_{t\rightarrow \infty} 1-\beta_{t+1}^2\le \lim_{t\rightarrow \infty} \alpha_{t+1}\cdot \frac{\eta(0)(1-\beta_0^2)}{\varepsilon\alpha_0}\cdot \prod_{i=0}^\tau (1-\alpha_i\beta_i)^2=0,
\end{equation}

as $\lim_{t\rightarrow\infty}\alpha_{t+1}=0$. This reveals $\lim_{t\rightarrow \infty}\beta_t^2=1$. Furthermore, since $\beta_t\ge 0$, we then have $\lim_{t\rightarrow \infty}\beta_t=1$.

Next, we are going to prove $\lim_{t\rightarrow\infty}\beta_t=1$ when $-1<\beta_0< 0$. According to Lemma <ref>, we just need prove that $\exists \tau>0$, s.t., $\beta_\tau\ge 0$.

For the sake of contradiction, suppose that $\beta_t< 0$ for all $t>0$, we then have $\lim_{t\rightarrow\infty}\alpha_t=0$. As a consequence, we obtain

\begin{equation}
\alpha_t+\beta_t-\alpha_t\beta_t^2< 0,\ \forall t\ge 0,
\end{equation}

and we know that $\exists t'>0$, such that

\begin{equation}
\label{contradiction}
\alpha_t<\frac{\varepsilon}{\eta(0)}\alpha_0(1-\beta_0^2),\ \forall t\ge t'.
\end{equation}

According to the iterations in <ref>, we can derive that

\begin{equation}
\notag
\begin{aligned}
\alpha_{t+1}(1-\beta_{t+1}^2)&=\frac{\xi_t\alpha_t}{1+\alpha_t^2-\alpha_t^2\beta_t^2}\left(1-\frac{(\alpha_t+\beta_t-\alpha_t\beta_t^2)^2}{1+\alpha_t^2-\alpha_t^2\beta_t^2}\right)\\
& =\xi_t \alpha_t(1-\beta_t^2)\left(1-\frac{\alpha_t(\alpha_t+\beta_t-\alpha_t\beta_t^2)}{1+\alpha_t^2-\alpha_t^2\beta_t^2}\right)^2,
\end{aligned}
\end{equation}

Since $1-\frac{\alpha_t(\alpha_t+\beta_t-\alpha_t\beta_t^2)}{1+\alpha_t^2-\alpha_t^2\beta_t^2}\ge1$, then for $t\ge t'$,

\begin{equation}
\alpha_t\ge\alpha_{t}(1-\beta_{t}^2)\ge \xi_{t-1}\alpha_{t-1}(1-\beta_{t-1}^2)\ge \ldots\ge \alpha_0(1-\beta_0^2)\prod_{i=0}^{t-1}\xi_i \ge\frac{\varepsilon}{\eta(0)}\alpha_0(1-\beta_0^2),
\end{equation}

which contradicts the fact in <ref>. Thus, $\exists \tau>0$, s.t. $\beta_\tau\ge 0$. Consider the dynamical system with an initial time $\tau$, we have $\lim_{t\rightarrow\infty} \beta_{t}=1$ according to Lemma <ref>.

To sum up, we have proven that $\lim_{t\rightarrow\infty}\beta_t=1$ when $\beta_0>-1$ and $\alpha_0>0$.

Theorem <ref> (Convergence in the Spherical Constrained Case)

Considering the discrete dynamics in <ref>, if $\forall i\in[N],c\in[C]$, $\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)>-1$, the learning rate $\eta(t)$ satisfies that $\frac{\eta(t)}{\|\vh_{i,c}(t)\|_2}$ is non-increasing,  $\frac{\eta(0)(1+\gamma)}{CN\|\vh_{i,c}(0)\|_2}\le \frac{1}{\|\vw_c\|_2}$,  $\lim_{t\rightarrow\infty}\frac{\eta(t+1)}{\eta(t)}=1$, and there exists a constant $\varepsilon>0$, s.t., $\eta(t)>\varepsilon$, then we have

\begin{equation}
\lim_{t\rightarrow\infty}\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|=0,
\end{equation}

and further if $\lim_{t\rightarrow\infty}\|\mH(t)\|<\infty$, then there exists a constant $\mu>0$, such that the error above shows exponential convergence:

\begin{equation}
\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|\le O(e^{-\mu t}).
\end{equation}

Moreover, if $\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)=-1$, then $\vh_{i,c}(t)=\vh_{i,c}(0)$.

Since $\mH(t+1)=\mH(t)+\frac{(1+\gamma)\eta(t)}{CN}\left(\frac{\partial \hat{\mH}}{\partial \mH}\big|_{\mH=\mH(t)}\right)^\top\mW(\mI_C\otimes \1_N^\top)$, then for $i\in[N], c\in[C]$,

\begin{equation}
\vh_{i,c}(t+1)=\vh_{i,c}(t)+\frac{(1+\gamma)\eta(t)}{CN\|\vh_{i,c}(t)\|_2}\left(\mI_p-\hat{\vh}_{i,c}(t)\hat{\vh}_{i,c}^\top(t)\right)\vw_{c}.
\end{equation}

According to Lemma <ref>, when $\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)>-1$ and $\frac{\eta(0)(1+\gamma)}{CN}<\frac{\|\vh_{i,c}(0)\|_2^2}{\|\vw_c\|_2}$, we have $\lim_{t\rightarrow \infty}\|\hat{\vh}_{i,c}-\hat{\vw}_c\|=0$, then $\lim_{t\rightarrow\infty}\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|=0$.

If further $\lim_{t\rightarrow\infty}\|\mH(t)\|_2<\infty$, let $L=\sup_{i,c,t}\|\vh_{i,c}(t)\|$. According to the proof of Lemma <ref>, $\forall i,c$, we have $\lim_{t\rightarrow\infty}\hat{\vw}_c^\top\hat{\vh}_{i,c}(t)=1$, then for a given constant $\delta>0$, $\exists \tau>0$, s.t., $\forall t> \tau$, $\hat{\vw}_c^\top\hat{\vh}_{i,c}(\tau)\ge \delta$. Consider $t>\tau$, we have

\begin{equation}
\notag
\begin{aligned}
=&\tfrac{\|\vh_{i,c}(t)\|_2^2}{\|\vh_{i,c}(t+1)\|_2^2}\left(1-\left(\hat{\vw}_c^\top\hat{\vh}_{i,c}(t)\right)^2\right)\left(1-\tfrac{(1+\gamma)\eta(t)\|\vw_c\|_2}{CN\|\vh_{i,c}(t)\|_2^2}\cdot \hat{\vw}_c^\top\hat{\vh}_{i,c}(t)\right)^2\\
\le &\left(1-\left(\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)\right)^2\right)\prod_{j=0}^t \left(1-\tfrac{(1+\gamma)\eta(j)\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}\cdot \hat{\vw}_c^\top\hat{\vh}_{i,c}(j)\right)^2\\
\le &\left(1-\left(\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)\right)^2\right)\prod_{j=0}^\tau \left(1-\tfrac{(1+\gamma)\eta(j)\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}\cdot \hat{\vw}_c^\top\hat{\vh}_{i,c}(j)\right)^2\prod_{j=\tau+1}^t \left(1-\tfrac{(1+\gamma)\varepsilon\delta\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}\right)^2\\
\le & \left(1-\left(\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)\right)^2\right)\prod_{j=0}^\tau \left(1-\tfrac{(1+\gamma)\eta(j)\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}\cdot \hat{\vw}_c^\top\hat{\vh}_{i,c}(j)\right)^2 \left(1-\tfrac{(1+\gamma)\varepsilon\delta\|\vw_c\|_2}{CNL^2}\right)^{2(t-\tau)}.
\end{aligned}
\end{equation}

where the first, the second, and the third inequalities are based on the facts that $\frac{\|\vh_{i,c}(t)\|_2^2}{\|\vh_{i,c}(t+1)\|_2^2}\le 1$, $1-\frac{(1+\gamma)\eta(j+1)\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}\cdot \hat{\vw}_c^\top\hat{\vh}_{i,c}(j)\le 1-\frac{(1+\gamma)\varepsilon\delta\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}$ for $t>\tau$, and $\|\vh_{i,c}(j)\|_2\le L$, respectively.

Let $c_1=\max_{i,c}(1-(\hat{\vw}_c^\top\hat{\vh}_{i,c}(0))^2)\prod_{j=0}^\tau \left(1-\frac{(1+\gamma)\eta(j)\|\vw_c\|_2}{CN\|\vh_{i,c}(j)\|_2^2}\cdot \hat{\vw}_c^\top\hat{\vh}_{i,c}(j)\right)^2 \left(1-\frac{(1+\gamma)\varepsilon\delta\|\vw_c\|_2}{CNL^2}\right)^{-2\tau}$, and $\mu=\min_c -2\log \left(1-\frac{(1+\gamma)\varepsilon\delta\|\vw_c\|_2}{CNL^2}\right)$, then $1-\hat{\vw}_c^\top\hat{\vh}_{i,c}(t+1)\le \frac{c_1e^{-\mu t}}{1+\delta}$.

Therefore, we have

\begin{equation}
\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|_2^2=2\sum_{i,c}\big(1-\hat{\vw}_c^\top\hat{\vh}_{i,c}(t+1)\big)\le \frac{2c_1CNe^{-\mu t}}{(1+\delta)},
\end{equation}

i.e., $\left\|\hat{\mH}(t)-\hat{\mW}(\mI_C\otimes \1_N^\top)\right\|=O(e^{-\mu t})$.

Moreover, if $\hat{\vw}_c^\top\hat{\vh}_{i,c}(0)=-1$, we have

\begin{equation}
\vh_{i,c}(t+1)=\vh_{i,c}(t)+\frac{(1+\gamma)\eta(t)}{CN\|\vh_{i,c}(t)\|_2}\left(\mI_p-\hat{\vh}_{i,c}(t)\hat{\vh}_{i,c}^\top(t)\right)\vw_{c}=\vh_{i,c}(t),
\end{equation}

thus $\vh_{i,c}(t)=\vh_{i,c}(0)$.

§.§ Proof of Theorem <ref>

Let $\vz(t)$ denote the row-first vectorization of $\begin{pmatrix}\boldsymbol{H}(t) & 0\\0 & \boldsymbol{W}(t)\end{pmatrix}$, $\boldsymbol B=\begin{pmatrix}0 & \boldsymbol M^\top\\\boldsymbol M & 0\end{pmatrix}$, and $\boldsymbol A=\begin{pmatrix}0 & \nabla_{\boldsymbol \Theta}\boldsymbol H ^\top\nabla_{\boldsymbol \Theta}\boldsymbol H\\\mathbf I_C & 0\end{pmatrix}$. Considering the eigendecomposition $\boldsymbol A=\boldsymbol U_{\boldsymbol A}\boldsymbol \Lambda_{\boldsymbol A} \boldsymbol U_{\boldsymbol A}^{-1}$ and $\boldsymbol B=\boldsymbol U_{\boldsymbol B}\boldsymbol \Lambda_{\boldsymbol B} \boldsymbol U_{\boldsymbol B}^{-1}$, we have

\begin{equation}
\mC\boldsymbol z(t)=\exp[(\boldsymbol \Lambda_{\boldsymbol A}\otimes \boldsymbol \Lambda_{\boldsymbol B})t]\mC\boldsymbol z(0),
\end{equation}

where $\mC=\boldsymbol U_{\boldsymbol B}^{-1}\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I$.

Since $\boldsymbol z(t)$ denote the row-first vectorization of $\begin{pmatrix}\boldsymbol{H}(t) & 0\\0 & \boldsymbol{W}(t)\end{pmatrix}$, $\boldsymbol A=\begin{pmatrix}0 & \nabla_{\boldsymbol \Theta}\boldsymbol H ^\top\nabla_{\boldsymbol \Theta}\boldsymbol H\\\mathbf I & 0\end{pmatrix}$, and $\boldsymbol B=\begin{pmatrix}0 & \boldsymbol M^\top\\\boldsymbol M & 0\end{pmatrix}$, we have

\boldsymbol z'(t)=(\boldsymbol A\otimes \boldsymbol{B}^\top) \boldsymbol z(t)=(\boldsymbol A\otimes \boldsymbol{B}) \boldsymbol z(t),

and then $\boldsymbol z(t)=e^{(\boldsymbol A\otimes \boldsymbol{B})t}\boldsymbol z(0)$. Considering the eigendecomposition $\boldsymbol A=\boldsymbol U_{\boldsymbol A}\boldsymbol \Lambda_{\boldsymbol A} \boldsymbol U_{\boldsymbol A}^{-1}$ and $\boldsymbol B=\boldsymbol U_{\boldsymbol B}\boldsymbol \Lambda_{\boldsymbol B} \boldsymbol U_{\boldsymbol B}^{-1}$, where $\boldsymbol \Lambda_{\boldsymbol A}=\text{diag}(\lambda_1^{\boldsymbol A},...,\lambda_n^{\boldsymbol A})$ and $\boldsymbol \Lambda_{\boldsymbol B}=\text{diag}(\lambda_1^{\boldsymbol B},...,\lambda_n^{\boldsymbol B})$, thus

\boldsymbol A\otimes \boldsymbol B= (\boldsymbol U_{\boldsymbol A}\boldsymbol \Lambda_{\boldsymbol A} \boldsymbol U_{\boldsymbol A}^{-1})\otimes \boldsymbol B=(\boldsymbol U_{\boldsymbol A}\otimes \mathbf I)(\boldsymbol \Lambda_{\boldsymbol A} \otimes \boldsymbol B)(\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I),

where $\otimes$ denotes Kronecker product. Moreover, we have

\begin{aligned}
\exp[(\boldsymbol A\otimes \boldsymbol{B})t]&=(\boldsymbol U_{\boldsymbol A}\otimes \mathbf I)\exp[(\boldsymbol \Lambda_{\boldsymbol A}\otimes \boldsymbol{B})t](\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I)\\
&=(\boldsymbol U_{\boldsymbol A}\boldsymbol U_{\boldsymbol B}\otimes \mathbf I)\exp[(\boldsymbol \Lambda_{\boldsymbol A}\otimes \boldsymbol \Lambda_{\boldsymbol B})t](\boldsymbol U_{\boldsymbol B}^{-1}\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I).\\
\end{aligned}

We know that $\exp[(\boldsymbol \Lambda_{\boldsymbol A}\otimes \boldsymbol \Lambda_{\boldsymbol B})t]$ is a diagonal matrix. Therefore, we have $(\boldsymbol U_{\boldsymbol B}^{-1}\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I)\boldsymbol z(t)=\exp[(\boldsymbol \Lambda_{\boldsymbol A}\otimes \boldsymbol \Lambda_{\boldsymbol B})t](\boldsymbol U_{\boldsymbol B}^{-1}\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I)\boldsymbol z(0)$, which leads to a concise closed-form dynamics on $(\boldsymbol U_{\boldsymbol B}^{-1}\boldsymbol U_{\boldsymbol A}^{-1}\otimes \mathbf I)\boldsymbol z(t)$.

§.§ Proof of Theorem <ref>

Theorem <ref>
Consider the continual gradient flow (<Ref>) in which the prototypes $\mW$ is fixed, we have the closed-form dynamics:

\begin{equation}
\mH(t)=e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\mH(0)+\frac{1-e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}}{\lambda} \mW\mM,
\end{equation}

which further indicates that $\left\|\mH(t)-\frac{1}{\lambda}\mW\mM\right\|=O\left(e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\right)$.

For the first order non-homogeneous linear difference equation in <ref>, the solution is

\begin{equation}
\begin{aligned}
\mH(t)=&e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\left(\mH(0)+\int_0^t\eta(s)e^{\lambda\int_0^s\eta(\tau)\mathrm{d}\tau}\mathrm{d}s \mW\mM\right)\\
=&e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\mH(0) + e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\int_0^t\frac{1}{\lambda}\mathrm{d}e^{\lambda\int_0^s\eta(\tau)\mathrm{d}\tau} \mW\mM\\
=&e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\mH(0) + e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\frac{e^{\lambda\int_0^s\eta(\tau)\mathrm{d}\tau}}{\lambda}\bigg|_0^t \mW\mM\\
=&e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}\mH(0)+\frac{1-e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}}{\lambda} \mW\mM,
\end{aligned}
\end{equation}

and then $\|\mH(t)-\frac{1}{\lambda}\mW\mM\|=\|e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau}(\mH(0)-\frac{\mW\mM}{\lambda})\|=O(e^{-\lambda\int_0^t\eta(\tau)\mathrm{d}\tau})$.

§.§ The Projections onto $\gE_1^+$, $\gE_1^-$, $\gE_2^+$, $\gE_2^-$ and $\gE_3$

Let $\gS$ denote the subspace $\gS=\{\mW: \mW\1_n=0, \mW\in\sR^{m\times n}\}$, then the projection of a point $\mA\in\sR^{m\times n}$ onto $\gS$ can be denoted as $\Pi_{\gS}\mA=\mA\left(\mI_n-\frac{1}{n}\1_n\1_n^\top\right)$.

Let $\mW=(\vw_1,\vw_2,...,\vw_n)\in \gS$ and $\mA=(\va_1,\va_2,...,\va_n)\in\sR^{m\times n}$, we have

\begin{equation}
\notag
\begin{aligned}
\|\mW-\mA\|_F^2&=\sum_{i=1}^n \|\vw_i-\va_i\|_2^2\ge \frac{1}{n}\left\|\sum_{i=1}^n\vw_i-\sum_{i=1}^n\va_i\right\|_2^2\\
\end{aligned}
\end{equation}

where we used the Cauchy-Schwarz inequality, and the equality holds if and only if $\vw_i-\va_i=\vw_n-\va_n$, $\forall i\in [n]$, and $\sum_{i=1}^n \vw_i=0$, i.e., $\mW=\mA\left(\mI_n-\frac{1}{n}\1_n\1_n^\top\right)$. Therefore, the projection of $\mA$ onto $\gS$ is $\Pi_{\gS}\mA=\argmin_{\mW\in \gS}\|\mW - \mA\|_F^2=\mA\left(\mI_n-\frac{1}{n}\1_n\1_n^\top\right).

Let $\gS$ denote the subspace $\{\mW: \mW(\mI_c\otimes \1_n)=0, \mW\in\sR^{m\times cn}\}$, then the projection of a point $\mA\in\sR^{m\times cn}$ onto $\gS$ can be denoted as $\Pi_{\gS}\mA=\mA\left(\mI_{cn}-\frac{1}{n}\mI_c\otimes \1_n\1_n^\top\right)$.

The proof is similar to Lemma <ref>. We can simply let $\mW=(\mW_1,...,\mW_n)\in\gS$ and $\mA=(\mA_1,...,\mA_n)\in\sR^{m\times cn}$, where $\mW_i,\mA_i\in\sR^{m\times c}$, then we have

\begin{equation}
\notag
\begin{aligned}
\|\mW-\mA\|_F^2&=\sum_{i=1}^n \|\mW_i-\mA_i\|_F^2\ge \frac{1}{n}\left\|\sum_{i=1}^n\mW_i-\sum_{i=1}^n\mA_i\right\|_2^2=\frac{1}{n}\|\mA(\mI_c\otimes \1_n)\|_2^2
\end{aligned}
\end{equation}

where the equality holds if and only if $\mW_i-\mA_i=\mW_n-\mA_n$, $\forall i\in [n]$, and $\sum_{i=1}^n \mW_i=0$, i.e., $\mW=\mA(\mI_{cn}-\tfrac{1}{n}\mI_c\otimes \1_n\1_n^\top)$. Therefore, the projection of $\mA$ onto $\gS$ is $\Pi_\gS\mA=\mA\left(\mI_{cn}-\frac{1}{n}\mI_c\otimes \1_n\1_n^\top\right)$.

For $\mH\in\sR^{p\times CN}$, $\mW\in\sR^{p\times C}$, the projection of $(\mH,\mW)$ onto $\gE_1^\eps$ is

\begin{equation}
\Pi_1^\eps (\mH,\mW)=(\tfrac{\eps}{\sqrt{N}}(\mP\otimes \1_N^\top), \mP),
\end{equation}

where $\mP=\tfrac{1}{2}(\tfrac{\eps}{\sqrt{N}}\mH(\mI_C\otimes \1_N)+\mW)(\mI_C-\tfrac{1}{C}\1_C\1_C^\top)$.

Let $\gS=\{\mZ:\mZ\1_C=0, \mZ\in\sR^{p\times C}\}$ and $\mH=\{\mH_1,...,\mH_N\}$ (where $\mH_i\in\sR^{p\times C}$), the minimizer of $\mZ\in\gS$ is

\begin{equation}
\begin{aligned}
&\argmin_{\mZ\in \gS}\|\tfrac{\eps}{\sqrt{N}}(\mZ\otimes \1_N^\top)-\mH\|_F^2+\|\mZ-\mW\|_F^2\\
=&\argmin_{\mZ\in \gS}\sum_{i=1}^N\|\tfrac{1}{\sqrt{N}}\mZ-\eps\mH_i\|_F^2+\|\mZ-\mW\|_F^2\\
=&\argmin_{\mZ\in \gS}\|\mZ\|_F^2-\tfrac{2\eps}{\sqrt{N}}\sum_{i=1}^N\langle\mZ,\mH_i\rangle+\|\mH\|_F^2+\|\mZ-\mW\|_F^2\\
=&\argmin_{\mZ\in \gS}\|\mZ-\tfrac{\eps}{\sqrt{N}}\sum_{i=1}^N\mH_i\|_F^2-\|\frac{1}{\sqrt{N}}\sum_{i=1}^N\mH_i\|_F^2+\|\mH\|_F^2+\|\mZ-\mW\|_F^2\\
=&\argmin_{\mZ\in \gS}\|\mZ-\tfrac{\eps}{\sqrt{N}}\mH(\mI_C\otimes \1_N)\|_F^2+\|\mZ-\mW\|_F^2+\|\mH\|_F^2-\|\tfrac{\eps}{\sqrt{N}}\mH(\mI_C\otimes \1_N)\|_F^2\\
=&\argmin_{\mZ\in \gS}\|\mZ-\tfrac{1}{2}(\tfrac{\eps}{\sqrt{N}}\mH(\mI_C\otimes \1_N)+\mW)\|_F^2\\
=&\tfrac{1}{2}(\tfrac{\eps}{\sqrt{N}}\mH(\mI_C\otimes \1_N)+\mW)(\mI_C-\tfrac{1}{C}\1_C\1_C^\top).
\end{aligned}
\end{equation}

Thus, $\Pi_1^\eps (\mH,\mW)=(\tfrac{\eps}{\sqrt{N}}(\mP\otimes \1_N^\top), \mP)$ with $\mP=\tfrac{1}{2}(\tfrac{\eps}{\sqrt{N}}\mH(\mI_C\otimes \1_N)+\mW)(\mI_C-\tfrac{1}{C}\1_C\1_C^\top)$.

For $\mH\in\sR^{p\times CN}$, $\mW\in\sR^{p\times C}$, the projection of $(\mH,\mW)$ onto $\gE_2^\eps$ is

\begin{equation}
\Pi_2^\eps (\mH,\mW)=(\tfrac{\eps}{\sqrt{N}}\vh\1_{CN}^\top, \vh\1_C^\top),
\end{equation}

where $\vh=\tfrac{1}{2C}(\tfrac{\eps}{\sqrt{N}}\mH\1_{CN}+\mW\1_C) $.

We have

\begin{equation}
\begin{aligned}
\end{aligned}
\end{equation}

For $\mH\in\sR^{p\times CN}$, $\mW\in\sR^{p\times C}$, the projection of $(\mH,\mW)$ onto $\gE_3$ is

\begin{equation}
\Pi_3(\mH,\mW)=(\mH(\mI_{CN}-\tfrac{1}{N}\mI_C\otimes \1_N\1_N^\top), 0).
\end{equation}

This can be easily derived by Lemma <ref>.