📂Matrix Algebra

A Trick for Reducing the Computational Cost of Tikhonov Regularization Using Singular Value Decomposition

Theorem

Optimal solution of ridge regression: $$ L \left( \beta \right) = \left\| Y - X \beta \right\|_{2}^{2} + \lambda \left\| \beta \right\|_{2}^{2} $$ When $\lambda$ is given as a constant, let the objective function $L$ of ridge regression be expressed as above. The optimal solution of ridge regression $\hat{\beta} = \argmin_{\beta} L \left( \beta \right)$ is as follows. $$ \hat{\beta} = \left( X^{\top} X + \lambda I \right)^{-1} X^{\top} Y $$ Here, $A^{\top}$ is the transpose of $A$, $I$ is the identity matrix, and $A^{-1}$ is the inverse matrix of $A$.

Given a matrix $A \in \mathbb{R}^{m \times n}$ and a vector $\mathbf{b} \in \mathbb{R}^{m}$, when finding the least squares solution of $A \mathbf{x} = \mathbf{b}$, applying Tikhonov regularization for $\lambda \ge 0$ leads to solving the following optimization problem. $$ \argmin_{\mathbf{x}} \left\| A \mathbf{x} - \mathbf{b} \right\|^2 + \lambda \left\| \mathbf{x} \right\|^2 $$ The solution to this problem is usually known in the closed form $\mathbf{x} = \left( A^{\top} A + \lambda I \right)^{-1} A^{\top} \mathbf{b}$, but using the singular value decomposition $A = U \Sigma V^{\top}$ of $A$, it can also be expressed by the following formula, which allows one to perform the matrix decomposition only once and compute while varying $\lambda$. $$ \mathbf{x}(\lambda) = V \left( \Sigma^{2} + \lambda I \right)^{-1} \Sigma U^{\top} \mathbf{b} $$

Explanation

If we compare the speed and compare the stability of implementing least squares via various matrix decompositions, we see a very clear trade-off as follows.

Method	Speed	Stability
LU	Fast	Unstable
QR	Medium	Stable
SVD	Slow	Very stable

In reality, looking only at the computational cost of QR and SVD themselves, SVD is slower, but for example when the matrix is very large and, for ridge regression, we repeat the computation while varying $\lambda$ to find the optimal value, the trick using SVD can reduce the computational cost. In the presented formula, $\left( \Sigma^{2} + \lambda I \right)$ is a diagonal matrix, so computing its inverse matrix is very easy, and $V$ and $U^{\top} \mathbf{b}$ can be computed only once and reused for all $\lambda$.

This trick not only reduces the computational cost but also secures the stability provided by SVD, so it can be a useful method the closer one gets to data of enormous size and high risk of ill-conditioning.

Of course, when the problem mentioned above does not exist, more often than not one cannot beat QR decomposition, or there is no need to. For instance, just because a task that took 1 second is reduced to 0.5 seconds, going out of one’s way to implement SVD-based ridge regression is excessive, and at the point where one writes unverified code oneself, the completeness of the implementation also becomes suspect. One must reasonably judge whether the situation calls for applying SVD, and even if the implementation succeeds, verification of whether there was actually a performance improvement is essential.

Proof

Let $A = U \Sigma V^{\top}$ be the singular value decomposition (SVD) of $A$. Here, $U$ and $V$ are orthogonal matrices, and $\Sigma$ is a diagonal matrix containing the singular values. Since $V$ is an orthogonal matrix, substituting $\mathbf{y} = V^{\top} \mathbf{x}$ gives the following. $$ \lambda \left\| \mathbf{x} \right\|^2 = \lambda \left\| V \mathbf{y} \right\|^2 = \lambda \mathbf{y}^{\top} V^{\top} V \mathbf{y} = \mathbf{y}^{\top} \lambda I \mathbf{y} $$ And $$ \begin{align*} \left\| A \mathbf{x} - \mathbf{b} \right\|^2 =& \left\| U \Sigma V^{\top} V \mathbf{y} - \mathbf{b} \right\|^2 \\ =& \left\| U \Sigma \mathbf{y} - \mathbf{b} \right\|^2 \\ =& \left( U \Sigma \mathbf{y} - \mathbf{b} \right)^{\top} \left( U \Sigma \mathbf{y} - \mathbf{b} \right) \\ =& \mathbf{y}^{\top} \Sigma^{\top} U^{\top} U \Sigma \mathbf{y} - 2 \mathbf{y}^{\top} \Sigma U^{\top} \mathbf{b} + \left\| \mathbf{b} \right\|^{2} \\ =& \mathbf{y}^{\top} \Sigma^{2} \mathbf{y} - 2 \mathbf{y}^{\top} \Sigma U^{\top} \mathbf{b} + \left\| \mathbf{b} \right\|^{2} \end{align*} $$ Therefore, the objective function with respect to $\mathbf{y}$ and its gradient can be expressed as follows. $$ \begin{align*} J \left( \mathbf{y} \right) =& \mathbf{y}^{\top} \left( \Sigma^{2} + \lambda I \right) \mathbf{y} - 2 \mathbf{y}^{\top} \Sigma U^{\top} \mathbf{b} + \left\| \mathbf{b} \right\|^{2} \\ \partial J \left( \mathbf{y} \right) / \partial \mathbf{y} =& 2 \left( \Sigma^{2} + \lambda I \right) \mathbf{y} - 2 \Sigma U^{\top} \mathbf{b} . \end{align*} $$ Finally, the optimal solution $\mathbf{x}(\lambda)$ for a given $\lambda$ is as follows. $$ \mathbf{x}(\lambda) = V \left( \Sigma^{2} + \lambda I \right)^{-1} \Sigma U^{\top} \mathbf{b} $$

■

Code

The following is the SVD-based Tikhonov regularization that the author actually uses in research.

using LinearAlgebra

function tikhonov_λ(A, b)
    U, s, V = svd(A)
    Uᵀb = U'b
    return λ -> V * (Diagonal(s ./ (s.^2 .+ λ)) * Uᵀb)
end

In the author’s case, there was a mistake of forgetting Uᵀb = U'b and plugging it directly into the anonymous function, which caused $U^{\top} \mathbf{b}$ to be computed every time $\lambda$ was changed, so that the computational cost fell behind the basic implementation (QR decomposition). Do not just trust the formulas; be sure to run an actual benchmark.