代码之家  ›  专栏  ›  技术社区  ›  Thomas Ahle

工作矩阵平方根

  •  4
  • Thomas Ahle  · 技术社区  · 15 年前

    我试着求矩阵的平方根。就是找到矩阵 B 所以 B*B=A . 我找到的方法都没有一个有效的结果。

    首先我发现这个公式 Wikipedia :

    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 应该收敛到 .

    然而,在python中实现算法(对逆矩阵使用numpy)给了我垃圾结果:

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

    我对矩阵平方根不太了解,但我想一定有比上面更好的算法?

    3 回复  |  直到 15 年前
        1
  •  8
  •   steabert    15 年前

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

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

    编辑 试试下面的,你会发现它是正确的:

    asmatrix(scipy.linalg.sqrtm(matrix('1,2;3,4')))**2
    
        2
  •  3
  •   David Heffernan    15 年前

    你的矩阵[12;34]不是正的,所以在实矩阵域中没有问题的解。

        3
  •  2
  •   peter karasev    15 年前

    你做矩阵平方根的目的是什么?我怀疑在实际应用中,矩阵可能是对称正定的(例如协方差),所以你不应该遇到复数。

    在这种情况下,您可以计算cholesky分解,如按比例LU分解,请参见此处: http://en.wikipedia.org/wiki/Cholesky_decomposition

    另一个实际的例子是如果你的矩阵是旋转的,那么你可以先用矩阵对数分解,然后在对数空间除以2,然后用矩阵指数返回旋转。。。在任何情况下,要求“通用矩阵平方根”听起来都很奇怪,您可能希望更深入地理解特定的应用程序。