I'm trying to take the square root of a matrix. That is find the matrix B
so B*B=A
. None of the methods I've found around gives a working result.
First I found this formula on Wikipedia:
Set Y_0 = A
and Z_0 = I
then the iteration:
Y_{k+1} = .5*(Y_k + Z_k^{-1}),
Z_{k+1} = .5*(Z_k + Y_k^{-1}).
Then Y
should converge to B
.
However implementing the algorithm in python (using numpy for inverse matrices), gave me rubbish results:
>>> def denbev(Y,Z,n):
if n == 0: return Y,Z
return denbev(.5*(Y+Z**-1), .5*(Z+Y**-1), n-1)
>>> denbev(matrix('1,2;3,4'), matrix('1,0;0,1'), 3)[0]**2
matrix([[ 1.31969074, 1.85986159],
[ 2.78979239, 4.10948313]])
>>> denbev(matrix('1,2;3,4'), matrix('1,0;0,1'), 100)[0]**2
matrix([[ 1.44409972, 1.79685675],
[ 2.69528512, 4.13938485]])
As you can see, iterating 100 times, gives worse results than iterating three times, and none of the results get within a 40% error margin.
Then I tried the scipy sqrtm method, but that was even worse:
>>> scipy.linalg.sqrtm(matrix('1,2;3,4'))**2
array([[ 0.09090909+0.51425948j, 0.60606061-0.34283965j],
[ 1.36363636-0.77138922j, 3.09090909+0.51425948j]])
>>> scipy.linalg.sqrtm(matrix('1,2;3,4')**2)
array([[ 1.56669890+0.j, 1.74077656+0.j],
[ 2.61116484+0.j, 4.17786374+0.j]])
I don't know a lot about matrix square rooting, but I figure there must be algorithms that perform better than the above?
(1) the square root of the matrix [1,2;3,4] should give something complex, as the eigenvalues of that matrix are negative. SO your solution can't be correct to begin with.
(2) linalg.sqrtm returns an array, NOT a matrix. Hence, using *
to multiply them is not a good idea. In your case, the solutions is thus correct, but you're not seeing it.
edit try the following, you'll see it's correct:
asmatrix(scipy.linalg.sqrtm(matrix('1,2;3,4')))**2
Your matrix [1 2; 3 4] isn't positive so there is no solution to the problem in the domain of real matrices.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With