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

13 views (last 30 days)
Mohammad Mehboudi on 6 Oct 2015
Commented: Bobby Cheng on 8 Oct 2015
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

Bobby Cheng on 6 Oct 2015
Edited: Bobby Cheng on 7 Oct 2015
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.
##### 3 CommentsShow 1 older commentHide 1 older comment
Steven Lord on 8 Oct 2015
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 on 8 Oct 2015
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.

Steven Lord on 6 Oct 2015
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 on 6 Oct 2015
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''.