# Regularized least squares with QR factorization

Consider the regularized least squares problem $$\min_x || b - A x ||^2 + \lambda^2 ||x||^2$$ which is equivalent to $$\min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{A \\ \lambda I} x \right|\right|^2$$ Suppose I have factored $$A$$ with a $$QR$$ factorization, so that $$A = QR$$ and therefore \begin{align} \min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{A \\ \lambda I} x \right|\right|^2 &= \min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{QR \\ \lambda I} x \right|\right|^2 \\ &= \min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{Q & 0 \\ 0 & I} \pmatrix{R \\ \lambda I} x \right|\right|^2 \\ &= \min_x \left|\left| \pmatrix{Q^T b \\ 0} - \pmatrix{R \\ \lambda I} x \right|\right|^2 \end{align} Now, I want to solve this regularized problem for many different $$\lambda$$ values, so my question is how to do this efficiently once I have already computed $$R$$ and $$Q^T b$$.

My thinking so far is to do another $$QR$$ factorization on the matrix $$\pmatrix{R \\ \lambda I} = W S$$ The matrix is almost in upper triangular form except for the $$\lambda I$$ term, so I initially thought one could use Givens rotations to zero out the $$\lambda I$$ part, but it seems this will not work. For example, in the $$N = 2$$ case, we would have $$\pmatrix{ R_{11} & R_{12} \\ 0 & R_{22} \\ \lambda & 0 \\ 0 & \lambda}$$ A first Givens transformation would zero out the $$\lambda$$ in the first column, however it would then destroy the sparse structure of the $$\lambda I$$ block. It would produce something like this: $$\pmatrix{ \tilde{R}_{11} & \tilde{R}_{12} \\ 0 & R_{22} \\ 0 & X \\ 0 & \lambda}$$ where the $$X$$ term used to be zero but isn't anymore.

So a Givens approach seems inefficient since it does not preserve sparse structure. This seems to indicate that an optimized Householder approach is better since Householder transformations do preserve sparse structure.

Just wondering if anyone has any clever insights before I try coding up a Householder approach.

EDIT:

The book "Numerical Methods in Matrix Computations" by Ake Bjorck has an Exercise 2.6.2 in which he claims it is possible to perform a QR factorization of the matrix $$\pmatrix{B \\ \mu D}$$ where $$B$$ is upper triangular and $$D$$ is diagonal, using Givens rotations in $$O(11 n)$$ flops. The only hint given (earlier in the chapter) is a reference to the following paper,

Eldén, L. (1984). An efficient algorithm for the regularization of ill-conditioned least squares problems with triangular Toeplitz matrix. SIAM journal on scientific and statistical computing, 5(1), 229-236.

This paper develops a Givens-based approach to QR factor the matrix $$\pmatrix{K \\ \mu L}$$ where both $$K$$ and $$L$$ are upper triangular Toeplitz matrices. His algorithm is a straightforward application of Givens rotations, which work in this case because both $$K$$ and $$L$$ have the same non-zero pattern. When applying a Givens transformation to vectors $$x$$ and $$y$$, the result $$(x',y')$$ has a non-zero pattern which is the union of the patterns of $$x$$ and $$y$$.

Elden's algorithm for upper triangular Toeplitz matrices is $$O(n^2)$$ but I don't see an immediate way to adapt his method to the case of an upper triangular matrix and diagonal matrix.

Unfortunately Bjorck's exercise 2.6.2 doesn't give any further hints.

• You're aware that this computation is easy using the SVD of A, right? – Brian Borchers Jul 7 '19 at 17:25
• Yes but I want to use a QR approach – vibe Jul 7 '19 at 18:55
• Have you tried using the method in the Elden paper and treating the diagonal matrix $D$as an upper triangular matrix with the same sparsity pattern as $B$? That would put it into the form assumed in the Elden paper. – Brian Borchers Jul 13 '19 at 4:35
• @Brian yes that would work but cost extra flops. Ake Bjorck sent me an email stating there is an error in his exercise 2.6.2 and his matrix $B$ is meant to be bidiagonal, instead of upper triangular. In this case, Givens transformations can be used to zero the $\mu D$ matrix, using 2 Givens rotations per row. So it seems one possible way forward is to bidiagonalize the $R$ factor after $QR$ has been computed. – vibe Jul 15 '19 at 19:46