In Machine Learning, why is Regularization called Regularization?

Shifting the eigenvalues of a Symmetric Matrix

Suppose we are asked to solve the linear inverse problem from above \(A\beta = y\), but that this time \(A\in \mathbb{R}^{N\times N}\) is a symmetric matrix. The spectral theorem for symmetric matrices tells us that \(A\) can be represented as

\[ A = QDQ^T \] where \(Q\) is an orthogonal matrix and \(D\) is a diagonal matrix. Moreover, the columns of \(Q\) are the eigenvectors of \(A\), and the diagonal elements of \(D\) are the corresponding (real) eigenvalues. Here’s a proof if you care for it.

This representation shows us exactly what the theoretical difficulty is in the inverse problem. Inverting \(Q\) is easy. Indeed orthogonal matrices are always invertible, with the inverse given by the transpose: \(Q^{-1} = Q^T\). The geometric significance of orthogonal matrices comes from the fact (basically their definition) that they preserve the inner product of vectors: If we denote the inner product by \(\langle \cdot,\cdot\rangle\) then for any \(x,y \in R^N\)

\[ \langle x, y \rangle = \langle Qx, Qy \rangle \] Thus (using the definition of the transpose) \(\langle x, y \rangle = \langle Q^TQx, y \rangle\). Since this holds for any \(y\), then \(x = Q^TQx\) for all \(x\). Since this holds for all \(x\) then \(Q^TQ = I\) and so \(Q^T\) is the inverse of \(Q\).

No, the difficulty is simply in inverting the diagonal matrix \(D = \text{diag}(\sigma_1, \sigma_2, ..., \sigma_N)\). If none of the eigenvalues \(\sigma_i\) are \(0\) then \(D^{-1} = \text{diag}(\sigma_1^{-1}, \sigma_2^{-1}, ..., \sigma_N^{-1})\). In this case there isn’t any direction that \(D\) (and hence \(A\)) squashes into \(0\). However, if some of the \(\sigma_i\)’s are \(0\) then we can not invert \(D\) and the problem will fail to satisfy at least one of the first 2 Hadamard conditions. Even if none of the \(\sigma_i\)’s are exactly \(0\), some may be numerically very close to \(0\) in comparison to the others. In that case the value of their reciprocals may be enormously large and may lead to numerical instability in the problem, violating the 3rd Hadamard condition. This would be a big problem in practice because computers hate mixing floating point numbers that are drastically different in size.

To address this issue we note that we can shift the eigenvalues of \(A\) by adding a multiple of an identity matrix:

\[ A \to A + \lambda I \] If \(\sigma\) is an eigenvalue of \(A\) then \(\sigma+\lambda\) is an eigenvalue of \(A+\lambda I\). This is because every vector is an eigenvector of \(I\): If \(Av = \sigma v\) then trivially \(\lambda Iv = \lambda v\), so adding gives \((A+\lambda I)v = (\sigma + \lambda)v\). This can also be seen form the representation above: \[ A + \lambda I = QDQ^T + \lambda QQ^T = Q(D + \lambda I)Q^T = Q\tilde{D}Q^T \] where \(\tilde{D}= \text{diag}(\sigma_1 + \lambda, \sigma_2 + \lambda, ..., \sigma_N + \lambda)\). Therefore, if we choose \(\lambda > \min\{\sigma_1, \sigma_2, ..., \sigma_N\}+ \delta\) for some \(\delta > 0\), then all the shifted eigenvalues satisfy \(\sigma_i + \lambda > \delta > 0\) and the new shifted problem

\[ (A+\lambda I)x = y \] will be solvable with solution given by \[ x = Q(D + \lambda I)^{-1}Q^Ty \] If \(\lambda\) is sufficiently large this inverse will exist and will be numerically stable (all of the eigenvalues will have been shifted away from 0).

Making the change \(A \to A + \lambda I\) regularized the problem into one that satisfied Hadamard’s conditions, which is fundamentally the point of regularization. The change we made was essentially to replace the \(A^{-1}\) with the approximation \((A+\lambda I)^{-1}\), but we could have used other approximations as well, for example partial sums of the Neumann Series Expansion of the \(A^{-1}\). Regardless, the general principle illustrated above is basically to replace one problem by an approximate problem that does not suffer the same existence/stability issues.

Regularization as a perturbation of an invertible matrix

Above we regularized the ill-posed problem \(A\beta = y\) by replacing it with the problem \((A+\lambda I)x = y\). Let’s go a bit deeper with this process.You may skip this section on your first read.

Dividing by \(\lambda\), the problem \((A+\lambda I)x = y\) is equivalent to the problem

