Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
1  
you are asking for output that is higher precision than the input –  george2079 Apr 17 at 11:31
2  
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 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 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 at 8:36

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

Your Answer

 
discard

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.