Lecture 21. 고윳값과 고유벡터 — 파이썬 실습 - 선형대수

L21 서술 파트에서 $A\vv{x} = \lambda\vv{x}$ 를 만족하는 특별한 방향을 찾았고, 그 열쇠가 $\det(A - \lambda I) = 0$ 이었다.

이 노트북에서는 방향을 한 바퀴 돌려 가며 $\vv{x}$ 와 $A\vv{x}$ 의 각도를 실제로 재서 "안 돌아가는 방향"을 눈으로 찾고, 특성다항식을 직접 세워 np.linalg.eig 와 대조한다. 마지막으로 그 특성다항식을 실제 계산에 쓰면 왜 안 되는지도 확인한다.

서술 파트의 내용	여기서 확인하는 방법
안 돌아가는 방향	$\vv{x}$ 와 $A\vv{x}$ 의 각도를 훑는다
$\det(A - \lambda I) = 0$	특성다항식을 세워 근을 구한다
고유벡터 $=$ 영공간	`null_space(A - lambda I)`
합은 대각합, 곱은 행렬식	무작위 행렬로 검산
네 가지 사례	단위원이 어디로 가는지 그린다
$P^2 = P \Rightarrow \lambda \in \{0, 1\}$	계산 없이 나오는 것을 확인
특성방정식은 실전용이 아니다	근이 계수에 얼마나 민감한지

0. 준비¶

import numpy as np
import plotly.graph_objects as go
import sympy as sp
from scipy.linalg import null_space

from linalg_viz import COLORS, arrow2d, layout2d, show_matrix, slider_figure

np.set_printoptions(precision=3, suppress=True)
rng = np.random.default_rng(21)
print("numpy", np.__version__, " sympy", sp.__version__)

numpy 2.5.2  sympy 1.14.0

1. 안 돌아가는 방향을 찾는다¶

$\vv{x}$ 를 한 바퀴 돌려 가며 $\vv{x}$ 와 $A\vv{x}$ 사이의 각도를 재 보자. 각도가 0이거나 180도가 되는 자리가 고유방향이다.

def 사잇각(x, y):
    """두 벡터의 사잇각을 도 단위로 돌려준다."""
    코사인 = (x @ y) / (np.linalg.norm(x) * np.linalg.norm(y))
    return np.degrees(np.arccos(np.clip(코사인, -1.0, 1.0)))

A = np.array([[2.0, 1.0], [1.0, 2.0]])
각도들 = np.linspace(0, 360, 721)
잰각 = np.array([사잇각(v, A @ v) for v in
                (np.array([np.cos(np.radians(t)), np.sin(np.radians(t))])
                 for t in 각도들)])

print(f"가장 작은 사잇각 : {잰각.min():.3f} 도")
print("그때의 방향들 :", np.round(각도들[잰각 < 1e-6], 1), "도")

가장 작은 사잇각 : 0.000 도
그때의 방향들 : [ 45. 135. 225. 315.] 도

go.Figure(
    data=[go.Scatter(x=각도들, y=잰각, mode="lines",
                     line=dict(color=COLORS["input"], width=3),
                     name="x 와 Ax 의 사잇각")],
    layout=go.Layout(
        title=dict(text="사잇각이 0 이 되는 방향이 고유방향이다"),
        xaxis=dict(title=dict(text="x 의 방향 (도)"), dtick=45),
        yaxis=dict(title=dict(text="사잇각 (도)")),
        height=380, margin=dict(l=70, r=20, t=60, b=50)))

45도와 225도에서 사잇각이 0이다. 곧 $(1, 1)$ 방향이 고유방향이고, 그 반대 방향도 같은 직선이니 같은 고유방향이다. 135도와 315도에서도 마찬가지로 0이 되어야 하는데, $(1, -1)$ 은 $\lambda = 1 > 0$ 이라 방향이 그대로이므로 역시 0이다.

