How to find determinant of large matrix

Question

I found some C++ code for finding the determinant of matrix, for 4x4 to 8x8. It works ok, but my project needs matrices that are 18x18 or more, and the code is too slow. The code is recursive, but is recursion the right concept to deal with an 18x18 matrix? How else can I find the determinant?

Answer 1

I assume you're using the naive method of expanding Laplace's formula. If you want to gain in speed, you can decompose your matrix M using LU decomposition (into two lower- and upper-diagonal matrices) which you can achieve with a modified Gauss-Jordan elimination in 2*n^3/3 FLOPS (edit: you can also use Crout's method, which is widely implemented) and then calculating the determinant as:

det(M) = det(L) * det(U), which for triangular matrices is just the product of the entries in their diagonal.

This process will still be faster than O(n!).

Answer 2

The laplace expansion algorithm is the simplest and is a recursive method. Sample code is available at the bottom of http://www.euclideanspace.com/maths/algebra/matrix/resources/code/index.htm#determinant.

Answer 3

Hi

Well, not many of us working in the field would regard 18x18 as a large matrix and almost any technique you choose should be fast enough on any modern computer. Nor would many of us tackle matrix questions with recursive algorithms, much more likely to use iterative ones -- but that could be a reflection of the fact that a lot of people working on matrix problems are scientists and engineers not computer scientists.

I suggest you look at Numerical Recipes in C++. Not necessarily the best code you'll find, but it is a text for studying and learning from. For better codes, BOOST has a good reputation and there's always BLAS and things like the Intel Maths Kernel Library or the AMD Core Maths Library. I think all of these have implementations of determinant-finding routines which will tackle an 18x18 matrix very quickly.

Regards

Mark

Answer 4

Since I can't comment, I wish to add this: the Cholesky decomposition (or its variant, LDL^T, L a unit lower triangular matrix and D a diagonal matrix) can be used to verify if a symmetric matrix is positive/negative definite: if it is positive definite, the elements of D are all positive, and the Cholesky decomposition will finish successfully without taking the square root of a negative number. If the matrix is negative definite, the elements of D are all negative, the matrix itself will not have a Cholesky decomposition, but the negative of it would.

"Calculating the determinant of a triangular matrix is simple: multiply the diagonal elements, as the cofactors of the off-diagonal terms are 0. Using an LU decomposition further simplifies this, as L is a unit, lower triangular matrix, i.e. its diagonal elements are all 1, in most implementations. Therefor, you often only have to calculate the determinant of U."

• You forgot here to take into account that all practical implementations of Gaussian elimination make use of (partial) pivoting for extra numerical stability; so your description is incomplete; one counts the number of row swaps done during the decomposition phase, and after multiplying together all the diagonal elements of U, this product should be negated if the number of swaps is odd.

As for code, NR is not free; I suggest taking a look at LAPACK/CLAPACK/LAPACK++ @ http://www.netlib.org/ instead. For reference, I can do no better than point you to "Matrix Computations" by Golub and Van Loan.