Рабочая матрица Квадратный корень

Question

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

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

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

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