Find inverse and determinant of a positive definite matrix

I need to find the inverse and the determinant of a positive definite matrix. The matrix typically has size 10000x10000.

What is the most efficient and reliable way to get the inverse? I need the inverse since it would be used numerous times in later calculations.

As of now, I am using cholesky to get the answer. Determinant is just the product of diagonals of the triangular matrix in this case.

Also, I found that inv(A) based on LU is slower and inaccurate. How come the default inverse in Matlab has so much error?

I created a test case to demonstrate both the speed and accuracy issue. I am using Matlab 2014a.

close all;
clear all;
clc;
format short;
rng(0,'twister');

size = 5000;
temp = rand(size);
A = temp*temp'; % A is positive definite

tic
B = inv(A);
fprintf('Time for inv(A)   = %.3f seconds.\n',toc);
err1 = max(max(abs(A*B - eye(size))));
err2 = max(max(abs(B*A - eye(size))));
err = (err1+err2)/2.0;
fprintf('Error for inv(A)  = %e\n',err);

tic
B = A\eye(size);
fprintf('Time for A\\eye    = %.3f seconds.\n',toc);
err1 = max(max(abs(A*B - eye(size))));
err2 = max(max(abs(B*A - eye(size))));
err = (err1+err2)/2.0;
fprintf('Error for A\\eye   = %e\n',err);

tic
L = chol(A,'lower');
U = L\eye(size);
B = L'\U;
fprintf('Time for chol(A)  = %.3f seconds.\n',toc);
err1 = max(max(abs(A*B - eye(size))));
err2 = max(max(abs(B*A - eye(size))));
err = (err1+err2)/2.0;
fprintf('Error for chol(A) = %e\n',err);

Answer 1

[John D'Errico]

Many questions. How come the default inverse has so much error? Etc...

Answer: The reason why you were given this homework assignment. (Assuming it is homework.) Double precision arithmetic has limitations. The inverse of a matrix can be a poorly defined thing. The determinant of a matrix is almost completely useless in floating point arithmetic. NEVER use a determinant to determine if a matrix is singular. I hope you will learn why after this assignment. By the way...

det(A)
ans =
   Inf

So expect the determinant to be garbage, because it will overflow the dynamic range of double precision arithmetic.

The best way to compute an inverse? Inv is entirely sufficient. If there was a hugely better way, they would have written inv to use it. You can also use pinv, useful sometimes when singularity is an issue, but pinv does NOT create a true inverse for a singular system, anymore than does any other method. And pinv will be slower by far. No inverse exists for a singular matrix, any more than you can compute the multiplicative inverse of 0.

In general however, the best way to compute an inverse is to not compute the inverse at all. It is relatively rare that you ever truly need to compute an inverse matrix. Most of the time when you think you need that inverse, it is because you need to solve a linear system of equations. Backslash does that for you.

So why is the inverse so "inaccurate"? Time for you to learn about the condition number of a matrix. In the case I ran, for a 5K square random matrix as you generate it I got

cond(A)
ans =
   3.5996e+11

This is a big number. It means that your inverse matrix is pretty close to being numerical garbage, although it is not truly numerically singular. I'd expect a 10Kx10K matrix formed this way to be closer yet to numerical singularity.

So I have no idea why you are doing what you are doing. If this is homework, then you will learn something, I hope. If not homework, then your questions indicate you would benefit from a class on linear algebra, especially numerical linear algebra. At least, crack a good text and do some reading.

[Santosh Tiwari]

Hello,

Thanks for your response.

This is not a homework assignment, but work I am doing for an actual project. The code was just a sample I wrote to demonstrate the issue.

I am NOT using the determinant to determine if a Matrix is singular. cholesky returns error if Matix is really badly conditioned/singular. We introduce regularization to improve conditioning in such a case.

I actually need the log of the determinant - which is actually a finite number.

log(det(A)) = log(L11*L22*L33*...) = log(L11) + log(L22) + ...

Since, the determinant in general will be a huge number, we take log of the left and right side and solve log likelihood estimate precisely because of this reason.

If lower triangular matrix L is available, then log of determinant is a finite number which can be computed without any overflow issues.

I know about the condition number and have implemented regularization to improve ill-conditioning of my matrix.

Whenever possible, I do not actually compute the inverse, but use back-substitution. In this case, I need the actual inverse which will be stored in some file, and used later thousands of times, so calculating the actual inverse is indeed desirable.

[Santosh Tiwari]

Sorry, if my question was not clear. Let me refine my question.

Let A be a positive definite matrix. To compute Ainv and log of its determinant, I have following code.

L = chol(A,'lower');
U = L\eye(sizeOfA);
Ainv = L'\U;
logDetA = sum(log(abs(diag(L))));

The above code computes the inverse as well as log(det(A)).

Can I do the above more reliably and more efficiently assuming matrix A has dimension of the order 10000.

Thanks :)