Ridge, Lasso, and Elastic Net in C++

July 2026 · 56 min read

The three solvers so far in this series, OLS, QR/SVD, and gradient descent, all minimize the same unpenalized objective $\lVert y - X_b\beta \rVert_2^2$. They can disagree on how they get there, but when the design is full rank they agree on the answer. That answer has a weakness that has nothing to do with numerical stability: when features are collinear, redundant, or simply numerous relative to the number of rows, the unique minimizer of squared error can have unstable, high-variance coefficients even though the fitted line predicts perfectly well. Regularization changes the objective itself, adding a penalty on the size of $\beta$, trading a small amount of bias for an often much larger reduction in variance. This post implements the three standard penalized estimators in standard-library C++: ridge regression as a closed-form modification of the normal equations, and lasso and elastic net via cyclic coordinate descent with soft-thresholding. All three keep the series' fit / predict contract, handle the intercept correctly (never penalized), and add a train-only standardization layer, warm-started coefficient paths, and deterministic cross-validation for choosing the penalty strength.

Why Penalize at All

The QR/SVD post already showed the mechanism: when two columns of $X_b$ are highly correlated, the design has a near-null direction, and the least-squares solution can move a long way along that direction while barely changing the fitted values. Formally, the variance of the OLS estimator is $\operatorname{Var}(\hat\beta) = \sigma^2 (X_b^\top X_b)^{-1}$, and a near-singular $X_b^\top X_b$ makes some entries of that inverse enormous. The fitted line is still fine; the individual coefficients are not trustworthy, and a slightly different sample can produce wildly different $\hat\beta$. Regularization addresses this directly by shrinking $\beta$ toward zero, which shrinks $\operatorname{Var}(\hat\beta)$ at the cost of introducing a small, known bias. The three penalties in this post make different trade-offs along that bias-variance line:

• Ridge ($\ell_2$ penalty $\lambda\lVert\beta\rVert_2^2$) shrinks every coefficient continuously and never sets one to exactly zero. It is a closed-form, always-invertible modification of the normal equations.
• Lasso ($\ell_1$ penalty $\lambda\lVert\beta\rVert_1$) shrinks and can zero out coefficients entirely, performing feature selection as a side effect of fitting. It has no closed form and needs an iterative solver.
• Elastic net blends both penalties. It selects features like lasso but handles groups of correlated features more gracefully, an issue lasso has that this post demonstrates directly.

Ridge: A Small, Always-Invertible Modification

Ridge regression minimizes squared error plus an $\ell_2$ penalty on the coefficients (never the intercept):

Differentiating with respect to $\beta$ (holding the intercept out of the penalty) and setting the gradient to zero gives a normal-equation system that looks almost identical to Part 1's, with one added term:

where $\Lambda$ is a $(d{+}1)\times(d{+}1)$ diagonal matrix that is zero in the intercept's row/column and $\lambda$ everywhere else, so the intercept is estimated but never shrunk. Adding $\lambda I$ to the feature block of $X_b^\top X_b$ has an immediate practical benefit: $X_b^\top X_b$ is only positive semi-definite (it can be singular under exact collinearity, which is why solveLinearSystem can throw on duplicated columns), but for any $\lambda > 0$ the sum $X_b^\top X_b + \Lambda$ is strictly positive definite and therefore always invertible, regardless of how collinear the features are. Ridge does not merely tolerate collinearity; it is one of the standard fixes for it.

The scaling convention matters when $\lambda$ values are compared with libraries or formulas elsewhere. This post consistently uses an unnormalized half-sum-of-squares term, $\tfrac12\lVert y-X\beta\rVert^2$, so the ridge system contains $\lambda I$, the elastic-net threshold is $\lambda\alpha$, and its denominator contains $\lambda(1-\alpha)$. Under the equally common normalized objective $\tfrac1{2n}\lVert y-X\beta\rVert^2+\lambda_{\mathrm{norm}}[\alpha\lVert\beta\rVert_1+\tfrac{1-\alpha}{2}\lVert\beta\rVert_2^2]$, multiplying the stationarity equation by $n$ produces $X^\top X+n\lambda_{\mathrm{norm}}I$, threshold $n\lambda_{\mathrm{norm}}\alpha$, and denominator $\sum_i x_{ij}^2+n\lambda_{\mathrm{norm}}(1-\alpha)$. The two conventions describe identical fits when $\lambda_{\mathrm{code}}=n\lambda_{\mathrm{norm}}$; copying the same numeric $\lambda$ between them would change the strength of regularization by a factor of $n$.

Standardization Is Not Optional Here

The penalty $\lambda\lVert\beta\rVert_2^2$ treats every coefficient identically, but a coefficient's size is only meaningful relative to the scale of its feature. A feature measured in kilometers and the same feature measured in millimeters would need $\beta$ values that differ by a factor of $10^6$ to represent the same relationship, and a fixed $\lambda$ would penalize them completely differently. Every regularized estimator in this post therefore reuses the train-only standardization from Part 3: center and scale every feature column to zero mean and unit population variance, with the same zero-variance guard (a constant column gets scale $1$ and centers to exactly zero). The one difference from Part 3 is the target: ridge and lasso only center $y$, they do not rescale it, because there is no coefficient attached to $y$ that a penalty could distort.

Centering both $X$ and $y$ has a second, load-bearing consequence: it removes the intercept from the optimization entirely. For any fixed $\tilde\beta$, the optimal standardized intercept is $\overline{\tilde y - \tilde X\tilde\beta} = \overline{\tilde y} - \overline{\tilde X}\tilde\beta = 0 - 0 = 0$, exactly, because both $\tilde X$ and $\tilde y$ have mean zero by construction. So the standardized-space fit never carries an intercept term at all; it only ever solves for the $d$ feature weights $\tilde\beta$. This is true for ridge and for the penalized objectives later in this post alike, since centering is a property of the data, not of the penalty.

Converting back to raw units mirrors Part 3's derivation but is simpler, because $y$ was never rescaled:

One consequence of standardizing to unit population variance is worth stating explicitly, because the rest of this post leans on it: for every non-constant column, $\frac1n\sum_i \tilde x_{ij}^2 = \sigma_j^2/\sigma_j^2 = 1$, so

Every standardized column has exactly the same sum of squares, $n$, regardless of the original units or spread of that feature. This single fact is what makes the coordinate-descent update later in this post collapse to a clean, uniform formula instead of a per-feature-dependent one.

Ridge as a Spectral Filter

The SVD post wrote the unpenalized least-squares solution as $\beta = \sum_k \frac{u_k^\top y}{\sigma_k}v_k$, a sum over singular directions weighted by $1/\sigma_k$. Substituting the SVD $\tilde X = U S V^\top$ into the ridge normal equations $(\tilde X^\top\tilde X + \lambda I)\tilde\beta = \tilde X^\top \tilde y$ and simplifying (using $V^\top V = I$ and $S$ diagonal) gives the ridge solution in exactly the same singular-direction form, with one extra factor:

At $\lambda = 0$ this is precisely the unpenalized $\sum_k \frac{u_k^\top \tilde y}{\sigma_k}v_k$. As $\lambda$ grows, every direction's contribution is scaled by the shrinkage factor $\sigma_k/(\sigma_k^2+\lambda)$, and that factor is not uniform across $k$: for a large singular value $\sigma_k \gg \sqrt\lambda$, the factor is close to $1/\sigma_k$ (barely shrunk), while for a small singular value $\sigma_k \ll \sqrt\lambda$, the factor is close to $\sigma_k/\lambda$ (shrunk almost proportionally to $\sigma_k$ itself, i.e., far more aggressively). Ridge therefore does not shrink every direction equally: it preferentially damps the directions the data constrains least (small $\sigma_k$, the near-null, high-variance directions responsible for the instability described above) while leaving the well-determined directions comparatively untouched.

Relative to the unpenalized contribution in singular direction $k$, Ridge retains the fraction

This exposes the selective damping directly. At $\lambda=4$, a well-determined direction with $\sigma=5$ retains $25/29\approx0.8621$ of its OLS contribution, while a weak direction with $\sigma=0.5$ retains only $0.25/4.25\approx0.0588$. The weak direction loses almost all of its contribution while the strong direction remains mostly intact.

A concrete coefficient calculation must preserve the implementation's unpenalized intercept. Reuse the one-feature data $x=[1,2,4]$, $y=[1,3,7]^\top$, whose OLS line is $\hat y=-1+2x$. The implementation centers $y$ and standardizes $x$, leaving a single standardized feature direction with $\sigma_{\tilde X}=\sqrt3$. Its retained fraction is $3/(3+\lambda)$, the raw slope is $6/(3+\lambda)$, and the intercept is recovered from the training means rather than penalized:

The slope approaches zero as $\lambda$ grows, but the intercept approaches $\bar y=11/3$ rather than zero. Its changing numeric value is the centering correction required to keep the fitted line anchored at $(\bar x,\bar y)$ while the slope shrinks, not evidence that the intercept was penalized. With several correlated features, the same singular-direction filtering produces a smooth, generally curved Ridge trace because each direction retains a different fraction.

RidgeRegressor Implementation

The struct matches the shape used throughout the series: a fitted flag, nFeatures, intercept, and coef. Fitting standardizes, solves the small $(d\times d)$ ridge normal equations with the exact solveLinearSystem from Part 1, and converts back:

struct Standardizer { vector<double> mu, sigma; double yMean; };