for 방향 in ([1.0, 1.0], [1.0, -1.0], [1.0, 0.0]):
    v = np.array(방향)
    Av = A @ v
    print(f"x = {v}   A x = {Av}   사잇각 {사잇각(v, Av):6.2f} 도", end="")
    print("   -> 고유방향" if 사잇각(v, Av) < 1e-4 else "")   # 반올림 오차만큼 여유

x = [1. 1.]   A x = [3. 3.]   사잇각   0.00 도   -> 고유방향
x = [ 1. -1.]   A x = [ 1. -1.]   사잇각   0.00 도   -> 고유방향
x = [1. 0.]   A x = [2. 1.]   사잇각  26.57 도

2. 특성다항식¶

$\det(A - \lambda I)$ 를 $\lambda$ 의 다항식으로 세운다. sympy 로 하면 그대로 보인다.

def 특성다항식(A):
    """det(A - lambda I) 를 sympy 다항식으로 돌려준다."""
    람 = sp.Symbol("lambda")
    M = sp.Matrix(np.asarray(A, dtype=float).tolist()).applyfunc(sp.nsimplify)
    return 람, sp.expand((M - 람 * sp.eye(M.shape[0])).det())

람, 식 = 특성다항식(A)
print("det(A - lambda I) =", 식)
print("인수분해           :", sp.factor(식))
print("근                 :", sp.solve(식, 람))
print()
print("np.linalg.eigvals  :", np.sort(np.linalg.eigvals(A).real))

det(A - lambda I) = lambda**2 - 4*lambda + 3
인수분해           : (lambda - 3)*(lambda - 1)
근                 : [1, 3]

np.linalg.eigvals  : [1. 3.]

$\lambda$ 를 실제로 넣어 행렬식이 0이 되는지도 확인하자.

print(f"{'lambda':>8}{'det(A - lambda I)':>20}{'rank':>7}{'dim N':>8}")
for l in (0.0, 1.0, 2.0, 3.0):
    M = A - l * np.eye(2)
    print(f"{l:>8.1f}{np.linalg.det(M):>20.4f}"
          f"{np.linalg.matrix_rank(M):>7}{null_space(M).shape[1]:>8}")

  lambda   det(A - lambda I)   rank   dim N
     0.0              3.0000      2       0
     1.0              0.0000      1       1
     2.0             -1.0000      2       0
     3.0              0.0000      1       1

3. 고유벡터는 영공간이다¶

$\lambda$ 를 알았으면 $(A - \lambda I)\vv{x} = \vv{0}$ 을 푸는 일만 남는다. L7에서 배운 영공간 구하기이다.

def 고유벡터(A, 람다):
    """N(A - lambda I) 의 기저를 돌려준다."""
    A = np.asarray(A, dtype=float)
    return null_space(A - 람다 * np.eye(A.shape[0]))

for l in (3.0, 1.0):
    V = 고유벡터(A, l)
    v = V[:, 0]
    print(f"lambda = {l:g}")
    print(f"  N(A - {l:g}I) 의 기저 : {v}")
    print(f"  정수로 보면          : {np.round(v / abs(v).min())}")
    print(f"  A v = {A @ v},   {l:g} v = {l * v}")
    print(f"  같은가 : {np.allclose(A @ v, l * v)}")
    print()

lambda = 3
  N(A - 3I) 의 기저 : [0.707 0.707]
  정수로 보면          : [1. 1.]
  A v = [2.121 2.121],   3 v = [2.121 2.121]
  같은가 : True

lambda = 1
  N(A - 1I) 의 기저 : [-0.707  0.707]
  정수로 보면          : [-1.  1.]
  A v = [-0.707  0.707],   1 v = [-0.707  0.707]
  같은가 : True

np.linalg.eig 는 이 둘을 한 번에 돌려준다.

