Cross-Validation vs Bayesian Regularization for Ridge Regression

We look at an example of time series data $y_t$. We fit the model:

where $\epsilon_t \overset{\text{i.i.d}}{\sim} N(0, \sigma^2)$, and

Estimation of $\beta_0, \dots, \beta_{n-1}$ will be done via ridge regression:

for a fixed $\lambda > 0$. The choice of $\lambda$ is done either via 5-fold Cross-Valication (for a particular selection of folds), or via the Bayesian method discussed in Lectures 12-13. The Bayesian method also yields an alternative estimate of $\beta_0, \dots, \beta_{n-1}$ (along with uncertainty quantification).

We work with the FRED data on weekly initial unemployment claims (from https://fred.stlouisfed.org/series/ICSA).

We use all available data before March 2020 so that the exceptional surge in unemployment claims during COVID-19 does not dominate the analysis.

We will apply our model to this data, and obtain estimates of the underlying trend. Note that, in the context of our model, trend means the fitted function:

Recall the following function for fitting ridge regressions.

Given $\lambda$, the ridge regression estimate is computed as follows.

It is clear that when $\lambda = \lambda_{\min} := 10^{-2}$, the fitted trend function nearly reproduces the data, while when $\lambda = \lambda_{\max} := 10^{13}$, the fitted trend function fits a line to the data. Below we see which value of $\lambda$ in the range $[\lambda_{\min}, \lambda_{\max}]$ will be picked by the various methods.

First, we use the Cross-Validation method that we used in lectures.

The $\lambda$ value selected by this method is on the higher side. So the fitted trend function will be very smooth.

Obviously this is a poor estimate of the underlying trend. This CV method uses consecutive data points as folds. The five folds used by CV can be visualized as follows.

If we fix a fold (for testing), and fit the model on the other 4 folds, the fitted function can be visualized (for a fixed $\lambda$) as follows.

Suppose the test fold is the 3rd fold. It should be clear from the above that if $\lambda$ is chosen small, then the fitted function looks quite off compared to the actual data for fold = 3. This is the reason why a large $\lambda$ is being chosen by CV.

Now we change the CV method by taking the five folds randomly. This was basically the CV method used in the last lab. We can fit this CV also with the same function ‘ridge_cv’ by permuting the rows of $X$ and $y$ before calling ‘ridge_cv’. Indeed, the original function makes consecutive folds in the row order it receives. Supplying a random permutation of the paired rows therefore creates random folds in calendar time. We do not alter ridge_cv (or solve_ridge_formula).

Note that it we should permute both X and y with the same indices. Shuffling only y, or constructing a new time basis after shuffling it, would change the regression problem. The original X, y, and calendar years below remain chronological for final fitting and plotting.

The $\lambda$ now chosen is quite small. This should lead to estimates that are much more noisy (compared to the previous case).

Perhaps, the trend estimate now is too smooth.

So, depending on how the folds are chosen, CV can give $\lambda$ that is quite high or quite low. Below, we plot the folds used in the second CV (random folds).

The distinct folds are hard to see now because they are randomly present all over the data.

Here is some intuition for why the two CV methods are behaving so differently. This is because the two CV schemes assess different prediction tasks. With consecutive folds, we remove a long stretch of observations, so the fitted curve must predict across a large gap or extrapolate beyond the training data. Flexible (overfitting) fits can perform poorly in these regions, leading CV to favor a large $\lambda$ (smooth fit). With random folds, most held-out observations have nearby training observations. A flexible fit can use these neighbors to track the data and predict the held-out values accurately, favoring a small $\lambda$. Thus, the preferred amount of smoothing depends strongly on whether we assess prediction across long gaps or interpolation between nearby observations.

Next let us try Bayesian regularization for this problem.

The $\lambda$-values along with the log posterior and posterior probabilities are printed below.

Here we are only using 25 $\lambda$ values (as in CV). The posterior probability is 1 exactly for one $\lambda$ (and zero for all other $\lambda$’s). This $\lambda$ value is the posterior maximizer. This value is small but not as small as the $\lambda$ selected by the random-fold CV. So the fitted function will still be noisy (tracking the original data)but not as noisy as the ridge estimate obtained by the random-fold CV method.

Below are the estimates of the trend functions corresponding to each of these $\lambda$ choices.

It is hard to spot the differences between the random CV estimate and the Bayesian estimate. One can zoom in to see that the Bayesian estimate is indeed smoother.

One can also run the full Bayesian method to get posterior samples of $\beta_0, \dots, \beta_{n-1}$, which also tracks uncertainty in $\lambda$. The relevant formulae here are given here. We work with the following prior:

We also treat $\sigma$ and $\lambda$ as unknown parameters and assign the prior:

where

As derived in the lecture, the posterior of $\beta, \sigma, \lambda$ then becomes:

and

and

The key here is the posterior for $\lambda$. We compute this posterior numerically below by taking a grid of values for $\log \lambda$. The computation proceeds more stably if we evaluate the posterior density of $\log \lambda$:

The code for doing all this computation is given below. Explanation for how the code relates to the above formulae is explained after the code.

We need to select a large number of $\lambda$ values, and evaluate the posterior at these values. Depending on the number of $\lambda$ values chosen, and the number of posterior samples we want to obtain, this can be computationally slower.

Here is an explanation of the above code, and how it relates to the posterior formulae.

Write

The quantities $X^TX$, $b$, and $c$ do not depend on $\lambda$, so the code computes them only once. For each value of $\lambda$, we need

Here $m_\lambda$ is the conditional posterior mean of $\beta$, and $q_\lambda/2$ is the rate parameter in the conditional Gamma distribution of $1/\sigma^2$.

Cholesky decomposition

Although these formulas involve a matrix inverse, the code does not construct that inverse. Instead, it uses a Cholesky decomposition:

The matrix $L_\lambda$ is lower triangular, meaning that all entries above its diagonal are zero. Its diagonal entries are positive: $$ L_\lambda=

$$

Such a factorization exists for every real symmetric positive definite matrix. In our problem, $A_\lambda$ has this property because $X$ has full column rank and, for every nonzero vector $v$,

The command L = np.linalg.cholesky(XTX + lam * J) computes this lower triangular factor.

The advantage of the factorization is that equations involving triangular matrices are easy to solve. For example, to solve $L_\lambda u=b$, we first obtain $u_1$ from the first equation:

The second equation then gives

We continue down the rows, finding one new unknown at each step. This is called forward substitution. Equations involving $L_\lambda^T$, which is upper triangular, can similarly be solved starting from the last row and working upward. This is called back substitution.

Computing the posterior mean

To compute $m_\lambda=A_\lambda^{-1}b$, we solve $A_\lambda m_\lambda=b$. Using the factorization, this becomes

We solve this through two triangular systems:

The command cho_solve((L, True), XTy) performs these two steps. The argument True specifies that the supplied factor is lower triangular. This computes the desired posterior mean without explicitly forming $A_\lambda^{-1}$.

Computing the determinant

The same Cholesky factor also gives the determinant. The determinant of a triangular matrix is the product of its diagonal entries, so

Consequently,

The code therefore obtains $-\frac12\log|A_\lambda|$ as -np.log(np.diag(L)).sum(). Working with a sum of logarithms avoids forming a product that might be extremely large or small.

Once $m_\lambda$ is available, the remaining quantity is

The code stores $m_\lambda$ in a column of beta_means and stores $q_\lambda/2$ in rates. These quantities will be needed again during posterior sampling, so saving them avoids repeating their calculation.

Computing the posterior grid weights

Let $x_g=\log\lambda_g$, for $g=1,\ldots,G$, be equally spaced grid points between $\log\lambda_{\min}$ and $\log\lambda_{\max}$.

The logarithm of the posterior density of $\log\lambda$, up to an additive constant that does not depend on $g$, is

These values are stored in logpost_loglambda.

To turn them into posterior grid weights, the code computes

Subtracting $h_{\max}$ does not change the normalized weights: it multiplies every unnormalized weight by the same constant. It improves numerical stability by making the largest exponential equal to one.

Because the grid is equally spaced in $\log\lambda$, its common spacing cancels when the weights are normalized. The resulting weights sum to one and approximate posterior masses associated with the grid points. They are stored in prob_loglambda.

If the grid spacing is $\Delta x$, then $w_g/\Delta x$ approximates the posterior density of $\log\lambda$ at $x_g$. These density values are stored in post_loglambda.

Sampling from the posterior

Posterior sampling follows the factorization $$ p(\beta,\sigma,\lambda\mid y)¶

p(\lambda\mid y), p(\sigma\mid\lambda,y), p(\beta\mid\sigma,\lambda,y). $$

First, the code samples grid indices using probabilities $w_1,\ldots,w_G$. A sampled index $g$ corresponds to $\lambda=\lambda_g$ and gives immediate access to its saved mean $m_{\lambda_g}$ and rate $r_g=q_{\lambda_g}/2$.

Next, conditional on this sampled $\lambda_g$, we draw

NumPy specifies a Gamma distribution using its shape and scale, rather than its shape and rate. Since scale is the reciprocal of rate, the appropriate call is np.random.gamma(shape, 1/rate). The code draws all the required values of $\sigma$ at once.

Finally, we must sample

To do this, generate a vector $z$ of independent standard normal variables and solve

Then set

To verify this construction, observe that

and hence

Therefore, the constructed vector has mean $m_\lambda$ and covariance $\sigma^2A_\lambda^{-1}$, as required.

The code computes $u$ using solve_triangular(L.T, z, lower=False). The transpose matters: we solve with $L_\lambda^T$, which is upper triangular.

Avoiding repeated calculations

Several posterior draws may select the same grid value of $\lambda$. The sampling loop therefore visits each distinct sampled grid index only once. For that index, it recomputes one Cholesky factor and uses it for all the corresponding draws.

The independent standard normal vectors are arranged as columns of a matrix, allowing the triangular solver to process them together. Each column uses its own sampled $\sigma$.

Thus, the code saves the coefficient means and Gamma rates from the grid calculation, while recomputing Cholesky factors only for the distinct grid values selected during sampling. This avoids storing a large collection of matrix factors and avoids a separate factorization for every posterior draw.

A Cholesky factorization requires order $n^3$ operations, whereas solving a triangular system for one vector requires order $n^2$ operations. Reusing the factors can therefore substantially reduce computation.

Constructing the posterior curves

For each sampled coefficient vector $\beta^{(s)}$, the corresponding mean curve is

The command muhats = X @ betahats computes all these curves together. Their average,

approximates the posterior mean $E[\mu\mid y]$.

The plot shows the observations in light gray, the sampled mean curves in light coral, and their average in black. The computation has a grid approximation from discretizing $\lambda$ and Monte Carlo error from using finitely many posterior draws.

To visualize these posterior samples better, we can zoom into certain shorter time intervals.

This notebook compares cross-validation and Bayesian regularization for estimating the trend in weekly U.S. initial unemployment claims before March 2020. Consecutive-fold CV selects a large $\lambda$ and produces a very smooth fit, while random-fold CV selects a much smaller $\lambda$ and follows more short-term variation. Bayesian tuning gives an intermediate amount of smoothing in this example. These differences illustrate how the choice of validation folds and modeling assumptions influences the estimated trend. Full Bayesian inference additionally accounts for uncertainty in $\lambda$, $\sigma$, and $\beta$, producing posterior samples of the trend. Cholesky decomposition and reuse of intermediate calculations make this computation more efficient.