An Algorithm for Bayesian Ridge Regression¶

Introduction¶

Let $\params = \left(\noise^2, \w\right)$ denote parameters for a linear regression model with weights $\w$ and normally distributed errors of variance $\noise^2$.

If $\X$ represents an $n\times p$ matrix of full rank of $p$ regressors and $n$ rows, then $\params$ specifies a probability distribution over possible target values $\y$:

Suppose we observe $\y$ and assume $\y$ is generated from a linear model of unknown parameters. A Bayesian approach to inference seeks to quantify our belief in the unknown parameters $\params$ given the observation.

Applying Bayes’ theorem, we can rewrite the probability as

where we refer to

• $\Prb\left(\params|\y\right)$ as the posterior distribution

• $\Prb\left(\y|\params\right)$ as the likelihood function

• $\Prb\left(\params\right)$ as the prior distribution

The prior distribution describes our belief of $\params$ before observing $\y$ and the posterior distribution describes our updated belief after observing $\y$.

Suppose the prior distribution can be expressed as

where $\priorf(\cdot,\hp)$ denotes a family of probability distributions parameterized by what we call a hyperparameter $\hp$.

Traditional approaches to Bayesian linear regression have used what are called conjugate priors. A family of priors $\priorf(\cdot,\hp)$ is conjugate if the posterior also belongs to the family

Conjugate priors are mathematically convenient as successive observations can be viewed as making successive updates to the parameters of a family of distributions, but requiring $\priorf$ to be conjugate is a strong assumption.

We’ll instead describe an algorithm where

1. Priors are selected to shrink $\w$, reflecting the prior hypothesis that $\w$ is not predictive, and be approximately noninformative for other parameters.

2. We fully integrate over hyperparameters so that no arbitrary choice of $\hp$ is required.

Let’s first consider what it means for a prior to be noninformative.

How to Pick Non-informative Priors?¶

Suppose $\y$ is data generated from a normal distribution of mean $0$ but unknown variance. Let $\sigma$ denote standard deviation and let $\like$ denote the likelihood function

Suppose we impose a uniform prior over $\sigma$ so that the posterior is

and the cumulative distribution function is

where $\N$ is some normalizing constant.

But now let’s suppose that instead of standard deviation we parameterize over variance. Making the substitution $u=\sigma^2$ into the cumulative distribution function gives us

Thus, we see that choosing a uniform prior over $\sigma$ is equivalent to choosing the improper prior

over variance. In general, if $u=\varphi\left(\params\right)$ is an alternative way of parameterizing the likelihood function where $\varphi$ is some monotonic, one-to-one, onto function. Then a uniform prior over $u$ is equivalent to choosing a prior

over $\params$.

So, when does it make sense to use a uniform prior if the choice is sensitive to parameterization?

Let’s consider how the likelihood is affected by changes in the observed value $\multt{\y}$.

Likelihood function for the standard deviation of normally distributed data with zero mean, $n=10$, and different values of $\multt{\y}$.¶

As we can see from the figure, as $\multt{\y}$ is increased the shape of the likelihood function changes: it becomes less peaked and more spread out.

Observe that we can rewrite the likelihood function as

where

Thus, in the parameter $\log \sigma$, likelihood has the form

where

And we say that the likelihood function is data translated in $\log \sigma$ because everything is known about the likelihood curve except the location and the value of the observation serves only to shift the curve.

Likelihood function for the log standard deviation of normally distributed data with zero mean, $n=10$, and different values of $\multt{\y}$.¶

When the likelihood function is data translated in a parameter, then it makes sense to use a uniform function for a noninformative prior. Before observing data, we have no reason to prefer one range of parameters $\left[a, b\right]$ over another range of the same width $\left[t + a, t + b\right]$ because they will be equivalently placed relative to the likelihood curve if the observed data translates by $t$.

Now, let’s return to picking our priors.

Picking Regularized Bayesian Linear Regression Priors¶

For the parameter $\noise$, we use the noninformative prior

which is equivalent to using a uniform prior over the parameter $\log \noise$. For $\w$, we want an informative prior that shrinks the weights, reflecting a prior belief that weights are non-predictive. Let $\hp$ denote a hyperparameter that controls the degree of shrinkage. Then we use the spherical normal distribution with covariance matrix $\left(\frac{\noise}{\lda}\right)^2\I$

Note that we haven’t described yet how $\hp$ parameterizes $\lda$ and we’ll also be integrating over $\hp$ so we additionally have a prior for $\hp$ (called a hyperprior) so that

Our goal is for the prior $\Pr(\hp)$ to be noninformative. So we want to know: In what parameterization, would $\Pr(\y|\hp)$ be data translated?

How to Parameterize the shrinkage prior?¶

Before thinking about how to parameterize $\lda$, let’s characterize the likelihood function for $\lda$.

Expanding the integrand, we have

where

and

