# Context

In connection to this question I am interested in orthogonalizing known matrices.

As a test case, let us consider the definite positive $$15 \times 15$$ matrix mat

mat={{1, Sqrt[Pi]/2, -Sqrt[Pi]/2, 0, 0, 0, (-3*Sqrt[Pi])/4, (3*Sqrt[Pi])/4, (-3*Sqrt[Pi])/4, (3*Sqrt[Pi])/4, -1, 1, -1, 1, -1},
{Sqrt[Pi]/2, 2, 0, (5*Sqrt[Pi])/4, -Sqrt[Pi]/4, (-3*Sqrt[Pi])/4, -1, 1, -1, 1, (-33*Sqrt[Pi])/8, (21*Sqrt[Pi])/8, (-9*Sqrt[Pi])/8,
(-3*Sqrt[Pi])/8, (15*Sqrt[Pi])/8}, {-Sqrt[Pi]/2, 0, 2, (3*Sqrt[Pi])/4, Sqrt[Pi]/4, (-5*Sqrt[Pi])/4, 1, -1, 1, -1, (-15*Sqrt[Pi])/8,
(3*Sqrt[Pi])/8, (9*Sqrt[Pi])/8, (-21*Sqrt[Pi])/8, (33*Sqrt[Pi])/8}, {0, (5*Sqrt[Pi])/4, (3*Sqrt[Pi])/4, 7, 1, -1, (39*Sqrt[Pi])/8,
Sqrt[Pi]/8, (-9*Sqrt[Pi])/8, (-15*Sqrt[Pi])/8, -10, 6, -2, -2, 6}, {0, -Sqrt[Pi]/4, Sqrt[Pi]/4, 1, 3, 1, (21*Sqrt[Pi])/8,
(11*Sqrt[Pi])/8, (-11*Sqrt[Pi])/8, (-21*Sqrt[Pi])/8, 2, -2, 2, -2, 2}, {0, (-3*Sqrt[Pi])/4, (-5*Sqrt[Pi])/4, -1, 1, 7,
(15*Sqrt[Pi])/8, (9*Sqrt[Pi])/8, -Sqrt[Pi]/8, (-39*Sqrt[Pi])/8, 6, -2, -2, 6, -10},
{(-3*Sqrt[Pi])/4, -1, 1, (39*Sqrt[Pi])/8, (21*Sqrt[Pi])/8, (15*Sqrt[Pi])/8, 36, 6, 0, -6, (423*Sqrt[Pi])/16, (45*Sqrt[Pi])/16,
(-33*Sqrt[Pi])/16, (-75*Sqrt[Pi])/16, (-105*Sqrt[Pi])/16}, {(3*Sqrt[Pi])/4, 1, -1, Sqrt[Pi]/8, (11*Sqrt[Pi])/8, (9*Sqrt[Pi])/8, 6,
12, 2, 0, (201*Sqrt[Pi])/16, (123*Sqrt[Pi])/16, (-31*Sqrt[Pi])/16, (-93*Sqrt[Pi])/16, (-135*Sqrt[Pi])/16},
{(-3*Sqrt[Pi])/4, -1, 1, (-9*Sqrt[Pi])/8, (-11*Sqrt[Pi])/8, -Sqrt[Pi]/8, 0, 2, 12, 6, (135*Sqrt[Pi])/16, (93*Sqrt[Pi])/16,
(31*Sqrt[Pi])/16, (-123*Sqrt[Pi])/16, (-201*Sqrt[Pi])/16}, {(3*Sqrt[Pi])/4, 1, -1, (-15*Sqrt[Pi])/8, (-21*Sqrt[Pi])/8,
(-39*Sqrt[Pi])/8, -6, 0, 6, 36, (105*Sqrt[Pi])/16, (75*Sqrt[Pi])/16, (33*Sqrt[Pi])/16, (-45*Sqrt[Pi])/16, (-423*Sqrt[Pi])/16},
{-1, (-33*Sqrt[Pi])/8, (-15*Sqrt[Pi])/8, -10, 2, 6, (423*Sqrt[Pi])/16, (201*Sqrt[Pi])/16, (135*Sqrt[Pi])/16, (105*Sqrt[Pi])/16, 249,
39, 9, -9, -39}, {1, (21*Sqrt[Pi])/8, (3*Sqrt[Pi])/8, 6, -2, -2, (45*Sqrt[Pi])/16, (123*Sqrt[Pi])/16, (93*Sqrt[Pi])/16,
(75*Sqrt[Pi])/16, 39, 69, 15, -3, -9}, {-1, (-9*Sqrt[Pi])/8, (9*Sqrt[Pi])/8, -2, 2, -2, (-33*Sqrt[Pi])/16, (-31*Sqrt[Pi])/16,
(31*Sqrt[Pi])/16, (33*Sqrt[Pi])/16, 9, 15, 41, 15, 9}, {1, (-3*Sqrt[Pi])/8, (-21*Sqrt[Pi])/8, -2, -2, 6, (-75*Sqrt[Pi])/16,
(-93*Sqrt[Pi])/16, (-123*Sqrt[Pi])/16, (-45*Sqrt[Pi])/16, -9, -3, 15, 69, 39},
{-1, (15*Sqrt[Pi])/8, (33*Sqrt[Pi])/8, 6, 2, -10, (-105*Sqrt[Pi])/16, (-135*Sqrt[Pi])/16, (-201*Sqrt[Pi])/16, (-423*Sqrt[Pi])/16,
-39, -9, 9, 39, 249}};


