Most efficient way to solve a system of linear equations

The best way to solve a system of linear equations of the form Ax = b is to do the following.

• decompose A into the format A = M1 * M2 (where M1 and M2 are triangular)
• Solve M1 * y = b for y using back substitution
• Solve M2 * x = y for x using back substitution

For square matrices, step 1 would use LU Decomposition.

For non square matrices, step 1 would use QR Decomposition.

If matrix A is positive definite and not sparse you’d use Cholesky Decomposition for the first step.

If you want to use eigen, you will have to first decompose it and then triangular solve it.

If this is still slow, thankfully, there are numerous linear algebra libraries available that can help improve the performance. The routine you should be looking for is dpotrs. Some libraries that have this implemented are as follows:

• Netlib’s lapack: Documentation and Download (free)
• Intel’s MKL: Documentation and Download (Free for non commercial use).
• AMD’s ACML: Download (free)
• PLASMA: Download (free, multi core optimized)
• MAGMA : Download (free, Implemented in CUDA, OpenCL)
• CULA: Download (freemium, Implemented in CUDA).

If you are using eigen in the overall project, you can interface the LAPACK routine you need as described here.