Spectral Dynamics Of Gradient Flow

Why the Spectrum Describes the Dynamics

Least-squares gradient flow changes a coefficient vector continuously in the direction that most rapidly decreases squared error. In the original coordinates, the coordinates are coupled: changing one coefficient can affect several fitted values, and the gradient in one coordinate can depend on all the others. The singular vectors of the design matrix separate this motion into independent modes. In each mode, the fitted component moves exponentially toward its least-squares target, and the rate is determined by a squared singular value, equivalently an eigenvalue of the design Gram matrix.

We use a fixed design matrix \(A\in\mathbb{R}^{n\times p}\), response vector \(y\in\mathbb{R}^n\), and coefficient vector \(\beta(t)\in\mathbb{R}^p\). The scaling of the loss matters: here it is \(1/(2n)\) times the sum of squared residuals. With this convention, the eigenvalue that determines a mode’s rate is the corresponding eigenvalue of \(A^\mathsf{T}A/n\), not of \(A^\mathsf{T}A\) without the factor \(1/n\).

Coordinate-Wise Solution

The singular-vector coordinates of the coefficients are \(z_j(t)=v_j^\mathsf{T}\beta(t)\). The corresponding fitted-response coordinate is \(u_j^\mathsf{T}A\beta(t)=s_jz_j(t)\). To find their evolution, differentiate the loss:

Taking the inner product of the gradient-flow equation with \(v_j\), and using \(Av_j=s_ju_j\) and \(A^\mathsf{T}u_j=s_jv_j\), gives a scalar differential equation for each \(j\le r\). Its solution also gives the fitted component directly.

To verify the solution, let \(z_j^\star=a_j/s_j\) and solve \(\dot z_j=-\lambda_j(z_j-z_j^\star)\). Multiplying by \(e^{\lambda_jt}\) yields \(\frac{d}{dt}[e^{\lambda_jt}(z_j-z_j^\star)]=0\). Therefore \(e^{\lambda_jt}(z_j(t)-z_j^\star)=z_j(0)-z_j^\star\), which rearranges to the stated expression. Multiplication by \(s_j\) gives the fitted-response formula.

Thus each fitted component is an exponential blend of its initial value and its least-squares target. The time constant is \(1/\lambda_j=n/s_j^2\): a large eigenvalue means rapid adjustment, while a small positive eigenvalue means slow adjustment. If \(v_k\) is in the nullspace of \(A\), then \(Av_k=0\) and \(v_k^\mathsf{T}\dot\beta(t)=0\). Those coefficient coordinates remain at their initial values. They affect the coefficient vector but not the fitted response.

Fitted Residuals and the Best Possible Fit

The full fitted vector has an especially useful form. Let \(P_A\) denote orthogonal projection onto the column space of \(A\). Within that space, the vectors \(u_1,\ldots,u_r\) form an orthonormal basis, so \(P_Ay=\sum_{j=1}^r a_ju_j\). The part of \(y\) orthogonal to that space cannot be produced by any coefficient vector.

For the residual identity, subtract \(y\) from the fitted-response expansion and decompose \(y=P_Ay+(I-P_A)y\). In mode \(j\), the fitted coordinate minus the response coordinate equals \[ e^{-\lambda_jt}s_jz_j(0)+(1-e^{-\lambda_jt})a_j-a_j =-e^{-\lambda_jt}(a_j-s_jz_j(0)). \] The orthogonal component of the residual is always \(-(I-P_A)y\), since every fitted vector \(A\beta(t)\) lies in the column space of \(A\). Squaring the norm adds the squared coordinates of these mutually orthogonal components, giving the loss identity.

This identity proves that the loss decreases to \(\|(I-P_A)y\|_2^2/(2n)\), the least-squares loss. It also shows why convergence need not mean exact interpolation: exact fit is possible precisely when \(y\) lies in the column space of \(A\). Each representable residual mode decays at its own exponential rate; the unrepresentable component does not decay.

Reading the Spectral Formula

The spectral solution supports several useful conclusions directly from the equations.

• Rates are determined by squared singular values. For the chosen loss normalization, the rate is \(s_j^2/n\). Describing a mode’s rate using \(s_j\) alone misses the square and the factor \(1/n\).
• Initialization affects the amplitude, not the rate. The initial fitted component determines \(a_j-s_jz_j(0)\), the amount of residual in that mode at time zero. Its subsequent decay rate remains \(\lambda_j\).
• Nullspace motion is absent from the fit. A coefficient coordinate in the nullspace is constant under this gradient flow because the gradient has no component there.
• Nonzero initialization can change an individual component’s direction of travel. If its initial fitted coordinate lies above its target, that coordinate decreases toward the target; if it lies below, it increases. In both cases its difference from the target decays exponentially.

The same rates can be seen in the prediction-space differential equation. Setting \(r(t)=A\beta(t)-y\), differentiation gives \(\dot r(t)=-AA^\mathsf{T}r(t)/n\). On the column space of \(A\), the eigenvectors are \(u_j\) and the associated eigenvalues are \(s_j^2/n\). On its orthogonal complement, \(AA^\mathsf{T}\) is zero, which accounts for the fixed residual \(-(I-P_A)y\).

Common Pitfall

A frequent error is to say that a singular value \(s_j\) is the rate of its mode, or that a mode with a small singular value converges faster because it has a smaller target coefficient \(a_j/s_j\). Neither claim follows. The differential equation is \(\dot z_j=-(s_j^2/n)(z_j-a_j/s_j)\): the target depends on \(1/s_j\), but the speed of approach is set by \(s_j^2/n\). The exponential factor is \(e^{-s_j^2t/n}\), so a smaller positive singular value gives slower decay of the residual in that mode. The target’s size and the rate are distinct features of the formula.