\[ (\epsilon A+I)x = \epsilon y \] where \(\epsilon := \frac{1}{\lambda}\). When \(\lambda\) is large \(\epsilon\) is small, and vice versa. Hence for large \(\lambda\) the matrix \(\epsilon A+I\) can be seen as a small perturbation from the identity matrix \(I\).

Now because \(I\) is invertible then for a small enough \(\epsilon\) (and hence for a large enough \(\lambda\)) the preturbed matrix \(\epsilon A+I\) is also invertible! Why? That’s a great question! The previous section gave one proof, but there are some much nicer ways to see why:

Proof: Consider the function \[ \det:\mathbb{R}^{N\times N} \to \mathbb{R} \] that maps a matrix to it’s determinant. Because the space \(\mathbb{R}^{N\times N}\) of \(N\times N\) matrices is just the Euclidean inner product space \(\mathbb{R}^{N^2}\) with some extra algebraic structure, and because \(\det(A)\) is a polynomial function of the elements of a matrix \(A\), \(\det\) is a continuous function on \(\mathbb{R}^{N\times N}\).

By Cramer’s Rule, a matrix \(A\) is invertible if and only if \(\det(A) \ne 0\). Because \(\det\) is a continuous function on \(\mathbb{R}^{N\times N}\) then the set of invertible matrices is an open subset of \(\mathbb{R}^{N\times N}\)! Hence for every invertible matrix \(L\in\mathbb{R}^{N\times N}\) and every arbitrary matrix \(A\in\mathbb{R}^{N\times N}\) there exists an \(\epsilon_0 > 0\) such that for all \(\epsilon < \epsilon_0\) the matrix \(\epsilon A+L\) is invertible. QED

There’s a version of the theorem in Banach spaces, but we don’t need it.

Noticed that the only thing we needed about the matrix \(I\) in the above proof was that it was invertible. Therefore, we never needed to restrict attention to just the identity matrix \(I\), but could have used any invertible matrix to regularize the problem:

\[ (A+\lambda L)x = y \]

where \(L\) is any convenient invertible matrix. Below, where we regularize OLS, we are not restricted to using the identity matrix \(I\) to regularize, but can use any invertible symmetric matrix (preferably one that is positive definite so that a minimizer continues to exist).

\(L^2\)-regularization of OLS

With this example we begin moving towards the statistical part of the post. One of the most widely known examples of regularization is what is often called \(L^2\)-regularization, or Tikhonov regularization of Ordinary Least Squares.

Suppose \(Y\) is a real valued random variable and \(\vec{X}\) is a random vector with values in \(\mathbb{R}^p\). Suppose that we have the conditional relationship

\[ Y \ \ | \ \ \vec{X} \sim \mathcal{N}(\ \langle \vec{X} , \beta\rangle \ ,\ \sigma^2) \]

where \(\mathcal{N}(\mu, \sigma^2)\) denotes the univariate normal distribution with mean \(\mu\) and variance \(\sigma^2\). Here (and everywhere else) the symbol \(\langle v,w\rangle\) represents the inner product of two vectors \(v\) and \(w\). This is the most natural probability model that leads to linear regression. In practice the parameters \(\beta\) and \(\sigma\) that specify this conditional distribution are unknown and it is desired that they be estimated from data.

In this canonical situation we assume that we have a data set \(\{ (Y_i, \vec{X}_i)\}_{i = 1}^N\) consisting of samples generated independently of one another from a fixed multivariate distribution for \((Y, \vec{X})\) (i.e. we assume our data was sampled IID). To fit the unknown parameters we use MLE. We may choose to use the conditional density \(p(Y|\vec{X})\) in the likelihood and this would make it conditional MLE. Or we may choose the unconditional multivariate density \(p(Y, \vec{X})\). However, if we assume that the marginal distribution of \(\vec{X}\), (i.e. \(p(\vec{X})\)), does not depend on either \(\beta\) or \(\sigma\) then because \(p(Y,\vec{X}) = p(Y|\vec{X})p(\vec{X})\) building the likelihood using either \(p(Y|\vec{X})\) or \(p(Y,\vec{X})\) will lead to the same maximization problem because they differ by a constant factor (constant in \(\beta\) and \(\sigma\) that is). So we will use the conditional density \(p(Y|\vec{X})\).

Because the data is assumed to be generated IID then the full likelihood of the data is

