0
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ you are asking for output that is higher precision than the input $\endgroup$
    – george2079
    Apr 17, 2014 at 11:31
  • 2
    $\begingroup$ 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. $\endgroup$ Apr 17, 2014 at 12:12
  • $\begingroup$ 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. $\endgroup$
    – celtschk
    Apr 17, 2014 at 14:26
  • $\begingroup$ 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. $\endgroup$
    – user18742
    Apr 21, 2014 at 8:36

1 Answer 1

3
$\begingroup$

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]];
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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