Matrix inversion or Cholesky?

Question

I am developing an algorithm which solves Ax= b, where A and b are known.

There are two ways to do this x= A^-1 b or using Cholesky. I know the matrix will always be square and positive definite although the det(A) maybe zero. In those rare cases I can just ignore it. But from a computation point of view and efficiency, is creating an inverse matrix too inefficient?

Answer 1

For large matrices, yes the inverse is very inefficient. However, special properties, such as the matrix being lower triangular make the inverse much simpler to compute.

In numerical analysis, the most typical solution to Ax=b is LU factorization of A(LUx=b), then solving Ly = b for y and Ux = y for x.

For a more stable approach, consider using QR decomposition, where Q has the special property of Q^T*Q = I, so Rx = Q^Tb where R is upper triangular(only one back solve as opposed to a forward and back solve with LU).

Other special properties, such as the matrix being symmetric(Cholesky) or banded(Gauss) make the certain solvers better than the other.

As always be aware of floating point errors in your calculations.

Might I add that iterative solvers are also a popular way to solve systems. Conguate gradient method is the most used and works well for sparse matrices. Jacobi and Gauss-Seidel are good for matrices that are diagonally dominate and sparse.

Answer 2

In general, you always want to use a solver; the actual solver should run about as fast as multiplying by an inverse. Not only is computing an inverse matrix inefficient compared to doing a decomposition, using an inverse matrix has precision problems that a decompose/solver approach avoids.

If you have a symmetric matrix, a Cholesky decomposition is a reasonable choice. The closely-related LDL decomposition has comparable precision, while also avoiding the need for square roots.

If your matrix is not symmetric, you can't use Cholesky or LDL decompositions -- use the LU decomposition method instead.