\[ \mathcal{L}(\beta, \sigma\ |\ \{ (Y_i, \vec{X}_i)\}) = \prod_{i = 1}^N p(Y_i|\vec{X}_i) = \prod_{i=1}^N\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\big(\frac{Y_i - \langle \vec{X}_i , \beta\rangle}{\sigma}\big)^2} = \frac{1}{(\sigma^22\pi)^{N/2}}e^{\sum_{i=1}^N-\frac{1}{2}\big(\frac{Y_i - \langle \vec{X}_i , \beta\rangle}{\sigma}\big)^2} \] Because the function \(f(x) := -\log(x)\) is decreasing we may instead minimize the negative of the log of this expression: \[ -\log(\mathcal{L}(\beta,\sigma\ |\ \{ (Y_i, \vec{X}_i)\})) = \frac{N}{2}\log(\sigma^22\pi) + \frac{1}{2\sigma^2}\sum_{i=1}^N\bigg(Y_i - \langle \vec{X}_i , \beta\rangle\bigg)^2 \] We first minimize with respect to \(\beta\) as this is necessary to do first before finding the minimizing value of \(\sigma\). To do this we need to minimize the only term that depends on \(\beta\), namely the sum of squares \(SSE(\beta) := \sum_{i=1}^N\bigg(Y_i - \langle \vec{X}_i , \beta\rangle\bigg)^2\) (hence Least Squares regression).

nice.paraboloid = function(x,y)
{
    return(x^2+0.5*y^2)
}
x = y = seq(from = -4, to = 4, by = 0.2)
z = outer(x, y, nice.paraboloid)

persp(x, y, z,
      main="Plot of a Non-degenerate 2D Paraboloid",
      theta = 30, phi = 15,
      col = "springgreen", shade = 0.5)

degenerate.paraboloid = function(x,y)
{
    return(x^2) #Does not change value as y changes
}
x = y = seq(from = -4, to = 4, by = 0.2)
z = outer(x, y, degenerate.paraboloid)

persp(x, y, z,
      main="Plot of a Degenerate 2D Paraboloid",
      theta = 30, phi = 15,
      col = "springgreen", shade = 0.5)

tricky.paraboloid = function(x,y)
{
    return(x^2+ 0.05*y^2) #Notice the coefficient of y is quite small
}
x = y = seq(from = -4, to = 4, by = 0.2)
z = outer(x, y, tricky.paraboloid)

persp(x, y, z,
      main="Plot of a Nearly-Degenerate 2D Paraboloid",
      theta = 30, phi = 15,
      col = "springgreen", shade = 0.5)

An Illustrative Example

Below we can see geometrically what regularization does. The sum of squares expression \((Y - X\beta)^T(Y - X\beta)\) is quadratic in \(\beta\), but may have a graph that is a degenerate paraboloid. This is what causes it to have multiple minimizers in OLS and what makes the matrix \(X^TX\) singular (more on this point in the next section). However the expression \(\lambda\langle \beta, \beta\rangle = \lambda\sum_{j = 1}^p\beta_j^2\) is a strictly positive-definite quadratic form. Its graph is a non-degenerate bowl shaped paraboloid.

Adding a non-degenerate paraboloid to something that is not bowl shaped makes the second graph more bowl shaped! Moreover it shifts the minimum of the 2nd graph towards the mimimum of the bowl. As an illustration, let’s take a look at an example where this is easy to see.

rm(list = ls())
bumpy.function = function(x,y)
{
    return(sin(x)+sin(y))
}

nice.paraboloid = function(x,y)
{
    return(0.15*(x^2 + y^2)) #lambda = 0.15 is used
}

x = y = seq(from = -4, to = 4, by = 0.2)
bumpy = outer(x, y, bumpy.function)
paraboloid = outer(x, y, nice.paraboloid)
bumpy.plus.paraboloid = bumpy + paraboloid

persp(x, y, bumpy,
      main="Graph of a bumpy function with multiple minima",
      theta = 30, phi = 15,
      col = "springgreen", shade = 0.5)

persp(x, y, paraboloid,
      main="Graph of a nice paraboloid",
      theta = 30, phi = 15,
      col = "springgreen", shade = 0.5)

persp(x, y, bumpy.plus.paraboloid,
      main="Graph of a regularized bumpy function = bumpy function + paraboloid",
      theta = 30, phi = 15,
      col = "springgreen", shade = 0.5)

Here we see that the most important geometric aspect of the regularizing term \(\lambda\langle\beta, \beta\rangle\) is the fact that it is strictly convex!! Although we will not dwell on it, it is impossible to overstate the theoretical importance of the previous sentence. As a matter of fact, geometrically speaking it’s clear that had we added any strictly convex function to the bumpy function we would have gotten something more bowl shaped. We will not go further into it here but you should know that convexity is one of those properties in mathematics out of which entire fields are created.

