The Normalized Maximum Likelihood for Regular Non-Smooth Models: Measure-Theoretic Foundations and Geometric Sampling

Abstract: The Normalized Maximum Likelihood (NML) codelength, or stochastic complexity, represents a principled criterion for universal coding. While recent coarea-based formulations provided a calculation method for smooth models, this framework collapses for the non-smooth estimators ubiquitous in modern machine learning (e.g., Lasso, Sparse SVMs). In this work, we provide a rigorous framework for computing the NML for regular path-differentiable Lipschitz (PDL) estimators. By applying classical geometric measure theory and bridging the coarea formula with conservative Jacobians, we prove that the stochastic complexity for non-smooth models is well-posed and theoretically consistent with the outputs of modern Automatic Differentiation. To compute this quantity exactly, we introduce the Propose-and-Project Metropolis-Hastings (PDL-PPMH) sampler, a geometric MCMC algorithm capable of traversing the non-differentiable level sets of the maximum likelihood estimator. We theoretically justify its components, including a stochastic tangent space proposal and a provably convergent non-smooth projection solver. We demonstrate the method's robustness by sampling from a high-dimensional Lasso posterior ($P=2000$), while simultaneously quantifying the computational scaling that governs the trade-off between exactness and mixing time. Crucially, we empirically demonstrate that our exact NML criterion provides a highly data-efficient alternative to cross-validation, achieving statistically indistinguishable predictive optima without requiring data splitting. Altogether, our work paves the way for the theoretical analysis of the NML codelength for regular non-smooth models.

• The paper extends Normalized Maximum Likelihood to path-differentiable Lipschitz estimators, proving well-posedness, uniqueness, minimax optimality, and consistency with pathwise automatic differentiation.
• It shows that regular non-smooth models retain the classical k/2 log N complexity penalty, while experiments confirm a slope of 2.5 for k = 5 and exclude singular models that require RLCT methods.
• The proposed Propose-and-Project Metropolis–Hastings sampler traverses non-differentiable MLE level sets with computational cost scaling as O((N+k)^3), remaining effectively independent of feature count and matching or outperforming cross-validation in small samples.

Overview

This paper extends the Normalized Maximum Likelihood (NML) framework of universal coding to a class of estimators that prior theory could not handle: path-differentiable Lipschitz (PDL) maps, which include the Lasso, Sparse SVMs, and other $L_1$-regularized estimators whose maximum likelihood estimators (MLEs) are non-smooth. Building on the coarea-based reformulation of stochastic complexity introduced by Suzuki et al. [Suzuki_2024], the authors replace the classical Fréchet Jacobian with conservative Jacobians from nonsmooth analysis [Bolte2021Conservative] and prove that the resulting NML integral is well-posed, unique, minimax optimal, and consistent with the output of modern pathwise Automatic Differentiation (AD). On the computational side, they introduce the Propose-and-Project Metropolis–Hastings (PDL-PPMH) sampler, a geometric MCMC algorithm that traverses non-differentiable level sets of the MLE exactly. Experiments on high-dimensional Lasso problems ($P = 2000$) validate both the theoretical scaling laws and a practical claim of considerable interest: that exact NML model selection can match or exceed cross-validation in data efficiency.

Extending the coarea framework to PDL estimators

The starting point is the classical NML universal model, whose normalization constant — the stochastic complexity — is an integral over the data space of the maximized likelihood. Suzuki et al. reformulated this via the coarea formula as an integral over parameter space, with the MLE's pushforward density expressed as an integral over its level sets weighted by the inverse Jacobian factor $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$. This construction requires smoothness.

The paper's central theoretical move is to show that smoothness can be dropped. Three results carry the argument:

