6. 고유공간과 고유분해

 

 

In [2]:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

 

6. 고유공간과 대각화¶

 

1) 대각화¶

 

(1) 대각화는 다음의 경우에 활용할 수 있다.¶

• 선형계의 발산의 판단 : 행렬A를 대각화 하는 경우, 대각행렬의 요소값을 통해 발산 여부를 판단할 수 있다.(|a| > 1이면 발산한다)
• 연산의 최적화 : 행렬 A를 연속적으로 적용하는 변환의 경우, 대각화 할 경우 대각행렬만 계속 곱해주면 A를 연속 적용하는것과 같다.

 

(2) 발산 판단 여부의 판단¶

 

- 대각행렬의 경우¶

In [3]:

A = np.diag([5,-3,0.8])
x = np.array(["X_1,X_2,X_3"])

 

• A와 x를 내적하면

In [4]:

Ax = np.matrix(["5 * x_1", "-3 * x_2", "0.8 * x_3"])
print(Ax.T)

 

[['5 * x_1']
 ['-3 * x_2']
 ['0.8 * x_3']]

 

• 꼴로 나타난다. 따라서, 차분방정식을 이용할경우 파라미터 A의 t제곱(t = infinity)에서 안정된 값이 구해지는데

In [5]:

A ** np.inf

Out[5]:

array([[inf,  0.,  0.],
       [ 0., inf,  0.],
       [ 0.,  0.,  0.]])

 

이 대각행렬은 a33 성분을 제외하고는 무한으로 발산하므로, 이 시스템은 t = infinity일때 발산한다.

 

- 대각행렬은 아니지만, 대각화가 가능한 경우¶

In [6]:

A = np.matrix([[5,1],[1,5]])
x_1 = np.array(["x1_1","x2_1"]) #단, x1_1은 변수 x1의 1차시 전 값을, x2_1은 변수 x2의 1차시 전 값을 의미한다. 

 

• 이 때, 이 행렬의 곱은

In [7]:

print(np.matrix([["X1","X2"]]))
print("=")
print(A)
print("*")
print(x_1)
print("=")
print(np.matrix([["(5 * x1_1) + (x2_1)"],["(x1_1) + (5 * x2_1)"]]))

 

[['X1' 'X2']]
=
[[5 1]
 [1 5]]
*
['x1_1' 'x2_1']
=
[['(5 * x1_1) + (x2_1)']
 ['(x1_1) + (5 * x2_1)']]

 

• 변수변환을 다음과 같이 실시한다.
• y1 = x1 + x2
• y2 = x1 - x2

In [8]:

C = np.matrix([[1,1],[1,-1]])
x = np.array(["X1,X2"]) #단, X1, X2는 변수 X1, X2의 현상태의 값을 의미한다.

Cx = np.matrix([["X1 + X2"],["X1 - X2"]])

In [9]:

print(np.matrix([["y1"],["y2"]]),"\n = \n", Cx)

 

[['y1']
 ['y2']] 
 = 
 [['X1 + X2']
 ['X1 - X2']]

 

• 변수변환과 원래 행렬을 동시에 정의하면

In [10]:

print(np.matrix([["y1"],["y2"]]))
print("=")
print(A,"\n * \n",C,"\n * \n",x_1)

 

[['y1']
 ['y2']]
=
[[5 1]
 [1 5]] 
 * 
 [[ 1  1]
 [ 1 -1]] 
 * 
 ['x1_1' 'x2_1']

 

• 한편, 우변을 하나로 정리해주면

In [11]:

print(np.matrix([["5 * X1_1 + X2_1 + X1_1 + 5 * X2_1"],["5 * X1_1 + X2_1 - X1_1 - 5 * X2_1"]]))
print("=")
print(np.matrix([["6X1 + 6X2"],["4X1 - 4X2"]]))

 

[['5 * X1_1 + X2_1 + X1_1 + 5 * X2_1']
 ['5 * X1_1 + X2_1 - X1_1 - 5 * X2_1']]
=
[['6X1 + 6X2']
 ['4X1 - 4X2']]

 

• 위 식을 다시 행렬꼴로 나타내면

In [12]:

K = np.matrix([[6,0],[0,4]])
print(np.matrix([["y1"],["y2"]]))
print("=")
print(np.matrix([["6X1 + 6X2"],["4X1 - 4X2"]]))
print("=")
print(K)
print("*")
print(C)
print("*")
print(x_1)

 