Observe that $\#1$ is independent of $\w$, so the integral with respect to $\w$ is equivalent to integrating over an unnormalized Gaussian.

Thus,

Next let’s consider the integral over $\noise$.

Put

Then we can rewrite

After making a change of variables, we can express the integral with respect to $\noise$ as a form of the Gamma function.

Thus, we can write

Let $\U$, $\bSigma$, and $\V$ denote the singular value decomposition of $\X$ so that

Let $\eig_1, \ldots \eig_p$ denote the non-zero diagonal entries of $\bSigma$. Put $\eigd = \multt{\bSigma}$. Then $\eigd$ is a diagonal matrix with entries $\eig_1^2, \ldots, \eig_p^2$ and

Put

Then we can rewrite $h$ and $g$ as

And

Here we adopt the terminology

Put

Then $\shrink$ is a monotonically increasing function that ranges from $0$ to $1$, and we can think of $\shrink$ as the average shrinkage factor across the eigenvectors.

Now, let’s make an approximation by replacing individual eigenvector shrinkages with the average:

Substituting the approximation into the likelihood then gives us

We see that the approximated likelihood can be expressed as a function of

and it follows that the likelihood is approximately data translated in $\log r(\lda)$.

Thus, we can achieve an approximately noninformative prior if we let $\hp$ represent the average shrinkage, put

and use the hyperprior

Making Predictions¶

The result of fitting a Bayesian model is the posterior distribution $\Pr(\params|\y)$. Let’s consider how we can use the distribution to make predictions given a new data point $\x'$.

We’ll start by computing the expected target value

And

To compute variance, we have

and

Starting with $\Ex{\noise^2|\y}$, recall that

And

Thus,

So,

And

It follows that

For $\Ex{\w\w^\top|\y}$, we have

Outline of Algorithm¶

We’ve seen that computing statistics about predictions involve integrating over the posterior distribution $\Prb(\hp|\y)$. We’ll briefly sketch out an algorithm for computing such integrals. We describe it only for computing the expected value of $\w$. Other expected values can be computed with straightforward modifications.

The proceedure SHRINK-INVERSE applies Newton’s algorithm for root-finding with $\shrink$ and $\shrink'$ to compute $\shrinki$.

To compute the hyperparameter posterior integration, we make use of an adaptive quadrature algorithm. Quadratures approximate an integral by a weighted sum at different points of evaluation.

In general, the more points of evaluation used, the more accurate the approximation will be. Adaptive quadrature algorithms compare the integral approximation at different levels of refinement to approximate the error and increase the number of points of evaluation until a desired tolerance is reached.

Experimental Results¶

To better understand how the algorithm works in practice, we’ll set up a small experimental data set, and we’ll compare a model fit with the Bayesian algorithm to Ordinary Least Squares (OLS), and a ridge regression model fit so as to minimize error on a Leave-one-out Cross-validation (LOOCV) of the data set.

Full source code for the experiment is available at github.com/rnburn/bbai/example/03-bayesian.py.

np.random.seed(0)

def generate_correlation_matrix(p, param):
    res = np.zeros(shape=(p, p))
    for s in range(p):
        for t in range(0, s+1):
            corr = param
            if s == t:
                corr = 1.0
            res[s, t] = corr
            res[t, s] = corr
    return res

def generate_design_matrix(n, K):
    mean = np.zeros(K.shape[0])
    return np.random.multivariate_normal(mean, K, size=n)

def generate_weights(p):
    return np.random.normal(size=p)

def generate_data_set(n, K, signal_noise_ratio):
    p = K.shape[0]
    X = generate_design_matrix(n, K)
    w = generate_weights(p)
    signal_var = np.dot(w, np.dot(K, w))
    w *= np.sqrt(signal_noise_ratio / signal_var)
    y = np.dot(X, w) + np.random.normal(size=n)
    return X, y, w

p = 10
n = 20
signal_noise_ratio = 1.0
K = generate_correlation_matrix(p, 0.5)
X, y, w_true = generate_data_set(n, K, signal_noise_ratio)

# OLS
from sklearn.linear_model import LinearRegression
model_ols = LinearRegression(fit_intercept=False)
model_ols.fit(X, y)

# Bayesian Linear Ridge Regression
from bbai.glm import BayesianRidgeRegression
model_bay = BayesianRidgeRegression(fit_intercept=False)
model_bay.fit(X, y)

# Ridge Regression fit to optimize LOOCV error
#
# For details refer to
#     https://arxiv.org/abs/2011.10218
from bbai.glm import RidgeRegression
model_rr = RidgeRegression(fit_intercept=False)
model_rr.fit(X, y)

def compute_prediction_error_variance(K, w_true, w):
    delta = w - w_true
    noise_variance = 1.0
    return noise_variance + np.dot(delta, np.dot(K, delta))

