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

enter image description here

shows catastrophic loss of accuracy.


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

enter image description here

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
  • 1
    $\begingroup$ Probably one has to increase the internal working precision... $\endgroup$ – Henrik Schumacher Jan 11 at 13:03
  • 1
    $\begingroup$ 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. $\endgroup$ – Roman Jan 11 at 13:56
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 Jan 11 at 17:11
  • $\begingroup$ @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? $\endgroup$ – chris Jan 11 at 17:13
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Jan 11 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"]

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"]
  • $\begingroup$ Thanks. The advantage being it will be faster? Do you know which version orthogonalise is using? $\endgroup$ – chris Mar 12 at 21:18
  • $\begingroup$ Faster and more accurate, when I checked. I do not think Orthogonalize[] uses this method, much less recognize a "weighted norm" problem. $\endgroup$ – J. M. will be back soon Mar 12 at 22:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.