Рабочая матрица квадратного корня

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 для обратных матриц) дала мне глупые результаты:

>>> 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 возвращает массив, а НЕ матрицу. Следовательно, использование * для их умножения не очень хорошая идея. В вашем случае решения, таким образом, верны, но вы этого не видите.

edit попробуйте следующее, вы увидите, что это правильно:

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

Answer 2

Ваша матрица [1 2;3 4] не является положительным, поэтому нет решения проблемы в области вещественных матриц.

Answer 3

Какова цель матричного квадратного корня, который вы делаете? Я подозреваю, что в практическом применении матрица действительно может быть симметричной положительно определенной (например, ковариация), поэтому вы не должны сталкиваться с комплексными числами.

В этом случае вы можете вычислить разложение Холецкого, например, масштабированное разложение LU, см. Здесь: http://en.wikipedia.org/wiki/Cholesky_decomposition

Другой практический пример: если ваши матрицы являются вращениями, то вы можете сначала разложить с помощью log log и просто разделить на 2 в пространстве журнала, а затем вернуться к вращению с показателем матрицы ... в любом случае это звучит странно, что вы спросите «общий матричный квадратный корень», вы, вероятно, хотите понять конкретное приложение более подробно.