Standardizer standardizeFit(const Matrix& X, const vector<double>& y,
                             Matrix& Xs, vector<double>& ys) {
    int n = (int)X.size(), d = (int)X[0].size();
    Standardizer s;
    s.mu.assign(d, 0.0); s.sigma.assign(d, 0.0);
    for (int i = 0; i < n; ++i) for (int j = 0; j < d; ++j) s.mu[j] += X[i][j];
    for (int j = 0; j < d; ++j) s.mu[j] /= n;

    vector<double> maxDev(d, 0.0);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < d; ++j) maxDev[j] = max(maxDev[j], fabs(X[i][j]-s.mu[j]));
    for (int j = 0; j < d; ++j) {
        if (maxDev[j] == 0.0) { s.sigma[j] = 1.0; continue; }
        double acc = 0.0;
        for (int i = 0; i < n; ++i) { double z = (X[i][j]-s.mu[j])/maxDev[j]; acc += z*z; }
        s.sigma[j] = maxDev[j] * sqrt(acc / n);
    }
    s.yMean = accumulate(y.begin(), y.end(), 0.0) / n;

    Xs.assign(n, vector<double>(d));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < d; ++j) Xs[i][j] = (X[i][j]-s.mu[j]) / s.sigma[j];
    ys.assign(n, 0.0);
    for (int i = 0; i < n; ++i) ys[i] = y[i] - s.yMean;
    return s;
}

void toRaw(const vector<double>& betaTilde, const Standardizer& s,
           double& intercept, vector<double>& coef) {
    int d = (int)betaTilde.size();
    coef.assign(d, 0.0);
    for (int j = 0; j < d; ++j) coef[j] = betaTilde[j] / s.sigma[j];
    intercept = s.yMean;
    for (int j = 0; j < d; ++j) intercept -= coef[j] * s.mu[j];
}

vector<double> ridgeFitStandardized(const Matrix& Xs, const vector<double>& ys, double lambda) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    Matrix A(d, vector<double>(d, 0.0));
    vector<double> b(d, 0.0);
    for (int i = 0; i < n; ++i)
        for (int a = 0; a < d; ++a) {
            b[a] += Xs[i][a] * ys[i];
            for (int c = 0; c < d; ++c) A[a][c] += Xs[i][a] * Xs[i][c];
        }
    for (int a = 0; a < d; ++a) A[a][a] += lambda;   // the entire ridge modification
    return solveLinearSystem(A, b);   // reuses Part 1's partial-pivoting solver
}

struct RidgeRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;

    void fit(const Matrix& X, const vector<double>& y, double lambda) {
        fitted = false;
        validateFitInput(X, y, "RidgeRegressor");
        if (!isfinite(lambda) || lambda < 0.0)
            throw invalid_argument("RidgeRegressor::fit: lambda must be finite and non-negative");
        Matrix Xs; vector<double> ys;
        Standardizer s = standardizeFit(X, y, Xs, ys);
        vector<double> betaTilde = ridgeFitStandardized(Xs, ys, lambda);
        toRaw(betaTilde, s, intercept, coef);
        nFeatures = (int)X[0].size();
        fitted = true;
    }
    // predict(x) / predict(X): identical shape to every other regressor in the series.
};

Ridge never allocates an $n\times(d{+}1)$ augmented matrix the way QR and SVD do: the feature block $\tilde X^\top\tilde X$ is only $d\times d$, and the intercept is recovered afterward from the means, exactly as it was in Part 1's unpenalized normal equations. The only new line, relative to a plain ridge-free normal-equation solve on standardized data, is A[a][a] += lambda.

Elastic Net: When the Penalty Needs an Iterative Solver

Lasso replaces the $\ell_2$ penalty with $\ell_1$, and elastic net keeps both, controlled by a mixing parameter $\alpha\in[0,1]$ (the code calls this l1Ratio, following the common convention that $\alpha=1$ is pure lasso and $\alpha=0$ is pure ridge):

The $\ell_1$ term is not differentiable at $\tilde\beta_j=0$, so there is no analog of the closed-form ridge solve; the standard fitting algorithm is cyclic coordinate descent: fix every coefficient except $\tilde\beta_j$, minimize the one-dimensional objective in $\tilde\beta_j$ exactly, then move to the next coordinate, repeating until nothing changes. Because $J$ is convex (strictly, whenever $\alpha<1$ or the design has full column rank), this converges to the global minimum regardless of the order coordinates are visited, though the speed of convergence does depend on that order and on how correlated the columns are, a point the coefficient path later in this post returns to.

Holding every coefficient but $\tilde\beta_j$ fixed, write $r_{ij} = \tilde y_i - \sum_{k\ne j}\tilde x_{ik}\tilde\beta_k$ for the partial residual with feature $j$'s contribution removed, and $z_j = \sum_i \tilde x_{ij}\,r_{ij}$ for that feature's correlation with the partial residual. The one-dimensional objective in $\tilde\beta_j$ is a quadratic plus an absolute value, and its minimizer is the classic soft-thresholding update:

$S$ shifts $z_j$ toward zero by $\gamma=\lambda\alpha$ and clips it at zero if it would cross, which is exactly what makes lasso ($\alpha>0$) capable of setting a coefficient to exactly zero, something the smooth ridge penalty can never do (ridge's shrinkage factor $\sigma_k/(\sigma_k^2+\lambda)$ is strictly positive for every finite $\lambda$). Here is where the standardization fact from earlier collapses the formula: because every standardized column satisfies $\sum_i\tilde x_{ij}^2=n$ exactly, the denominator is the same constant, $n+\lambda(1-\alpha)$, for every feature, and only the numerator $z_j$ varies across coordinates. Without standardization, the denominator would be a different, data-dependent number for every feature, and comparing penalty strength across features would not be meaningful.

Coordinate Descent by Hand

A concrete pass makes the update mechanical rather than abstract. Six rows carry two informative features and one weak, largely uninformative one: raw $x_1=[1,2,3,4,5,6]$, $x_2=[2,1,4,3,6,5]$, $x_3=[0.3,-0.4,0.1,-0.2,0.5,-0.3]$, and $y=[3.2,2.9,7.8,7.6,12.4,11.9]$ (so $x_1$ and $x_2$ are correlated with each other at $\rho=0.83$, and $x_3$ is nearly unrelated to $y$). Standardizing gives $\tilde x_1 = [-1.4639,-0.8783,-0.2928,0.2928,0.8783,1.4639]$, $\tilde x_2=[-0.8783,-1.4639,0.2928,-0.2928,1.4639,0.8783]$, $\tilde x_3=[0.9186,-1.2247,0.3062,-0.6124,1.5309,-0.9186]$, and $\tilde y = [-4.4333,-4.7333,0.1667,-0.0333,4.7667,4.2667]$; every column's sum of squares is exactly $6=n$, as the standardization fact guarantees. Starting from $\tilde\beta=(0,0,0)$ and using $\lambda=6,\ \alpha=1$ (pure lasso, so $\gamma=6$ and the denominator is just $n=6$ for every coordinate):

1. Coordinate 1. With the other two coefficients at zero, $z_1 = \tilde x_1^\top\tilde y = 21.0209$. Soft-thresholding: $S(21.0209,6)=15.0209$, so $\tilde\beta_1 \leftarrow 15.0209/6 = 2.5035$.
2. Coordinate 2. The partial residual now excludes $\tilde\beta_1\tilde x_1$: $z_2 = \tilde x_2^\top(\tilde y - \tilde\beta_1\tilde x_1) = 9.1605$. $S(9.1605,6)=3.1605$, so $\tilde\beta_2\leftarrow 3.1605/6=0.5268$.
3. Coordinate 3. $z_3 = \tilde x_3^\top(\tilde y - \tilde\beta_1\tilde x_1 - \tilde\beta_2\tilde x_2) = 5.1045$. Since $|z_3|=5.1045 < \gamma=6$, soft-thresholding clips it to exactly zero: $\tilde\beta_3\leftarrow 0.0000$. One sweep already finds the weak feature uninformative enough to remove.
4. Sweep 2, coordinate 1. $z_1 = 18.4021$ (smaller than before, because the fit already explains part of $y$); $\tilde\beta_1\leftarrow (18.4021-6)/6=2.0670$.
5. Sweep 2, coordinate 2. $z_2=11.3304$; $\tilde\beta_2\leftarrow(11.3304-6)/6=0.8884$.
6. Sweep 2, coordinate 3. $z_3=3.8972$, still below $\gamma=6$; $\tilde\beta_3$ stays at exactly $0.0000$.

After these two sweeps, $\tilde\beta=(2.0670,\ 0.8884,\ 0.0000)$. Continuing the same cycle (verified by script, not shown step by step) converges after 62 sweeps to $\tilde\beta\approx(1.1111,\ 1.6804,\ 0.0000)$: coordinates 1 and 2 keep trading weight back and forth, slowly, because they are correlated ($\rho=0.83$) and the objective is nearly flat along the direction that keeps $\tilde\beta_1+\tilde\beta_2$ constant, while coordinate 3 never leaves zero once it arrives there, since its $z_3$ only shrinks further as the fit improves. The animation below replays exactly these six update steps.

ElasticNetRegressor Implementation

Rather than writing a separate Lasso struct, the implementation follows the same convention real libraries use: one ElasticNetRegressor whose alpha (l1Ratio) parameter recovers lasso at $\alpha=1$. A dedicated closed-form RidgeRegressor still exists for the reasons above (it is cheaper and exact), but as a consistency check, $\alpha=0$ in the coordinate-descent code below must converge to the same answer as the closed-form ridge solve, since both are minimizing the identical convex quadratic; the complete program later confirms this bit-for-bit.

double softThreshold(double z, double gamma) {
    if (z > gamma) return z - gamma;
    if (z < -gamma) return z + gamma;
    return 0.0;
}

// Coordinate descent on already-standardized data. `beta` is both the
// warm-start input and the output (modified in place).
int coordinateDescentFit(const Matrix& Xs, const vector<double>& ys, double lambda, double alpha,
                          vector<double>& beta, int maxIter, double tol,
                          bool* didConverge = nullptr) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    if (didConverge) *didConverge = false;
    vector<double> colSq(d, 0.0);
    for (int j = 0; j < d; ++j) {
        double s = 0.0; for (int i = 0; i < n; ++i) s += Xs[i][j]*Xs[i][j]; colSq[j] = s;
    }  // == n for every non-constant column, by construction

    vector<double> resid(n);
    for (int i = 0; i < n; ++i) {
        double pred = 0.0;
        for (int j = 0; j < d; ++j) pred += Xs[i][j] * beta[j];
        resid[i] = ys[i] - pred;
    }

    int iter = 0;
    for (; iter < maxIter; ++iter) {
        double maxDelta = 0.0;
        for (int j = 0; j < d; ++j) {
            if (colSq[j] == 0.0) continue;  // constant column: coefficient stays at 0
            double zj = 0.0;
            for (int i = 0; i < n; ++i) zj += Xs[i][j] * (resid[i] + Xs[i][j] * beta[j]);
            double newBeta = softThreshold(zj, lambda * alpha) / (colSq[j] + lambda * (1.0 - alpha));
            double delta = newBeta - beta[j];
            if (delta != 0.0) for (int i = 0; i < n; ++i) resid[i] -= Xs[i][j] * delta;
            maxDelta = max(maxDelta, fabs(delta));
            beta[j] = newBeta;
        }
        if (maxDelta < tol) {
            ++iter;
            if (didConverge) *didConverge = true;
            break;
        }
    }
    return iter;   // reported as iterationsRun; no throw if maxIter is hit
}

