# 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 Apr 17, 2014 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. Apr 17, 2014 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. Apr 17, 2014 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. Apr 21, 2014 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]];