41
$\begingroup$

I am having some issues with the derivation of the solution for ridge regression.

I know the regression solution without the regularization term:

$$\beta = (X^TX)^{-1}X^Ty.$$

But after adding the L2 term $\lambda\|\beta\|_2^2$ to the cost function, how come the solution becomes

$$\beta = (X^TX + \lambda I)^{-1}X^Ty.$$

$\endgroup$
25
$\begingroup$

It suffices to modify the loss function by adding the penalty. In matrix terms, the initial quadratic loss function becomes $$ (Y - X\beta)^{T}(Y-X\beta) + \lambda \beta^T\beta.$$ Deriving with respect to $\beta$ leads to the normal equation $$ X^{T}Y = \left(X^{T}X + \lambda I\right)\beta $$ which leads to the Ridge estimator.

$\endgroup$
  • 1
    $\begingroup$ How come the derivative of $\lambda \beta^T \beta$ is equal to $\lambda I \beta$ $\endgroup$ – user34790 Sep 4 '13 at 16:21
  • 4
    $\begingroup$ @user34790 It's not. It's equal to $2\lambda\beta$. But the 2 cancels with similar 2s on the other terms. Of course, the $I$ factor is like a factor of 1 in "regular" algebra, you can multiply it anywhere you like without changing anything. $\endgroup$ – Bill Sep 4 '13 at 18:42
  • 4
    $\begingroup$ @bill: here you need the $I$ to get a matrix of the correct dimension so the addition works with $X^TX$: $\lambda$ is just a scalar $\endgroup$ – Henry Aug 4 '15 at 17:23
49
$\begingroup$

Let's build on what we know, which is that whenever the $n\times p$ model matrix is $X$, the response $n$-vector is $y$, and the parameter $p$-vector is $\beta$, the objective function

$$f(\beta) = (y - X\beta)^\prime(y - X\beta)$$

(which is the sum of squares of residuals) is minimized when $\beta$ solves the Normal equations

$$(X^\prime X)\beta = X^\prime y.$$

Ridge regression adds another term to the objective function (usually after standardizing all variables in order to put them on a common footing), asking to minimize

$$(y - X\beta)^\prime(y - X\beta) + \lambda \beta^\prime \beta$$

for some non-negative constant $\lambda$. It is the sum of squares of the residuals plus a multiple of the sum of squares of the coefficients themselves (making it obvious that it has a global minimum). Because $\lambda\ge 0$, it has a positive square root $\nu^2 = \lambda$.

Consider the matrix $X$ augmented with rows corresponding to $\nu$ times the $p\times p$ identity matrix $I$:

$$X_{*} = \pmatrix{X \\ \nu I}$$

When the vector $y$ is similarly extended with $p$ zeros at the end to $y_{*}$, the matrix product in the objective function adds $p$ additional terms of the form $(0 - \nu \beta_i)^2 = \lambda \beta_i^2$ to the original objective. Therefore

$$(y_{*} - X_{*}\beta)^\prime(y_{*} - X_{*}\beta) = (y - X\beta)^\prime(y - X\beta) + \lambda \beta^\prime \beta.$$

From the form of the left hand expression it is immediate that the Normal equations are

$$(X_{*}^\prime X_{*})\beta = X_{*}^\prime y_{*}.$$

Because we adjoined zeros to the end of $y$, the right hand side is the same as $X^\prime y$. On the left hand side $\nu^2 I=\lambda I$ is added to the original $X^\prime X$. Therefore the new Normal equations simplify to

$$(X^\prime X + \lambda I)\beta = X^\prime y.$$


Besides being conceptually economical--no new manipulations are needed to derive this result--it also is computationally economical: your software for doing ordinary least squares will also do ridge regression without any change whatsoever. (It nevertheless can be helpful in large problems to use software designed for this purpose, because it will exploit the special structure of $X_{*}$ to obtain results efficiently for a densely spaced interval of $\lambda$, enabling you to explore how the answers vary with $\lambda$.)

