2010-11-21 2 views
3

행렬의 제곱근을 취하려고합니다. 그 행렬은 B이므로 B*B=A입니다. 내가 찾은 방법 중 어느 것도 작동 결과를 얻지 못합니다.행렬 제곱근

먼저 나는 Wikipedia에이 공식을 발견

설정 Y_0 = AZ_0 = I을 다음 반복 :

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

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

그런 다음 YB으로 수렴해야한다.

그러나, (역 행렬에 대한 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]]) 

,하지만 난 위보다 더 나은 수행이 있어야 알고리즘을 파악 ?

답변

7

(1) 행렬 [1,2,3,4]의 제곱근은 그 행렬의 고유 값이 음수이므로 복잡한 것을 제공해야합니다. 그래서 당신의 솔루션은 처음부터 정확할 수 없습니다.

(2) linalg.sqrtm은 배열이 아니라 배열을 반환합니다. 따라서 *을 사용하여 곱하면 좋은 생각이 아닙니다. 귀하의 경우, 솔루션은 정확하지만 당신은 그것을 보지 못합니다.

편집가 다음을 시도, 당신은 맞습니다 볼 수 있습니다 :

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

감사합니다. 완벽하게 작동합니다. 단지'sqrtm (matrix ('1,2,3,4') ** 2)'이 작동하지 않는 이유는 아직도 모르겠다. –

+0

그것은 효과가 있습니다 ... 솔루션을 제곱 해보십시오. http://www.mathworks.com/help/techdoc/ref/sqrtm을 참조하십시오.htaccess Admin Home English Language Content – steabert

+1

맞아, 거기에 솔루션을 많이하고 sqrtm은 단순히 나보다 다른 하나를 선택합니다. –

3

행렬 [1 2; 3 4]는 양의 값이 아니므로 실제 행렬의 영역에서 문제에 대한 해결책이 없습니다.

2

당신이하고있는 매트릭스 제곱근의 목적은 무엇인가? 실용적인 응용 프로그램이 매트릭스가 실제로 대칭 긍정 (예 : 공분산) 될 수 있으므로 복잡한 숫자가 발생하지 않아야합니다. 당신이 확장 LU 인수 분해처럼 콜레 분해를 계산할 수있는 경우

여기를 참조 : http://en.wikipedia.org/wiki/Cholesky_decomposition

또 다른 실제적인 예를 들어 당신의 행렬은 회전이 있다면, 먼저 매트릭스 로그로 분해 할 수 있고 단지로 분할이다 2 로그 공간에서 다시 행렬 지수로 회전 시키십시오 ... "제네릭 행렬 제곱근"을 요구하는 이상한 소리가 들리면 특정 응용 프로그램을 더 깊이 이해하고 싶을 것입니다.