7
$\begingroup$

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

$\endgroup$
  • $\begingroup$ This may be obvious to you, but just because you have a decomposition of $A$ doesn't mean that you can accurately compute its determinant. In general, computing the determinant can not be done accurately for large matrices (and neither is it typically necessary). $\endgroup$ – Wolfgang Bangerth Jan 28 '14 at 4:05
  • $\begingroup$ This is a very interesting point. Could you suggest an accurate method for large matrices? $\endgroup$ – yannick Jan 28 '14 at 10:15
11
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This answer is great because it is general. I will probably accept it because my question was too general also. For my particular case, I rephrased the problem in scicomp.stackexchange.com/q/10638/733 which I would appreciate you have a look at. $\endgroup$ – yannick Jan 28 '14 at 10:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.