값, 벡터 = np.linalg.eig(A)
print("고윳값 :", 값)
print(show_matrix(벡터.real, "고유벡터 (열마다 하나)"))
print()
for k in range(2):
    print(f"A v{k + 1} = {A @ 벡터[:, k]},   lambda{k + 1} v{k + 1} ="
          f" {값[k] * 벡터[:, k]},   같은가 {np.allclose(A @ 벡터[:, k], 값[k] * 벡터[:, k])}")

고윳값 : [3.+0.j 1.+0.j]
고유벡터 (열마다 하나)
[   0.707   -0.707 ]
[   0.707    0.707 ]

A v1 = [2.121+0.j 2.121+0.j],   lambda1 v1 = [2.121+0.j 2.121+0.j],   같은가 True
A v2 = [-0.707+0.j  0.707+0.j],   lambda2 v2 = [-0.707+0.j  0.707+0.j],   같은가 True

4. 검산 도구 두 개¶

고윳값의 합은 대각합이고 곱은 행렬식이다. 무작위 행렬로 확인하자.

print(f"{'n':>4}{'합':>12}{'trace':>12}{'곱':>14}{'det':>14}")
for _ in range(6):
    n = int(rng.integers(2, 6))
    M = rng.integers(-4, 5, (n, n)).astype(float)
    값들 = np.linalg.eigvals(M)
    print(f"{n:>4}{값들.sum().real:>12.4f}{np.trace(M):>12.4f}"
          f"{np.prod(값들).real:>14.4f}{np.linalg.det(M):>14.4f}")

   n           합       trace             곱           det
   3      6.0000      6.0000        4.0000        4.0000
   4      6.0000      6.0000      630.0000      630.0000
   2      6.0000      6.0000       11.0000       11.0000
   5     -6.0000     -6.0000      990.0000      990.0000
   2     -7.0000     -7.0000       21.0000       21.0000
   2     -2.0000     -2.0000        4.0000        4.0000

$2 \times 2$ 에서는 이 둘만으로 특성방정식을 풀지 않고 고윳값을 알 수 있다. 합이 $s$ 이고 곱이 $p$ 이면 두 수는 $t^2 - st + p = 0$ 의 근이다.

s, p = np.trace(A), np.linalg.det(A)
판별식 = s ** 2 - 4 * p
두근 = ((s + np.sqrt(판별식)) / 2, (s - np.sqrt(판별식)) / 2)

print(f"합 s = {s:g},  곱 p = {p:g}")
print(f"t^2 - {s:g}t + {p:g} = 0 의 근 : {두근}")
print("eigvals 와 같은가 :", np.allclose(np.sort(두근), np.sort(np.linalg.eigvals(A).real)))

합 s = 4,  곱 p = 3
t^2 - 4t + 3 = 0 의 근 : (np.float64(3.0), np.float64(0.9999999999999998))
eigvals 와 같은가 : True

5. 네 가지 사례¶

서술 파트의 네 행렬을 차례로 넣어 보자.

사례 = {
    "투영 (y = x 위로)": np.array([[0.5, 0.5], [0.5, 0.5]]),
    "삼각 [3 1 ; 0 2]": np.array([[3.0, 1.0], [0.0, 2.0]]),
    "회전 90 도":        np.array([[0.0, -1.0], [1.0, 0.0]]),
    "결함 [3 1 ; 0 3]": np.array([[3.0, 1.0], [0.0, 3.0]]),
}

print(f"{'':>18}{'고윳값':>26}{'독립 고유벡터':>14}")
for 이름, M in 사례.items():
    값들 = np.linalg.eigvals(M)
    개수 = np.linalg.matrix_rank(np.linalg.eig(M)[1])
    보기 = ", ".join(f"{v:.2f}" if abs(v.imag) < 1e-12 else f"{v:.2f}"
                    for v in 값들)
    print(f"{이름:>18}{보기:>26}{개수:>14}")

                                         고윳값       독립 고유벡터
     투영 (y = x 위로)    1.00+0.00j, 0.00+0.00j             2
    삼각 [3 1 ; 0 2]    3.00+0.00j, 2.00+0.00j             2
           회전 90 도    0.00+1.00j, 0.00-1.00j             2
    결함 [3 1 ; 0 3]    3.00+0.00j, 3.00+0.00j             1

