Catastrophic loss of accuracy in Orthogonalize
Map[Precision, tt2, {2}]
shows the Gram-Schmidt process loses about 33 digits of precision on the OP's example. This precision loss estimate is the result of Mathematica estimating the propagated error and is normally an upper bound on the error. When the error is greater than the point-estimate of the result, the result is returned as an arbitrary-precision zero, which might be viewed as similar to underflow. Starting with 20-digit precision then seems to be insufficient.
However, there is a workaround. It seems if the precision of the first argument of Orthogonalize
is not MachinePrecision
nor Infinity
, then the precision of the argument sets the working precision. (Presumably this is done by setting $MaxPrecision
and $MinPrecision
.) Since arbitrary precision numbers carry extra guard digits, the actual rounding error is usually less than error estimate calculated by the arbitrary-precision software. Fixing the precision via $MaxPrecision
and $MinPrecision
prevents the calculation of the propagated error, and we won't get the underflow-like zeros mentioned above. That seems to be enough to produce an accurate result starting with 20-digit precision. This produces an accurate result:
tt1 = Orthogonalize[N[IdentityMatrix[Length[mat]], 20], N[#1.mat.#2, 20] &]
Then tt1.mat.Transpose[tt1]
agrees with the identity matrix to ~16 digits, and tt1
agrees with tt2
to ~6 digits, which is not that bad considering that the precision of tt2
is only ~7.5 digits.
At least for the OP's special case, there is a simpler way to calculate the lower triangular matrix tt
:
tt = LowerTriangularize[Inverse[Transpose[CholeskyDecomposition[N[mat, 20]]]]];
Norm[tt.mat.Transpose[tt] - IdentityMatrix[Length[mat]], "Frobenius"]
0
To avoid the cleanup needed to be done with LowerTriangularize[]
, here is an equivalent method (using undocumented functionality):
tt = CholeskyDecomposition[N[mat, 20]]; chk = False;
LinearAlgebra`LAPACK`TRTRI["U", "N", tt, chk];
tt = Transpose[tt];
An alternative algorithm, if CholeskyDecomposition[]
is unsuitable (i.e. if mat
is ill-conditioned), uses the eigendecomposition:
{Λ, v} = Eigensystem[N[mat, 20]];
tt = First[QRDecomposition[DiagonalMatrix[Sqrt[Λ]].v]].DiagonalMatrix[1/Sqrt[Λ]].v;
Norm[tt.mat.Transpose[tt] - IdentityMatrix[Length[mat]], "Frobenius"]
0