Clever ways to update LU factorization for ridge regression

Question

Ridge regression can be posed as minimizing the following objective function (over $x$):

$$\frac{1}{2} \lVert Ax - b \lVert_2^2 ~+ \frac{\lambda}{2} \lVert x \lVert_2^2 $$

Which has a closed form solution:

$$x = (A^TA + \lambda I)^{-1} A^T b $$

Let's say we want to solve this problem for several different values of $b$. Then its a good idea to cache the LU decomposition of $A^TA + \lambda I$. In MATLAB I this translates to:

[L,U] = factor(A'*A + lambda*I)
x1 = L \ U \ A'*b1
x2 = L \ U \ A'*b2
% etc...

Here's my question: if $\lambda$ is updated after several solves, is there a simple/efficient way of updating $L$ and $U$? Or do I have to completely redo the LU factorization?

Answer 1

What is the size of your $A$ matrix? Is $A$ sparse? Does $A$ have some other special structure? How many values of $\lambda$ do you want to try?

Normally, you'd use the Cholesky factorization of $A^{T}A+\lambda I$ rather than the LU factorization since $A^{T}A+\lambda I$ is symmetric and positive definite. However, updating the Cholesky or LU factorization of a matrix after a full rank diagonal update is not easy in practice (low rank updates are a different matter.)

One alternative is to do the computation of the regularized solutions using the SVD of $A$. There's a simple formula for the Tikhonov regularized solution in terms of the SVD so that you only have to compute the SVD of $A$ once. This is a reasonable way to go for small sized problems (e.g. $A$ has no more than a few thousand rows/columns.) You could also use the eigenvalue decomposition of $A^{T}A$ to do this instead of the SVD of $A$.

For large problems, especially where $A$ is sparse, iterative methods are typically used. You can solve the problem with explicit regularization as a linear least squares problem (hint: work from largest to smallest values of the parameter $\lambda$ and use the previous solution as the starting point for each new value of $\lambda$.) A more common approach is to use an iteration such as the Landweber iteration or CGLS that will give you an implicitly regularized solution when you stop the iterations before full convergence.