Generalization Error Estimation for Primal--Dual Algorithms in Non-Smooth Regression

Abstract: This paper studies trajectory-wise estimation of generalization error for primal--dual algorithms in non-smooth regression. Motivating examples include (\ell_1)-penalized least absolute deviations regression and square-root Lasso regression, where the data-fitting loss is non-differentiable and existing risk estimators for gradient-type optimization paths do not apply directly. We develop a general recursive framework that includes the Chambolle--Pock algorithm and related primal--dual splitting methods. We estimate risk by correcting each in-sample fitted value with a weighted combination of past dual iterates. The ideal weights are Stein derivative contractions and depend on the design covariance. We construct replacement weights from observable derivative contractions of the fitted-signal trajectory, yielding a covariance-free, data-driven correction. For high-dimensional Gaussian designs and fixed finite iteration horizon, we prove finite-sample guarantees for both estimators. For square-root ridge, we further establish a matched-Gaussian universality result beyond Gaussian designs. Numerical experiments show that the proposed estimators accurately track the out-of-sample risk along finite optimization paths.

• The paper develops trajectory-wise leave-one-out estimators for non-smooth primal–dual algorithms, correcting training predictions with weighted past dual iterates rather than primal iterates.
• Under Gaussian designs and fixed iteration horizons, the covariance-dependent estimator achieves an expected error rate of O(n^-1/2), while its observable covariance-free counterpart achieves O(n^-1/4).
• The results support asymptotically oracle early stopping and establish matched-Gaussian universality for square-root ridge, while leaving growing horizons and nonlinear proximal maps as open problems.

Overview

The paper develops trajectory-wise estimators of generalization error for primal–dual optimization algorithms applied to non-smooth regression problems, with the Chambolle–Pock recursion as the central example (2608.13870). The motivating settings are $\ell_1$-penalized least absolute deviations (LAD) regression and square-root Lasso, where the data-fitting loss is non-differentiable and existing risk-estimation methods for gradient-type paths do not apply directly. The target is the conditional prediction risk $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$ of each iterate along a finite trajectory of length $T$, which supports early stopping without rerunning the algorithm.

The construction follows a leave-one-out (LOO) logic: since $b^{t,-i}$ is independent of $(x_i,y_i)$, the LOO empirical risk approximates $r_t$. Gaussian Stein's formula converts the difference between full-data and deleted fitted values into a derivative contraction, yielding a corrected training error of the form

$$\widetilde r_t = \frac1n \sum_{i=1}^n \ell\Bigl(y_i,\; x_i^\top b^t - \sum_{s=1}^t w_{t+1,s}\, u_i^s\Bigr),$$

where the correction subtracts a weighted combination of past dual iterates $u^s$. This structural feature distinguishes the paper from prior work on gradient descent [bellec2024uncertainty] and proximal SGD [tan2024estimating], where corrections were built from the primal trajectory; here the correction weights live in a strictly lower-triangular matrix $W$ defined by derivative contractions of the coupled primal–dual recursion.

Algorithmic framework

The authors consider penalized regression $\min_b f_y(Xb) + g(b)$ with possibly non-differentiable $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$0 and non-smooth $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$1. The Chambolle–Pock updates involve proximal maps of $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$2 and $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$3; for LAD the dual proximal map is coordinatewise clipping to $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$4, and for square-root loss it is Euclidean projection onto the unit ball. After rewriting in residual form with normalized residual vectors $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$5 and dual-gradient vectors $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$6, the recursion is abstracted into a general history-dependent iteration

$$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$7

with Lipschitz update maps. This class contains Chambolle–Pock, proximal-gradient methods for smooth losses, and linearized ADMM. Jacobian blocks of these maps, taken as modified almost-everywhere derivatives so that a Lipschitz chain rule holds, drive all subsequent constructions.

Covariance-dependent estimator

Under i.i.d. Gaussian designs with covariance $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$8 satisfying eigenvalue bounds, proportional asymptotics $$r_t = \mathbb{E}[\ell(y_{\rm new}, x_{\rm new}^\top b^t)\mid (X,y)]$$9, bounded signal, and either 1-Lipschitz test losses or centered squared error, the main theorem establishes

$T$0

for every fixed iterate $T$1. The proof decomposes the total error into five residuals: LOO-to-risk approximation, empirical fluctuation of centered deleted losses (handled via a leave-two-out decoupling identity), second-order Stein comparison, contraction of derivatives into the weight matrix $T$2, and replacement of $T$3 by its data-driven counterpart. Notably, the Lipschitz-loss case requires only first moments of the noise, whereas the squared-error case uses second moments.

Covariance-free estimator

The population weight matrix $T$4 depends on the unknown design covariance $T$5. To remove this dependence, the authors define two observable matrices built from empirical derivative contractions of the trajectory: $T$6, lower triangular with deterministic diagonal entries $T$7, and $T$8. The covariance-free correction is then $T$9, computable by forward substitution in $b^{t,-i}$0 after propagating Hutchinson probe vectors through the linearized recursion at cost $b^{t,-i}$1. Because $b^{t,-i}$2's diagonal is deterministic, invertibility holds unconditionally — an improvement over the corresponding matrix in the SGD analysis of Tan and Bellec, which was invertible only on a high-probability event.