[['y1']
 ['y2']]
=
[['6X1 + 6X2']
 ['4X1 - 4X2']]
=
[[6 0]
 [0 4]]
*
[[ 1  1]
 [ 1 -1]]
*
['x1_1' 'x2_1']

 

이 때, 대각행렬은 차분방정식에서 t = infinity일때 안정화되므로, 이 계는 궁극적으로 발산한다.

In [13]:

# - 한편
print(np.matrix([["y1"],["y2"]]))
print("=")
print(C)
print("*")
print(x)
# 이므로

 

[['y1']
 ['y2']]
=
[[ 1  1]
 [ 1 -1]]
*
['X1,X2']

In [14]:

print("C행렬의 역 : \n",np.linalg.inv(C)) #을 이용하여 표현하면
print("\n")
print(x)
print("=")
print(np.linalg.inv(C))
print("*")
print(np.matrix([["y1"],["y2"]]))

#이다.

 

C행렬의 역 : 
 [[ 0.5  0.5]
 [ 0.5 -0.5]]


['X1,X2']
=
[[ 0.5  0.5]
 [ 0.5 -0.5]]
*
[['y1']
 ['y2']]

In [15]:

# 앞서 

print(np.matrix([["y1"],["y2"]]))
print("=")
print(C)
print("*")
print(x)
# 이므로, 이를 앞에서 구한 역행렬과 내적하면

 

[['y1']
 ['y2']]
=
[[ 1  1]
 [ 1 -1]]
*
['X1,X2']

In [16]:

print(x)
print("=")
print(np.linalg.inv(C))
print("*")
print(K)
print("*")
print(C)
print("*")
print(x_1)

 

['X1,X2']
=
[[ 0.5  0.5]
 [ 0.5 -0.5]]
*
[[6 0]
 [0 4]]
*
[[ 1  1]
 [ 1 -1]]
*
['x1_1' 'x2_1']

 

• 이 때, C와 C^-1 사이에 대각행렬 diag(6,4)가 끼어있는 형태가 나오고, 이 형식을 맨 처음에 언급한

In [17]:

print(np.matrix([["X1","X2"]]))
print("=")
print(A)
print("*")
print(x_1)

#와 비교하면

 

[['X1' 'X2']]
=
[[5 1]
 [1 5]]
*
['x1_1' 'x2_1']

In [18]:

print(A)
print("=")
print(np.linalg.inv(C))
print("*")
print(K)
print("*")
print(C)

#임이 명백하다.

 

[[5 1]
 [1 5]]
=
[[ 0.5  0.5]
 [ 0.5 -0.5]]
*
[[6 0]
 [0 4]]
*
[[ 1  1]
 [ 1 -1]]

 

• A = inv(C) diag(ㅅ) C 꼴로 분해하였는데, 이 것을 대각화라고 한다.
• 이제 관건은, 앞서 대각화를 가능하게 했던 행렬

In [19]:

print(C)

#를 찾는 것이다.

 

[[ 1  1]
 [ 1 -1]]

 

2) 고윳값과 고유벡터¶

 

(1) 고유벡터의 기하학적 이해¶

 

(2) 고윳값과 고유벡터의 도출¶

 

- 특성방정식¶

 

• 어떤 행렬이 특이행렬이 되도록 만들어주는 고윳값을 구하는 방정식
• 즉, det(A - aI) = 0이 되도록 하는 a를 구한다.

In [40]:

A = np.matrix([[6,-3,5],[-1,4,-5],[-3,3,-4]])

 

• 3X3 행렬이므로, 사루스의 법칙을 이용하여 행렬식을 구하면

det(A - aI) = {(6-a)(4-a)(-4-a)]} + {-3 -5 -3} + {5 -1 3} - {5 (4-a) -3} - {(6-a) -5 3} - {-3 -1 (4-a)} = a^3-6a^2+11a-6

• 이 때 해는

In [41]:

a = np.matrix([3,2,1])
print(a)

 

[[3 2 1]]

 

이다

 

• 이제 고윳값에 맞는 고유벡터를 구하면
• 각 고윳값을 대각성분에서 뺀 행렬의 영공간을 구한다.

In [60]:

A_eig_1 = (A - a[0,0] * np.identity(3))
print(A_eig_1)

 

[[ 3. -3.  5.]
 [-1.  1. -5.]
 [-3.  3. -7.]]

 

이 행렬을 행연산하면 다음과 같은 행렬이 된다.