struct ElasticNetRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    int iterationsRun = 0;
    bool converged = false;

    void fit(const Matrix& X, const vector<double>& y, double lambda, double alpha,
             int maxIter = 20000, double tol = 1e-9) {
        fitted = false;
        validateFitInput(X, y, "ElasticNetRegressor");
        if (!isfinite(lambda) || lambda < 0.0)
            throw invalid_argument("ElasticNetRegressor::fit: lambda must be finite and non-negative");
        if (!isfinite(alpha) || alpha < 0.0 || alpha > 1.0)
            throw invalid_argument("ElasticNetRegressor::fit: alpha (l1Ratio) must be in [0, 1]");

        Matrix Xs; vector<double> ys;
        Standardizer s = standardizeFit(X, y, Xs, ys);
        vector<double> betaTilde(X[0].size(), 0.0);
        bool didConverge = false;
        iterationsRun = coordinateDescentFit(
            Xs, ys, lambda, alpha, betaTilde, maxIter, tol, &didConverge);
        converged = didConverge;
        toRaw(betaTilde, s, intercept, coef);
        nFeatures = (int)X[0].size();
        fitted = true;
    }
    // predict(x) / predict(X): identical shape to every other regressor in the series.
};

Two cheap checks confirm the coordinate-descent solver reaches the true minimum rather than merely stopping. First, the elastic-net objective $J(\tilde\beta)=\tfrac12\lVert\tilde y-\tilde X\tilde\beta\rVert_2^2+\lambda\left[\tfrac{1-\alpha}{2}\lVert\tilde\beta\rVert_2^2+\alpha\lVert\tilde\beta\rVert_1\right]$ is non-increasing across sweeps; the complete program prints it after the first three sweeps of the $\lambda=5$ lasso fit ($J_1=69.5524$, $J_2=69.5376$, $J_3=69.5375$) down to the converged $J_{\text{final}}=69.5294$. Second, at a minimum every coordinate satisfies the KKT (subgradient stationarity) condition: $\tilde x_j^\top(\tilde y-\tilde X\tilde\beta)=\lambda\alpha\operatorname{sign}(\tilde\beta_j)+\lambda(1-\alpha)\tilde\beta_j$ for an active coefficient, and $\lvert\tilde x_j^\top(\tilde y-\tilde X\tilde\beta)\rvert\le\lambda\alpha$ for one held at zero. The largest violation of that condition at the CV-selected $\lambda=5$ solution is $1.97\times10^{-7}$, certifying a coordinate-wise optimum to solver precision.

Coefficient Paths and Warm Starts

A single $\lambda$ is rarely the right unit of analysis; what a reader (or a cross-validation loop) usually wants is the whole path of solutions across a range of $\lambda$. Fitting each $\lambda$ from a cold start ($\tilde\beta=0$) works but wastes effort: adjacent $\lambda$ values on a reasonably fine grid have similar solutions, so initializing coordinate descent from the previous $\lambda$'s answer (a warm start) needs far fewer sweeps to reconverge. The natural traversal order is descending $\lambda$: start at the most-penalized (simplest, fastest-converging) end of the path and relax the penalty gradually.

vector<vector<double>> elasticNetPath(const Matrix& X, const vector<double>& y,
                                       const vector<double>& lambdasDesc, double alpha,
                                       vector<int>* itersOut = nullptr) {
    Matrix Xs; vector<double> ys;
    Standardizer s = standardizeFit(X, y, Xs, ys);
    vector<double> beta(X[0].size(), 0.0);      // carried across the whole loop: the warm start
    vector<vector<double>> path;
    for (double lambda : lambdasDesc) {
        int iters = coordinateDescentFit(Xs, ys, lambda, alpha, beta, 20000, 1e-9);
        if (itersOut) itersOut->push_back(iters);
        double intercept; vector<double> coef;
        toRaw(beta, s, intercept, coef);
        vector<double> row{intercept};
        for (double c : coef) row.push_back(c);
        path.push_back(row);
    }
    return path;
}

On the six-feature dataset used in Scenario B below, fitting a ten-point lasso path with a fresh cold start at every $\lambda$ costs $20{,}132$ total coordinate sweeps; warm-starting the identical path costs $16{,}372$, about $19\%$ fewer. That saving is smaller than the roughly eight-fold reduction a less collinear toy example shows, and the reason is instructive: the slow-converging direction in this dataset is the near-duplicate pair of features (introduced in Scenario B), and that direction is slow to resolve at every $\lambda$, warm-started or not, because $\lambda$ mostly regularizes the other four coordinates rather than breaking the near-degeneracy between the duplicated pair. Warm-starting still helps, just less dramatically when the bottleneck is collinearity rather than the distance between successive $\lambda$ values.

The path has a natural upper endpoint as well. Starting from $\tilde\beta=0$, the soft-thresholding update leaves coordinate $j$ at zero for as long as $\lvert\tilde x_j^\top\tilde y\rvert\le\lambda\alpha$, so every coefficient is exactly zero once $\lambda\ge\lambda_{\max}=\max_j\lvert\tilde x_j^\top\tilde y\rvert$ (taking $\alpha=1$). Any larger $\lambda$ returns the same empty model, which is why a principled lasso grid begins at $\lambda_{\max}$ and descends from there. The complete program computes $\lambda_{\max}\approx675.67$ for Scenario B and confirms the fit at $\lambda_{\max}$ keeps zero coefficients nonzero; the cross-validation grid below tops out at $\lambda=50$, well inside the path, so even its most-penalized point already holds three features active.

Deterministic K-Fold Cross-Validation

Choosing $\lambda$ from the training data alone (e.g., minimizing training error) always prefers $\lambda\to 0$, since less regularization can only fit the training set better. $\lambda$ has to be chosen by estimated out-of-sample error, and $k$-fold cross-validation is the standard way to estimate that without setting aside a separate validation set permanently. The fold assignment is seeded exactly like the shuffles elsewhere in this series, so the same seed always reproduces the same folds:

vector<vector<int>> kFoldIndices(int n, int k, unsigned seed) {
    vector<int> order(n);
    iota(order.begin(), order.end(), 0);
    mt19937 rng(seed);
    shuffle(order.begin(), order.end(), rng);
    vector<vector<int>> folds(k);
    for (int i = 0; i < n; ++i) folds[i % k].push_back(order[i]);
    return folds;
}

