Perturbation of Cholesky decomposition for matrix inversion

Question

I am looking for a computationally cheap way to compute $x$ such that $$(L L^T + \mu^2 I)x = y$$ where $L \in \mathbb{R}^{n \times n}$ is a lower triangular definite positive matrix (with some very small eigenvalues), $y \in \mathbb{R}^n$ and $\mu \in \mathbb{R}$ are known. If it is necessary, I can assume that $$\mu \ll 1$$ But $\mu^2$ is larger than the smallest eigenvalue of $LL^T$.

Basically, I would like to make the most of my knowledge of the Cholesky decomposition $L L^T$. Eventually, I hope to be able to compute $x$ in $\mathcal{O}(n^2)$. Approximate approaches are also welcomed.

I have seen here that this does not seem to be doable in a more general situation, but I hope the smallness of $\mu$ may help...

Any idea, reference or warning? Thanks for your help.

Answer 1

Unfortunately, although it is easy to update the Cholesky factorization after a low-rank update, updating the Cholesky factorization after a full rank update like this isn't any easier than recomputing the Cholesky factorization.

An alternative to consider would be to compute the eigenvalue decomposition of your matrix,

$A=LL^{T}=UDU^{T}$.

Adding $\mu^{2}I$ to $A$ simply adds $\mu^{2}$ to the eigenvalues of $A$, so

$A+\mu^{2}I=UVU^{T}$

where $V=D+\mu^{2}I$. Then, to solve $(A+\mu^{2}I)x=y$, we let

$UVU^{T}x=y$

or

$x=UV^{-1}U^{T}y$.

For each value of $\mu$ you're interested in, this takes only $O(n^{2})$ additional work (three matrix-vector multiplications and one of them is just a diagonal scaling matrix. Be careful not to multiply the matrices!) In comparison with computing the Cholesky factorization of $A$, computing the eigenvalue decomposition is also $O(n^{3})$, but it's typically an order of magnitude slower than the Cholesky factorization. On an amortized basis this will pay for itself if you have enough (say dozens to hundreds) of values of $\mu$ to consider.

Answer 2

You can rephrase your question as an optimization problem and solve numerically with gradient descent. You would probably converge quite fast, because your initial guess will be the solution you already have for $\mu^2=0$.

More formally, we like to find $x$ that minimizes:

$f(x) = ||(LL^T+μ^2I)x -y||^2$, where $L, μ^2, y$ are given.

To solve it, start with $f(x_0) = f(A^{-1}y)$ and update $f(x)$ with gradient descent steps until convergence. Each evaluation of the gradient should cost you $O(d^2)$.

Since your initial guess is good, other iterative approaches will probably work well. See for example chapter 2, here