Wrong output from det(Array{Complex})

First of all, you should understand that accurate computation of the determinant is a notoriously difficult problem. See, for example,

In particular, that paper says that the condition number of det(A) for an n×n matrix A is on the order of norm(A, Inf)^n. For your matrix, that gives about 2.44×10⁶. So, I would naively expect that computing determinants of your matrix could easily lose 6 digits to roundoff errors, which means that for Float16 arithmetic (which keeps 3 digits) I would expect to lose everything, and even Float32 arithmetic (which keeps about 7 digits) is problematic.

det(A) works by performing the LU factorization of A and then multiplying the diagonal entries of U. This is efficient for large n, O(n³), but it incurs roundoff errors even for integer-valued matrix entries.

For a 3×3 matrix, you can instead use the direct permutation formula for the determinant, which is exact for integer values that aren’t too big. This is built-in to the StaticArrays.jl package (it is a good idea to use StaticArrays anyway, for efficiency, if you are working with fixed size ≤ 4×4 matrices):

julia> using StaticArrays

julia> A = @SMatrix Complex{Float32}[20-50im -10+20im -10+30im; -10.0+20.0im 26.0-52.0im -16.0+32.0im; -10.0+30.0im -16.0+32.0im 26.0-62.0im];

julia> det(A)
0.0f0 + 0.0f0im

Even with SMatrix, Float16 won’t work because the individual terms in the determinant formula overflow the maximum Float16 value. But it is really inadvisable to do computations with Float16 anyway — Float32 or Float64 are much faster. Float16 is mainly useful for low-precision storage of the result.

All this being said, you should really think carefully about whether you need to compute the determinant at all. It is rarely used nowadays in numerical computations — there is often a better way to do whatever you wanted the determinant for.

7 Likes