2

I need to invert a p x p symmetric banded hessian matrix H, which has 7 diagonals. p may be very high (=1000 or 10000).

H^{-1} can be considered as banded, and thus, I do not need to compute the complete inverse matrix, but rather its approximation. (It could be assumed to have 11 or 13 diagonals for example.) I am looking for a method which does not imply parallelization.

Is there any possibility to build such an algorithm with R, in linear time? 

double-beep
  • 5,031
  • 17
  • 33
  • 41
0spirit0
  • 83
  • 4
  • If something like this is available, I expect it in package Matrix. E.g., there: https://stat.ethz.ch/R-manual/R-devel/library/Matrix/html/solve-methods.html. Of course, matrix inversion should usually be avoided. – Roland Aug 10 '16 at 12:21

1 Answers1

1

There is no linear time algorithm for this, to the best of my knowledge. But you're not totally without hope:

  1. Your matrix is not really that large, so using a relatively optimised implementation might be reasonably fast for p < 10K. For example a dense LU decomposition requires at most O(p^3), with p = 1000, that would probably take less than a second. In practice, an implementation for sparse matrices should achieve much better performance, taking advantage of the sparsity;
  2. Do you really, really, really need to compute the inverse? Very often explicitly computing the inverse can be replaced by solving an equivalent linear system; with some methods such as iterative solvers (e.g. conjugate gradient) solving the linear system is significantly more efficient because the sparsity pattern of the source matrix is preserved, leading to reduced work; when computing the inverse, even if you do know it's OK to approximate with a banded matrix, there will still be substantial amount of fill-in (added nonzero values)

Putting it all together I'd suggest you try out the R matrix package for your matrix. Try all available signatures and ensure you have a high performance BLAS implementation installed. Also try to rewrite your call to compute the inverse:

# e.g. rewrite...
A_inverse = solve(A)
x = y * A_inverse
# ... as
x = solve(A, y)

This may be more subtle for your purposes, but there is a very high chance you should be able to do it, as suggested in the package docs: 

 solve(a, b, ...) ## *the* two-argument version, almost always preferred to
 solve(a)         ## the *rarely* needed one-argument version

If all else fails, you may have to try more efficient implementations available in: Matlab, Suite Sparse, PetSC, Eigen or Intel MKL.

paul-g
  • 3,797
  • 2
  • 21
  • 36