err_variance_ols = compute_prediction_error_variance(K, w_true, model_ols.coef_)
err_variance_bay = compute_prediction_error_variance(K, w_true, model_bay.weight_mean_vector_)
err_variance_rr = compute_prediction_error_variance(K, w_true, model_rr.coef_)

err_ols = y - np.dot(X, model_ols.coef_)
s_ols = np.dot(err_ols, err_ols) / (n - p)
print("noise_variance_ols =", s_ols)
print("noise_variance_bay =", model_bay.noise_variance_mean_)

Appendix A: Comparison with Other Bayesian Algorithms¶

Here, we’ll give a brief comparison of the algorithm presented to scikit-learn’s algorithm for Bayesian ridge regression and describe the advantages.

Scikit-learn’s algorithm makes use of conjugate priors and because of that is restricted to use the Gamma prior which requires four hyperparameters chosen arbitrarily to be small values. Additionally, it requires initial values for parameters $\alpha$ and $\lambda$ that are then updated from the data.

The algorithm’s performance can be sensitive to the choice of values for these parameters, and scikit-learn’s documentation provides a curve fitting example where the defaults perform poorly.

# Author: Yoshihiro Uchida <nimbus1after2a1sun7shower@gmail.com>

import numpy as np
import matplotlib.pyplot as plt

from sklearn.linear_model import BayesianRidge


def func(x):
    return np.sin(2 * np.pi * x)


# #############################################################################
# Generate sinusoidal data with noise
size = 25
rng = np.random.RandomState(1234)
x_train = rng.uniform(0.0, 1.0, size)
y_train = func(x_train) + rng.normal(scale=0.1, size=size)
x_test = np.linspace(0.0, 1.0, 100)


# #############################################################################
# Fit by cubic polynomial
n_order = 3
X_train = np.vander(x_train, n_order + 1, increasing=True)
X_test = np.vander(x_test, n_order + 1, increasing=True)

# #############################################################################
# Plot the true and predicted curves with log marginal likelihood (L)
reg = BayesianRidge(tol=1e-6, fit_intercept=False, compute_score=True)
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
for i, ax in enumerate(axes):
    # Bayesian ridge regression with different initial value pairs
    if i == 0:
        init = [1 / np.var(y_train), 1.0]  # Default values
    elif i == 1:
        init = [1.0, 1e-3]
        reg.set_params(alpha_init=init[0], lambda_init=init[1])
    reg.fit(X_train, y_train)
    ymean, ystd = reg.predict(X_test, return_std=True)

    ax.plot(x_test, func(x_test), color="blue", label="sin($2\\pi x$)")
    ax.scatter(x_train, y_train, s=50, alpha=0.5, label="observation")
    ax.plot(x_test, ymean, color="red", label="predict mean")
    ax.fill_between(
        x_test, ymean - ystd, ymean + ystd, color="pink", alpha=0.5, label="predict std"
    )
    ax.set_ylim(-1.3, 1.3)
    ax.legend()
    title = "$\\alpha$_init$={:.2f},\\ \\lambda$_init$={}$".format(init[0], init[1])
    if i == 0:
        title += " (Default)"
    ax.set_title(title, fontsize=12)
    text = "$\\alpha={:.1f}$\n$\\lambda={:.3f}$\n$L={:.1f}$".format(
        reg.alpha_, reg.lambda_, reg.scores_[-1]
    )
    ax.text(0.05, -1.0, text, fontsize=12)

plt.tight_layout()
plt.show()

Shows predictions for sklearn’s BayesianRidge model on a curve fitting example. The left shows the predictions with default parameters (performs poorly), the right shows the predictions after initial parameters have been tweaked (performs better).¶

In comparison, the algorithm we presented requires no initial parameters; and because the hyperparameter is integrated over, poor performing values contribute little to the posterior probability mass.

We can see that our approach handles the curve-fitting example without requiring any tweaking.

# #############################################################################
# Plot the true and predicted curves for bbai's BayesianRidgeRegression model
from bbai.glm import BayesianRidgeRegression
reg = BayesianRidgeRegression(fit_intercept=False)
fig, ax = plt.subplots(1, 1, figsize=(4, 4))

# Note: there are no parameters to tweak
reg.fit(X_train, y_train)
ymean, ystd = reg.predict(X_test, return_std=True)

ax.plot(x_test, func(x_test), color="blue", label="sin($2\\pi x$)")
ax.scatter(x_train, y_train, s=50, alpha=0.5, label="observation")
ax.plot(x_test, ymean, color="red", label="predict mean")
ax.fill_between(
    x_test, ymean - ystd, ymean + ystd, color="pink", alpha=0.5, label="predict std"
)
ax.set_ylim(-1.3, 1.3)
ax.legend()

plt.tight_layout()
plt.show()

Shows predictions from our Bayesian ridge regression algorithm on the curve fitting example. In comparison to sklearn, our algorithm perform well without requiring any tweaking.¶

The full curve-fitting comparison example is available here.