Lecture 22. 대각화와 A의 거듭제곱 — 파이썬 실습 - 선형대수

L22 서술 파트의 결론은 $A = S\Lambda S^{-1}$ 이었고, 그 뜻은 번역 → 계산 → 역번역이었다.

이 노트북에서는 $S$ 와 $\Lambda$ 를 직접 만들어 $AS = S\Lambda$ 를 확인하고, $A^k$ 를 두 가지 방법으로 구해 대조하며, 고윳값의 위치가 운명을 어떻게 가르는지 잰다. 마지막으로 피보나치 비가 황금비로 가는 것을 확인한다.

서술 파트의 내용	여기서 확인하는 방법
$AS = S\Lambda$	고유벡터를 열에 세우고 곱해 본다
세 단계	$S^{-1}\vv{x}$ , $\Lambda\vv{c}$ , $S(\cdot)$ 를 따로 찍는다
$A^k = S\Lambda^k S^{-1}$	직접 곱한 것과 대조
$\lvert\lambda\rvert < 1$ 이면 $A^k \to 0$	세 경우를 나란히
성분 크기와 무관	성분에 20이 있는데도 0으로
지배 고윳값	나머지 항이 얼마나 빨리 죽는지
대각화 실패	$3I$ 와 $[3\;1;0\;3]$ 대조
피보나치	비가 황금비로 수렴

0. 준비¶

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, show_matrix

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(22)
print("numpy", np.__version__)

numpy 2.5.2

1. $S$ 와 $\Lambda$ 를 만든다¶

고유벡터를 열에 세우고 고윳값을 같은 순서로 대각에 놓는다.

def 대각화(A):
    """A = S Lambda S^-1 의 세 조각을 돌려준다. 대각화가 안 되면 알려 준다."""
    A = np.asarray(A, dtype=float)
    값, S = np.linalg.eig(A)
    if np.linalg.matrix_rank(S) < A.shape[0]:
        raise ValueError("독립인 고유벡터가 모자라 대각화할 수 없다")
    return S, np.diag(값), np.linalg.inv(S)

A = np.array([[2.0, 1.0], [1.0, 2.0]])
S, L, S역 = 대각화(A)

print(show_matrix(S.real, "S  (열이 고유벡터)"))
print(show_matrix(L.real, "Lambda  (대각이 고윳값)"))
print()
print("A S 와 S Lambda 가 같은가 :", np.allclose(A @ S, S @ L))
print("S Lambda S^-1 이 A 인가   :", np.allclose(S @ L @ S역, A))
print("S^-1 A S 가 Lambda 인가   :", np.allclose(S역 @ A @ S, L))

S  (열이 고유벡터)
[   0.707   -0.707 ]
[   0.707    0.707 ]
Lambda  (대각이 고윳값)
[  3   0 ]
[  0   1 ]

A S 와 S Lambda 가 같은가 : True
S Lambda S^-1 이 A 인가   : True
S^-1 A S 가 Lambda 인가   : True

순서를 어긋나게 놓으면 등식이 깨진다. 고유벡터의 순서를 바꾸면 고윳값도 함께 바꿔야 한다.

S뒤집 = S[:, ::-1]                       # 고유벡터만 순서를 바꾼다
print("고유벡터만 바꿨을 때  A S = S L 인가 :", np.allclose(A @ S뒤집, S뒤집 @ L))

L뒤집 = np.diag(np.diag(L)[::-1])        # 고윳값도 함께 바꾼다
print("둘 다 바꿨을 때        A S = S L 인가 :", np.allclose(A @ S뒤집, S뒤집 @ L뒤집))

고유벡터만 바꿨을 때  A S = S L 인가 : False
둘 다 바꿨을 때        A S = S L 인가 : True

2. 세 단계를 따로 본다¶

$A\vv{x}$ 를 한 번에 계산하지 말고 번역 · 계산 · 역번역으로 나눠 찍어 보자.

x = np.array([3.0, 1.0])

c = S역 @ x                  # 1단계 : 고유좌표로 번역
Lc = L @ c                   # 2단계 : 성분마다 lambda 곱하기
되돌림 = S @ Lc              # 3단계 : 표준좌표로 역번역