Another beauty of this way of looking at things is how it can help us understand ridge regression. When we want to really understand regression, it almost always helps to think of it geometrically: the columns of $X$ constitute $p$ vectors in a real vector space of dimension $n$. By adjoining $\nu I$ to $X$, thereby prolonging them from $n$-vectors to $n+p$-vectors, we are embedding $\mathbb{R}^n$ into a larger space $\mathbb{R}^{n+p}$ by including $p$ "imaginary", mutually orthogonal directions. The first column of $X$ is given a small imaginary component of size $\nu$, thereby lengthening it and moving it out of the space generated by the original $p$ columns. The second, third, ..., $p^\text{th}$ columns are similarly lengthened and moved out of the original space by the same amount $\nu$--but all in different new directions. Consequently, any collinearity present in the original columns will immediately be resolved. Moreover, the larger $\nu$ becomes, the more these new vectors approach the individual $p$ imaginary directions: they become more and more orthonormal. Consequently, the solution of the Normal equations will immediately become possible and it will rapidly become numerically stable as $\nu$ increases from $0$.

This description of the process suggests some novel and creative approaches to addressing the problems Ridge Regression was designed to handle. For instance, using any means whatsoever (such as the variance decomposition described by Belsley, Kuh, and Welsch in their 1980 book on Regression Diagnostics, Chapter 3), you might be able to identify subgroups of nearly collinear columns of $X$, where each subgroup is nearly orthogonal to any other. You only need adjoin as many rows to $X$ (and zeros to $y$) as there are elements in the largest group, dedicating one new "imaginary" dimension for displacing each element of a group away from its siblings: you don't need $p$ imaginary dimensions to do this.

$\endgroup$
  • 2
    $\begingroup$ The last author of the book is Welsch, not Welsh. $\endgroup$ – Mark L. Stone Mar 25 '16 at 17:25
  • 1
    $\begingroup$ Whoa, this just blew my mind. Is there any discussion about what happens when this is generalized outside linear models, i.e. to glm's? The penalty shouldn't line up to be the same as ridge regression...but this interpretation implies that it would still be a potential useful estimator! $\endgroup$ – Cliff AB Mar 25 '16 at 19:14
  • 2
    $\begingroup$ @Cliff That's a very interesting suggestion. Since, however, GLM estimates depend in a more complicated way on $X$ and their estimators cannot usually be factored in the form $$\hat\beta = g(X)\cdot h(y)$$ as they are for OLS (where $g(X)=(X^\prime X)^{-1}X^\prime$ and $h(y)=y$), it may be difficult to establish a useful relationship between imposing a penalty function and modifying the columns of $X$. In particular, it's unclear how the values in $y$ would need to be augmented in order to make this work. $\endgroup$ – whuber Mar 25 '16 at 20:45
  • 1
    $\begingroup$ Yes, it would take some thought to try establish what the penalty is, but I'm not so concerned about that. The idea of what $y_*$ to use is not generally easy either...except perhaps in the case of logistic regression, where we could add two $y_*$'s; one of 0's and one of 1's. This augmentation would then be a more general version of the "+2 binomial estimator" (there's a more proper name for this estimator that I'm blanking on, which is basically when you are estimating $p$ from a binomial distribution using the posterior mean as the estimate with a uniform prior on $p$). $\endgroup$ – Cliff AB Mar 25 '16 at 21:15
  • $\begingroup$ @Mark Thank you for the correction. You can tell I was going from memory... :-). $\endgroup$ – whuber Apr 12 '18 at 17:59
20
$\begingroup$

The derivation includes matrix calculus, which can be quite tedious. We would like solve the following problem: \begin{equation} \min_\beta (Y-\beta^T X)^T(Y-\beta^T X)+\lambda \beta^T \beta \end{equation}

Now note that \begin{equation} \frac{\partial (Y-\beta^T X)^T (Y-\beta^T X)}{\partial \beta}=-2X^T(Y-\beta^T X) \end{equation} and \begin{equation} \frac{\partial \lambda \beta^T \beta}{\partial \beta}=2\lambda\beta. \end{equation} Together we get to the first order condition \begin{equation} X^TY = X^TX\beta + \lambda\beta. \end{equation} Isolating $\beta$ yields the solution: \begin{equation} \beta = (X^TX+ \lambda I )^{-1}X^T Y. \end{equation}

