Is there a big error in the method matlab uses to finds sqrtm?

By definition a pure [density] matrix "A" is a square complex matrix with an eigenvalue equal to one and all others equal to zero. Such a matrix has the unique property that its sqrtm (and all other matrix powers) is equal to itself, i.e. sqrtm(A)=A holds.

The problem is that when I define some specific pure matrices in matlab and find their sqrtm, their sqrtm are not equal to them, with a considerable error. Here is my specific example:

theta=0.589048622548086;

e0=[-1i*cos(theta/2),0,0,sin(theta/2)];

%%%it is clear this is a normalized vector, i.e. e0*e0'=1;

A=e0'*e0; 

%%%A here is pure. why? Because A^2=(e0'*e0)*(e0'*e0)=e0'*(e0*e0')*e0=e0'*e0=A, hence taking matrix square root of both sides gives A=sqrtm(A) which verifies A is pure.

Now do this numerically in matlab (look at sqrtm(A)-A for example, or alternatively see how trace(sqrtm(A)) is different from trace(A)=1), and you will find out that matlab makes a big [to me] error.

I would be happy if there is a way to ask matlab do this more precisely, or there is another precise method which can be used.

Thanks

Answer 1

[Bobby Cheng]

This is an interesting example. It goes to the heart of the issue when computes with floating point arithmetic. First consider this

 >> theta=0.589048622548086;
 >> e0=[-1i*cos(theta/2),0,0,sin(theta/2)];
 >> A=e0'*e0;
 >> eig(A)
 ans =
             0
             0
    1.3878e-17
    1.0000e+00

Clearly there is some mistake here. Matrix A should only have rank 1. But MATLAB thinks that it has rank 2. And it is the small nonzero eigenvalues that is causing the problem.

 >> sqrt(1.3878e-17)
 ans =
    3.7253e-09

Has EIG gone dotty? Hardly, it is trying its very best. Let double check the result using symbolic toolbox.

 >> sA = sym(A,'f')
  sA =
 [    8248205863506131/9007199254740992, 0, 0, 5004131788810439i/18014398509481984]
 [                                    0, 0, 0,                                   0]
 [                                    0, 0, 0,                                   0]
 [ -5004131788810439i/18014398509481984, 0, 0,    379496695617431/4503599627370496]
 >> double(eig(sA))
 ans =
             0
             0
    7.7070e-18
    1.0000e+00

Better but still rank 2. So one can conclude that the problem exists in Symbolic Toolbox too. Maybe, but this is not the root cause. What I did not say was that I told Symbolic Toolbox to treat A as exact, and that assumption was wrong. Matrix A was not the "A" that you think it is. It is a very good approximation of the real thing. To see this, use symbolic calculation to form A.

 >> theta= sym(0.589048622548086)
 theta =
 (3*pi)/16
 >> e0=[-1i*cos(theta/2),0,0,sin(theta/2)];
 >> A=e0'*e0  
 A =
 [                  cos((3*pi)/32)^2, 0, 0, cos((3*pi)/32)*sin((3*pi)/32)*1i]
 [                                 0, 0, 0,                                0]
 [                                 0, 0, 0,                                0]
 [ -cos((3*pi)/32)*sin((3*pi)/32)*1i, 0, 0,                 sin((3*pi)/32)^2] 
 >> double(eig(A))
 ans =
      0
      0
      0
      1

Yes, rank 1 as expected. But wait, can use this to form a better A? Maybe, but no guarantee. On my computer,

 >> eig(double(A))
 ans =
             0
             0
    1.3878e-17
    1.0000e+00

This is a true limitation of floating point arithmetic. By the time SQRTM sees A, it is too late.

[Mohammad Mehboudi]

Thanks a lot Bobby for your answer.

I need to make a comment though. It is true that even symbolically this problem exists yet, but even if one could overcome this issue symbolically the main problem is still there, which is probably the method matlab uses to find sqrtm, (or mpower in general). I say this because I solve the same problem with mathematica (I define everything numerically in mathematica), but mathematica gives the exact solution, i.e. eigenvalues are exact and I can find the fidelity=trace(sqrtm(A)*A*sqrtm(A))=1 in mathematica, while matlab gives NaN for fidelity.

I was always thinking that matlab is more trustable with numbers!

[Steven Lord]

If you want to avoid floating-point arithmetic entirely, start off with a symbolic theta.

 theta = 3*sym('pi')/16;
 e0 = [-1i*cos(theta/2), 0, 0 sin(theta/2)];
 A = e0'*e0;
 T = sqrtm(A);
 isAlways(T == A)

By taking floating-point arithmetic out of the picture, and using symbolic calculations or variable-precision arithmetic (which is what I suspect Mathematica is using under the covers, even though it looks like it's doing computations in double precision) you'll see that all the elements in T are exactly equal to the corresponding elements in A.

[Bobby Cheng]

Thank you for putting your trust in MATLAB. I am not sure what Mathematica does. But for numerical computing, there are limitations. The algorithm B = SQRTM(A) tries to compute B such that B*B-A is small.

 >> format short e
 >> theta=0.589048622548086;
 >> e0=[-1i*cos(theta/2),0,0,sin(theta/2)];
 >> A=e0'*e0;
 >> B = sqrtm(A)
 Warning: Matrix is singular and may not have a square root. 
 > In sqrtm (line 92) 
 B =
   9.1573e-01 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 2.7779e-01i
   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i
   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i
   0.0000e+00 - 2.7779e-01i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   8.4265e-02 + 0.0000e+00i
>> B*B-A
 ans =
   2.2204e-16 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 5.5511e-17i
   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i
   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i
   0.0000e+00 - 5.5511e-17i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i   0.0000e+00 + 0.0000e+00i

For special matrices like yours, we may need specialized algorithms.

Answer 2

[Steven Lord]

Which release are you using? We've upgraded the algorithms in SQRTM, LOGM, and EXPM in release R2015b as indicated in the Release Notes. But note that calling SQRTM on your matrix issues a warning that your matrix is singular and may not have a square root.

[Mohammad Mehboudi]

I was using R2015a until yesterday, and upgraded to R20015b today, with the hope not to deal the same problem, but the problem is still there.

You are right about Singularity, for that, one can do as follows:

[B,res]=sqrtm(A);

Then there will be no warning. The sqrtm of A is now returned in B. But apart from the warning, on paper we know that this matrix has a sqrtm, which is equal to itself indeed. The question is how can matlab reproduce it for us?

Another thing is that, you may wonder why a very small error (order of -8 or something) is important? because In the rest of my calculations I need to find such a function:

fidelty=trace(sqrtm(sqrtm(A)*A*sqrtm(A)));

On paper we know this is equal to one for this case, but matlab gives ''NaN''.