# How to get more accurate orthonormal eigenvectors?

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.

-
you are asking for output that is higher precision than the input – george2079 Apr 17 '14 at 11:31
Basically you cannot have more precision in your eigenvectors than you have in your matrix M. So either you increase the precision of your matrix values "at source" or you artificially increase the precision, for example using SetPrecision‌​. But using SetPrecision to increase precision is arbitrarily padding your values with (base 2) zeros, in reality the calculated eigenvectors will be no more precise. Browse the Related posts for more information. – MikeLimaOscar Apr 17 '14 at 12:12
Also note that a higher precision at the source may mean significant higher calculation time, since the processor floating point capabilities cannot be used any more. – celtschk Apr 17 '14 at 14:26
Thanks for your comments/answers. In order to increase precision SetPrecision adds noise to the original data so I don't think I want to use it. – user18742 Apr 21 '14 at 8:36

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];


 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]];

-