I am having some issues with the derivation of the solution for ridge regression.
I know the regression solution without the regularization term is given by:
$$\beta = (X^\top X)^{-1}X^\top y.$$
But after adding the L2 term $\lambda\|\beta\|_2^2$ to the cost function, how does the solution become:
$$\beta = (X^\top X + \lambda I)^{-1}X^\top y?$$
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.
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 orthogonal. 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.
A method of performing this was proposed by Rolf Sundberg some 30 years ago. See Continuum Regression and Ridge Regression. Journal of the Royal Statistical Society. Series B (Methodological) Vol. 55, No. 3 (1993), pp. 653-659.
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}
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
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$
There are a few important things that are missing in the answers given.
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.
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