In [81]:

A_eig_1_reduced = np.matrix([[1,-1,0],[0,0,1],[0,0,0]])
print(A_eig_1_reduced)

 

[[ 1 -1  0]
 [ 0  0  1]
 [ 0  0  0]]

 

위 행줄임 행렬을 영공간의 Span으로 나타내면

In [82]:

print(np.matrix(["x1","x2","x3"]).T)
print("=")
print(np.matrix([1,1,0]).T)
print("*")
print("x2")
print("단, x2는 0 이외의 모든 실수이다.")

 

[['x1']
 ['x2']
 ['x3']]
=
[[1]
 [1]
 [0]]
*
x2
단, x2는 0 이외의 모든 실수이다.

In [73]:

A_eig_2 = (A - a[0,1] * np.identity(3))
print(A_eig_2)

 

[[ 4. -3.  5.]
 [-1.  2. -5.]
 [-3.  3. -6.]]

 

이 행렬을 행연산하면 다음과 같은 행렬이 된다.

In [76]:

A_eig_2_reduced = np.matrix([[1,0,-1],[0,1,-3],[0,0,0]])
print(A_eig_2_reduced)

 

[[ 1  0 -1]
 [ 0  1 -3]
 [ 0  0  0]]

 

위 행줄임 행렬을 영공간의 Span으로 나타내면

In [79]:

print(np.matrix(["x1","x2","x3"]).T)
print("=")
print(np.matrix([1,3,1]).T)
print("*")
print("x3")
print("단, x3는 0 이외의 모든 실수이다.")

 

[['x1']
 ['x2']
 ['x3']]
=
[[1]
 [3]
 [1]]
*
x3
단, x3는 0 이외의 모든 실수이다.

In [85]:

A_eig_3 = (A - a[0,2] * np.identity(3))
print(A_eig_3)

 

[[ 5. -3.  5.]
 [-1.  3. -5.]
 [-3.  3. -5.]]

 

이 행렬을 행연산하면 다음과 같은 행렬이 된다.

In [88]:

A_eig_3_reduced = np.matrix([[1,0,0],[0,1,(-5/3)],[0,0,0]])
print(A_eig_3_reduced)

 

[[ 1.          0.          0.        ]
 [ 0.          1.         -1.66666667]
 [ 0.          0.          0.        ]]

 

위 행줄임 행렬을 영공간의 Span으로 나타내면

In [90]:

print(np.matrix(["x1","x2","x3"]).T)
print("=")
print(np.matrix([0,5,3]).T)
print("*")
print("x3")
print("단, x3는 0 이외의 모든 실수이다.")

 

[['x1']
 ['x2']
 ['x3']]
=
[[0]
 [5]
 [3]]
*
x3
단, x3는 0 이외의 모든 실수이다.

In [99]:

eigen_value = np.diag([3,2,1])
eigen_vector = np.matrix([[1,1,0],[1,3,1],[0,5,3]]).T

In [100]:

eigen_value

Out[100]:

array([[3, 0, 0],
       [0, 2, 0],
       [0, 0, 1]])

In [109]:

eigen_vector

Out[109]:

matrix([[1, 1, 0],
        [1, 3, 5],
        [0, 1, 3]])

 

이다

 

3) 요르단 표준형¶

 

(1) 개요¶

- 정석적인 방법으로 대각화가 이루어지지 않는 경우 고려할 수 있는 대안이다.¶

- 주로, 고유공간을 구할 때 원공간의 차원 만큼 고유공간의 기저가 생성되지 않을때 채택할 수 있다.¶

• 고윳값이 n중해를 가지고, 이 때 n개의 선형독립한 고유벡터가 생성되지 않을경우 기저가 하나 모자르게 된다.

- 이 때, 요르단 분해는 유사 대각행렬의 역할을 하며, 어떤 정칙행렬 P에 대하여¶

• Ax = y일때, (P^(-1) J P) * x = y로 분해할 수 있다.

 

- 요르단 표준형은 다음의 형태를 띈다.¶

In [156]:

upper = np.vstack([np.hstack([a,np.zeros((4,3)),np.zeros((4,3))]),np.hstack([np.zeros((3,4)),b,np.zeros((3,3))])])
J = np.vstack([upper,np.hstack([np.zeros((3,7)),c])])
print(J)

 

[[3. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 3. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 3. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 3. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 3. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 3. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 3. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 5. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 5. 1.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 5.]]

 