For each candidate $\lambda$, every fold in turn is held out as validation while the model is refit (from scratch, standardized on the training folds only, exactly as a real train/validation split would be) on the remaining $k-1$ folds. The reported curve pools squared errors across all held-out predictions before dividing by the total validation-row count, so every row receives equal weight even when fold sizes differ by one. With equal-sized folds this is exactly the arithmetic mean of the fold MSE values. The full CV loop for ridge (elastic net's is identical except for the model being refit) is:

vector<double> cvCurveRidge(const Matrix& X, const vector<double>& y,
                            const vector<double>& lambdas, int k, unsigned seed) {
    vector<vector<int>> folds = kFoldIndices((int)X.size(), k, seed);
    vector<double> curve(lambdas.size(), 0.0);
    for (size_t li = 0; li < lambdas.size(); ++li) {
        double total = 0.0; int count = 0;
        for (int f = 0; f < k; ++f) {
            vector<int> trainIdx;
            for (int g = 0; g < k; ++g) if (g != f) for (int i : folds[g]) trainIdx.push_back(i);
            RidgeRegressor model;
            model.fit(rowsOf(X, trainIdx), valsOf(y, trainIdx), lambdas[li]);
            auto pred = model.predict(rowsOf(X, folds[f]));
            auto yva = valsOf(y, folds[f]);
            for (size_t i = 0; i < yva.size(); ++i) { double e = yva[i]-pred[i]; total += e*e; }
            count += (int)yva.size();
        }
        curve[li] = total / count;
    }
    return curve;
}

The single CV minimum is not the only defensible choice of $\lambda$. Because each fold contributes its own held-out error, the cross-validation estimate carries a standard error across folds, and the one-standard-error rule selects the largest (most regularized, simplest) $\lambda$ whose CV error still lies within one SE of the best. On Scenario B's lasso curve the minimum is $0.39667$ at $\lambda=5$ with a fold-level standard error of $0.04041$, and every $\lambda$ up to $20$ falls inside that one-SE band, so the rule prefers $\lambda=20$: a much sparser model at negligible cost in estimated error. The complete program reports both the raw minimum and this one-SE choice.

Validation and Errors

• Row count mismatch, ragged rows, empty input, or non-finite entries throw std::invalid_argument, exactly as in every earlier post.
• RidgeRegressor::fit rejects a negative or non-finite lambda. It never throws for a collinear design, because $\lambda>0$ makes the modified normal equations positive definite regardless.
• ElasticNetRegressor::fit rejects a negative or non-finite lambda and rejects an alpha outside $[0,1]$. Non-convergence within maxIter is not an exception, following the same convention as gradient descent's epochsRun: iterationsRun and converged are exposed so the caller can decide whether to raise maxIter or accept the current iterate.
• A non-finite coefficient after conversion back to raw units (which would indicate an internal solver failure, not a user error) throws std::runtime_error.
• Calling predict before fit throws std::logic_error; a feature-count mismatch or non-finite feature throws std::invalid_argument.

Complete Program and Output

The program below runs two scenarios. Scenario A refits the exact shared two-feature dataset from Part 1 (seed 7, $n=240$, no collinearity) with a vanishingly small $\lambda$, confirming that both ridge and lasso reduce to the known OLS answer when the penalty is switched off. Scenario B builds a six-feature dataset ($n=200$) with $x_1,x_2$ as true signal (coefficients $2.0,-1.5$), $x_3$ a near-duplicate of $x_1$ (correlation $0.9965$), and $x_4,x_5,x_6$ pure noise (true coefficient $0$), then compares OLS, cross-validated ridge, cross-validated lasso, and cross-validated elastic net, prints both a ridge and a lasso coefficient path, and closes with the $\alpha=0$-equals-ridge cross-check.

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;
using Matrix = vector<vector<double>>;

// ============================================================
// Shared validation and metrics (same conventions as Parts 1-3).
// ============================================================
static void validateFitInput(const Matrix& X, const vector<double>& y, const string& who) {
    if (X.size() != y.size())
        throw invalid_argument(who + "::fit: X has " + to_string(X.size()) +
                               " rows but y has " + to_string(y.size()) + " entries");
    int n = (int)X.size();
    if (n == 0) throw invalid_argument(who + "::fit: empty training set");
    int d = (int)X[0].size();
    if (d == 0) throw invalid_argument(who + "::fit: zero feature columns");
    for (int i = 0; i < n; ++i) {
        if ((int)X[i].size() != d)
            throw invalid_argument(who + "::fit: row " + to_string(i) +
                                   " has " + to_string(X[i].size()) +
                                   " columns, expected " + to_string(d));
        for (int j = 0; j < d; ++j)
            if (!isfinite(X[i][j]))
                throw invalid_argument(who + "::fit: non-finite feature in row " + to_string(i));
        if (!isfinite(y[i]))
            throw invalid_argument(who + "::fit: non-finite target at row " + to_string(i));
    }
    if (n < d + 1)
        throw invalid_argument(who + "::fit: need at least " + to_string(d + 1) +
                               " rows for " + to_string(d) + " features plus intercept, got " +
                               to_string(n));
}

static void validateMetricInputs(const vector<double>& yTrue, const vector<double>& yPred, const string& name) {
    if (yTrue.size() != yPred.size()) throw invalid_argument(name + ": size mismatch");
    if (yTrue.empty()) throw invalid_argument(name + ": vectors must be non-empty");
    for (size_t i = 0; i < yTrue.size(); ++i)
        if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
            throw invalid_argument(name + ": non-finite value at index " + to_string(i));
}

static double mse(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "mse");
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) { double e = yTrue[i] - yPred[i]; s += e * e; }
    return s / yTrue.size();
}
static double rmse(const vector<double>& yTrue, const vector<double>& yPred) { return sqrt(mse(yTrue, yPred)); }
static double r2(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "r2");
    double meanY = accumulate(yTrue.begin(), yTrue.end(), 0.0) / yTrue.size();
    double ssRes = 0.0, ssTot = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        ssRes += (yTrue[i]-yPred[i])*(yTrue[i]-yPred[i]);
        ssTot += (yTrue[i]-meanY)*(yTrue[i]-meanY);
    }
    if (ssTot == 0.0) throw invalid_argument("r2: target values are constant, so R^2 is undefined");
    return 1.0 - ssRes / ssTot;
}

// ============================================================
// Partial-pivoting Gaussian solver (identical to Part 1).
// ============================================================
static vector<double> solveLinearSystem(Matrix A, vector<double> b) {
    int n = (int)A.size();
    if (n == 0) throw invalid_argument("solveLinearSystem: A must be non-empty");
    if (b.size() != A.size()) throw invalid_argument("solveLinearSystem: A and b dimensions do not match");

    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("solveLinearSystem: A must be square; row " +
                                   to_string(i) + " has " + to_string(A[i].size()) +
                                   " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v)) throw invalid_argument("solveLinearSystem: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
        if (!isfinite(b[i])) throw invalid_argument("solveLinearSystem: b contains a non-finite value");
    }
    const double tol = numeric_limits<double>::epsilon() * max(1.0, matrixScale) * n * 100.0;

    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i) if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        if (best <= tol)
            throw runtime_error("solveLinearSystem: matrix is singular or near-singular "
                                 "(no usable pivot in column " + to_string(k) + "); "
                                 "check for duplicated/collinear features");
        if (piv != k) { swap(A[k], A[piv]); swap(b[k], b[piv]); }
        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            if (factor == 0.0) continue;
            for (int j = k; j < n; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }
    vector<double> beta(n, 0.0);
    for (int k = n - 1; k >= 0; --k) {
        double s = b[k];
        for (int j = k + 1; j < n; ++j) s -= A[k][j] * beta[j];
        beta[k] = s / A[k][k];
    }
    return beta;
}

// ============================================================
// Direct OLS via the exact Part-1 normal equations (used as the
// unpenalized baseline; independent of the penalized solvers below).
// ============================================================
struct OLSFit { double intercept; vector<double> coef; };

static OLSFit fitOLSDirect(const Matrix& X, const vector<double>& y) {
    validateFitInput(X, y, "fitOLSDirect");
    int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;
    Matrix A(p, vector<double>(p, 0.0));
    vector<double> b(p, 0.0);
    for (int i = 0; i < n; ++i) {
        vector<double> row(p);
        row[0] = 1.0;
        for (int j = 0; j < d; ++j) row[j+1] = X[i][j];
        for (int a = 0; a < p; ++a) {
            b[a] += row[a] * y[i];
            for (int c = 0; c < p; ++c) A[a][c] += row[a] * row[c];
        }
    }
    vector<double> beta = solveLinearSystem(A, b);
    OLSFit out;
    out.intercept = beta[0];
    out.coef.assign(beta.begin()+1, beta.end());
    return out;
}

// ============================================================
// Train-only standardization: center y, standardize X columns to
// zero mean / unit population variance (constant columns -> scale 1).
// Because every standardized column has population variance 1,
// sum_i xTilde_ij^2 == n exactly for every non-constant column.
// ============================================================
struct Standardizer { vector<double> mu, sigma; double yMean; };

static Standardizer standardizeFit(const Matrix& X, const vector<double>& y, Matrix& Xs, vector<double>& ys) {
    int n = (int)X.size(), d = (int)X[0].size();
    Standardizer s;
    s.mu.assign(d, 0.0);
    s.sigma.assign(d, 0.0);
    for (int i = 0; i < n; ++i) for (int j = 0; j < d; ++j) s.mu[j] += X[i][j];
    for (int j = 0; j < d; ++j) s.mu[j] /= n;

    vector<double> maxDev(d, 0.0);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < d; ++j) maxDev[j] = max(maxDev[j], fabs(X[i][j] - s.mu[j]));
    for (int j = 0; j < d; ++j) {
        if (maxDev[j] == 0.0) { s.sigma[j] = 1.0; continue; }
        double acc = 0.0;
        for (int i = 0; i < n; ++i) { double z = (X[i][j] - s.mu[j]) / maxDev[j]; acc += z * z; }
        s.sigma[j] = maxDev[j] * sqrt(acc / n);
    }
    s.yMean = accumulate(y.begin(), y.end(), 0.0) / n;

    Xs.assign(n, vector<double>(d));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < d; ++j) Xs[i][j] = (X[i][j] - s.mu[j]) / s.sigma[j];
    ys.assign(n, 0.0);
    for (int i = 0; i < n; ++i) ys[i] = y[i] - s.yMean;
    return s;
}

static void toRaw(const vector<double>& betaTilde, const Standardizer& s,
                   double& intercept, vector<double>& coef) {
    int d = (int)betaTilde.size();
    coef.assign(d, 0.0);
    for (int j = 0; j < d; ++j) coef[j] = betaTilde[j] / s.sigma[j];
    intercept = s.yMean;
    for (int j = 0; j < d; ++j) intercept -= coef[j] * s.mu[j];
}