If I use arbitrary precision to Orthogonalize it and subtract the result with 20 and 40 digits

tt = Orthogonalize[IdentityMatrix[Length[mat]], N[#1.mat.#2, 20] &];
tt2 = Orthogonalize[IdentityMatrix[Length[mat]],
N[#1.mat.#2, 40] &];


then the difference

  tt - tt2


shows catastrophic loss of accuracy.

# Question

Any idea what is going wrong there?

FYI, the conditioning of the matrix is not that bad.

mat // N // Eigenvalues // MinMax // First[#]/Last[#] &


2.3 10^-7

Of course I can push the number of digits, say to 40 and then works but it seems a bit extreme for such a simple matrix?

tt = Orthogonalize[IdentityMatrix[Length[mat]], N[#1.mat.#2, 40] &];
tt2 = Orthogonalize[IdentityMatrix[Length[mat]],
N[#1.mat.#2, 50] &]; tt - tt2


In the more realistic case of interest to me, the matrix is $$66\times 66$$ and working with 500 digits seems to do the trick, but it seems a bit excessive! (In that case the conditioning of the matrix is $$10^{21}$$).

Note that it seems about 430 digits have been lost!

tt = Orthogonalize[IdentityMatrix[Length[mat]], N[#1.mat.#2, 500] &];
tt2 = Orthogonalize[IdentityMatrix[Length[mat]], N[#1.mat.#2, 510] &];
tt - tt2 // Flatten // Abs // Max

1. 10^67
• Probably one has to increase the internal working precision... Jan 11, 2019 at 13:03
• Using native precision seems to do better: t = Orthogonalize[IdentityMatrix[Length[mat]], N[#1.mat.#2] &]; t - tt2 // Abs // Max gives 6.90505*10^-11. It's just the case of 20 digits that goes wrong. Looking at the rows of tt you see that the precision keeps getting worse (and the accuracy too), while for native precision it doesn't. Jan 11, 2019 at 13:56
• Map[Precision, tt2, {2}] shows the Gram-Schmidt process loses about 33 digits of precision. The same loss is observed for higher precision input. Also, tt.mat.Transpose[tt] is "almost" equal to the identity matrix, except for the last four rows/columns, where the the loss of precision yields zero. -- I expect if one analyzes Gram-Schmidt, the precision lost is as expected. Note Orthogonalize[N[IdentityMatrix[Length[mat]], 20], N[#1.mat.#2, 20] &] does better: it keeps the working precision at 20, which seems to help because the extra guard digits are not lost, I suppose. Jan 11, 2019 at 17:11
• @MichaelE2 thanks. I wondered about wrapping N[#,20]& around the identity but naively I assumed it would be arbitrary precision. In what sense doesn't this answer my question though? Jan 11, 2019 at 17:13
• Okay, but I'm not familiar with the numerics of Gram-Schmidt, and it would take too long to figure out how to explain "what is going wrong here." I'm not surprised that the iterative step of "subtract projection and rescale" might lead to significant errors, and I just guess that the rescaling magnifies the error in some steps in this case. Jan 11, 2019 at 17:23

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;
LinearAlgebraLAPACKTRTRI["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

• Thanks. The advantage being it will be faster? Do you know which version orthogonalise is using? Mar 12, 2019 at 21:18
• Faster and more accurate, when I checked. I do not think Orthogonalize[] uses this method, much less recognize a "weighted norm" problem. Mar 12, 2019 at 22:19