- 이 때, 특정한 규칙에 따라 요르단 행렬을 블록화 할 수 있는데,¶

• 동일한 고윳값을 공유하면서(위에서 3과 5)
• 그 비스듬한 요소에는 1이 들어간다.
• 그리고, 구분점이 되는 지점에선 비스듬한 요소가 1 대신 0이 들어간다.

 

- 예를들어, 위 요르단 표준형은 다음과 같은 두개의 요르단 셀로 분해할 수 있다.¶

In [138]:

a = np.matrix([[3,1,0,0],[0,3,1,0],[0,0,3,1],[0,0,0,3]])
b = np.matrix([[3,1,0],[0,3,1],[0,0,3]])
c = np.matrix([[5,1,0],[0,5,1],[0,0,5]])
print("첫 번째 요르단 셀")
print(a)
print("두 번째 요르단 셀")
print(b)
print("세 번째 요르단 셀")
print(c)

 

첫 번째 요르단 셀
[[3 1 0 0]
 [0 3 1 0]
 [0 0 3 1]
 [0 0 0 3]]
두 번째 요르단 셀
[[3 1 0]
 [0 3 1]
 [0 0 3]]
세 번째 요르단 셀
[[5 1 0]
 [0 5 1]
 [0 0 5]]

 

(2) 요르단 표준형의 성질¶

 

- 각 셀의 행(열)수는 그 고윳값의 중복해의 갯수를 의미한다.¶

In [149]:

print(c)

 

[[5 1 0]
 [0 5 1]
 [0 0 5]]

 

• 에서, 고윳값 5의 중복해는 3개가 존재한다는것을 의미한다.

In [151]:

print(a)
print(b)

 

[[3 1 0 0]
 [0 3 1 0]
 [0 0 3 1]
 [0 0 0 3]]
[[3 1 0]
 [0 3 1]
 [0 0 3]]

 

• 에서, 고윳값 3의 중복해는 총 7개가 존재한다는것을 의미한다

 

- 각 고윳값에 대응하는 셀의 수는 그 고유공간의 선형독립인 기저의 수를 의미한다.¶

In [153]:

print(c)

 

[[5 1 0]
 [0 5 1]
 [0 0 5]]

 

• 에서, 고윳값 5의 경우 선형독립인 고유기저는 하나만 존재한다는 것을 의미한다.
• 즉, 3개의 중복해에 하나의 선형독립 고유벡터가 존재한다.

In [154]:

print(a)
print(b)

 

[[3 1 0 0]
 [0 3 1 0]
 [0 0 3 1]
 [0 0 0 3]]
[[3 1 0]
 [0 3 1]
 [0 0 3]]

 

• 란 특이 케이스에서, 고윳값 3의 경우, 선형독립인 고유기저는 두 개가 존재한다는 것을 의미한다.
• 즉 7개의 중복해에 두 개의 선형독립 고유벡터가 존재한다.

 

- 요르단 표준형의 거듭제곱¶

• 요르단 표준형의 한 셀은 다음과 같이 분리가 가능하다.

In [161]:

z = np.matrix([[0,1,0,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]])
print(3 * np.identity(4) + z)

 

[[3. 1. 0. 0.]
 [0. 3. 1. 0.]
 [0. 0. 3. 1.]
 [0. 0. 0. 3.]]

 

• 이 때, z는 다음과 같다.

In [162]:

print(z)

 

[[0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]
 [0 0 0 0]]

 

• z를 거듭제곱하면 신기한 형태를 관찰할 수 있다.

In [165]:

print("2제곱")
print(z ** 2)
print("3제곱")
print(z ** 3)
print("4제곱")
print(z ** 4)

 

2제곱
[[0 0 1 0]
 [0 0 0 1]
 [0 0 0 0]
 [0 0 0 0]]