// ============================================================
// Ridge: closed form on standardized data via the Part-1 solver.
// (X~^T X~ + lambda*I) beta~ = X~^T y~   (intercept excluded: centering
// already makes the standardized-space intercept exactly zero).
// ============================================================
static vector<double> ridgeFitStandardized(const Matrix& Xs, const vector<double>& ys, double lambda) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    Matrix A(d, vector<double>(d, 0.0));
    vector<double> b(d, 0.0);
    for (int i = 0; i < n; ++i)
        for (int a = 0; a < d; ++a) {
            b[a] += Xs[i][a] * ys[i];
            for (int c = 0; c < d; ++c) A[a][c] += Xs[i][a] * Xs[i][c];
        }
    for (int a = 0; a < d; ++a) A[a][a] += lambda;
    return solveLinearSystem(A, b);
}

struct RidgeRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;

    void fit(const Matrix& X, const vector<double>& y, double lambda) {
        fitted = false;
        validateFitInput(X, y, "RidgeRegressor");
        if (!isfinite(lambda) || lambda < 0.0)
            throw invalid_argument("RidgeRegressor::fit: lambda must be finite and non-negative");
        Matrix Xs; vector<double> ys;
        Standardizer s = standardizeFit(X, y, Xs, ys);
        vector<double> betaTilde = ridgeFitStandardized(Xs, ys, lambda);
        toRaw(betaTilde, s, intercept, coef);
        for (double c : coef) if (!isfinite(c))
            throw runtime_error("RidgeRegressor::fit: conversion produced a non-finite coefficient");
        nFeatures = (int)X[0].size();
        fitted = true;
    }
    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("RidgeRegressor::predict: model is not fitted");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("RidgeRegressor::predict: expected " + to_string(nFeatures) +
                                   " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) {
            if (!isfinite(x[j])) throw invalid_argument("RidgeRegressor::predict: non-finite feature at index " + to_string(j));
            out += coef[j] * x[j];
        }
        if (!isfinite(out)) throw runtime_error("RidgeRegressor::predict: prediction became non-finite");
        return out;
    }
    vector<double> predict(const Matrix& X) const {
        vector<double> out; out.reserve(X.size());
        for (const auto& row : X) out.push_back(predict(row));
        return out;
    }
};

static vector<vector<double>> ridgePath(const Matrix& X, const vector<double>& y, const vector<double>& lambdas) {
    Matrix Xs; vector<double> ys;
    Standardizer s = standardizeFit(X, y, Xs, ys);
    vector<vector<double>> path;
    for (double lambda : lambdas) {
        vector<double> betaTilde = ridgeFitStandardized(Xs, ys, lambda);
        double intercept; vector<double> coef;
        toRaw(betaTilde, s, intercept, coef);
        vector<double> row; row.push_back(intercept);
        for (double c : coef) row.push_back(c);
        path.push_back(row);
    }
    return path;
}

// ============================================================
// Elastic Net (Lasso is the alpha=1 special case): cyclic coordinate
// descent with soft-thresholding on standardized data.
//   beta~_j <- S(z_j, lambda*alpha) / (colSq_j + lambda*(1-alpha))
//   z_j = x~_j . (y~ - X~ beta~ + x~_j * beta~_j)     [partial residual]
// colSq_j == n exactly for every standardized column.
// ============================================================
static inline double softThreshold(double z, double gamma) {
    if (z > gamma) return z - gamma;
    if (z < -gamma) return z + gamma;
    return 0.0;
}

// Coordinate descent on already-standardized data. beta is both the
// warm-start input and the result (modified in place).
static int coordinateDescentFit(const Matrix& Xs, const vector<double>& ys, double lambda, double alpha,
                                 vector<double>& beta, int maxIter, double tol,
                                 bool* didConverge = nullptr) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    if (didConverge) *didConverge = false;
    vector<double> colSq(d, 0.0);
    for (int j = 0; j < d; ++j) { double s = 0.0; for (int i = 0; i < n; ++i) s += Xs[i][j]*Xs[i][j]; colSq[j] = s; }

    vector<double> resid(n);
    for (int i = 0; i < n; ++i) {
        double pred = 0.0;
        for (int j = 0; j < d; ++j) pred += Xs[i][j] * beta[j];
        resid[i] = ys[i] - pred;
    }

    int iter = 0;
    for (; iter < maxIter; ++iter) {
        double maxDelta = 0.0;
        for (int j = 0; j < d; ++j) {
            if (colSq[j] == 0.0) continue;  // constant column, coefficient stays at its warm-start value (0)
            double zj = 0.0;
            for (int i = 0; i < n; ++i) zj += Xs[i][j] * (resid[i] + Xs[i][j] * beta[j]);
            double newBeta = softThreshold(zj, lambda * alpha) / (colSq[j] + lambda * (1.0 - alpha));
            if (!isfinite(newBeta))
                throw runtime_error("coordinateDescentFit: coefficient became non-finite; check lambda/alpha");
            double delta = newBeta - beta[j];
            if (delta != 0.0) for (int i = 0; i < n; ++i) resid[i] -= Xs[i][j] * delta;
            maxDelta = max(maxDelta, fabs(delta));
            beta[j] = newBeta;
        }
        if (maxDelta < tol) {
            ++iter;
            if (didConverge) *didConverge = true;
            break;
        }
    }
    return iter;
}

struct ElasticNetRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    int iterationsRun = 0;
    bool converged = false;

    void fit(const Matrix& X, const vector<double>& y, double lambda, double alpha,
             int maxIter = 20000, double tol = 1e-9) {
        fitted = false;
        validateFitInput(X, y, "ElasticNetRegressor");
        if (!isfinite(lambda) || lambda < 0.0)
            throw invalid_argument("ElasticNetRegressor::fit: lambda must be finite and non-negative");
        if (!isfinite(alpha) || alpha < 0.0 || alpha > 1.0)
            throw invalid_argument("ElasticNetRegressor::fit: alpha (l1Ratio) must be in [0, 1]");
        if (maxIter <= 0) throw invalid_argument("ElasticNetRegressor::fit: maxIter must be positive");
        if (!isfinite(tol) || tol <= 0.0) throw invalid_argument("ElasticNetRegressor::fit: tol must be finite and positive");

        Matrix Xs; vector<double> ys;
        Standardizer s = standardizeFit(X, y, Xs, ys);
        vector<double> betaTilde(X[0].size(), 0.0);
        bool didConverge = false;
        int iters = coordinateDescentFit(
            Xs, ys, lambda, alpha, betaTilde, maxIter, tol, &didConverge);
        iterationsRun = iters;
        converged = didConverge;
        toRaw(betaTilde, s, intercept, coef);
        for (double c : coef) if (!isfinite(c))
            throw runtime_error("ElasticNetRegressor::fit: conversion produced a non-finite coefficient");
        nFeatures = (int)X[0].size();
        fitted = true;
    }
    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("ElasticNetRegressor::predict: model is not fitted");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("ElasticNetRegressor::predict: expected " + to_string(nFeatures) +
                                   " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) {
            if (!isfinite(x[j])) throw invalid_argument("ElasticNetRegressor::predict: non-finite feature at index " + to_string(j));
            out += coef[j] * x[j];
        }
        if (!isfinite(out)) throw runtime_error("ElasticNetRegressor::predict: prediction became non-finite");
        return out;
    }
    vector<double> predict(const Matrix& X) const {
        vector<double> out; out.reserve(X.size());
        for (const auto& row : X) out.push_back(predict(row));
        return out;
    }
};

// Pathwise fit across a descending sequence of lambdas, warm-starting
// each fit from the previous (larger-lambda) solution.
static vector<vector<double>> elasticNetPath(const Matrix& X, const vector<double>& y,
                                              const vector<double>& lambdasDesc, double alpha,
                                              vector<int>* itersOut = nullptr,
                                              int maxIter = 20000, double tol = 1e-9) {
    Matrix Xs; vector<double> ys;
    Standardizer s = standardizeFit(X, y, Xs, ys);
    int d = (int)X[0].size();
    vector<double> beta(d, 0.0);
    vector<vector<double>> path;
    for (double lambda : lambdasDesc) {
        int iters = coordinateDescentFit(Xs, ys, lambda, alpha, beta, maxIter, tol);
        if (itersOut) itersOut->push_back(iters);
        double intercept; vector<double> coef;
        toRaw(beta, s, intercept, coef);
        vector<double> row; row.push_back(intercept);
        for (double c : coef) row.push_back(c);
        path.push_back(row);
    }
    return path;
}

// ============================================================
// Elastic-net solver diagnostics: objective value, lambda_max, and the
// KKT (subgradient stationarity) violation of a solution. All operate on
// already-standardized data, matching the coordinate-descent solver above.
// ============================================================
static double enetObjective(const Matrix& Xs, const vector<double>& ys,
                            const vector<double>& beta, double lambda, double alpha) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    double sse = 0.0;
    for (int i = 0; i < n; ++i) {
        double pred = 0.0;
        for (int j = 0; j < d; ++j) pred += Xs[i][j] * beta[j];
        double e = ys[i] - pred; sse += e * e;
    }
    double l2 = 0.0, l1 = 0.0;
    for (int j = 0; j < d; ++j) { l2 += beta[j] * beta[j]; l1 += fabs(beta[j]); }
    return 0.5 * sse + lambda * (0.5 * (1.0 - alpha) * l2 + alpha * l1);
}

// Smallest lambda that leaves every coefficient at zero (pure-lasso scale):
// from beta=0 the update keeps coord j at zero iff |x~_j . y~| <= lambda*alpha,
// so lambda_max = max_j |x~_j . y~| (taking alpha = 1).
static double lambdaMaxLasso(const Matrix& Xs, const vector<double>& ys) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    double m = 0.0;
    for (int j = 0; j < d; ++j) {
        double z = 0.0;
        for (int i = 0; i < n; ++i) z += Xs[i][j] * ys[i];
        m = max(m, fabs(z));
    }
    return m;
}