print("x            =", x)
print("c = S^-1 x   =", c.real, "   <- 고유기저에서의 좌표")
print("Lambda c     =", Lc.real, "   <- 성분마다 3배, 1배")
print("S (Lambda c) =", 되돌림.real)
print("A x          =", A @ x)
print("같은가 :", np.allclose(되돌림.real, A @ x))

x            = [3. 1.]
c = S^-1 x   = [ 2.8284 -1.4142]    <- 고유기저에서의 좌표
Lambda c     = [ 8.4853 -1.4142]    <- 성분마다 3배, 1배
S (Lambda c) = [7. 5.]
A x          = [7. 5.]
같은가 : True

$\vv{c}$ 가 정말 고유기저에서의 좌표인지 확인하자. $\vv{x} = c_1\vv{s}_1 + c_2\vv{s}_2$ 여야 한다.

조립 = c[0] * S[:, 0] + c[1] * S[:, 1]
print("c1 s1 + c2 s2 =", 조립.real)
print("x             =", x)
print("같은가 :", np.allclose(조립.real, x))

c1 s1 + c2 s2 = [3. 1.]
x             = [3. 1.]
같은가 : True

3. $A^k$ 는 $\Lambda^k$ 로 끝난다¶

def 거듭제곱(A, k):
    """A^k = S Lambda^k S^-1. k 가 아무리 커도 계산량이 같다."""
    S, L, S역 = 대각화(A)
    return (S @ np.diag(np.diag(L) ** k) @ S역).real

for k in (2, 5, 10):
    직접 = np.linalg.matrix_power(A, k)
    공식 = 거듭제곱(A, k)
    print(f"k = {k:>2} :  최대 차이 {np.abs(직접 - 공식).max():.2e}")

