Is there a big error in the method matlab uses to finds sqrtm?
5 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(2 个)
Bobby Cheng
2015-10-6
编辑:Bobby Cheng
2015-10-7
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 个评论
Steven Lord
2015-10-8
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
2015-10-8
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
2015-10-6
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!