회전행렬만 고윳값이 복소수이고, 결함 행렬만 독립 고유벡터가 하나뿐이다. 단위원이 어디로 가는지 슬라이더로 넘겨 가며 보자.

각 = np.linspace(0, 2 * np.pi, 200)
원 = np.array([np.cos(각), np.sin(각)])

고유선 = {
    "투영 (y = x 위로)": [([1, 1], 1.0), ([1, -1], 0.0)],
    "삼각 [3 1 ; 0 2]": [([1, 0], 3.0), ([1, -1], 2.0)],
    "회전 90 도": [],
    "결함 [3 1 ; 0 3]": [([1, 0], 3.0)],
}

프레임, 라벨 = [], []
for 이름, M in 사례.items():
    상 = M @ 원
    자료 = [
        go.Scatter(x=원[0], y=원[1], mode="lines", name="단위원",
                   line=dict(color=COLORS["nullspace"], dash="dash")),
        go.Scatter(x=상[0], y=상[1], mode="lines", name="그 상",
                   line=dict(color=COLORS["colspace"], width=3)),
    ]
    방향들 = 고유선[이름] + [([0, 0], 0.0)] * (2 - len(고유선[이름]))
    for (방향, 람다값), 색 in zip(방향들, (COLORS["input"], COLORS["second"])):
        v = np.asarray(방향, dtype=float)
        진짜 = float(np.linalg.norm(v)) > 0
        if 진짜:
            v = v / np.linalg.norm(v) * 2.2
        자료 += arrow2d([0, 0], v, 색, f"lambda = {람다값:g}" if 진짜 else "",
                      legend=진짜)
    프레임.append(자료)
    라벨.append(이름)

slider_figure(프레임, 라벨, layout2d("단위원이 어디로 가는가", extent=3.6),
              initial=0)

6. $P^2 = P$ 하나로 고윳값이 정해진다¶

서술 파트에서 계산을 한 줄도 하지 않고 투영행렬의 고윳값이 0과 1임을 보였다. 같은 논법이 다른 행렬에도 통한다.

def 무작위투영(m, r, 씨앗):
    """랭크 r 짜리 투영행렬을 만든다."""
    지역 = np.random.default_rng(씨앗)
    B = 지역.normal(size=(m, r))
    return B @ np.linalg.inv(B.T @ B) @ B.T


for m, r in ((5, 2), (6, 3), (7, 1)):
    P = 무작위투영(m, r, m * 10 + r)
    값들 = np.sort(np.linalg.eigvals(P).real)
    print(f"{m} x {m}, 랭크 {r} :  P^2 = P {np.allclose(P @ P, P)},"
          f"  고윳값 {np.round(값들, 6)}")

5 x 5, 랭크 2 :  P^2 = P True,  고윳값 [-0. -0.  0.  1.  1.]
6 x 6, 랭크 3 :  P^2 = P True,  고윳값 [-0. -0.  0.  1.  1.  1.]
7 x 7, 랭크 1 :  P^2 = P True,  고윳값 [-0. -0. -0.  0.  0.  0.  1.]

1이 랭크만큼, 0이 나머지만큼 나온다. $\lambda^2 = \lambda$ 이니 다른 값은 나올 수 없다.

같은 방법으로 다른 행렬도 다룰 수 있다. $Q^{\mathsf{T}}Q = I$ 인 직교행렬이라면 $\|Q\vv{x}\| = \|\vv{x}\|$ 이므로 $|\lambda| = 1$ 이어야 한다.

각도 = np.pi / 7
Q = np.array([[np.cos(각도), -np.sin(각도)], [np.sin(각도), np.cos(각도)]])
값들 = np.linalg.eigvals(Q)

