Search code examples
matlabmathmatrixdeterminants

Determinant is showing infinity instead of zero! Why?

This is my matlab code I wrote for a problem I got as homework. after multiplication of A and its transpose the resulting square matrix should have determinant zero according all classmates as their codes (different one) gave them so. Why is my code not giving the determinant of c and d to be infinity

A = rand(500,1500);

b = rand(500,1);

c = (A.')*A;

detc = det(c);

cinv = inv((A.')*A);

d = A*(A.');

detd = det(d);

dinv = inv(A*(A.'));

x1 = (inv((A.')*A))*((A.')*b);

x2 = A.'*((inv(A*(A.')))*b);
Solution

  • This behavior is explained in the Limitations section of the det's documentation and exemplified in the Find Determinant of Singular Matrix subsection where it is stated:

    The determinant of A is quite large despite the fact that A is singular. In fact, the determinant of A should be exactly zero! The inaccuracy of d is due to an aggregation of round-off errors in the MATLAB® implementation of the LU decomposition, which det uses to calculate the determinant.

    That said, in this instance, you can produce your desired result by using the m-code implementation given on that same page but sorting the diagonal elements of U in an ascending matter. Consider the sample script:

    clc();
clear();

A = rand(500,1500);
b = rand(500,1);
c = (A.')*A;

[L,U] = lu(c);
% Since det(L) is always (+/-)1, it doesn't impact anything
diagU = diag(U);
detU1 = prod(diagU);
detU2 = prod(sort(diagU,'descend'));
detU3 = prod(sort(diagU,'ascend'));

fprintf('Minimum: %+9.5e\n',min(abs(diagU)));
fprintf('Maximum: %+9.5e\n',max(abs(diagU)));
fprintf('Determinant:\n');
fprintf('\tNo Sort:          %g\n'  ,detU1);
fprintf('\tDescending Sort:  %g\n'  ,detU2);
fprintf('\tAscending Sort:   %g\n\n',detU3);

    This produces the output:

    Minimum: +1.53111e-13
Maximum: +1.72592e+02
Determinant:
    No Sort:          Inf
    Descending Sort:  Inf
    Ascending Sort:   0

    Notice that the direction of the sort matters, and that no-sorting gives Inf since a true 0 doesn't exist on the diagonal. The descending sort sees the largest values multiplied first, and apparently, they exceed realmax and are never multiplied by a true 0, which would generate a NaN. The ascending sort clumps together all of the near-zero diagonal values with very few large negative values (in truth, a more robust method would sort based on magnitude, but that was not done here), and their multiplication generates a true 0 (meaning that the value falls below the smallest denormalized number available in IEEE-754 arithmetic) that produces the "correct" result.

    All that written, and as others have implied, I'll quote original Matlab developer and Mathworks co-founder Cleve Moler:

    [The determinant] is useful in theoretical considerations and hand calculations, but does not provide a sound basis for robust numerical software.