2

I am interested in solving a large number of small linear systems of equations, $Ax=b$, with $A$ either $2\times2$ or $3\times3$. Assuming none of these systems are actually singular, is there anything to deter me from forming $A^{-1}$ explicitly using Cramer's rule and computing $x=A^{-1}b$? The only things that comes to mind is possible loss of precision in the calculation of $\mathrm{det}A$. Can I expect any improvement over this naive approach by using a canned solver, e.g., LAPACK's dgesv? Namely, would doing so reap numerical benefits that justify the overhead of calling out to LAPACK?

(I am open to methodologies in between these two extremes as well, e.g., hard-coding Gaussian elimination).

Endulum
  • 735
  • 3
  • 9
  • 2
    Did you check the answer to this question? – nicoguaro Oct 16 '18 at 17:18
  • Yes, I did. One answer (GertVdE's) does not actually address numerics. The other (Geoff's) does, but only in a cursory, qualitative way. I am hoping for more detail. – Endulum Oct 16 '18 at 19:16
  • mind sharing a bit more detail about your problem? if you're doing e.g. computer vision and use GPUs, here's a thought: group your many matrices into batches and use specialized code to process them all at once. this has been done. a couple of references to get the ball rolling: Lemaitre et al. - batch Cholesky, Dong - batch LU – GoHokies Oct 16 '18 at 19:42
  • Intel's MKL (v2018) supports batch operations through vectorized compact routines, see here. another optimization technique (again, MKL-specific) is the use of MKL_DIRECT_CALL's. AFAIK, what this does is to disable error checking. The Intel folks report decent DGEMM speedups over the standard (non-direct?) MKL lib functions. YMMV. – GoHokies Oct 16 '18 at 19:56
  • 1
    anyway, to quote from here: "libraries compiled for specific sizes, either via templating or just-in-time (JIT) compilation, provide a large performance boost for small (matrix) sizes; for instance, Intel’s libxsmm3 for extra-small matrix-multiply or batched BLAS for sets of small matrices. Thus, users with such small matrices are encouraged to use special purpose interfaces rather than try to optimize overheads in a general purpose interface." – GoHokies Oct 16 '18 at 20:03

1 Answers1

2

Cramer's rule using a more stable determinant:

https://hal.archives-ouvertes.fr/hal-01500199/document

Implementation of the algorithm is consistently 2.5 times slower than LU factorization of matlab regardless of size of matrix.

R zu
  • 163
  • 9
  • 2
    To be precise, according to the data in the paper you linked, their Cramer's rule is faster than Matlab LU (whatever it's calling internally) for $N\lessapprox180$. However, I would be much more interested in comparison with an efficient LAPACK exactly for the matrix sizes of interest. – Anton Menshov Oct 16 '18 at 17:53