print("고윳값 :", 값들)
print("절댓값 :", np.abs(값들), " -> 전부 1")
print("반사행렬도 :", np.abs(np.linalg.eigvals(np.array([[1.0, 0], [0, -1]]))))

고윳값 : [0.901+0.434j 0.901-0.434j]
절댓값 : [1. 1.]  -> 전부 1
반사행렬도 : [1. 1.]

7. 고유벡터의 스케일과 부호¶

고유벡터는 방향이라 상수배해도 여전히 고유벡터이다. 라이브러리가 주는 부호에 기대면 안 된다.

v = 고유벡터(A, 3.0)[:, 0]
for c in (1.0, 2.5, -1.0, 100.0):
    w = c * v
    print(f"c = {c:>6.1f} :  A w = {A @ w},   3 w = {3 * w},"
          f"   고유벡터인가 {np.allclose(A @ w, 3 * w)}")

c =    1.0 :  A w = [2.121 2.121],   3 w = [2.121 2.121],   고유벡터인가 True
c =    2.5 :  A w = [5.303 5.303],   3 w = [5.303 5.303],   고유벡터인가 True
c =   -1.0 :  A w = [-2.121 -2.121],   3 w = [-2.121 -2.121],   고유벡터인가 True
c =  100.0 :  A w = [212.132 212.132],   3 w = [212.132 212.132],   고유벡터인가 True

_, V1 = np.linalg.eig(A)
_, V2 = np.linalg.eig(A + 0.0)          # 같은 행렬을 다시
print("두 번 부른 결과가 같은가 :", np.allclose(V1, V2))
print()
print("길이는 1 로 맞춰 준다 :", np.linalg.norm(V1, axis=0))
print("부호는 정해 주지 않는다. 비교할 때는 부호를 맞추거나 방향만 비교할 것.")
print("  V1[:,0] 과 -V1[:,0] 이 같은 공간인가 :",
      np.linalg.matrix_rank(np.column_stack([V1[:, 0], -V1[:, 0]])) == 1)

두 번 부른 결과가 같은가 : True

길이는 1 로 맞춰 준다 : [1. 1.]
부호는 정해 주지 않는다. 비교할 때는 부호를 맞추거나 방향만 비교할 것.
  V1[:,0] 과 -V1[:,0] 이 같은 공간인가 : True

8. 특성방정식은 실전용이 아니다¶

$2 \times 2$ 에서는 특성방정식이 편하다. 그러나 차수가 커지면 계수를 거쳐 가는 것 자체가 위험해진다. 고윳값이 $1, 2, \dots, 15$ 인 행렬로 확인하자.

참값 = np.arange(1.0, 16.0)
계수 = np.poly(참값)                       # 특성다항식의 계수
대각 = np.diag(참값)

근 = np.sort(np.roots(계수).real)
직접 = np.sort(np.linalg.eigvals(대각).real)

print("다항식의 근으로   : 최대 오차", f"{np.abs(근 - 참값).max():.2e}")
print("eigvals 로        : 최대 오차", f"{np.abs(직접 - 참값).max():.2e}")
print()
print("이 크기에서는 둘 다 잘 나온다. 문제는 다른 데 있다.")

다항식의 근으로   : 최대 오차 9.24e-06
eigvals 로        : 최대 오차 0.00e+00

이 크기에서는 둘 다 잘 나온다. 문제는 다른 데 있다.

문제는 계수가 조금만 틀려도 근이 크게 움직인다는 것이다. 계수 하나를 10^-10 만큼만 건드려 보자. 행렬 쪽을 같은 만큼 건드린 것과 비교한다.

흔든계수 = 계수.copy()
흔든계수[1] *= (1 + 1e-10)                 # 계수 하나를 1e-10 만큼만
흔든근 = np.sort(np.roots(흔든계수).real)

