Estimating implicit regularization in deep learning

In Figure 3, $L(\theta )=\frac{1}{2}(\theta -2{)}^2$ and $R(\theta ;\lambda )=\lambda {\theta }^2$ with $\lambda \ge 0$. Their derivatives give the residual $(\widehat{\theta }-2)+2\lambda \widehat{\theta }$. For $\widehat{\theta }\ne 0$, minimizing its square has the closed-form solution:

“Estimate $\lambda$” animates from the current guess to this exact least-squares solution; it is not a simulation of an iterative optimizer.

For $0<\widehat{\theta }\le 2$, the gradients can cancel exactly. Outside that interval, nonnegative ridge cannot explain the observed weight as an optimum: its best fit leaves a residual. At $\widehat{\theta }=0$, every $\lambda$ supplies the same zero penalty gradient, so $\lambda$ is unidentifiable. A fitted coefficient is an explanation within the chosen family, not evidence that this mechanism generated the weights.

This is a simplified version of the paper’s explicit elastic-net recovery experiment. Here we estimate one ${\ell }_2$ coefficient in a small, single-input network. The paper jointly estimates smoothed ${\ell }_1$ and ${\ell }_2$ coefficients in a 200-unit ReLU network with 10 inputs, across repeated datasets and penalty settings.

We use mean squared error $L=\frac{1}{n}{\sum }_i(f_{\theta }(x_i)-y_i{)}^2$ and $R=\lambda \Vert \theta {\Vert }^2$, penalizing all parameters, including biases (3 × width + 1 parameters). Full-batch gradient descent runs for up to 20,000 steps. The learning rate is 0.04 × 16 / width (0.04, 0.02, or 0.01) to keep wider networks stable. The data and initialization are reproducible for each seed.

$$\widehat{\lambda }=max(0,\frac{-{\theta }^{\top }\nabla L}{2\Vert \theta {\Vert }^2})$$

This is a nonnegative least-squares projection of $-\nabla L$ onto $2\theta$. “Gradient explained” is $1-\frac{\Vert \nabla L+2\widehat{\lambda }\theta {\Vert }^2}{\Vert \nabla L{\Vert }^2}$. It is a fit diagnostic, not predictive $R^2$. Near-zero weights make the coefficient unidentifiable; a near-zero loss gradient makes the explained fraction uninformative. We display a dash for unavailable diagnostics.

For ReLU, we use a zero derivative at a kink. Tanh uses its smooth derivative, $1-{\tanh }^2(z)$. Finite training time and nonsmoothness can leave a residual even when the loss appears stable. This experiment demonstrates explicit-penalty recovery, not proof that an arbitrary implicit bias is ridge regularization.

Let $H=\frac{X^{\top }X}{n}=Vdiag(h_i)V^{\top }$ and use the loss $L(\theta )=\frac{\Vert X\theta -y{\Vert }^2}{2n}$. Full-batch gradient descent starts at zero with step size $\eta$ satisfying $0<\eta h_i<1$. At step $t\ge 1$, its weights minimize $L(\theta )+{\theta }^{\top }{\Lambda }^{(t)}\theta$, where

$${\Lambda }^{(t)}=\frac{1}{2}Vdiag\left[\frac{h_i(1-\eta h_i{)}^t}{1-(1-\eta h_i{)}^t}\right]V^{\top }$$

The eigenvectors give the penalty’s directions, and its eigenvalues decrease over training. The matrix is generally not diagonal in the original input coordinates. The factor ½ reflects the half-MSE loss convention used here; the paper writes its theoretical penalty with the loss $\Vert X\theta -y{\Vert }^2/n$.

For any estimated symmetric matrix $\widehat{\Lambda }$, the explicit solution is ${\widehat{\theta }}^R=(H+2\widehat{\Lambda }{)}^{-1}\frac{X^{\top }y}{n}$, provided $H+2\widehat{\Lambda }$ is positive definite. We solve this system independently to test weight recovery. The estimator receives weights and loss gradients, not the theoretical matrix or its eigenvectors.

Absolute matrix error is the Frobenius norm of the difference between the estimated and theoretical matrices. Zero means exact agreement. Gradient matching can be exact even when the matrix is not identified; agreement with theory and recovery of the weights test different claims.

These comparisons correspond to Figure 4 of the paper: full-matrix recovery in panel B, diagonal estimation in panel C, and a scalar ridge summary in panel E. The curvature ratio is the Hessian condition number: the largest curvature divided by the smallest.

We use 160 observations and 10 inputs; the paper uses 1000 observations. We construct X with Hessian eigenvalues ranging from 2 to $2/\kappa$, use coefficient seed 7 and Gaussian noise with standard deviation 0.3, and run full-batch GD with $\eta =0.2$. Changing $\kappa$ rebuilds the design. All endpoints use the same X and the same selected training step; only the true coefficients and noise are resampled.

Diagonal: fit ten signed coefficients in the original input coordinates, ${\widehat{\Lambda }}_{ii}=-\frac{{\partial }_{i}L}{2{\widehat{\theta }}_i}$. These are not coefficients in the Hessian eigenbasis. A zero weight makes its coefficient unidentifiable. Even with an exact gradient match, a unique minimum requires $H+2\widehat{\Lambda }$ to be positive definite; the display checks this condition.

Full symmetric: fit all 55 independent entries by minimizing ${\sum }_k\Vert 2\widehat{\Lambda }{\widehat{\theta }}_k+\nabla L_k{\Vert }^2$ over symmetric matrices. If the system is underdetermined, we select the solution with smallest Frobenius norm. Ten linearly independent endpoints identify the matrix; fewer leave unobserved directions. No theoretical eigenvectors are supplied. The endpoints share a stopping time; the paper also studies estimation across different stopping times.

Scalar ridge: fit $\widehat{\Lambda }=\widehat{\lambda }I$, using $\widehat{\lambda }=max(0,\frac{-{\widehat{\theta }}^{\top }\nabla L}{2\Vert \widehat{\theta }{\Vert }^2})$. The right plot always shows this separate scalar summary, so it has the same meaning under every matrix selection. Its decline need not imply that scalar ridge recovers the full matrix or the stopped weights.

The matrix heatmaps compare against the independent theoretical formula above. Random initialization changes the starting weights for every endpoint while keeping the data fixed; the zero-start theoretical comparison is hidden in that case.