11
$\begingroup$

Given a positive definite symmetric matrix, what is the fastest algorithm for computing the inverse matrix and its determinant? For problems I am interested in, the matrix dimension is 30 or less.

  1. High accuracy and speed is really necessary. (millions matrices are performed)
  2. The determinant is necessary.In each calculation, only one element of the iverse matrix is required. Thanks!
| cite | improve this question | |
$\endgroup$
  • $\begingroup$ Do you have to invert millions of such matrices? Otherwise, speed should not be an issue. $\endgroup$ – Wolfgang Bangerth May 12 '15 at 12:37
  • $\begingroup$ I edited your title and question for clarity. If I made any errors, please let me know. $\endgroup$ – Geoff Oxberry May 13 '15 at 0:11
  • $\begingroup$ @Wolfgang Bangerth Yes, speed should be considered. $\endgroup$ – Orders Jun 11 '16 at 0:44
  • 2
    $\begingroup$ Do you know which element of the inverse matrix is needed? Or can it be a random entry? $\endgroup$ – Memming Jun 11 '16 at 1:43
  • 2
    $\begingroup$ @Orders Your comment and edit seem contradicting: do you need one element of the inverse, or all of them? $\endgroup$ – Federico Poloni Jun 11 '16 at 11:14
13
$\begingroup$

For problems I am interested in, the matrix dimension is 30 or less.

As WolfgangBangerth notes, unless you have a large number of these matrices (millions, billions), performance of matrix inversion typically isn't an issue.

Given a positive definite symmetric matrix, what is the fastest algorithm for computing the inverse matrix and its determinant?

If speed is an issue, you should answer the following questions:

  • Do you really need the whole inverse? (Many applications don't need to form an explicit inverse.)
  • Do you really need the determinant? (Determinants are uncommon, but certainly not unheard of in computational science.)
  • Do you need either to high accuracy? (Low accuracy algorithms tend to be faster.)
  • Would a probabilistic approximation suffice? (Probabilistic algorithms tend to be faster.)

The standard response to your problem of inverting a small, positive definite matrix and calculating its determinant would be Cholesky decomposition. If $A = LL^{T}$, then $\det(A) = \prod_{i=1}^{n}l_{ii}^{2}$, and $\det(A^{-1}) = \prod_{i=1}^{n}l_{ii}^{-2}$.

Assuming $A$ is $n$ by $n$, the Cholesky decomposition can be computed in around $n^{3}/3$ flops, which is about half the cost of an LU decomposition. However, such an algorithm would not be considered "fast". A randomized LU decomposition might be a faster algorithm worth considering if (1) you really do have to factor a large number of matrices, (2) the factorization is really the limiting step in your application, and (3) any error incurred in using a randomized algorithm is acceptable. Your matrices are probably too small for sparse algorithms to be worthwhile, so the only other opportunities for faster algorithms would require additional matrix structure (e.g., banded), or exploiting problem structure (e.g., maybe you can cleverly restructure your algorithm so that you no longer need to calculate a matrix inverse or its determinant). Efficient determinant algorithms are roughly the cost of solving a linear system, to within a constant factor, so the same arguments used for linear systems apply to calculating determinants as well.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Just a brief note: if $B = A^{-1}$, to compute a single element $b_{ij}$ one should compute only the $j$th column of $B$. Once the Cholesky factorisation is computed, this is done by forward and back substitution with respect to a rhs vector of all zeros and with just a single one in row $j$. Since the computation can be interrupted as soon as $b_{ij}$ has been computed, best case is for $b_{nn} = l_{nn} ^{-2}$ worst case for $b_{11}$ where one has to compute full back and forward substitutions. $\endgroup$ – Stefano M May 13 '15 at 20:40
  • 1
    $\begingroup$ @StefanoM Even better, you can permute your matrix before the beginning of the computation so that you are always in the best case. $\endgroup$ – Federico Poloni Jun 11 '16 at 11:15
0
$\begingroup$

Briefly, referencing the Julia documentation on linear algebra subroutines, they note that the Bunch-Kaufman factorization method is more appropriate for symmetric matrices.(old source from NASA) It may go without saying that positive definite matrices are a subset of symmetric matrices, so while Bunch-Kaufman factorization is an improvement, it isn't necessarily the best they can do.

A 2015 matlab user submission entitled "Fast and Accurate Symmetric Positive Definite Matrix Inverse Using Cholesky Decomposition" clearly suggests the Cholesky decomposition, and the RFast package shares that opinion, but another stack exchange conversation suggests that the best method is really application dependent - so if you have the time to benchmark various methods, that would give you a more definitive answer.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Bunch-Kaufman is a pivoting strategy for the LDL^T decomposition, not a separate factorization. PD matrices are a much simpler case though (there are never off-diagonal pivots), so this strategy it can be greatly simplified: replaced by just a max diagonal entry, or even omitted altogether. $\endgroup$ – Federico Poloni Jul 29 at 6:37

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.