흔든대각 = 대각.copy()
흔든대각[0, 0] *= (1 + 1e-10)              # 행렬 성분을 같은 만큼
흔든값 = np.sort(np.linalg.eigvals(흔든대각).real)

print(f"계수를 1e-10 흔들었을 때 근이 움직인 거리   : "
      f"{np.abs(흔든근 - 근).max():.3e}")
print(f"행렬을 1e-10 흔들었을 때 고윳값이 움직인 거리 : "
      f"{np.abs(흔든값 - 직접).max():.3e}")
print()
증폭 = np.abs(흔든근 - 근).max() / max(np.abs(흔든값 - 직접).max(), 1e-300)
print(f"계수 쪽이 {증폭:.0e} 배 더 민감하다.")

계수를 1e-10 흔들었을 때 근이 움직인 거리   : 6.479e-02
행렬을 1e-10 흔들었을 때 고윳값이 움직인 거리 : 1.000e-10

계수 쪽이 6e+08 배 더 민감하다.

print(f"{'참값':>6}{'흔든 뒤의 근':>16}{'움직인 거리':>14}")
for k in range(15):
    print(f"{참값[k]:>6.0f}{흔든근[k]:>16.4f}{abs(흔든근[k] - 참값[k]):>14.4f}")

    참값         흔든 뒤의 근        움직인 거리
     1          1.0000        0.0000
     2          2.0000        0.0000
     3          3.0000        0.0000
     4          4.0000        0.0000
     5          5.0000        0.0000
     6          6.0000        0.0000
     7          7.0003        0.0003
     8          7.9979        0.0021
     9          9.0096        0.0096
    10          9.9730        0.0270
    11         11.0538        0.0538
    12         11.9352        0.0648
    13         13.0482        0.0482
    14         13.9780        0.0220
    15         15.0040        0.0040

큰 근일수록 크게 흔들린다. 계수는 근들의 곱과 합으로 만들어지므로, 큰 근 쪽의 정보가 계수 안에서 서로 상쇄되어 묻히기 때문이다.

그래서 실제 소프트웨어는 특성방정식을 세우지 않는다. numpy.linalg.eig 는 행렬을 직접 다루는 반복법( $QR$ 알고리즘)을 쓴다. L19에서 만난 태도가 또 나온다. 옳은 공식과 쓰는 방법은 다르다.

마치며...¶

서술 파트의 내용	이 노트북의 코드
안 돌아가는 방향	`사잇각(v, A @ v)` 을 한 바퀴 훑는다
특성다항식	`특성다항식(A)` — sympy 로 그대로
고유벡터 $=$ 영공간	`고유벡터(A, l)` $=$ `null_space(A - lI)`
검산	합과 대각합, 곱과 행렬식
네 사례	단위원의 상을 슬라이더로
$P^2 = P$	랭크만큼의 1과 나머지 0
스케일 자유	상수배해도 고유벡터
특성방정식의 위험	근이 계수에 극도로 민감하다

더 해 볼 것¶

1절의 사잇각 곡선을 다른 행렬로 그려 보자. 회전행렬이면 어떤 모양인가? 결함 행렬이면?
고유벡터(A, l) 에 고윳값이 아닌 $\lambda$ 를 넣으면 무엇이 나오는가? 왜 그런가?
대칭행렬을 여러 개 만들어 고유벡터끼리 내적해 보자. 무엇이 눈에 띄는가? (이 관찰은 L25에서 정리로 승격된다.)
8절에서 고윳값을 $1, 2, \dots, n$ 으로 두고 $n$ 을 키워 보자. 몇 차부터 다항식의 근이 못 쓸 정도가 되는가?

다음 강의에서는 고유벡터를 기저로 삼는다. 좌표계를 갈아끼우는 순간 행렬 곱셈이 수 곱셈이 되고, $A^{100}$ 이 우스워진다.

Lecture 21. 고윳값과 고유벡터 — 파이썬 실습