full rank update to cholesky decomposition

Question

Let $A$ be a real, symmetric, positive definite matrix. It has at least 500 rows, possibly much more. I compute its Cholesky decomposition, which allows me to calculate

• $det(A)$
• $A^{-1}X$ for some vector/matrix $X$
• $A^{-1}$ (I am aware of the warnings here)

Unfortunately, this operation needs to be performed a lot of times (say, a million), on matrices which are closely related but different from $A$. Formally, $A$ and any other similar matrix $A'$ differ by a full-rank update $E$ such that

$$A' = A+E \quad \frac{\|E\|}{\|A\|}\ll 1 \quad rank(E)=dim(E)$$

For now, I am computing the Cholesky decomposition of every matrix at every step. In order to speed up my calculation, I would like to use the information that two consecutive matrices are closely related. How can I do that?

I am looking for something $O(N^2)$ instead of the $O(N^3)$ of the Cholesky decomposition. Currently the matrices are dense, but I am also interested in solutions for banded sparse matrices.

EDIT Since this question was general and gave a great general answer, I will mark it as accepted. However, as far as my more specific problem is concerned, see full rank update to cholesky decomposition for multivariate normal distribution

Answer 1

In general, there is no shortcut other than completely re-factoring the matrix. There have been a few similar questions on this SE that cover the topic in more depth than I can:

Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation?

LU Decom of PSD Matrix + Diagonal Matrix

Perturbation of Cholesky decomposition for matrix inversion

Diagonal update of a symmetric positive definite matrix

Computational complexity and implementation of UDU Modified Cholesky Rank 1 Update

Long story short, you can do it in $O(\textrm{rank}(E)\cdot n^2)$. If $E$ is rank-1 or of some rank asymptotically lower than the dimension of the matrix you're trying to update, then you can do better than $n^3$, otherwise you're out of luck.

For banded matrices, things are a little nicer; you can do a Cholesky factorization in $O(\textrm{bandwidth}^2\cdot n)$, so provided the bandwidth is appreciably smaller than $n$, the cost of doing a single factorization isn't so punitive.

Is your implementation serial or parallel? If it's serial, you may want to investigate using a package like Elemental to parallelize the decomposition of each individual matrix. Sometimes, the matrices are too small to get much advantage out of parallelism, but if you have 1000-dimensional dense matrices I suspect that won't be the case.

You've said that the matrices are related; are you given a whole mess of matrices $\{A_i\}$ right at the start of the program, or do you update them sequentially? If it's the former, you may be able to use MPI to get several processors working on the matrices in parallel.