ElasticNet

Coordinate-descent elastic net regression

from _api_doc_utils import *

1 Where it fits

Group: Regression

ElasticNet estimates a penalized linear model with a convex combination of L1 and L2 penalties. The current delegated coordinate-descent implementation centers the outcome, but not the feature columns. With \(\bar y=n^{-1}\sum_i y_i\) and \(y_c=y-\bar y\mathbf1\), it solves

\[ \begin{aligned} \min_{\beta}\quad &\frac{1}{2n}\|y_c-X\beta\|_2^2 \\ &+\lambda\rho\|\beta\|_1 +\frac{\lambda(1-\rho)}{2}\|\beta\|_2^2. \end{aligned} \]

and predicts \(\hat y=\bar y\mathbf1+X\hat\beta\). Thus the reported intercept is \(\bar y\), not \(\bar y-\bar x'\hat\beta\). This equals the usual jointly optimized unpenalized-intercept formulation when feature columns are centered; otherwise it is a different criterion. This behavior is retained for compatibility and should be accounted for when preparing features.

2 Criterion and solver

Here \(\lambda\) is penalty and \(\rho\) is l1_ratio. Coordinate updates use soft thresholding. The wrapper independently recomputes the final duality gap rather than treating the existence of returned coefficients as convergence. Let

\[ \begin{aligned} r&=y_c-X\hat\beta, \\ a&=X'r-n\lambda(1-\rho)\hat\beta, \\ L&=n\lambda\rho. \end{aligned} \]

and let \(c=\min\{1,L/\|a\|_\infty\}\), with \(c=1\) when \(\|a\|_\infty\leq L\). The reported gap is

\[ \begin{aligned} G={}&\frac12(1+c^2)\|r\|_2^2 +L\|\hat\beta\|_1 \\ &-c\,r'y_c +\frac{n\lambda(1-\rho)}{2}(1+c^2)\|\hat\beta\|_2^2. \end{aligned} \]

The fit is accepted only when \(G\) is finite and \(G\leq\tau\|y_c\|_2^2\), where \(\tau\) is tolerance. Otherwise fit() raises ValueError and clears any previous fitted state. A successful summary exposes duality_gap, duality_gap_tolerance, converged, iterations, termination_reason, and the final primal objective. Inputs and all hyperparameters are validated before fitting.

The class does not standardize \(X\). Because both the L1 and L2 penalties act on raw coefficients, users must scale features explicitly when a common penalty across columns is intended.

3 Implementation walkthrough and delegation boundary

Coordinate descent itself is delegated to linfa-elasticnet; the surrounding fit contract and convergence audit are package-owned.

  1. The wrapper clears prior state, validates the dense arrays and every hyperparameter, clones the arrays into a Linfa dataset, and configures Linfa with the requested penalty, mixing ratio, intercept flag, tolerance, and iteration cap. It performs no feature centering or scaling before delegation.
  2. Linfa runs coordinate descent and returns its hyperplane, intercept, and number of sweeps. The package does not expose or modify individual coordinate updates, active-set rules, or Linfa’s internal stopping decision.
  3. Crabbymetrics independently reconstructs \(y_c\) when an intercept is requested, computes the residual and elastic-net dual certificate shown above, and compares the gap with tolerance * dot(y_centered, y_centered). This second check, not merely Linfa returning a model, decides whether the public fit succeeds.
  4. It separately recomputes the public primal objective from Linfa predictions and coefficients. A finite accepted gap produces package-standard FitDiagnostics; an excessive or nonfinite gap is labeled as iteration-budget failure and the wrapper raises without installing the returned model.
  5. Prediction delegates to the accepted Linfa object. summary() extracts its intercept and hyperplane but combines them with package-owned gap and termination fields. The pairs bootstrap repeats this entire delegated-fit-plus-native-audit pipeline and aborts at the first failed replicate.

This boundary is why the page documents the observed intercept convention rather than assuming the textbook jointly centered problem. The wrapper can validate the returned solution against the criterion it exposes, but it cannot change the coordinate-descent mechanics without replacing the delegated solver.

4 Inference

The summary deliberately returns no analytic covariance or standard errors. L1 selection makes naive inverse-Hessian inference inappropriate, and the implementation does not provide debiasing, selective inference, or cross-validated penalty selection. The pairs bootstrap refits the same fixed \((\lambda,\rho)\) in every resample and returns raw intercept and coefficient draws. Those draws can describe algorithmic and sampling stability, but the class does not turn them into confidence intervals and does not account for tuning uncertainty.

5 Performance and numerical behavior

One complete coordinate sweep is \(O(np)\) for dense input, so runtime is approximately \(O(Inp)\) for \(I\) iterations. The solver uses \(O(n+p)\) working memory in addition to the dense \(O(np)\) input retained by the estimator. The final gap check costs another \(O(np)\). Sparse coefficients do not reduce the stored dense design. Highly correlated or poorly scaled columns can slow coordinate descent and make the selected support unstable. At \(\rho=0\) this solver is a coordinate-descent ridge fit with a different penalty normalization from the package’s dedicated Ridge class; the same numeric penalty therefore does not imply the same objective.

6 Python API

Constructor: cm.ElasticNet

Use ElasticNet(penalty, l1_ratio, tolerance, max_iterations), then fit(x, y), predict(x), summary(), and optionally bootstrap(B, seed=None). The summary reports point estimates, convergence diagnostics, and the final duality gap, and marks analytic inference unavailable. Bootstrap coefficient draws are stability diagnostics, not automatic confidence intervals; any nonconverged replicate aborts the bootstrap.

print(inspect.signature(cm.ElasticNet))
(penalty=1.0, l1_ratio=0.5, tolerance=0.0001, max_iterations=1000)
cls = cm.ElasticNet
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
bootstrap(self, /, n_bootstrap, seed=None)
fit(self, /, x, y)
predict(self, /, x)
summary(self, /)

7 Minimal example

rng = np.random.default_rng(4)
x = rng.normal(size=(180, 8))
y = 0.4 + x[:, :3] @ np.array([1.0, -0.8, 0.5]) + rng.normal(scale=0.5, size=180)
model = cm.ElasticNet(penalty=0.05, l1_ratio=0.7)
model.fit(x, y)
fit = model.summary()
print({key: fit[key] for key in ['converged', 'iterations', 'duality_gap', 'duality_gap_tolerance']})
print(fit['coef'])
print(model.predict(x[:3]))
{'converged': True, 'iterations': 5, 'duality_gap': 0.000140643302657395, 'duality_gap_tolerance': 0.04063251586480301}
[ 0.90319681 -0.79760879  0.45662966 -0.01442286 -0.01857245 -0.
  0.         -0.        ]
[ 0.6312702  -1.2665473  -2.43483936]

8 summary() contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

rng = np.random.default_rng(104)
x = rng.normal(size=(90, 5))
y = 0.4 + x[:, :2] @ np.array([1, -0.8]) + rng.normal(size=90) * 0.3
model = cm.ElasticNet(penalty=0.05, l1_ratio=0.7)
model.fit(x, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))
summary() key shape
intercept ()
coef (5,)
penalty ()
l1_ratio ()
duality_gap ()
duality_gap_tolerance ()
converged ()
iterations ()
termination_reason ()
objective ()
inference_available ()
intercept_se ()
coef_se ()