The guarantee is

$b^{t,-i}$3

on the event $b^{t,-i}$4 where $b^{t,-i}$5 is controlled (probability at least $b^{t,-i}$6), giving overall consistency of $b^{t,-i}$7 for $b^{t,-i}$8 at rate $b^{t,-i}$9. The authors state plainly that they do not expect the $(x_i,y_i)$0 rate to be optimal; it arises from the weighted stability bound $(x_i,y_i)$1 governing the average squared change in corrected fitted values. A corollary shows that minimizing $(x_i,y_i)$2 over the first $(x_i,y_i)$3 iterates yields an early-stopping rule whose risk exceeds the oracle-optimal risk on the path by $(x_i,y_i)$4 for fixed $(x_i,y_i)$5.

Universality beyond Gaussian designs

For square-root ridge with squared test loss, the paper proves that both the true risk and both estimators have identical limits under a non-Gaussian design $(x_i,y_i)$6 (sub-Gaussian entries) and its covariance-matched Gaussian counterpart $(x_i,y_i)$7:

$(x_i,y_i)$8

Consequently both estimators remain consistent under such non-Gaussian designs in this setting. The proof exploits the linearity of the ridge proximal map: iterates admit a finite Krylov representation $(x_i,y_i)$9 with polynomials driven by $r_t$0, and all fixed-time quantities reduce to finitely many "admissible primitive scalars" (trace functionals, signal alignments, quadratic forms). Matched-Gaussian universality of these scalars follows from Chatterjee's Lindeberg principle plus conditional variance bounds for noise quadratic forms. Two restrictions are explicit: the theorem requires ridge penalties (linear proximal map), and the estimator universality requires a no-boundary condition ensuring the limiting trajectory avoids the nondifferentiability kink of the square-root projection $r_t$1; true-risk universality needs only continuity and holds without this condition.

Numerical experiments

Simulations use $r_t$2, $r_t$3, $r_t$4, Toeplitz covariance $r_t$5, sparse signals with $r_t$6, and Chambolle–Pock constants $r_t$7. Eight procedures combine LAD or square-root loss with ridge, Lasso, MCP, or SCAD penalties; MCP and SCAD are included only as empirical checks using single-valued proximal maps, with no convergence guarantee claimed for these non-convex objectives. Across all panels, both $r_t$8 and $r_t$9 track the true test risk and the location of its minimum along the computed path, including trajectories that have not stabilized by iteration 300. For square-root ridge, Gaussian and covariance-matched uniform designs produce nearly overlapping curves for the true risk and both estimators, matching the universality theorem empirically; analogous comparisons for other loss–penalty combinations suggest, but do not prove, broader universality.

Limitations and open questions

Several limitations are conceded explicitly. The main risk bounds are proved under Gaussian designs and hold for fixed iteration horizon $$\widetilde r_t = \frac1n \sum_{i=1}^n \ell\Bigl(y_i,\; x_i^\top b^t - \sum_{s=1}^t w_{t+1,s}\, u_i^s\Bigr),$$0 as $$\widetilde r_t = \frac1n \sum_{i=1}^n \ell\Bigl(y_i,\; x_i^\top b^t - \sum_{s=1}^t w_{t+1,s}\, u_i^s\Bigr),$$1 grows; extension to growing horizons is open. The covariance-free rate $$\widetilde r_t = \frac1n \sum_{i=1}^n \ell\Bigl(y_i,\; x_i^\top b^t - \sum_{s=1}^t w_{t+1,s}\, u_i^s\Bigr),$$2 is acknowledged as likely suboptimal. The universality theorem is confined to square-root ridge with squared loss because the proof relies on the spectral structure of the linear ridge proximal map; whether universality extends to nonlinear proximal maps such as square-root Lasso, or to broader design classes suggested by the uniform-design simulations, remains unestablished. The framework also assumes convex penalties for its guarantees, so the MCP/SCAD experiments carry no theoretical backing.

Conclusion

This paper extends pathwise generalization-error estimation from gradient-type algorithms to coupled primal–dual recursions for non-smooth regression. Its two estimators — one requiring the design covariance, one fully data-driven via an always-invertible triangular system — achieve $$\widetilde r_t = \frac1n \sum_{i=1}^n \ell\Bigl(y_i,\; x_i^\top b^t - \sum_{s=1}^t w_{t+1,s}\, u_i^s\Bigr),$$3 and $$\widetilde r_t = \frac1n \sum_{i=1}^n \ell\Bigl(y_i,\; x_i^\top b^t - \sum_{s=1}^t w_{t+1,s}\, u_i^s\Bigr),$$4 accuracy respectively per fixed iterate, support asymptotically oracle early stopping over finite horizons, and, for square-root ridge, enjoy matched-Gaussian universality. The results provide a statistically grounded basis for monitoring and tuning iterative solvers in settings where smooth-loss theory does not apply.