How to get more accurate orthonormal eigenvectors?

Question

I have a matrix M with real entries in MachinePrecision. I compute its eigenvectors and construct matrix U:

  ev = Eigenvectors[M];    
  U = Transpose @ ev;

U is supposed to be an orthogonal matrix and thus Transpose[U].U == IdentityMatrix[ Length @ U].
Mathematica gives the latter equality with the accuracy of 10^(-16):

Transpose[U].U - IdentityMatrix[ Length @ U]// Abs// Max == 8.8131 * 10^(-16).

Is it possible to make this difference smaller, let's say, order of 10^(-32)? I need to get eigenvectors which are orthonormal with accuracy of order of 10(-32) or higher.

Answer 1

a machine precision symmetric matrix:

 n = 5;
 m = Nest[ 
     Join[#, {Table[#[[i, (Length[#] + 1)]], {i, Length[#]}]~Join~
            RandomReal[{-10, 10}, {n - Length[#]}]}] & , {} , n];

convert the matrix to an exact rational form, do the eignevector extraction and convert back to high precision floats:

 ev = N[Normalize /@ Eigenvectors[(Rationalize[#, 10^-16] & /@ m)], 32];

your test:

 U = Transpose@ev;
 (Transpose[U].U - IdentityMatrix[Length@U] ) // Abs // Max

0.*10^-32

This is going to quickly get real slow for larger n..

edit

after working through it the hard way, I'll note the @MikeLimaOscars approach works better (faster), ie.

 ev = Eigenvectors[SetPrecision[m, 32]];