#4957 determinant and invert of a float matrix are catastrophically wrong

[David Scherfgen]

determinant of a float matrix returns the wrong sign and 130x the magnitude; invert misses its own defining identity by 6627

Maxima 5.49post, SBCL 2.6.7.

A Hilbert matrix is symmetric positive definite, so its determinant cannot be negative. X.A = I needs no external reference at all.

(%i1) display2d:false$
(%i2) load(linearalgebra)$
(%i3) H: float(genmatrix(lambda([i,j],1/(i+j-1)),8,8))$
(%i4) determinant(H);
(%o4) -3.5571424818265393e-31
(%i5) block([fpprec:40], bfloat(determinant(matrixmap(rationalize,H))));
(%o5) 2.737050121755728820366909340980048809227b-33
(%i6) determinant_by_lu(H,'floatfield);
(%o6) 2.737050139851956e-33
(%i7) mat_cond(H,1);
(%o7) 3.3872793842156906e10
(%i8) mx(R):=lmax(map(abs,flatten(args(R))))$
(%i9) mx(invert(H).H-ident(8));
(%o9) 6627.625557230902
(%i10) mx(invert_by_gausselim(H).H-ident(8));
(%o10) 2.9802322387695313e-8

Expected: (%o4) = +2.7370501217557288e-33; (%o9) should be around 1e-6 or better, as (%o10) is on the identical matrix.

(%o5) is the exact rational determinant of exactly the doubles passed in, so no perturbation slack is involved. (%o7) shows the conditioning permits about 7e-6 relative error, and three other Maxima routes achieve 1e-8; determinant gets zero correct digits.

Solving with the inverse shows the effect: invert(H).b for b = H.(1,...,1) gives components from -29619.6 to +36556.5 instead of all ones, while invert_by_gausselim gives 1.0000000 to 7 digits.

[David Billinghurst]

A cross-check with octave using double-precision floating point.

octave:12> H=hilb(8);
octave:13> [D, RCOND] = det(H)
D = 2.7371e-33
RCOND = 2.9522e-11
octave:14> 1/RCOND
ans = 3.3873e+10

The Hilbert matrices are ill-conditioned, but these results show that floating-point calculations can achieve results similar to the bfloat calcs above. RCOND is the reciprocal of the condition number. maxima can do better.