We study matrix completion via deep matrix factorization (a.k.a. deep linear neural networks) as a simplified testbed to examine how network depth influences training dynamics. Despite the simplicity and importance of the problem, prior theory largely focuses on shallow (depth-2) models and does not fully explain the implicit low-rank bias observed in deeper networks. We identify coupled dynamics as a key mechanism behind this bias and show that it intensifies with increasing depth. Focusing on gradient flow under block-diagonal observations, we prove: (a) networks of depth ≥3 exhibit coupling unless initialized diagonally, and (b) convergence to rank-1 occurs if and only if the dynamics is coupled -- resolving an open question by Menon (2024) for a family of initializations. We also revisit the loss of plasticity phenomenon in matrix completion (Kleinman et al., 2024), where pre-training on few observations and resuming with more degrades performance. We show that deep models avoid plasticity loss due to their low-rank bias, whereas depth-2 networks pre-trained under decoupled dynamics fail to converge to low-rank, even when resumed training (with additional data) satisfies the coupling condition -- shedding light on the mechanism behind this phenomenon.
Figures & tables
Figure 1: ( a ) Examples of bipartite graphs corresponding to observation patterns of MD (disconnected) and MC (connected). ( b-c ) Training results showing effective rank ( cf . Roy and Vetterli (2007) ) for completing rank-1 matrices MD and MC , respectively. The rank-1 ground truth matrices were generated via uv⊤ , where u,v∈R2 with entries sampled i.i.d. from a standard normal distribution. We initialized each layer’s entries by sampling from a Gaussian distribution with mean zero and standard deviation α , chosen to ensure the initial scale of the product matrix WL:1(0) is approximately invariant to depth L . Each result shows an average of 300 independent random trials.
Figure 2: Singular values of WL:1(∞) (numerically obtained from Theorem 3.3 ) against initialization scale αL , for the diagonal observation task where s=1 . Solid lines represent the largest singular value σ1 ; dashed lines denote the other (identical) singular values σr for r≥2 . For finite m , these results illustrate that both greater depth L and a smaller initial scale α strengthen the low-rank bias, in contrast to the L=2 case. Conversely, a very large m ( m=1010 ), approximating an αId (rank- d ) initialization, leads to decoupled dynamics and a full-rank solution, independent of both L and α .
Figure 3: Experiments use a 100×100 rank-5 ground-truth matrix. Pre-training utilizes 2000 randomly sampled entries ( Ωpre ; ∣Ωpre∣=2000 ), while post-training adds 1000 more, forming Ωpost ( Ωpre⊂Ωpost ; ∣Ωpost∣=3000 ). The top row of panels displays effective rank, and the bottom row shows reconstruction error, both measured at convergence. The leftmost panels depict training on Ωpre , and the rightmost on Ωpost , both starting from random Gaussian initialization. The middle panels show warm-start training on Ωpost , initialized from converged pre-trained models with Ωpre .
Appendix figures & tables30 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 4: The left panel shows the averaged effective rank of all possible connected patterns as a function of the initial scale αL . The right panel displays the averaged effective rank of all possible disconnected patterns.
Figure 5: Singular values of WL:1(∞) (numerically obtained from Theorem 3.3 ) against initialization scale αL for the block-diagonal observation task. Solid lines represent the largest singular value σ1 ; dashed lines denote the identical singular values σi for i∈{2,…,n} . Note that σj for j∈{n+1,…,d} are all zero. For finite m , these results show that both greater depth L and a smaller initial scale α strengthen the low-rank bias, in contrast to the L=2 case. Conversely, when m is extremely large (e.g., m=1010 ), approximating an αId rank d initialization, the dynamics decouple and cannot achieve the minimal low-rank solution, regardless of L or α .
Figure 6: Numerical conditions identical to those in Figure 2 , except with ground truth value w∗=10 and dimension d=10 where the block size is s=1 .
Figure 7: Numerical conditions identical to those in Figure 2 , except with ground truth value w∗=0.1 and dimension d=10 where the block size is s=1 .
Figure 8: Numerical conditions identical to those in Figure 2 , except with ground truth value w∗=1 and dimension d=3 where the block size is s=1 .
Figure 9: Gradient descent experiments conducted under conditions identical to those in Figure 8 .
Figure 10: (a) Effective rank for the initialization scheme in ( 7 ). The x-axis denotes the parameter m , which controls the initial rank characteristics of the model, while the y-axis represents the corresponding effective rank after convergence. (b) Effective rank distributions for Gaussian initialization. The results are from 1000 independent trials, sorted by their converged effective rank. The x-axis denotes the sorted trial index (from lowest to highest converged rank), and the y-axis represents the corresponding effective rank after convergence.
Figure 11: Limiting stable rank (y-axis) as a function of αL (x-axis) under noisy diagonal observations. Dashed lines indicate the noiseless baseline, and solid lines indicate the noisy case. The noise standard deviation is set to σ=0.01 (leftmost), σ=0.05 (middle), and σ=0.1 (rightmost). The depth dependent low-rank bias persists and follows a trend similar to the noiseless setting.
Figure 12: Limiting stable rank (y-axis) as a function of the initialization scale (x-axis) under non-equal diagonal observations. The low-rank bias induced by coupled training dynamics persists and closely matches the behavior in the equal-diagonal setting described in Theorem 3.3 and Figure 2 .
Figure 13: Final stable rank as a function of depth. Each panel corresponds to a different initialization scale. Results are obtained using SGD.
Figure 14: Final stable rank as a function of depth. Each panel corresponds to a different initialization scale. Results are obtained using GD with momentum.
Figure 15: Final stable rank as a function of depth. Each panel corresponds to a different initialization scale. Results are obtained using Adam.
Figure 16: Final stable rank as a function of depth. Each panel corresponds to a different initialization scale. Results are obtained using RMSProp.
Figure 17: Final stable rank as a function of depth. Each panel corresponds to a different initialization scale. Results are obtained using Adagrad.
Figure 18: We train CIFAR-10 with ResNet models ranging from 18 to 101 layers using SGD with momentum, averaging results over five independent runs with 95% confidence intervals. The leftmost plot reports the training accuracy, the middle plot the test accuracy, and the rightmost plot the average effective rank. As depth increases, the average effective rank decreases.
Figure 19: The results for CIFAR-100 with ResNet-18 to 101, under the same conditions as in Figure 18 .
Figure 20: The results for CIFAR-10 with VGG-11 to 19, under the same conditions as in Figure 18 .
Figure 21: The results for CIFAR-100 with VGG-11 to 19, under the same conditions as in Figure 18 .
Figure 22: We train CIFAR-10 with ResNet models ranging from 18 to 101 layers using Adam, averaging results over five independent runs with 95% confidence intervals. The leftmost plot reports the training accuracy, the middle plot the test accuracy, and the rightmost plot the average effective rank. As depth increases, the average effective rank decreases.
Figure 23: The results for CIFAR-100 with ResNet-18 to 101, under the same conditions as in Figure 22 .
Figure 24: The results for CIFAR-10 with VGG-11 to 19, under the same conditions as in Figure 22 .
Figure 25: The results for CIFAR-100 with VGG-11 to 19, under the same conditions as in Figure 22 .
Figure 26: We train CIFAR-10 with ResNet models ranging from 18 to 101 layers using RMSProp, averaging results over five independent runs with 95% confidence intervals. The leftmost plot reports the training accuracy, the middle plot the test accuracy, and the rightmost plot the average effective rank. As depth increases, the average effective rank decreases.
Figure 27: The results for CIFAR-100 with ResNet-18 to 101, under the same conditions as in Figure 26 .
Figure 28: The results for CIFAR-10 with VGG-11 to 19, under the same conditions as in Figure 26 .
Figure 29: The results for CIFAR-100 with VGG-11 to 19, under the same conditions as in Figure 26 .
Figure 30: Left: training loss (log scale). Right: average stable rank across all layers except the last. Solid lines correspond to Gaussian initialization and dashed lines to identity-based initialization. Gaussian initialization (corresponding to coupled training dynamics) converges to noticeably lower rank than identity-based initialization (corresponding to more decoupled training dynamics).
Figure 31: Experimental results for a 2×2 rank-1 ground-truth matrix W∗ with w11∗=w22∗=1 and w12∗=0.5 (implying w21∗=2 for rank-1 structure). Models, initialized according to ( 7 ), are first pre-trained on diagonal entries. After achieving zero-loss convergence in pre-training, the off-diagonal element w12∗ is introduced, and models are subsequently trained on combined diagonal and off-diagonal observations. The plots display: (Left and Middle) effective rank under different settings; (Right) converged value of w21(∞) . Key observations: (1) Warm-starting with a model that converged to a high-rank solution during pre-training tends to maintain this high rank, even when presented with the same subsequent observations as a cold-started model. (2) In the theoretically analyzed m=∞ case, w21(∞)<0 is observed, which correlates with the highest effective rank.
Figure 32: Experimental conditions identical to those in Figure 31 , except with ground truth value w12∗=10 . The model have to predict w21∗ as 0.1
Figure 33: Experimental conditions identical to those in Figure 31 , except with ground truth value w11∗=1,w22∗=2 , and w12∗=0.5 . The model have to predict w21∗ as 4.
Implicit bias studies how training algorithms and model structure favor particular solutions, representations, and generalization behaviors without explicit penalties.
Neural collapse (NC) describes the structured geometry that emerges in the features and weights of trained classifiers. Recent theory suggests NC can be suboptimal in deep architectures, attributing this to an explicit low-rank bias from L2 regularization. We study the deep unconstrained feature model (UFM)-equivalent to a deep linear network with orthogonal inputs-trained without regularization, to isolate how gradient descent and depth alone shape NC. We show that depth induces an implicit low-rank bias: low-rank matrices propagate norm more efficiently through successive multiplications, promoting low-rank alternatives to NC. These alternatives, we argue, correspond to softmax codes: max-margin solutions previously found in width-bottlenecked networks. Analyzing training dynamics under spectral initialization, we identify an early-time repulsion among singular values that drives low-rank emergence, and characterize how depth shrinks NC's basin of attraction. Finally, we show that some effects act in the opposite direction: for randomly initialized networks, increasing width biases training toward higher-rank solutions. Our results provide the first asymptotic and dynamic characterization of implicit bias in deep UFMs trained with unregularized multiclass cross-entropy.
Connall Garrod, Jonathan P. Keating, Christos Thrampoulidis
Mathematical Institute, University of Oxford · Department of Electrical and Computer Engineering, University of British Columbia
This paper studies the stability of low-rank implicit regularization in perturbed deep matrix factorization, where the target matrix is corrupted by a noise matrix. We first derive sufficient spectral conditions under which gradient descent exhibits a low-rank phase in the noiseless setting. These conditions show how the target spectrum, initialization, and step size jointly determine the existence of a nonempty low-rank interval. We then analyze the perturbed gradient descent dynamics, proving convergence guarantees and quantifying how the perturbation affects iteration complexity and eigenvalue recovery. Finally, we show that the low-rank phase persists under perturbation, with explicit dependence on the perturbation size. Numerical experiments support the theoretical findings.
Jingzhe Wang, Hung-Hsu Chou
Department of Informatics and Networked Systems, University of Pittsburgh · Department of Mathematics, University of Pittsburgh
Weight decay is ubiquitous in training deep neural network architectures. Its empirical success is often attributed to capacity control; nonetheless, our theoretical understanding of its effect on the loss landscape and the set of minimizers remains limited. In this paper, we show that ℓ2-regularized deep matrix factorization/deep linear network training problems with squared-error loss admit a unique end-to-end minimizer for all target matrices subject to factorization, except for a set of Lebesgue measure zero formed by the depth and the regularization parameter. This observation reveals fundamental properties of the loss landscape of regularized deep matrix factorization problems: the Hessian spectrum is constant across all minimizers of the regularized deep scalar factorization problem with squared-error loss. Moreover, we show that, in regularized deep matrix factorization problems with squared-error loss, if the target matrix does not belong to the Lebesgue measure-zero set, then the Frobenius norm of each layer is constant across all minimizers. This, in turn, yields a global lower bound on the trace of the Hessian evaluated at any minimizer of the regularized deep matrix factorization problem. Furthermore, we establish a critical threshold for the regularization parameter above which the unique end-to-end minimizer collapses to zero.