工作矩阵平方根

Question

我试图取矩阵的平方根.这是找到矩阵B等等B*B=A.我找到的所有方法都没有给出一个有效的结果.

首先我在维基百科上找到了这个公式:

设置Y_0 = A和Z_0 = I随后迭代:

Y_{k+1} = .5*(Y_k + Z_k^{-1}),

Z_{k+1} = .5*(Z_k + Y_k^{-1}).

然后Y应该收敛B.

然而,在python中实现算法(使用numpy for inverse matrices),给了我垃圾结果:

>>> 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]])

正如您所看到的,迭代100次,得到的结果比迭代三次更糟,并且没有任何结果在40%的误差范围内.

然后我尝试了scipy sqrtm方法,但更糟糕的是:

>>> 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]])

我不太了解矩阵平方根,但我认为必须有比上述更好的算法？

Answer 1

(1)矩阵[1,2; 3,4]的平方根应该给出一些复杂的东西,因为该矩阵的特征值是负的.所以你的解决方案一开始就不正确.

(2)linalg.sqrtm返回一个数组,而不是一个矩阵.因此,使用*它们乘以它们并不是一个好主意.在您的情况下,解决方案是正确的,但您没有看到它.

编辑尝试以下,你会看到它是正确的:

asmatrix(scipy.linalg.sqrtm(matrix('1,2;3,4')))**2