The Hessian matrix and more complex models

In the OLS problem above the Maximum Likelihood estimator turned out to be the one that minimized the sum of squares expression \(SSE(\beta) := (Y - X\beta)^T(Y - X\beta)\). This expression can be expanded in matrix notation as

\[ Y^TY - Y^TX\beta - \beta^TX^TY + \beta^TX^TX\beta \]

We see that there is a quadratic term (\(\beta^TX^TX\beta\)) and the rest are terms with powers of \(\beta\) less than 2. The matrix of second derivatives of this expression (known as the Hessian matrix) is therefore just the matrix of coefficients of this quadratic term: \(X^TX\). This Hessian is exactly what was the star of the show in the OLS problem!

The Hessian of a function at a point tells us the convexity of the function at the point. If the Hessian is positive definite, then near the minimizing point the function is bowl shaped. If the Hessian is negative definite then near a maximizing point the function is shaped like an upside down bowl.

Moreover, The Hessian is always a symmetric matrix by the equality of cross-derivatives. So the previous point is really a statement about Hessian’s eigenvalues. In short, if the Hessian has all positive eigenvalues then it is positive definite and the function is bowl shaped near its minimizer.

The OLS Hessian matrix \(X^TX\) is symmetric and non-negative definite, so the graph of the sum of squares \(SSE(\beta) := (Y - X\beta)^T(Y - X\beta)\) can only fail to be bowl shaped near the minimizer \(\hat{\beta}\) if the matrix \(X^TX\) has eigenvalues that are equal 0 (or positive but close to 0 in the case numerical instability). In which case, the graph of \(SSE(\beta)\) is a degenerate non-hyperbolic paraboloid and there are multiple minimizing solutions \(\beta\).

Thus regularizing the matrix \(X^TX\) is just regularizing the Hessian matrix of the function \(SSE(\beta)\) we want to minimize!! This fact is what allows us to take the idea beyond OLS.

Indeed if \(f(\beta)\) is any 2nd order differentiable cost function in any machine learning model then by the linearity of the derivative

\[ \text{Hessian}(f + \lambda \langle \beta, \beta\rangle) = \text{Hessian}(f) + \lambda\text{Hessian}(\langle \beta, \beta\rangle) = \text{Hessian}(f) + \lambda I \]

There we used the fact that \(\text{Hessian}(\langle \beta, \beta\rangle) = \text{Hessian}(\sum_i^p \beta_i^2) = I\). Because \(f\) is almost arbitrary we see that we can apply \(L^2\)-regularization to a very large family of problems, with the goal being to regularize the Hessian of \(f\). As an example let’s look at some other kinds of regression problems. For OLS we assumed the conditional distribution \(Y \ \ | \ \ \vec{X} \sim \mathcal{N}(\ \langle \vec{X} , \beta\rangle \ ,\ \sigma^2)\), but we may have chosen a different conditional distribution.

If \(Y\) takes on only values in the set \(\{0,1\}\) then it is a Bernoulli random variable. This is the case in logistic regression where the conditional distribution of \(Y\) given \(\vec{X}\) is \[ Y \ \ | \ \ \vec{X} \sim \mathcal{B}(\ p = \phi(\langle \vec{X} , \beta\rangle) \ ) \]

where \(\mathcal{B(p)}\) is a Bernoulli distribution with probability of a positive event equal to \(p\), \(p = \text{E}[Y|\vec{X}] = \phi(\langle \vec{X} , \beta\rangle)\) is the probability of a positive event given \(\vec{X}\), and \(\phi(t) = \frac{e^t}{1+e^t}\) is the standard logit. In this case the conditional density \(p(Y|\vec{X})\) can be written as

\[ p(Y|\vec{X}) = p^Y\big(1 - p)\big)^{1-Y} =\phi(\langle \vec{X} , \beta\rangle)^Y\big(1 - \phi(\langle \vec{X} , \beta\rangle)\big)^{1-Y} \] So the negative log-likelihood is given by

\[ -\mathcal{l}(\beta) = -\sum_{i=1}^N Y_i\log(\phi(\langle \vec{X}_i , \beta\rangle)) + (1-Y_i)\log(1-\phi(\langle \vec{X}_i , \beta\rangle)) \]

This function may or may not look bowl shaped (i.e. strictly convex) near the \(\hat{\beta}\) that minimizes it. In case it doesn’t we can make it so by adding \(\lambda \langle \beta, \beta\rangle = \lambda\sum_{j = 1}^p\beta_j^2\) for some sufficiently large \(\lambda\) and minimizing this new problem. The same applies to generalized linear models, neural networks, etc.