// Maximum KKT stationarity violation of an elastic-net solution. For each j let
// g_j = x~_j . (y~ - X~ beta). Optimality requires g_j = lambda*alpha*sign(beta_j)
// + lambda*(1-alpha)*beta_j when beta_j != 0, and |g_j| <= lambda*alpha otherwise.
// A value near zero certifies a coordinate-wise optimum.
static double enetKktViolation(const Matrix& Xs, const vector<double>& ys,
                               const vector<double>& beta, double lambda, double alpha) {
    int n = (int)Xs.size(), d = (int)Xs[0].size();
    double worst = 0.0;
    for (int j = 0; j < d; ++j) {
        double g = 0.0;
        for (int i = 0; i < n; ++i) {
            double pred = 0.0;
            for (int k = 0; k < d; ++k) pred += Xs[i][k] * beta[k];
            g += Xs[i][j] * (ys[i] - pred);
        }
        double v;
        if (beta[j] != 0.0) {
            double target = lambda * alpha * (beta[j] > 0.0 ? 1.0 : -1.0) + lambda * (1.0 - alpha) * beta[j];
            v = fabs(g - target);
        } else {
            v = max(0.0, fabs(g) - lambda * alpha);
        }
        worst = max(worst, v);
    }
    return worst;
}

// ============================================================
// Deterministic k-fold cross-validation.
// ============================================================
static vector<vector<int>> kFoldIndices(int n, int k, unsigned seed) {
    vector<int> order(n);
    iota(order.begin(), order.end(), 0);
    mt19937 rng(seed);
    shuffle(order.begin(), order.end(), rng);
    vector<vector<int>> folds(k);
    for (int i = 0; i < n; ++i) folds[i % k].push_back(order[i]);
    return folds;
}

static Matrix rowsOf(const Matrix& X, const vector<int>& idx) { Matrix out; out.reserve(idx.size()); for (int i : idx) out.push_back(X[i]); return out; }
static vector<double> valsOf(const vector<double>& y, const vector<int>& idx) { vector<double> out; out.reserve(idx.size()); for (int i : idx) out.push_back(y[i]); return out; }

static vector<double> cvCurveRidge(const Matrix& X, const vector<double>& y, const vector<double>& lambdas, int k, unsigned seed) {
    int n = (int)X.size();
    vector<vector<int>> folds = kFoldIndices(n, k, seed);
    vector<double> curve(lambdas.size(), 0.0);
    for (size_t li = 0; li < lambdas.size(); ++li) {
        double total = 0.0; int totalCount = 0;
        for (int f = 0; f < k; ++f) {
            vector<int> valIdx = folds[f], trainIdx;
            for (int g = 0; g < k; ++g) if (g != f) for (int i : folds[g]) trainIdx.push_back(i);
            RidgeRegressor model;
            model.fit(rowsOf(X, trainIdx), valsOf(y, trainIdx), lambdas[li]);
            Matrix Xva = rowsOf(X, valIdx); vector<double> yva = valsOf(y, valIdx);
            auto pred = model.predict(Xva);
            for (size_t i = 0; i < yva.size(); ++i) { double e = yva[i]-pred[i]; total += e*e; }
            totalCount += (int)yva.size();
        }
        curve[li] = total / totalCount;
    }
    return curve;
}

static vector<double> cvCurveEnet(const Matrix& X, const vector<double>& y, const vector<double>& lambdas, double alpha, int k, unsigned seed) {
    int n = (int)X.size();
    vector<vector<int>> folds = kFoldIndices(n, k, seed);
    vector<double> curve(lambdas.size(), 0.0);
    for (size_t li = 0; li < lambdas.size(); ++li) {
        double total = 0.0; int totalCount = 0;
        for (int f = 0; f < k; ++f) {
            vector<int> valIdx = folds[f], trainIdx;
            for (int g = 0; g < k; ++g) if (g != f) for (int i : folds[g]) trainIdx.push_back(i);
            ElasticNetRegressor model;
            model.fit(rowsOf(X, trainIdx), valsOf(y, trainIdx), lambdas[li], alpha);
            Matrix Xva = rowsOf(X, valIdx); vector<double> yva = valsOf(y, valIdx);
            auto pred = model.predict(Xva);
            for (size_t i = 0; i < yva.size(); ++i) { double e = yva[i]-pred[i]; total += e*e; }
            totalCount += (int)yva.size();
        }
        curve[li] = total / totalCount;
    }
    return curve;
}

// Per-fold validation MSE for elastic net (same folds and seed as cvCurveEnet),
// used only for the one-standard-error rule; the pooled cvCurveEnet curve above is
// left untouched so its printed numbers do not change.
static vector<vector<double>> cvFoldMseEnet(const Matrix& X, const vector<double>& y, const vector<double>& lambdas, double alpha, int k, unsigned seed) {
    int n = (int)X.size();
    vector<vector<int>> folds = kFoldIndices(n, k, seed);
    vector<vector<double>> foldMse(lambdas.size());
    for (size_t li = 0; li < lambdas.size(); ++li) {
        for (int f = 0; f < k; ++f) {
            vector<int> valIdx = folds[f], trainIdx;
            for (int g = 0; g < k; ++g) if (g != f) for (int i : folds[g]) trainIdx.push_back(i);
            ElasticNetRegressor model;
            model.fit(rowsOf(X, trainIdx), valsOf(y, trainIdx), lambdas[li], alpha);
            Matrix Xva = rowsOf(X, valIdx); vector<double> yva = valsOf(y, valIdx);
            auto pred = model.predict(Xva);
            double sderr = 0.0; for (size_t i = 0; i < yva.size(); ++i) { double e = yva[i]-pred[i]; sderr += e*e; }
            foldMse[li].push_back(sderr / yva.size());
        }
    }
    return foldMse;
}

static int argminIdx(const vector<double>& v) { return (int)(min_element(v.begin(), v.end()) - v.begin()); }

static void printVec(const char* label, const vector<double>& v) {
    printf("%s=[", label);
    for (size_t i = 0; i < v.size(); ++i) printf("%s%.4f", i ? ", " : "", v[i]);
    printf("]");
}