3제곱
[[0 0 0 1]
 [0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]
4제곱
[[0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]

 

• 즉, 유사 대각인 1 성분들이 점점 우상단으로 밀려나다가, 모두 밀려나면 0행렬이 된다.

 

• J = aI + z에서. J^2을 하면(단, a는 해당 셀의 고윳값)
• J^2 = (aI + z)^2 = a^2I + 2aZ + z^2

In [174]:

# 요르단 표준형의 거듭제곱

print("a^2I")
print(3**2 * np.identity(4))
print("2a*z")
print(2*3*z)
print("z^2")
print(z**2)
print("선형결합")
print(3**2 * np.identity(4) + 2*3*z + z**2)

 

a^2I
[[9. 0. 0. 0.]
 [0. 9. 0. 0.]
 [0. 0. 9. 0.]
 [0. 0. 0. 9.]]
2a*z
[[0 6 0 0]
 [0 0 6 0]
 [0 0 0 6]
 [0 0 0 0]]
z^2
[[0 0 1 0]
 [0 0 0 1]
 [0 0 0 0]
 [0 0 0 0]]
선형결합
[[9. 6. 1. 0.]
 [0. 9. 6. 1.]
 [0. 0. 9. 6.]
 [0. 0. 0. 9.]]

 

• J = aI + z에서. J^3을 하면(단, a는 해당 셀의 고윳값)
• J^3 = (aI + z)^3 = (a^3I) + (3 a^2 z) + (3 a z^2) + (z^3)

In [183]:

print("a^3I")
print(3**3 * np.identity(4))
print("3* a^2 * z")
print(3 * (3 ** 2) * z)
print("3 * a * z^2")
print(3 * 3 * (z ** 2))
print("z^3")
print(z ** 3)
print("선형결합")
print(3**3 * np.identity(4) + 3 * (3 ** 2) * z + 3 * 3 * (z ** 2) + z ** 3)

 

a^3I
[[27.  0.  0.  0.]
 [ 0. 27.  0.  0.]
 [ 0.  0. 27.  0.]
 [ 0.  0.  0. 27.]]
3* a^2 * z
[[ 0 27  0  0]
 [ 0  0 27  0]
 [ 0  0  0 27]
 [ 0  0  0  0]]
3 * a * z^2
[[0 0 9 0]
 [0 0 0 9]
 [0 0 0 0]
 [0 0 0 0]]
z^3
[[0 0 0 1]
 [0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]
선형결합
[[27. 27.  9.  1.]
 [ 0. 27. 27.  9.]
 [ 0.  0. 27. 27.]
 [ 0.  0.  0. 27.]]

 

• J = aI + z에서. J^4을 하면(단, a는 해당 셀의 고윳값)
• J^4 = (aI + z)^4 = (a^4I) + (4 a^3 z) + (6 a^2 z^2) + (4 a z^3) + (z^4)

In [189]:

print("a^4I")
print(3 ** 4 * np.identity(4))
print("4 * a^3 * z")
print(4 * (3 ** 3) * z)
print("6 * a^2 * z^2")
print(6 * (3 ** 2) * (z ** 2))
print("4 * a * z^3")
print( 4 * 3 * (z ** 3))
print("z ^4(영행렬)")
print( z ** 4)
print("선형결합")
print(3 ** 4 * np.identity(4) + 4 * (3 ** 3) * z + 6 * (3 ** 2) * (z ** 2) + 4 * 3 * (z ** 3) + z ** 4)

 

a^4I
[[81.  0.  0.  0.]
 [ 0. 81.  0.  0.]
 [ 0.  0. 81.  0.]
 [ 0.  0.  0. 81.]]
4 * a^3 * z
[[  0 108   0   0]
 [  0   0 108   0]
 [  0   0   0 108]
 [  0   0   0   0]]
6 * a^2 * z^2
[[ 0  0 54  0]
 [ 0  0  0 54]
 [ 0  0  0  0]
 [ 0  0  0  0]]
4 * a * z^3
[[ 0  0  0 12]
 [ 0  0  0  0]
 [ 0  0  0  0]
 [ 0  0  0  0]]
z ^4(영행렬)
[[0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]
선형결합
[[ 81. 108.  54.  12.]
 [  0.  81. 108.  54.]
 [  0.   0.  81. 108.]
 [  0.   0.   0.  81.]]

 

이상을 보면, 어떤 패턴을 확인할 수 있다.

• 거듭제곱 행렬의 첫번째 대각성분은 고윳값 ^ (거듭제곱수) 이다.
• 거듭제곱 행렬의 두번째 대각성분은 (거듭제곱수) 고윳값 ^ (거듭제곱수 - 1) z행렬
• 거듭제곱 행렬의 세번째 대각성분은 (거듭제곱수)C(2) 고윳값 ^ (거듭제곱수 - 2) z행렬^2 (단, (a)C(b) 는 a개에서 b개를 뽑는 조합)
• 거듭제곱 행렬의 네번째 대각성분은 (거듭제곱수)C(3) 고윳값^(거듭제곱수 - 3) z행렬^3

 

위를 일반화하면 다음과 같다.

• 거듭제곱 행렬의 첫번째 대각성분은 고윳값 ^ (거듭제곱수)이다.
• 거듭제곱 행렬의 두번째 대각성분은 (거듭제곱수) 고윳값 ^ (거듭제곱수 -1) z행렬
• 거듭제곱 행렬의 s번째 대각성분은 (거듭제곱수)C(s) 고윳값 ^ (거듭제곱수 - s) z행렬^s

 

위를 다른 형태로 변환하면 더 흥미로운 형식이 된다.

• f(a) = a^t. 즉 고윳값을 t제곱하는 함수를 정의하자.
• d(s,f(a)) 를 f(a)를 s번 미분하는 함수라고 정의할 때, 일반화한 요르단 행렬의 모습은

In [208]:

row_1 = np.array(["f(a)","d(1,f(a))","(1/2) * d(2,f(a))","(1/3!) * d(3,f(a)","...","1/(s-1) * d(s-1,f(a))"])
row_2 = np.array(["0","f(a)","d(1,f(a))","(1/2) * d(2,f(a))","...","1/(s-2) * d(s-2,f(a))"])
row_3 = np.array(["0","0","f(a)","(1/2) * d(2,f(a)),","...","1/(s-3) * d(s-3,f(a))"])
row_4 = np.array(["0","0","0","f(a)","...","1/(s-4) * d(s-4,f(a))"])

print(np.matrix([row_1,row_2,row_3,row_4]))

 

[['f(a)' 'd(1,f(a))' '(1/2) * d(2,f(a))' '(1/3!) * d(3,f(a)' '...'
  '1/(s-1) * d(s-1,f(a))']
 ['0' 'f(a)' 'd(1,f(a))' '(1/2) * d(2,f(a))' '...'
  '1/(s-2) * d(s-2,f(a))']
 ['0' '0' 'f(a)' '(1/2) * d(2,f(a)),' '...' '1/(s-3) * d(s-3,f(a))']
 ['0' '0' '0' 'f(a)' '...' '1/(s-4) * d(s-4,f(a))']]

 

각 행(row)만 봤을때, row의 각 요소들은 테일러 전개에서 각 항의 Parameter를 따서 입력한것과 같은 꼴이다. 즉, 이 상삼각행렬의 각 요소들은 테일러전개의 Parameter를 값으로 갖는다.

 

3) 총 정리 예제¶

 

(1) 어떤 자기회귀(AR)모델이 다음과 같이 정의되었다.¶

• y(t) = -0.5y(t-1) + 0.34y(t-2) + 0.08y(t-3) + 2u(t)
• 초기조건 y(0) = 0.78, y(-1) = 0.8, y(-2) = 1.5

In [90]:

def init_x(t):
    if t == 0:
        y_1 = 0.8
        y_2 = 1.5
        y_3 = (0.78+0.5*(0.8) - 0.34*(1.5))/0.08
        return ([y_1,y_2,y_3])

In [91]:

A = np.matrix([[-0.5,0.34,0.08],[1,0,0],[0,1,0]])
x = np.matrix(init_x(0)).T
print(A)

 

[[-0.5   0.34  0.08]
 [ 1.    0.    0.  ]
 [ 0.    1.    0.  ]]

 

• 위 방정식의 경우, A와 x를 내적하면

In [92]:

np.dot(A,x) # 이것이 t=0일때의 초깃값 y_1, y_2, y_3이다.

Out[92]:

matrix([[0.78],
        [0.8 ],
        [1.5 ]])

In [98]:

# 행렬 A를 고유분해하면

A_eig = np.linalg.eig(A)
print("행렬 A의 고윳값 : ","\n",A_eig[0])
print("행렬 A의 고유벡터 : ","\n",A_eig[1])

 

행렬 A의 고윳값 :  
 [-0.8  0.5 -0.2]
행렬 A의 고유벡터 :  
 [[-0.447039   -0.21821789  0.03919309]
 [ 0.55879876 -0.43643578 -0.19596545]
 [-0.69849845 -0.87287156  0.97982725]]

In [117]:

#이 계의 고윳값 대각행렬의 안정값은
np.diag(A_eig[0]) ** np.inf

Out[117]:

array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

In [118]:

# 이므로, 모두 0이기 때문에 궁극적으로 발산하지 않고 수렴한다.