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

13 views (last 30 days)

Show older comments

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 Comments

### Answers (2)

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 Comments

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

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!