// ============================================================
// main
// ============================================================
int main() {
    // === Scenario A: sanity check against Part 1's OLS on the series' shared dataset ===
    printf("=== Scenario A: sanity check on the shared 2-feature dataset (seed 7) ===\n");
    {
        mt19937 rng(7);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        normal_distribution<double> noise(0.0, 0.6);
        const int N = 240;
        Matrix X; vector<double> y;
        for (int i = 0; i < N; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            X.push_back({x1, x2});
            y.push_back(3.0 + 2.0*x1 - 1.5*x2 + noise(rng));
        }
        vector<int> order(N); iota(order.begin(), order.end(), 0);
        shuffle(order.begin(), order.end(), rng);
        Matrix Xtr, Xte; vector<double> ytr, yte;
        int nTrain = (int)(0.8 * N);
        for (int k = 0; k < N; ++k) {
            int i = order[k];
            if (k < nTrain) { Xtr.push_back(X[i]); ytr.push_back(y[i]); }
            else            { Xte.push_back(X[i]); yte.push_back(y[i]); }
        }
        OLSFit ols = fitOLSDirect(Xtr, ytr);
        printf("[OLS]      intercept=%.4f  coef=[%.4f, %.4f]\n", ols.intercept, ols.coef[0], ols.coef[1]);

        RidgeRegressor ridge; ridge.fit(Xtr, ytr, 1e-6);
        printf("[Ridge]    intercept=%.4f  coef=[%.4f, %.4f]  (lambda=1e-6)\n", ridge.intercept, ridge.coef[0], ridge.coef[1]);

        ElasticNetRegressor lasso; lasso.fit(Xtr, ytr, 1e-6, 1.0);
        printf("[Lasso]    intercept=%.4f  coef=[%.4f, %.4f]  (lambda=1e-6, iters=%d)\n",
               lasso.intercept, lasso.coef[0], lasso.coef[1], lasso.iterationsRun);

        auto predRidge = ridge.predict(Xte);
        auto predLasso = lasso.predict(Xte);
        printf("[Ridge]    test RMSE=%.4f  test R^2=%.4f\n", rmse(yte, predRidge), r2(yte, predRidge));
        printf("[Lasso]    test RMSE=%.4f  test R^2=%.4f\n", rmse(yte, predLasso), r2(yte, predLasso));
    }

    // === Scenario B: collinear duplicate + irrelevant noise features ===
    printf("\n=== Scenario B: collinear duplicate + 3 irrelevant features (n=200, d=6) ===\n");
    Matrix Xb; vector<double> yb;
    {
        mt19937 rng(2026);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        normal_distribution<double> dupNoise(0.0, 0.15);
        normal_distribution<double> irrelevant(0.0, 1.0);
        normal_distribution<double> yNoise(0.0, 0.6);
        const int N = 200;
        for (int i = 0; i < N; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            double x3 = x1 + dupNoise(rng);
            double x4 = irrelevant(rng), x5 = irrelevant(rng), x6 = irrelevant(rng);
            Xb.push_back({x1, x2, x3, x4, x5, x6});
            yb.push_back(3.0 + 2.0*x1 - 1.5*x2 + yNoise(rng));
        }
    }
    printf("corr(x1,x3) = %.6f\n", [&]{
        double m1=0,m3=0; int n=(int)Xb.size();
        for (auto&r:Xb){m1+=r[0];m3+=r[2];} m1/=n; m3/=n;
        double cov=0,v1=0,v3=0;
        for (auto&r:Xb){cov+=(r[0]-m1)*(r[2]-m3); v1+=(r[0]-m1)*(r[0]-m1); v3+=(r[2]-m3)*(r[2]-m3);}
        return cov/sqrt(v1*v3);
    }());

    OLSFit olsB = fitOLSDirect(Xb, yb);
    printf("[OLS]        intercept=%.4f  ", olsB.intercept); printVec("coef", olsB.coef); printf("\n");

    vector<double> ridgeLambdas = {0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100, 300};
    vector<double> enetLambdas  = {0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50};

    auto ridgeCv = cvCurveRidge(Xb, yb, ridgeLambdas, 5, 123);
    auto lassoCv = cvCurveEnet(Xb, yb, enetLambdas, 1.0, 5, 123);
    auto enetCv  = cvCurveEnet(Xb, yb, enetLambdas, 0.5, 5, 123);

    printf("\n-- 5-fold CV mean-squared error (Ridge) --\n");
    for (size_t i = 0; i < ridgeLambdas.size(); ++i) printf("  lambda=%-8g cvMSE=%.5f\n", ridgeLambdas[i], ridgeCv[i]);
    int ridgeBest = argminIdx(ridgeCv);
    printf("-- 5-fold CV mean-squared error (Lasso, alpha=1.0) --\n");
    for (size_t i = 0; i < enetLambdas.size(); ++i) printf("  lambda=%-8g cvMSE=%.5f\n", enetLambdas[i], lassoCv[i]);
    int lassoBest = argminIdx(lassoCv);
    printf("-- 5-fold CV mean-squared error (Elastic Net, alpha=0.5) --\n");
    for (size_t i = 0; i < enetLambdas.size(); ++i) printf("  lambda=%-8g cvMSE=%.5f\n", enetLambdas[i], enetCv[i]);
    int enetBest = argminIdx(enetCv);

    printf("\nCV-selected lambda: Ridge=%g  Lasso=%g  ElasticNet(0.5)=%g\n",
           ridgeLambdas[ridgeBest], enetLambdas[lassoBest], enetLambdas[enetBest]);

    RidgeRegressor ridgeB; ridgeB.fit(Xb, yb, ridgeLambdas[ridgeBest]);
    ElasticNetRegressor lassoB; lassoB.fit(Xb, yb, enetLambdas[lassoBest], 1.0);
    ElasticNetRegressor enetB; enetB.fit(Xb, yb, enetLambdas[enetBest], 0.5);

    printf("\n[Ridge, CV]      intercept=%.4f  ", ridgeB.intercept); printVec("coef", ridgeB.coef); printf("\n");
    printf("[Lasso, CV]      intercept=%.4f  ", lassoB.intercept); printVec("coef", lassoB.coef);
    printf("  (iters=%d, converged=%s)\n", lassoB.iterationsRun, lassoB.converged ? "true" : "false");
    printf("[ElasticNet, CV] intercept=%.4f  ", enetB.intercept); printVec("coef", enetB.coef);
    printf("  (iters=%d, converged=%s)\n", enetB.iterationsRun, enetB.converged ? "true" : "false");

    auto predOls  = [&]{ vector<double> p; for (auto&r:Xb) { double v=olsB.intercept; for(int j=0;j<6;++j) v+=olsB.coef[j]*r[j]; p.push_back(v); } return p; }();
    printf("\n[OLS]            train RMSE=%.4f  train R^2=%.4f\n", rmse(yb, predOls), r2(yb, predOls));
    printf("[Ridge, CV]      train RMSE=%.4f  train R^2=%.4f\n", rmse(yb, ridgeB.predict(Xb)), r2(yb, ridgeB.predict(Xb)));
    printf("[Lasso, CV]      train RMSE=%.4f  train R^2=%.4f\n", rmse(yb, lassoB.predict(Xb)), r2(yb, lassoB.predict(Xb)));
    printf("[ElasticNet, CV] train RMSE=%.4f  train R^2=%.4f\n", rmse(yb, enetB.predict(Xb)), r2(yb, enetB.predict(Xb)));

    // Held-out clean test set from the same true model (independent RNG, so the
    // n=200 training data above is unchanged); measures recovery, not training fit.
    printf("\n-- Held-out test metrics on an independent clean test set (n=100) --\n");
    Matrix Xtest; vector<double> ytest;
    {
        mt19937 rng(9001);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        normal_distribution<double> dupNoise(0.0, 0.15);
        normal_distribution<double> irrelevant(0.0, 1.0);
        normal_distribution<double> yNoise(0.0, 0.6);
        for (int i = 0; i < 100; ++i) {
            double x1 = ux(rng), x2 = ux(rng);
            double x3 = x1 + dupNoise(rng);
            double x4 = irrelevant(rng), x5 = irrelevant(rng), x6 = irrelevant(rng);
            Xtest.push_back({x1, x2, x3, x4, x5, x6});
            ytest.push_back(3.0 + 2.0*x1 - 1.5*x2 + yNoise(rng));
        }
    }
    auto predOlsTe = [&]{ vector<double> p; for (auto&r:Xtest){ double v=olsB.intercept; for(int j=0;j<6;++j) v+=olsB.coef[j]*r[j]; p.push_back(v);} return p; }();
    printf("[OLS]            test RMSE=%.4f  test R^2=%.4f\n", rmse(ytest, predOlsTe), r2(ytest, predOlsTe));
    printf("[Ridge, CV]      test RMSE=%.4f  test R^2=%.4f\n", rmse(ytest, ridgeB.predict(Xtest)), r2(ytest, ridgeB.predict(Xtest)));
    printf("[Lasso, CV]      test RMSE=%.4f  test R^2=%.4f\n", rmse(ytest, lassoB.predict(Xtest)), r2(ytest, lassoB.predict(Xtest)));
    printf("[ElasticNet, CV] test RMSE=%.4f  test R^2=%.4f\n", rmse(ytest, enetB.predict(Xtest)), r2(ytest, enetB.predict(Xtest)));

    // Coefficient paths: ridge (shrink-and-group) vs lasso (shrink-to-zero).
    printf("\n-- Ridge coefficient path (shrink toward zero, grouping on x1/x3) --\n");
    auto rPath = ridgePath(Xb, yb, ridgeLambdas);
    for (size_t i = 0; i < ridgeLambdas.size(); ++i) {
        printf("  lambda=%-8g intercept=%.4f  coef=[%.4f, %.4f, %.4f, %.4f, %.4f, %.4f]\n",
               ridgeLambdas[i], rPath[i][0], rPath[i][1], rPath[i][2], rPath[i][3], rPath[i][4], rPath[i][5], rPath[i][6]);
    }

    printf("\n-- Lasso coefficient path, warm-started high-to-low lambda (shrink to exact zero) --\n");
    vector<double> lassoDesc = enetLambdas; sort(lassoDesc.rbegin(), lassoDesc.rend());
    vector<int> pathIters;
    auto lPath = elasticNetPath(Xb, yb, lassoDesc, 1.0, &pathIters);
    for (size_t i = 0; i < lassoDesc.size(); ++i) {
        printf("  lambda=%-8g intercept=%.4f  coef=[%.4f, %.4f, %.4f, %.4f, %.4f, %.4f]  iters=%d\n",
               lassoDesc[i], lPath[i][0], lPath[i][1], lPath[i][2], lPath[i][3], lPath[i][4], lPath[i][5], lPath[i][6], pathIters[i]);
    }
    int totalWarm = 0; for (int it : pathIters) totalWarm += it;
    vector<int> coldIters;
    for (double lam : lassoDesc) {
        Matrix Xs; vector<double> ys; Standardizer s = standardizeFit(Xb, yb, Xs, ys);
        vector<double> beta(6, 0.0);
        coldIters.push_back(coordinateDescentFit(Xs, ys, lam, 1.0, beta, 20000, 1e-9));
    }
    int totalCold = 0; for (int it : coldIters) totalCold += it;
    printf("Total iterations across the path: warm-started=%d, cold-started=%d\n", totalWarm, totalCold);

    // Solver diagnostics on the lasso path: lambda_max, objective history, KKT, one-SE rule.
    printf("\n-- Solver diagnostics (lasso, alpha=1) --\n");
    {
        Matrix Xs; vector<double> ys;
        standardizeFit(Xb, yb, Xs, ys);

        double lmax = lambdaMaxLasso(Xs, ys);
        vector<double> bzero(6, 0.0);
        coordinateDescentFit(Xs, ys, lmax, 1.0, bzero, 20000, 1e-9);
        int nzAtMax = 0; for (double c : bzero) if (c != 0.0) ++nzAtMax;
        printf("lambda_max = max_j |x~_j . y~| = %.4f; the fit at lambda_max keeps %d of 6 coefficients nonzero\n", lmax, nzAtMax);

        printf("objective J(lambda=5) by sweep: ");
        for (int m = 1; m <= 3; ++m) {
            vector<double> bm(6, 0.0);
            coordinateDescentFit(Xs, ys, 5.0, 1.0, bm, m, 1e-9);
            printf("J%d=%.4f  ", m, enetObjective(Xs, ys, bm, 5.0, 1.0));
        }
        vector<double> bfin(6, 0.0);
        coordinateDescentFit(Xs, ys, 5.0, 1.0, bfin, 20000, 1e-9);
        printf("J_final=%.4f (monotone non-increasing)\n", enetObjective(Xs, ys, bfin, 5.0, 1.0));
        printf("KKT max violation at CV-selected lambda=5: %.2e (certifies a coordinate-wise optimum)\n",
               enetKktViolation(Xs, ys, bfin, 5.0, 1.0));
    }
    {
        vector<vector<double>> foldMse = cvFoldMseEnet(Xb, yb, enetLambdas, 1.0, 5, 123);
        double mean = 0.0; for (double v : foldMse[lassoBest]) mean += v; mean /= (double)foldMse[lassoBest].size();
        double var = 0.0; for (double v : foldMse[lassoBest]) var += (v-mean)*(v-mean); var /= (double)(foldMse[lassoBest].size()-1);
        double se = sqrt(var) / sqrt((double)foldMse[lassoBest].size());
        double thresh = lassoCv[lassoBest] + se;
        double lam1se = enetLambdas[lassoBest];
        for (size_t i = 0; i < enetLambdas.size(); ++i) if (lassoCv[i] <= thresh && enetLambdas[i] > lam1se) lam1se = enetLambdas[i];
        printf("one-SE rule: min cvMSE %.5f at lambda=%g, fold SE %.5f; largest lambda within one SE = %g\n",
               lassoCv[lassoBest], enetLambdas[lassoBest], se, lam1se);
    }

    // Cross-check: ElasticNet(alpha=0) must agree with the closed-form Ridge at the same lambda.
    printf("\n-- Cross-check: ElasticNet(alpha=0) vs closed-form Ridge, lambda=10 --\n");
    ElasticNetRegressor enetAsRidge; enetAsRidge.fit(Xb, yb, 10.0, 0.0, 20000, 1e-12);
    RidgeRegressor ridgeDirect; ridgeDirect.fit(Xb, yb, 10.0);
    printf("[ElasticNet a=0] intercept=%.6f  ", enetAsRidge.intercept); printVec("coef", enetAsRidge.coef);
    printf("  (iters=%d)\n", enetAsRidge.iterationsRun);
    printf("[Ridge closed]   intercept=%.6f  ", ridgeDirect.intercept); printVec("coef", ridgeDirect.coef); printf("\n");

    // Error-handling demonstrations.
    printf("\n-- Error handling --\n");
    try { RidgeRegressor bad; bad.fit(Xb, yb, -1.0); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }
    try { ElasticNetRegressor bad; bad.fit(Xb, yb, 1.0, 1.5); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }
    try { RidgeRegressor m; m.fit(Xb, yb, 1.0); m.predict(vector<double>{1.0, 2.0}); }
    catch (const invalid_argument& e) { printf("caught: %s\n", e.what()); }

    return 0;
}