A10 = 거듭제곱(A, 10)
print()
print("A^10 의 성분 :", np.round(A10).astype(int).ravel())
print("(3^10 + 1)/2 =", (3 ** 10 + 1) // 2, "  <- 대각")
print("(3^10 - 1)/2 =", (3 ** 10 - 1) // 2, "  <- 대각 밖")
print("맞는가 :", np.allclose(A10, [[(3 ** 10 + 1) / 2, (3 ** 10 - 1) / 2],
                                   [(3 ** 10 - 1) / 2, (3 ** 10 + 1) / 2]]))

k =  2 :  최대 차이 0.00e+00
k =  5 :  최대 차이 0.00e+00
k = 10 :  최대 차이 0.00e+00

A^10 의 성분 : [29525 29524 29524 29525]
(3^10 + 1)/2 = 29525   <- 대각
(3^10 - 1)/2 = 29524   <- 대각 밖
맞는가 : True

$k$ 가 커져도 계산량은 그대로이다. 마르코프 행렬처럼 고윳값이 1과 0.5인 경우를 보자. $k$ 를 십억으로 두어도 곧바로 답이 나온다.

M = np.array([[0.8, 0.3],
              [0.2, 0.7]])                # 열의 합이 1
print("고윳값 :", np.linalg.eigvals(M))

for k in (10, 100, 10 ** 9):
    print(f"k = {k:>12,} :")
    print(show_matrix(거듭제곱(M, k), f"  M^{k}"))

고윳값 : [1. +0.j 0.5+0.j]
k =           10 :
  M^10
[    0.6   0.599 ]
[    0.4   0.401 ]
k =          100 :
  M^100
[  0.6   0.6 ]
[  0.4   0.4 ]
k = 1,000,000,000 :
  M^1000000000
[  0.6   0.6 ]
[  0.4   0.4 ]

$\lambda = 1$ 은 그대로 남고 $\lambda = 0.5$ 는 죽으므로, $M^k$ 가 랭크 1 행렬로 수렴한다. 그 극한의 열이 $\lambda = 1$ 의 고유벡터이다. 다음다음 강의의 주제이다.

극한 = 거듭제곱(M, 10 ** 9)
값, 벡터 = np.linalg.eig(M)
하나 = 벡터[:, np.argmin(np.abs(값 - 1.0))].real

print("M^(10^9) 의 랭크 :", np.linalg.matrix_rank(극한, tol=1e-8))
print("그 열            :", 극한[:, 0])
print("lambda=1 의 고유벡터 (합이 1 이 되게) :", 하나 / 하나.sum())
print("같은가 :", np.allclose(극한[:, 0], 하나 / 하나.sum()))

M^(10^9) 의 랭크 : 1
그 열            : [0.6 0.4]
lambda=1 의 고유벡터 (합이 1 이 되게) : [0.6 0.4]
같은가 : True

4. 단위원이 운명을 가른다¶

경우 = {
    "모두 안쪽": np.array([[0.8, 0.3], [0.0, 0.5]]),
    "원 위":     np.array([[np.cos(0.5), -np.sin(0.5)],
                          [np.sin(0.5), np.cos(0.5)]]),
    "바깥쪽":    np.array([[1.2, 0.3], [0.0, 1.1]]),
}

print(f"{'':>12}{'고윳값의 절댓값':>22}{'|A^30 x|':>14}")
x0 = np.array([1.0, 1.0])
for 이름, M in 경우.items():
    크기 = np.abs(np.linalg.eigvals(M))
    끝 = np.linalg.norm(np.linalg.matrix_power(M, 30) @ x0)
    print(f"{이름:>12}{np.round(크기, 4)!s:>22}{끝:>14.4e}")

                          고윳값의 절댓값      |A^30 x|
       모두 안쪽             [0.8 0.5]    2.4759e-03
         원 위               [1. 1.]    1.4142e+00
         바깥쪽             [1.2 1.1]    8.9733e+02

$A^k \to 0$ 인 것은 성분이 작아서가 아니다. 성분에 20이 있어도 고윳값이 작으면 0으로 간다.

큰성분 = np.array([[0.5, 20.0],
                 [0.0, 0.5]])

print(show_matrix(큰성분, "성분에 20 이 있다"))
print("고윳값 :", np.linalg.eigvals(큰성분))
print()
print(f"{'k':>5}{'|A^k| 의 최대 성분':>22}")
for k in (1, 5, 10, 20, 50, 100):
    print(f"{k:>5}{np.abs(np.linalg.matrix_power(큰성분, k)).max():>22.4e}")

성분에 20 이 있다
[  0.5    20 ]
[    0   0.5 ]
고윳값 : [0.5+0.j 0.5+0.j]

    k         |A^k| 의 최대 성분
    1            2.0000e+01
    5            6.2500e+00
   10            3.9062e-01
   20            7.6294e-04
   50            1.7764e-12
  100            3.1554e-27

처음에는 오히려 커진다. 20이라는 성분이 한동안 힘을 쓰기 때문이다. 그러나 결국 $0.5^k$ 가 이겨서 0으로 간다. 긴 눈으로 보면 고윳값만 남는다.

걸음 = np.arange(0, 61)
자료 = []
for 이름, M in list(경우.items()) + [("성분이 큰 경우", 큰성분)]:
    크기 = [np.linalg.norm(np.linalg.matrix_power(M, int(k)) @ x0) for k in 걸음]
    자료.append(go.Scatter(x=걸음, y=크기, mode="lines", name=이름,
                          line=dict(width=3)))

go.Figure(data=자료, layout=go.Layout(
    title=dict(text="|A^k x| 의 운명"),
    xaxis=dict(title=dict(text="k")),
    yaxis=dict(title=dict(text="|A^k x|"), type="log"),
    height=440, margin=dict(l=70, r=20, t=60, b=50)))

5. 지배 고윳값만 살아남는다¶

$A^k\vv{x} = \sum c_i\lambda_i^k\vv{s}_i$ 에서 가장 큰 항만 남는다. 남는 방향이 정말 지배 고유벡터인지 각도로 재 보자.

값, 벡터 = np.linalg.eig(A)
으뜸 = 벡터[:, np.argmax(np.abs(값))].real

print("지배 고윳값 :", 값[np.argmax(np.abs(값))].real)
print("그 고유벡터 :", 으뜸)
print()
print(f"{'k':>4}{'A^k x 의 방향':>26}{'으뜸과의 사잇각 (도)':>22}")
v = np.array([3.0, 1.0])
for k in range(0, 13, 2):
    w = np.linalg.matrix_power(A, k) @ v
    단위 = w / np.linalg.norm(w)
    코사인 = abs(단위 @ 으뜸) / np.linalg.norm(으뜸)
    각 = np.degrees(np.arccos(np.clip(코사인, -1, 1)))
    print(f"{k:>4}{np.round(단위, 4)!s:>26}{각:>22.6f}")

지배 고윳값 : 3.0
그 고유벡터 : [0.7071 0.7071]

   k                A^k x 의 방향          으뜸과의 사잇각 (도)
   0           [0.9487 0.3162]             26.565051
   2           [0.7452 0.6668]              3.179830
   4           [0.7115 0.7027]              0.353673
   6           [0.7076 0.7066]              0.039298
   8           [0.7072 0.7071]              0.004366
  10           [0.7071 0.7071]              0.000485
  12           [0.7071 0.7071]              0.000054

몇 번만 곱해도 방향이 지배 고유벡터에 붙는다. 이것이 고윳값을 구하는 가장 단순한 수치 방법(거듭제곱법)의 원리이다.

6. 대각화가 되는가¶

서술 파트의 대조를 코드로 확인한다. 고윳값이 같은데도 운명이 갈린다.

셋 = {
    "3 I":            3.0 * np.eye(2),
    "[3 1 ; 0 3]":    np.array([[3.0, 1.0], [0.0, 3.0]]),
    "[2 1 ; 1 2]":    A,
}

print(f"{'':>14}{'고윳값':>18}{'고유공간 차원':>14}{'대각화':>10}")
for 이름, M in 셋.items():
    값들 = np.linalg.eigvals(M).real
    _, V = np.linalg.eig(M)
    차원 = np.linalg.matrix_rank(V)
    가능 = "가능" if 차원 == M.shape[0] else "불가능"
    print(f"{이름:>14}{np.round(값들, 2)!s:>18}{차원:>14}{가능:>10}")

                             고윳값       고유공간 차원       대각화
           3 I           [3. 3.]             2        가능
   [3 1 ; 0 3]           [3. 3.]             1       불가능
   [2 1 ; 1 2]           [3. 1.]             2        가능

B = np.array([[3.0, 1.0], [0.0, 3.0]])
print(show_matrix(B - 3 * np.eye(2), "B - 3I"))
print("랭크 :", np.linalg.matrix_rank(B - 3 * np.eye(2)))
print("영공간의 차원 :", 2 - np.linalg.matrix_rank(B - 3 * np.eye(2)),
      " -> 고유벡터가 하나뿐")
print()
try:
    대각화(B)
except ValueError as 오류:
    print("대각화(B) ->", 오류)

B - 3I
[  0   1 ]
[  0   0 ]
랭크 : 1
영공간의 차원 : 1  -> 고유벡터가 하나뿐

대각화(B) -> 독립인 고유벡터가 모자라 대각화할 수 없다

그래도 $B^k$ 는 계산할 수 있다. 손으로 몇 번 곱해 보면 규칙이 보인다.

print(f"{'k':>4}{'대각 성분':>12}{'3^k':>10}{'오른쪽 위':>12}{'k 3^(k-1)':>12}")
for k in range(1, 7):
    Bk = np.linalg.matrix_power(B, k)
    print(f"{k:>4}{Bk[0, 0]:>12.0f}{3 ** k:>10}{Bk[0, 1]:>12.0f}{k * 3 ** (k - 1):>12}")

   k       대각 성분       3^k       오른쪽 위   k 3^(k-1)
   1           3         3           1           1
   2           9         9           6           6
   3          27        27          27          27
   4          81        81         108         108
   5         243       243         405         405
   6         729       729        1458        1458

$3^k$ 옆에 $k\,3^{k-1}$ 이라는 항이 붙는다. 대각화가 안 되면 이런 $k$ 가 붙은 항이 생기는데, 그 정체는 L28에서 밝힌다.

7. 피보나치와 황금비¶

피보 = np.array([[1.0, 1.0], [1.0, 0.0]])
값들 = np.linalg.eigvals(피보)
황금 = (1 + np.sqrt(5)) / 2

print("고윳값       :", np.sort(값들)[::-1])
print("황금비       :", 황금)
print("1 - 황금비   :", 1 - 황금)
print()
print("trace =", np.trace(피보), " (= 두 고윳값의 합)")
print("det   =", np.linalg.det(피보), " (= 두 고윳값의 곱)")

고윳값       : [ 1.618+0.j -0.618+0.j]
황금비       : 1.618033988749895
1 - 황금비   : -0.6180339887498949

trace = 1.0  (= 두 고윳값의 합)
det   = -1.0  (= 두 고윳값의 곱)

u0 = np.array([1.0, 0.0])                  # (F1, F0)
print(f"{'k':>4}{'F(k)':>12}{'F(k+1)/F(k)':>18}{'황금비와의 차이':>18}")
for k in range(1, 21):
    u = 거듭제곱(피보, k) @ u0
    F다음, F = u
    print(f"{k:>4}{F:>12.0f}{F다음 / F:>18.10f}{abs(F다음 / F - 황금):>18.2e}")

   k        F(k)       F(k+1)/F(k)          황금비와의 차이
   1           1      1.0000000000          6.18e-01
   2           1      2.0000000000          3.82e-01
   3           2      1.5000000000          1.18e-01
   4           3      1.6666666667          4.86e-02
   5           5      1.6000000000          1.80e-02
   6           8      1.6250000000          6.97e-03
   7          13      1.6153846154          2.65e-03
   8          21      1.6190476190          1.01e-03
   9          34      1.6176470588          3.87e-04
  10          55      1.6181818182          1.48e-04
  11          89      1.6179775281          5.65e-05
  12         144      1.6180555556          2.16e-05
  13         233      1.6180257511          8.24e-06
  14         377      1.6180371353          3.15e-06
  15         610      1.6180327869          1.20e-06
  16         987      1.6180344478          4.59e-07
  17        1597      1.6180338134          1.75e-07
  18        2584      1.6180340557          6.70e-08
  19        4181      1.6180339632          2.56e-08
  20        6765      1.6180339985          9.77e-09

비네 공식도 확인하자. 정수만 나오는 수열의 닫힌 꼴에 $\sqrt5$ 가 세 번 들어 있다.

def 비네(k):
    """F(k) 의 닫힌 꼴."""
    뿌리 = np.sqrt(5.0)
    return ((1 + 뿌리) ** k - (1 - 뿌리) ** k) / (2 ** k * 뿌리)


print(f"{'k':>4}{'비네 공식':>18}{'점화식':>12}")
F = [0, 1]
while len(F) < 25:
    F.append(F[-1] + F[-2])
for k in (1, 5, 10, 20, 24):
    print(f"{k:>4}{비네(k):>18.6f}{F[k]:>12}")

   k             비네 공식         점화식
   1          1.000000           1
   5          5.000000           5
  10         55.000000          55
  20       6765.000000        6765
  24      46368.000000       46368

마치며...¶

서술 파트의 내용	이 노트북의 코드
$AS = S\Lambda$	`대각화(A)` 와 세 가지 검산
순서를 맞춰야 한다	고유벡터만 바꾸면 등식이 깨진다
세 단계	`S역 @ x`, `L @ c`, `S @ Lc` 를 따로
$A^k = S\Lambda^k S^{-1}$	`거듭제곱(A, k)` — $k = 10^9$ 도 즉시
단위원	세 경우의 $\lvert A^k\vv{x}\rvert$
성분 크기와 무관	성분에 20이 있어도 0으로
지배 고윳값	방향이 몇 걸음 만에 붙는다
대각화 실패	고유공간 차원으로 판정, $B^k$ 에 $k\,3^{k-1}$
피보나치	비가 황금비로, 비네 공식

더 해 볼 것¶

거듭제곱 을 고윳값이 복소수인 회전행렬에 써 보자. 결과가 실수로 돌아오는가? 왜 그런가?
3절의 마르코프 행렬에서 시작 벡터를 바꿔 보자. $M^k\vv{x}$ 의 극한이 달라지는가?
6절의 $B^k$ 에서 오른쪽 위 성분을 $k$ 에 대해 그려 보자. 언제 최대가 되는가?
피보나치 행렬의 $\vv{s}_1$ 을 구해 보자. 성분의 비가 황금비인가?

다음 강의에서는 대각화를 무기 삼아 미분방정식을 푼다. $\lambda^k$ 자리에 $e^{\lambda t}$ 가 들어가면 그대로 통한다.

Lecture 22. 대각화와 A의 거듭제곱 — 파이썬 실습