Arbeiten Matrix Square root

Ich versuche, nehmen Sie die Quadratwurzel einer matrix. Das die matrix B so B*B=A. Keine der Methoden, die ich gefunden habe um gibt ein funktionierendes Ergebnis.

Zuerst fand ich diese Formel auf Wikipedia:

Set Y_0 = A und Z_0 = I dann die iteration:

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

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

Dann Y konvergieren sollte, um B.

Jedoch die Umsetzung des Algorithmus in python (mit numpy für inverse Matrizen), gab mir den Müll Ergebnisse:

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

Wie Sie sehen können, Durchlaufen 100 mal, gibt schlimmer Ergebnisse als die Iteration von drei mal, und keines der Ergebnisse erhalten Sie innerhalb einer 40% Fehlerquote.

Dann habe ich versucht die scipy sqrtm-Methode, aber das war noch schlimmer:

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

Ich weiß nicht viel über die matrix quadratisch wühlen, aber ich glaube, es muss algorithmen, die besser sind als die oben genannten?