Why the Ridge Regression is NOT scale-invariant?

Question

In the Element of Statistical Learning, Chapter 3, we know that the linear regression is scale-invariant since the scale matrix for coefficient will be canceled eventually, but the Ridge regression doesn't have it? Since the form of Ridge coefficient has the closed-form $$ \beta = (X^{T}X + \lambda I)^{-1}X^{T}Y, $$ I don't see why the scale-invariance doesn't hold in here? Can anyone suggest a prove of it?

Answer 1

The intuition here is that there's a sleight-of-hand happening when you use the same symbol $X$ for both the original data and the rescaled data. It's misleading because the rescaling $\tilde{X}= XD$ is not the same as the original $X$, so we should make that explicit and write down how we're rescaling.

We can demonstrate this by considering two cases, first with the original units in $X$ and second the case where we use a rescaled matrix $\tilde{X}= XD$ where $D$ is a diagonal matrix that has all positive entries on the diagonal. If $X$ has shape $n \times p$ then $D$ has shape $p \times p$. (You can actually use any $D_{ii} \neq 0$ but "rescaling" is almost always meant to be restricted to multiplication by a positive scalar.)

In the first case, we have $$\beta(X) = (X^TX + \lambda I)^{-1}X^T y$$ which is just as written in the question.

In the second case, we apply the rescaling to $X$ and we have $$\begin{aligned} \beta(\tilde{X}) &= (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T y\\ &= (DX^TXD + \lambda I)^{-1}D X^Ty \\ &= (D(X^\top X + \lambda D^{-2})D)^{-1}DX^Ty \\ &= D^{-1}(X^T X + \lambda D^{-2})^{-1}X^Ty \end{aligned}$$

(remembering that $D$ is diagonal, so $D^T = D$).

From this we can conclude that the coefficients $\beta_X$ and $\beta_\tilde{X}$ are only the same if $D=I$.

The final line shows that the rescaling two effects on the coefficients.

1. It has a multiplicative effect on the coefficients, just as we would intuitively expect based on what happens when we rescale in the OLS case.
2. The last line makes explicit that the change in scale is "absorbed" in $\lambda$, and that the change in scale is gives $\beta(\tilde{X})_i$ penalized inversely to the square of the rescaling $D_{ii}$. (Thanks to Firebug for this helpful suggestion.)

Answer 2

Write it in terms of the cost function

$$ (y - X\beta)^2 + \sum_i \lambda \beta_i^2 $$

As you can see, each of the model parameters $\beta_i$ has the same penalty $\lambda$ (the $\lambda I$ part). If we want it to have the same degree of penalization for each parameter, we need them to have the same scale. This can be achieved either by scaling the data or by using different values of $\lambda$ per parameter, but scaling $\lambda$ is equivalent of scaling the data.

Answer 3

Note that being tractable doesn't mean not scale invariant. For example, PCA is tractable but not scale invariant.

Now let's look at the solution.

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

We have just one penalisation parameter, $\lambda$ which is added to $X^T X$ to form the penalty. Regardless of the units of each variable (columns of $X$) we have the same penalty. If I had $x_1$ in meters but then converted to km with the same penalty I wouldn't achieve the same fit to the data.

Lets consider a simple numerical example, two variables with some R code. The main thing to consider is that, under OLS, if I multiply $x_1$ by a factor $k$ then the new $\beta_1$ would be divided by $k$. In this numerical example I multiply $x_1$ by $100$, but we can see that the resulting $\beta_1$ is only divided by $90$. However, under OLS I recover this scaling factor of $100$.

set.seed(134221)
x1 <- runif(10)
x2 <- runif(10)
eps <- rnorm(10)*0.1
y <- 2 - 0.3*x1 + 0.9*x2 + eps
## fix lambda = 0.1
X0 <- cbind(1, x1, x2)
b0 <- solve(t(X0)%*%X0 + diag(0.1, 3))%*%t(X0)%*%y
b0
## now change x1 <- 100*x1
X1 <- cbind(1, 100*x1, x2)
## now change x1 <- 100*x1
b1 <- solve(t(X1)%*%X1 + diag(0.1, 3))%*%t(X1)%*%y
b1
## under OLS if we have x1 := k*x1
## then beta1 := beta1/k
b0[2]/b1[2]
[1] 90.42805

## OLS regression
fit0 <- summary(lm(y ~ X0))
fit1 <- summary(lm(y ~ X1))
fit0$coefficients[2]/fit1$coefficients[2]
[1] 100