Compiled with g++ -O2 -std=c++17 -Wall -Wextra -pedantic, this prints:

Reading the Two Scenarios

Scenario A confirms correctness before anything interesting happens: with $\lambda$ effectively switched off, ridge and lasso reproduce Part 1's OLS coefficients ($2.9639$, $1.9777$, $-1.5173$) to four decimal places, and lasso does it in only $6$ coordinate sweeps because there is no collinearity to slow it down. This is the same kind of sanity check the QR/SVD post ran against OLS: different algorithm, same answer, when the underlying objective is the same.

Scenario B is where the three penalties diverge. OLS splits the near-duplicate pair unevenly and unstably: $1.8355$ on $x_1$ and $0.1811$ on $x_3$, summing to about $2.02$, close to the combined true effect of $2.0$, but the individual split is an artifact of exactly which noise realization landed in the training data (the same identifiability argument the SVD post made for exact rank deficiency applies here in a softer form for near rank deficiency). Cross-validated ridge ($\lambda=0.3$) gives $1.5747$ and $0.4392$, only slightly more even than OLS at this small a $\lambda$; the coefficient path table above shows the grouping effect much more clearly at larger $\lambda$: by $\lambda=10$, the split is $1.0289$ and $0.9285$, nearly equal, and it stays nearly equal all the way to $\lambda=300$. This is ridge's signature behavior under collinearity: correlated features are shrunk together, toward a shared value, rather than arbitrarily apart. Cross-validated lasso ($\lambda=5$) keeps the noise features at exactly $0.0000$ and puts most of the weight on $x_1$ ($1.8292$) while leaving $x_3$ small but nonzero ($0.1713$); the lasso path table shows this pattern holds across the whole range of $\lambda$, lasso does not group the duplicated pair evenly the way ridge does, it just picks a dominant one. Cross-validated elastic net ($\lambda=0.5,\ \alpha=0.5$) lands in between, $1.6078$ and $0.4056$, closer to ridge's even split than to lasso's lopsided one, while still shrinking the noise coefficients ($0.0044,\ 0.0018,\ -0.0152$) more than ridge does at a comparable $\lambda$.

Worth noting honestly: the $\lambda$ that cross-validation selects for prediction accuracy is not the same $\lambda$ that produces the cleanest, most interpretable coefficient pattern. The CV curves are fairly flat across a wide range of $\lambda$ for this dataset (ridge's cvMSE only moves from $0.4052$ to $0.4043$ between $\lambda=0.01$ and its $\lambda=0.3$ optimum), because the noise level ($\sigma=0.6$) is large enough that heavy shrinkage costs almost nothing in validation error, but also does not buy much. A slightly larger $\lambda$ than the CV optimum, chosen deliberately from the path table rather than from the CV curve, sacrifices negligible test accuracy for a substantially sparser or more evenly grouped model, a legitimate reason to look at the whole path rather than only the single CV-selected point. The closing cross-check, ElasticNetRegressor with $\alpha=0$ against the closed-form RidgeRegressor at the same $\lambda=10$, matches to six decimal places ($2.951789$ intercept, identical coefficients down to the last printed digit), confirming the iterative solver and the closed form are minimizing the same objective.

The held-out numbers keep that comparison honest. Evaluated on an independent $100$-row clean test set drawn from the same model, cross-validated lasso posts the best test RMSE ($0.6178$), just ahead of OLS ($0.6197$), with ridge ($0.6268$) and elastic net ($0.6254$) a little behind at their CV-selected $\lambda$. The gaps are small because the noise floor ($\sigma=0.6$) dominates, but lasso's edge is exactly the payoff of setting the three irrelevant coefficients to zero instead of fitting them to noise. The solver diagnostics printed alongside, $\lambda_{\max}\approx675.67$, a monotone objective history, a KKT violation of $1.97\times10^{-7}$, and a one-SE $\lambda$ of $20$, together certify that each reported fit is a real optimum of its stated objective and quantify how much simpler a model the one-standard-error rule would accept.

When to Use Each

Complexity and Memory

Ridge's cost is fixed per fit regardless of how the data behaves. Coordinate descent's cost is data-dependent: well-separated, weakly-correlated features converge in a handful of sweeps (Scenario A's $6$), while strongly collinear ones can need thousands (Scenario B's lasso path, $1300$-$2000$ sweeps per $\lambda$), because the per-sweep contraction rate for cyclic coordinate descent degrades as the condition number of the correlated block grows. Warm-starting the path amortizes this cost across $\lambda$ values but does not eliminate the underlying slow direction.

Practical Diagnostics

• Comparison against the closed form. Fitting ElasticNetRegressor with $\alpha=0$ and comparing it to RidgeRegressor at the same $\lambda$, as the complete program does, is a strong regression test: any mismatch beyond floating-point tolerance means the coordinate-descent implementation has a bug, not that the two methods legitimately disagree.
• The CV curve's shape, not just its minimum. A curve that is flat across a wide range (Scenario B's ridge and lasso curves both are) means the exact $\lambda$ matters little for prediction, and a larger, more interpretable $\lambda$ from that flat region is a reasonable, low-cost choice.
• iterationsRun tracked against maxIter. A path where several $\lambda$ values report converged=false is a signal to either raise maxIter or inspect the design for the kind of near-collinearity that produces a slow-converging coordinate pair, exactly as Scenario B's $x_1$/$x_3$ pair does.
• The coefficient path, read for grouping versus selection. Two correlated features that move together as $\lambda$ increases indicate ridge- or elastic-net-like grouping; one moving to zero while the other stays nonzero indicates lasso-like arbitrary selection. Either can be the "correct" answer depending on whether the two features are substantively redundant or only accidentally correlated in this sample.

Limitations and Edge Cases

• Lasso's arbitrary tie-breaking. When two columns are (near-)identical, lasso's choice of which one to keep is sensitive to floating-point order of operations and the exact noise realization; it is not a meaningful signal about which feature is "truly" more important. Elastic net's added $\ell_2$ term is the standard fix when this ambiguity is undesirable.
• Coordinate descent's convergence rate depends on collinearity, not just on $\lambda$. As shown above, warm-starting does not fix a fundamentally slow-converging correlated pair; only a larger $\lambda$ (more $\ell_2$ curvature) or a different algorithm (a coordinate system aligned with the correlated pair, or an accelerated solver) resolves that.
• Cross-validation variance. With only $5$ folds and $n=200$, the CV curve itself is a noisy estimate of true out-of-sample error; the CV-selected $\lambda$ can shift meaningfully under a different fold seed, especially in a flat region of the curve like the one this dataset produces.
• None of these three penalties address outliers. Ridge, lasso, and elastic net all keep the squared-error term as the fit-quality measure; a single large outlier still pulls the fit exactly as hard as it would under plain OLS, regardless of how much $\ell_1$ or $\ell_2$ penalty is layered on top. Penalizing the coefficients and robustifying the loss are two independent design choices.

Next Step

Every solver in this series so far, penalized or not, has minimized some form of squared error, which is the right loss when noise is roughly symmetric and light-tailed but is exactly the wrong loss when a handful of extreme observations should not be allowed to dominate the fit. The next post reuses the weighted-least-squares idea underneath this post's ridge normal equations, a fixed diagonal penalty $\Lambda$ folded into $X_b^\top X_b$ turns out to generalize naturally to a diagonal reweighting $X_b^\top W X_b$ that changes every iteration, but replaces the loss itself: Huber loss for outlier-resistant fitting, and the pinball loss for quantile regression, both fit by iteratively reweighted least squares rather than coordinate descent.