Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
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];

your test:

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


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


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

 ev = Eigenvectors[SetPrecision[m, 32]];
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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