Most efficient way to solve a system of linear equations

c++ matrix linear-algebra eigen
aesir picture aesir · Nov 24, 2013 · Viewed 29.1k times · Source

I have an (n x n) symmetric matrix A and a (n x 1) vector B. Basically, I just need to solve Ax = b for x. The issue is that A is going to potentially be massive. So I'm looking for the most efficient algorithm for solving linear equations in C++. I looked over the Eigen library. Apparently it has an SVD method, but I've been told it's slow. Solving x=inverse(A)*b also seems like it would be suboptimal. Is uBLAS faster? Are there any more efficient methods? Thanks.

Edit: matrix A is positive definite and not sparse.

Answer

Pavan Yalamanchili picture Pavan Yalamanchili · Nov 25, 2013

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:

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

Related questions

Eigen efficient type for dense symmetric matrix

Does Eigen have efficient type for store dense, fixed-size, symmetric matrix? (hey, they are ubiquitous!) I.e. for N=9, it should store only (1+9)*9/2==45 elements and it has appropriate operations. For instance there should be efficient addition of two symmetric matrices, …

Read More
What are the most widely used C++ vector/matrix math/linear algebra libraries, and their cost and benefit tradeoffs?

It seems that many projects slowly come upon a need to do matrix math, and fall into the trap of first building some vector classes and slowly adding in functionality until they get caught building a half-assed custom linear algebra …

Read More
Is there around a straightforward way to invert a triangular (upper or lower) matrix?

I'm trying to implement some basic linear algebra operations and one of these operations is the inversion of a triangular (upper and/or lower) matrix. Is there an easy and stable algorithm to do that? Thank you.

Read More