• Rectifiability: By Federer's theorem, level sets of any Lipschitz map are $(N-k)$-rectifiable for almost all parameter values; via Kirszbraun extension and intersection with a Borel domain, the same holds for the MLE's level sets. This licenses integration against Hausdorff measure.
• Coarea formula with conservative Jacobians (Theorem 3): For a PDL map with full-rank conservative selections $G_x \in D(x)$, the identity $$\int u(x)\,j(x)\,dL(x) = \int_{\mathbb{R}^k}\int_{m^{-1}(z)} u\,dH\,dz$$ holds. The proof rests on Rademacher's theorem: the Clarke Jacobian collapses to the singleton $\{\nabla m(x)\}$ almost everywhere, so the conservative factor equals the classical one Lebesgue-a.e., and the substitution under the integral is exact.
• Generalized NML estimator PDF (Theorem 4): Under Lipschitz continuity, path-differentiability, and a local full-rank condition on all elements of the Clarke Jacobian, the density $$p[\#\mu_{\theta_0}](\theta'') = \int_{m^{-1}(\theta'')} p[x|\theta_0]/j(x)\,dH$$ is well-posed, unique, and reduces to the classical formula at differentiability points.

A consequence with practical weight is Theorem 5: because pathwise AD systems return selections from the Clarke Jacobian, and these coincide with the Fréchet derivative almost everywhere, evaluating the Jacobian factor with standard backpropagation yields integrands that agree with the measure-theoretic definition almost everywhere. The framework therefore provides "mathematical clearance" for computing non-smooth stochastic complexity with off-the-shelf AD tooling. The authors are careful to note a residual caveat: at individual kinks, different AD implementations may select matrices with different singular-value products, so pointwise values of $j(G)$ need not be unique — though this ambiguity lives on a measure-zero set and does not affect the integral.

Information-theoretic guarantees

Two results establish that non-smoothness carries no information-theoretic penalty for regular models. First, extending Shtarkov's characterization, the generalized NML distribution remains the unique minimax-optimal universal code, with regret $\log C_N(\Theta)$; since the non-differentiable set has Lebesgue measure zero by Rademacher's theorem, the coarea normalization equals the classical Lebesgue normalization exactly.

Second, the paper derives the asymptotic expansion under a Uniform Surjectivity Constraint Qualification (USCQ) and compactness:

$P = 2000$0

where $P = 2000$1 is the finite union of lower-dimensional manifolds of non-differentiability. The proof strategy splits the parameter domain into a smooth interior handled by Laplace approximation and an $P = 2000$2-tube around the singular set with $P = 2000$3. Because semi-algebraic sets admit finite $P = 2000$4 stratification, the tube volume scales as $P = 2000$5, which vanishes relative to the leading term. The implication is sharp: regular non-smooth models inherit the classical $P = 2000$6 penalty without divergent curvature artifacts, in explicit contrast to Singular Learning Theory, where rank-deficient Jacobians force resolution of the Real Log Canonical Threshold (RLCT). This delineates the paper's scope precisely — overparameterized ReLU networks violate the USCQ on positive-measure sets and fall strictly outside the framework.

The PDL-PPMH sampler

Computing the inner level-set integral requires sampling from densities proportional to $P = 2000$7 on non-smooth manifolds, where ambient methods suffer from the curse of dimensionality (hitting probability of thickened level sets scaling as $P = 2000$8). The proposed solution is a propose-and-project scheme on the level set itself:

• Proposal: A Stochastic Jacobian Oracle (SJO) based on gradient sampling draws multiple Fréchet derivatives from an $P = 2000$9-ball and uses the kernel of a randomly selected matrix as the tangent space. This design responds directly to a diagnosed failure mode: the deterministic AD selection map is provably discontinuous exactly where the Clarke subdifferential is not a singleton (Theorem 8), which would destabilize the proposal mechanism.
• Projection: Solving a non-smooth constrained least-squares problem onto the level set, with convergence guaranteed by prox-regularity, the Kurdyka–Łojasiewicz property of semi-algebraic constraints, strong metric regularity of the KKT map, and superlinear local convergence of semismooth Newton methods.
• Acceptance: A Metropolis–Hastings ratio corrected by Radon–Nikodym derivatives of the projection map, computed stably via BD-regularity of the generalized KKT matrix.

The main ergodicity result (Theorem 9) establishes strong Feller continuity of the expected kernel and geometric ergodicity via the Weak Harris Theorem, using compactness of the Clarke subdifferential to bound the variance introduced by randomization. Two honest qualifications accompany this guarantee. First, practical iterative projections terminate at tolerance $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$0, breaking detailed balance; perturbation analysis bounds the stationary distribution error as $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$1, making bias linearly controllable but not zero. Second, geometric ergodicity does not imply fast mixing — the spectral gap in high dimensions can be practically demanding, positioning the method as a rigorous baseline rather than a fast general-purpose tool.

Computational scaling

An active-set exploitation reduces determinant computations from naive $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$2 to $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$3, decoupling cost from ambient dimension $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$4. Empirically, wall-clock time per MCMC step is flat as $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$5 ranges from 100 to 2000 at fixed $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$6, while aligning closely with the theoretical $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$7 curve as $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$8 varies from 100 to 1000. The bottleneck is therefore sample size, not feature count — favorable for "large $J_K(x) = \sqrt{\det(\nabla m(x)\nabla m(x)^T)}$9, small $(N-k)$0" regimes but limiting for large-$(N-k)$1 problems.

Numerical validation

Experiments use sparse linear regression with Toeplitz-correlated designs ($(N-k)$2), true sparsity $(N-k)$3, SNR 3.0, and 50,000 samples per chain across 120 regularization parameters. Key findings:

• Scaling laws: Dimension invariance in $(N-k)$4 and cubic scaling in $(N-k)$5 match theory.
• Sampler health: Despite ridge-like posterior geometry from correlated features, autocorrelation decays to near-zero within a few lags, and acceptance rates remain healthy in sparse regimes, degrading only as $(N-k)$6.
• Asymptotics: The computed complexity against $(N-k)$7 is linear with empirical slope matching the theoretical $(N-k)$8, and the residual converges to a stable constant — direct empirical confirmation that no divergent non-smooth penalty arises.
• Model consistency: Truth recovery probability approaches 1.0 at the optimal $(N-k)$9; the MSE minimum coincides with the NML optimum and peak truth recovery. In dense regimes ($G_x \in D(x)$0 approaching $G_x \in D(x)$1), high-frequency volatility reflects genuine ill-conditioning of $G_x \in D(x)$2 rather than algorithmic failure.

Comparison against baselines

Against BIC, asymptotic NML, and 5-fold cross-validation, two claims stand out. Standard criteria apply penalties that are step functions of the active set size $G_x \in D(x)$3, whereas exact PDL-NML yields a continuous geometric volume capturing subspace ill-conditioning; at small $G_x \in D(x)$4 the two diverge substantially, converging only as $G_x \in D(x)$5 grows. More consequentially, in the moderate regime ($G_x \in D(x)$6), the NML-selected model achieves predictive error statistically indistinguishable from cross-validation ($G_x \in D(x)$7, Welch's t-test) without any data splitting, and at $G_x \in D(x)$8 it identifies a model with statistically significantly lower MSE than 5-fold CV ($G_x \in D(x)$9), because CV sacrifices 20% of observations while NML uses all data analytically. This is the paper's strongest practical claim: exact stochastic complexity is not merely a theoretical construct but a data-efficient alternative to cross-validation precisely where data splitting is most harmful.

Limitations and open questions

The paper is candid about boundaries. The framework requires the USCQ and excludes singular models — overparameterized networks with positive-measure rank deficiency — whose stochastic complexity demands RLCT computation via algebraic geometry. Exactness comes at $$\int u(x)\,j(x)\,dL(x) = \int_{\mathbb{R}^k}\int_{m^{-1}(z)} u\,dH\,dz$$0 per-step cost, restricting applicability to moderate sample sizes; mixing-time constants in high dimensions remain unquantified beyond geometric ergodicity. The uniqueness of the scalar Jacobian factor at kinks depends on a restrictive singular-value condition not generally satisfied by pathwise AD. Inexact projections introduce solver-tolerance-proportional bias. Open questions left explicitly include bridging pseudo-inverse-based coarea formulas for rank-deficient stratified spaces with the RLCT, and developing scalable approximations validated against this exact baseline.

Conclusion

This work supplies the missing measure-theoretic foundation for NML computation on regular non-smooth estimators, proving well-posedness, minimax optimality, AD-consistency, and classical asymptotic penalties for PDL models, and delivering an exact geometric MCMC reference implementation. Its empirical demonstration that exact NML matches or beats cross-validation without data splitting gives the theory immediate practical relevance for high-dimensional sparse inference, while the clearly marked boundary at singular models defines the precise frontier for subsequent work.