ElasticNet

Coordinate-descent elastic net regression

1 Where it fits

Group: Regression

ElasticNet estimates a penalized linear model with a convex combination of L1 and L2 penalties. With an intercept it centers both predictors and outcome: \(X_c=X-\mathbf1\bar x'\) and \(y_c=y-\bar y\mathbf1\). It solves

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

The reconstructed intercept is \(\hat b_0=\bar y-\bar x'\hat\beta\), and predictions are \(\hat b_0+X\hat\beta\). This corrects the old outcome-only centering behavior: shifted-design predictions can change relative to earlier releases. The public constructor always fits an intercept; there is no fit_intercept keyword.

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_c\hat\beta, \\ a&=X_c'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, centers feature columns when fitting an intercept, and delegates the mixed/L1 problem to Linfa. Linfa centers the outcome. No feature standardization is performed.
  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. The primal objective and gap use the same centered problem. The pure-L2 or zero-penalty endpoint instead uses augmented QR, with slope penalty \(n\lambda\), and reports a zero duality gap for the direct solve.
  5. The wrapper owns the accepted coefficients and reconstructed original-scale intercept. Prediction validates feature width and finiteness before multiplying; empty batches are allowed. The pairs bootstrap repeats the corrected fitting pipeline and aborts at the first failed replicate.

Centering changes the wrapper’s statistical contract without replacing Linfa’s coordinate-descent mechanics. Prediction does not retain a borrowed view of the training data.

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. Centering and the final gap check add dense \(O(np)\) work. Sparse coefficients do not reduce stored design size. Highly correlated or poorly scaled columns can slow coordinate descent. The pure-L2 endpoint uses QR; its equivalent dedicated Ridge penalty is \(n\lambda\), not \(\lambda\).

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.

Code 
print(inspect.signature(cm.ElasticNet))
(penalty=1.0, l1_ratio=0.5, tolerance=0.0001, max_iterations=1000)
Code 
cls = cm.ElasticNet
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
ElasticNet(penalty=1.0, l1_ratio=0.5, tolerance=0.0001, max_iterations=1000)
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.00017123025401133418, 'duality_gap_tolerance': 0.04063251586480301}
[ 0.91385898 -0.8087009   0.45861765 -0.03020589 -0.0299712  -0.
  0.00393059  0.        ]
[ 0.78632815 -1.16934817 -2.29436094]

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.

Code 
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 ()