$\endgroup$
9
$\begingroup$

I have recently stumbled upon the same question in the context of P-Splines and as the concept is the same I want to give a more detailed answer on the derivation of the ridge estimator.

We start with a penalized criterion function that differs from the classic OLS-criterion function by its penalization term in the last summand:

$Criterion_{Ridge} = \sum_{i=1}^{n}(y_i-x_i^T\beta)^2 + \lambda \sum_{j=1}^p\beta_j^2$

where

  • $p=$ the amount of covariables used in the model
  • $x_i^T\beta = $ your standard linear predictor
  • the first summand respresents the MSE (squared divergence of the prediction from the actual value) that we want to minimize as usual
  • the second summand represents the penalization we apply on the coefficients. Here we are in the Ridge-context which implies a Euclidian Distance Measure and therefore the degree of 2 in the penalization term. In the case of a Lasso-Penalization we would apply a degree of 1 and yield a totally different estimator.

We can rewrite this criterion in matrix-notation and further break it down:

$Criterion_{Ridge} = (y-X\beta)^T(y-X\beta) + \lambda\beta^T\beta$

$ = y^Ty - \beta^TX^Ty - y^TX\beta+ \beta^Tx^TX\beta + \lambda\beta^T\beta$

$ = y^Ty - \beta^TX^Ty - \beta^TX^Ty + \beta^TX^TX\beta + \beta^T\lambda I\beta$ with $I$ being the identity matrix

$ = y^Ty - 2\beta^TX^Ty + \beta^T(X^TX + \lambda I)\beta$

Now we search for the $\beta$ that minimizes our criterion. Amongst others we make use of the matrix differentiation rule $\frac{\partial x^TAx}{\partial x} = (A+A^T)x \overset{\text{A symmetric}}{=} 2Ax$ which we can apply here as $(X^TX + \lambda I) \in \mathbb{R}^{n \times n}$:

$\frac{\partial Criterion_{Ridge} }{\partial\beta} = -2X^Ty + 2(X^TX + \lambda I)\beta \overset{!}{=}0$

$(X^TX + \lambda I)\beta = X^Ty$

$\overset{\text{et voilà}}{\Rightarrow} \hat\beta = (X^TX + \lambda I)^{-1} X^Ty$

$\endgroup$
  • $\begingroup$ @Jahn, can you please explain how $$y^TX\beta$$ became $$\beta ^TX^Ty$$? I think you just applied transpose on it, right. But, you cannot just apply transpose on one term without applying it on all equation. What am I missing here? $\endgroup$ – theateist Jun 17 '18 at 23:31
  • 1
    $\begingroup$ @theateist A transposed scalar is the same scalar. $\endgroup$ – Konstantin Jun 20 '18 at 14:40
2
$\begingroup$

There are a few important things that are missing in the answers given.

  1. The solution for $\beta$ is derived from the first-order necessary condition: $\frac{\partial f_{ridge}(\beta, \lambda)}{\partial \beta} = 0$ which yields $\beta = (X^TX+ \lambda I )^{-1}X^T Y$. But is this sufficient? That is, the solution is a global minimum only if $f_{ridge}(\beta, \lambda)$ is strictly convex. This can be shown to be true.

  2. Another way to look at the problem is to see the equivalence between $f_{ridge}(\beta, \lambda)$ and $f_{OLS}(\beta) = (Y-\beta^T X)^T(Y-\beta^T X)$ constrained to $||\beta||^2_2 \leq t$. OLS stands for Ordinary Least Squares. From this perspective $f_{ridge}(\beta, \lambda)$ is just the Lagrangian function used to find the global minima of the convex objective function $f_{OLS}(\beta)$ constrained with convex function $||\beta||^2_2$.

A good explanation of these points and the derivation of $\beta$ can be found in these fine lecture notes: http://math.bu.edu/people/cgineste/classes/ma575/p/w14_1.pdf

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.