Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Lecture 28. 유사 행렬과 조르당 표준형 — 파이썬 실습

Similar Matrices and Jordan Form — 실습

L28 서술 파트의 결론은 씁쓸했다. 조르당 형은 아름답지만 수치적으로는 존재하지 않는다.

이 노트북에서는 결함 행렬을 손으로 조르당 형까지 옮겨 보고, teλtte^{\lambda t} 가 어디서 나오는지 확인한 다음, 섭동을 슬라이더로 키워 가며 고윳값이 ϵ\sqrt\epsilon 만큼 갈라지는 것을 본다. 마지막으로 numpy 가 결함을 어떻게 놓치는지 확인한다.

서술 파트의 내용여기서 확인하는 방법
유사이면 고윳값이 같다무작위 300쌍
보존되는 것과 아닌 것대각합·랭크는 같고 대칭성은 깨진다
대수 vs 기하 중복도랭크로 센다
3I3I vs [[3,1],[0,3]][[3,1],[0,3]]갈리는 것은 랭크
일반화 고유벡터(A3I)v2=v1(A-3I)\vv{v}_2 = \vv{v}_1 을 직접 푼다
M1AM=JM^{-1}AM = J손으로 만든 MM 으로 확인
teλtte^{\lambda t}N2=0N^2 = 0 이라 급수가 끊긴다
ϵ\sqrt\epsilon 법칙슬라이더와 로그-로그 기울기
컴퓨터에는 안 보인다고유벡터 행렬의 조건수
대안슈어 분해는 늘 안정하다

0. 준비

import numpy as np
import plotly.graph_objects as go
from scipy.linalg import expm, schur

from linalg_viz import COLORS, layout2d, show_matrix, slider_figure

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

1. 유사이면 고윳값이 같다

어긋 = 0
for _ in range(300):
    n = int(rng.integers(2, 6))
    X = rng.normal(size=(n, n))
    M = rng.normal(size=(n, n))
    if abs(np.linalg.det(M)) < 1e-6:
        continue
    a = np.sort_complex(np.linalg.eigvals(X))
    b = np.sort_complex(np.linalg.eigvals(np.linalg.inv(M) @ X @ M))
    if not np.allclose(a, b, atol=1e-7):
        어긋 += 1
print("무작위 300쌍 중 고윳값이 다른 경우 :", 어긋, "건")
무작위 300쌍 중 고윳값이 다른 경우 : 0 건

무엇이 보존되고 무엇이 안 되는지 눈으로 보자. 대칭행렬을 유사변환하면 대칭성이 깨진다.

S = np.array([[2.0, 1.0], [1.0, 2.0]])              # 대칭
M = np.array([[1.0, 2.0], [0.0, 1.0]])              # 직교가 아니다
T = np.linalg.inv(M) @ S @ M
print(show_matrix(S, "S  (대칭)"))
print(show_matrix(T, "M^-1 S M"))
print(f"{'':>16}{'S':>14}{'M^-1 S M':>14}")
for 이름, 값 in (("고윳값", lambda X: np.sort(np.linalg.eigvals(X).real)),
                ("대각합", np.trace),
                ("행렬식", np.linalg.det),
                ("랭크", np.linalg.matrix_rank)):
    print(f"{이름:>16}{str(np.round(값(S), 4)):>14}{str(np.round(값(T), 4)):>14}")
print(f"{'대칭인가':>16}{str(np.allclose(S, S.T)):>14}{str(np.allclose(T, T.T)):>14}")
S  (대칭)
[  2   1 ]
[  1   2 ]
M^-1 S M
[   0   -3 ]
[   1    4 ]
                             S      M^-1 S M
             고윳값       [1. 3.]       [1. 3.]
             대각합           4.0           4.0
             행렬식           3.0           3.0
              랭크             2             2
            대칭인가          True         False

2. 두 가지 중복도

대수적 중복도는 특성다항식의 근으로서 몇 번인지이고, 기하적 중복도N(AλI)\Nul(A - \lambda I) 의 차원이다.

def 중복도(A, lam, 눈감아=1e-8):
    """(대수적, 기하적) 중복도를 돌려준다."""
    A = np.asarray(A, dtype=float)
    n = A.shape[0]
    대수 = int(np.sum(np.abs(np.linalg.eigvals(A) - lam) < 눈감아))
    기하 = n - np.linalg.matrix_rank(A - lam * np.eye(n), tol=눈감아)
    return 대수, 기하


def 결함인가(A, 눈감아=1e-8):
    """어떤 고윳값에서든 기하 < 대수 이면 결함 행렬이다."""
    값 = np.linalg.eigvals(A)
    본것 = []
    for l in 값:
        if abs(l.imag) > 눈감아 or any(abs(l.real - m) < 눈감아 for m in 본것):
            continue
        본것.append(l.real)
        대수, 기하 = 중복도(A, l.real, 눈감아)
        if 기하 < 대수:
            return True
    return False
경우 = (("3I", 3.0 * np.eye(2)),
        ("[[3,1],[0,3]]", np.array([[3.0, 1.0], [0.0, 3.0]])),
        ("[[5,4],[-1,1]]", np.array([[5.0, 4.0], [-1.0, 1.0]])),
        ("[[2,1],[1,2]]", np.array([[2.0, 1.0], [1.0, 2.0]])))
print(f"{'':>18}{'고윳값':>18}{'대수':>6}{'기하':>6}{'A-lI 랭크':>12}{'결함':>8}")
for 이름, A0 in 경우:
    값 = np.linalg.eigvals(A0).real
    l = 값[0]
    대수, 기하 = 중복도(A0, l)
    랭크 = np.linalg.matrix_rank(A0 - l * np.eye(2), tol=1e-8)
    print(f"{이름:>18}{str(np.round(np.sort(값), 3)):>18}{대수:>6}{기하:>6}"
          f"{랭크:>12}{str(결함인가(A0)):>8}")
                                 고윳값    대수    기하     A-lI 랭크      결함
                3I           [3. 3.]     2     2           0   False
     [[3,1],[0,3]]           [3. 3.]     2     1           1    True
    [[5,4],[-1,1]]           [3. 3.]     0     1           1   False
     [[2,1],[1,2]]           [1. 3.]     1     1           1   False

3I3I[[3,1],[0,3]][[3,1],[0,3]] 은 특성다항식이 같은데 운명이 갈린다. 갈리는 것은 고윳값이 아니라 AλIA - \lambda I 의 랭크이다.

그런데 세 번째 줄이 이상하다. [[5,4],[1,1]][[5,4],[-1,1]](λ3)2(\lambda-3)^2 이라 대수적 중복도가 2여야 하는데 0 이라고 나온다. 함수가 틀린 것이 아니다.

A0 = np.array([[5.0, 4.0], [-1.0, 1.0]])
값 = np.linalg.eigvals(A0)
print("numpy 가 준 고윳값 :")
for v in 값:
    print(f"   {v:.17g}")
print()
print("3 에서 떨어진 거리 :", [f"{d:.4e}" for d in np.abs(값 - 3.0)])
print()
print("det A =", f"{np.linalg.det(A0):.20g}", "   (정확한 값은 9)")
print("판별식 36 - 4 det =", f"{36 - 4*np.linalg.det(A0):.4e}",
      "   (정확한 값은 0)")
numpy 가 준 고윳값 :
   3+2.9802322387695312e-08j
   3-2.9802322387695312e-08j

3 에서 떨어진 거리 : ['2.9802e-08', '2.9802e-08']

det A = 8.9999999999999982236    (정확한 값은 9)
판별식 36 - 4 det = 7.1054e-15    (정확한 값은 0)
print("눈감아 주는 폭을 sqrt(eps) 보다 크게 잡으면 제대로 센다.")
print(f"{'눈감아':>10}{'대수':>6}{'기하':>6}{'결함':>8}")
for 눈 in (1e-10, 1e-8, 1e-6, 1e-4):
    대수, 기하 = 중복도(A0, 3.0, 눈)
    print(f"{눈:>10.0e}{대수:>6}{기하:>6}{str(결함인가(A0, 눈)):>8}")
print()
print("-> 무엇을 '같다' 고 볼지 사람이 정해 주어야 한다.")
print("   조르당 형을 컴퓨터에 맡길 수 없는 이유가 이것이다.")
눈감아 주는 폭을 sqrt(eps) 보다 크게 잡으면 제대로 센다.
       눈감아    대수    기하      결함
     1e-10     0     1   False
     1e-08     0     1   False
     1e-06     2     1    True
     1e-04     2     1    True

-> 무엇을 '같다' 고 볼지 사람이 정해 주어야 한다.
   조르당 형을 컴퓨터에 맡길 수 없는 이유가 이것이다.

3. 앵커를 조르당 형으로

A=[5411]A = \begin{bmatrix} 5 & 4 \\ -1 & 1\end{bmatrix} 을 손으로 옮겨 보자.

A = np.array([[5.0, 4.0], [-1.0, 1.0]])
print("대각합 :", np.trace(A), "  행렬식 :", np.linalg.det(A),
      "  -> (lambda - 3)^2")
print("고윳값 :", np.linalg.eigvals(A))
print()
print(show_matrix(A - 3 * np.eye(2), "A - 3I"))
print("랭크 :", np.linalg.matrix_rank(A - 3 * np.eye(2)),
      " -> 영공간이 1차원. 고유벡터가 하나뿐")
대각합 : 6.0   행렬식 : 8.999999999999998   -> (lambda - 3)^2
고윳값 : [3.+0.j 3.-0.j]

A - 3I
[   2    4 ]
[  -1   -2 ]
랭크 : 1  -> 영공간이 1차원. 고유벡터가 하나뿐
v1 = np.array([2.0, -1.0])
print("(A-3I) v1 =", (A - 3 * np.eye(2)) @ v1, "  -> v1 은 고유벡터")

# (A - 3I) v2 = v1 을 푼다. 특이행렬이라 lstsq 를 쓴다.
v2, *_ = np.linalg.lstsq(A - 3 * np.eye(2), v1, rcond=None)
print("lstsq 가 준 v2 :", np.round(v2, 6))
v2 = np.array([1.0, 0.0])                            # 서술 파트에서 고른 것
print("서술 파트의 v2 :", v2)
print("(A-3I) v2 =", (A - 3 * np.eye(2)) @ v2, " = v1 인가 :",
      np.allclose((A - 3 * np.eye(2)) @ v2, v1))
(A-3I) v1 = [0. 0.]   -> v1 은 고유벡터
lstsq 가 준 v2 : [0.2 0.4]
서술 파트의 v2 : [1. 0.]
(A-3I) v2 = [ 2. -1.]  = v1 인가 : True
M = np.column_stack([v1, v2])                        # 고유벡터가 먼저
print(show_matrix(M, "M = [v1  v2]"))
print("det M =", np.linalg.det(M))
J = np.linalg.inv(M) @ A @ M
print(show_matrix(J, "M^-1 A M"))
print("[[3,1],[0,3]] 인가 :", np.allclose(J, [[3, 1], [0, 3]]))
print()
print("v2 에 v1 을 더해도 같은 결과가 나오는가 :")
M2 = np.column_stack([v1, v2 + 7 * v1])
print("  ", np.allclose(np.linalg.inv(M2) @ A @ M2, [[3, 1], [0, 3]]))
M = [v1  v2]
[   2    1 ]
[  -1    0 ]
det M = 1.0
M^-1 A M
[  3   1 ]
[  0   3 ]
[[3,1],[0,3]] 인가 : True

v2 에 v1 을 더해도 같은 결과가 나오는가 :
   True
print("서술 파트의 열별 읽기 확인")
print("  A v1 =", A @ v1, " = 3 v1 =", 3 * v1)
print("  A v2 =", A @ v2, " = 3 v2 + v1 =", 3 * v2 + v1)
print()
print("두 계수 (3,0) 과 (1,3) 이 J 의 두 열이다 :")
print(show_matrix(J, "J"))
서술 파트의 열별 읽기 확인
  A v1 = [ 6. -3.]  = 3 v1 = [ 6. -3.]
  A v2 = [ 5. -1.]  = 3 v2 + v1 = [ 5. -1.]

두 계수 (3,0) 과 (1,3) 이 J 의 두 열이다 :
J
[  3   1 ]
[  0   3 ]

4. teλtte^{\lambda t} 는 어디서 오는가

J=λI+NJ = \lambda I + N 이고 N2=0N^2 = 0 이라 급수가 두 항에서 끊긴다.

lam = 3.0
N = np.array([[0.0, 1.0], [0.0, 0.0]])
print(show_matrix(N @ N, "N^2  <- 사라진다"))
print()
J = lam * np.eye(2) + N
for t in (0.0, 0.5, 1.0, 2.0):
    손 = np.exp(lam * t) * np.array([[1.0, t], [0.0, 1.0]])
    print(f"t={t:>4} :  expm(Jt) 와 e^(3t)[[1,t],[0,1]] 의 최대 차이 "
          f"{np.abs(expm(J * t) - 손).max():.2e}")
N^2  <- 사라진다
[  0   0 ]
[  0   0 ]

t= 0.0 :  expm(Jt) 와 e^(3t)[[1,t],[0,1]] 의 최대 차이 0.00e+00
t= 0.5 :  expm(Jt) 와 e^(3t)[[1,t],[0,1]] 의 최대 차이 0.00e+00
t= 1.0 :  expm(Jt) 와 e^(3t)[[1,t],[0,1]] 의 최대 차이 3.88e-12
t= 2.0 :  expm(Jt) 와 e^(3t)[[1,t],[0,1]] 의 최대 차이 0.00e+00
# 크기 3 짜리 블록이면 t^2 까지
J3 = 2.0 * np.eye(3) + np.diag([1.0, 1.0], 1)
print(show_matrix(J3, "J_3(2)"))
N3 = J3 - 2 * np.eye(3)
print("N^2 =\n", N3 @ N3, "\nN^3 =\n", N3 @ N3 @ N3, "  <- 여기서 끊긴다")
t = 1.3
손 = np.exp(2 * t) * np.array([[1, t, t*t/2], [0, 1, t], [0, 0, 1]])
print(f"\nt={t} : expm 과 손계산의 최대 차이 {np.abs(expm(J3*t) - 손).max():.2e}")
J_3(2)
[  2   1   0 ]
[  0   2   1 ]
[  0   0   2 ]
N^2 =
 [[0. 0. 1.]
 [0. 0. 0.]
 [0. 0. 0.]] 
N^3 =
 [[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]   <- 여기서 끊긴다

t=1.3 : expm 과 손계산의 최대 차이 2.90e-12

미분방정식의 해를 직접 보자. dudt=Ju\frac{d\vv{u}}{dt} = J\vv{u} 에서 첫 성분에 tt 가 붙는다.

시각 = np.linspace(0, 2.5, 300)
u0 = np.array([0.0, 1.0])
경로 = np.array([expm(np.array([[3.0, 1.0], [0.0, 3.0]]) * t) @ u0 for t in 시각])
공식1 = np.exp(3 * 시각) * (u0[0] + 시각 * u0[1])
print("u1 = e^(3t)(u1(0) + t u2(0)) 와의 최대 차이 :",
      np.abs(경로[:, 0] - 공식1).max())

go.Figure(
    data=[go.Scatter(x=시각, y=경로[:, 0], mode="lines", name="u1 — t e^(3t) 가 섞였다",
                     line=dict(color=COLORS["output"], width=4)),
          go.Scatter(x=시각, y=경로[:, 1], mode="lines", name="u2 — 순수한 e^(3t)",
                     line=dict(color=COLORS["input"], width=3, dash="dash")),
          go.Scatter(x=시각, y=시각 * np.exp(3 * 시각), mode="lines",
                     name="t e^(3t) 만", line=dict(color="#999999", width=2, dash="dot"))],
    layout=go.Layout(title=dict(text="대각 위의 1 이 t 를 만든다"),
                     xaxis=dict(title=dict(text="t")),
                     yaxis=dict(title=dict(text="값"), type="log"),
                     height=420, margin=dict(l=70, r=20, t=60, b=50)))
u1 = e^(3t)(u1(0) + t u2(0)) 와의 최대 차이 : 3.2315483622369356e-11
Loading...

5. 아주 작게 건드리면 무너진다

Bϵ=[31ϵ3]λ=3±ϵB_\epsilon = \begin{bmatrix} 3 & 1 \\ \epsilon & 3 \end{bmatrix} \qquad\Longrightarrow\qquad \lambda = 3 \pm \sqrt\epsilon
B = np.array([[3.0, 1.0], [0.0, 3.0]])
print(f"{'eps':>10}{'고윳값':>34}{'이동':>14}{'sqrt(eps)':>14}{'결함인가':>10}")
for eps in (0.0, 1e-12, 1e-8, 1e-4, 1e-2):
    Bp = B + np.array([[0.0, 0.0], [eps, 0.0]])
    값 = np.linalg.eigvals(Bp)
    이동 = float(np.abs(값 - 3).max())
    보기 = ", ".join(f"{v.real:.8f}" for v in np.sort_complex(값))
    print(f"{eps:>10.0e}{보기:>34}{이동:>14.2e}"
          f"{np.sqrt(eps):>14.2e}{str(결함인가(Bp)):>10}")
       eps                               고윳값            이동     sqrt(eps)      결함인가
     0e+00            3.00000000, 3.00000000      0.00e+00      0.00e+00      True
     1e-12            2.99999900, 3.00000100      1.00e-06      1.00e-06     False
     1e-08            2.99990000, 3.00010000      1.00e-04      1.00e-04     False
     1e-04            2.99000000, 3.01000000      1.00e-02      1.00e-02     False
     1e-02            2.90000000, 3.10000000      1.00e-01      1.00e-01     False

ϵ\epsilon10-12 만 되어도 결함이 사라진다. 고윳값이 서로 달라지므로 대각화가 가능해진다. 결함 행렬은 ϵ=0\epsilon = 0 이라는 한 점에만 있다.

정상 = np.array([[2.0, 1.0], [1.0, 2.0]])            # 고윳값 1, 3
엡실론 = 10.0 ** np.arange(-14, -1.4, 0.5)
결함이동, 정상이동 = [], []
for e in 엡실론:
    Bp = B + np.array([[0.0, 0.0], [e, 0.0]])
    결함이동.append(float(np.abs(np.linalg.eigvals(Bp) - 3.0).max()))
    Np = 정상 + np.array([[0.0, 0.0], [e, 0.0]])
    정상이동.append(float(np.abs(np.sort(np.linalg.eigvals(Np).real)
                              - np.sort(np.linalg.eigvals(정상).real)).max()))
결함이동, 정상이동 = np.array(결함이동), np.array(정상이동)

기울기 = lambda x, y: float(np.polyfit(np.log(x), np.log(y), 1)[0])
print(f"결함 행렬의 기울기       : {기울기(엡실론, 결함이동):.3f}   (이론값 0.5)")
print(f"대각화 가능한 것의 기울기 : {기울기(엡실론, 정상이동):.3f}   (이론값 1)")
print()
print(f"{'eps':>10}{'결함':>14}{'정상':>14}{'몇 배':>14}")
for e, a, b in zip(엡실론, 결함이동, 정상이동):
    if abs(np.log10(e) - round(np.log10(e))) < 1e-9 and round(np.log10(e)) % 4 == 0:
        print(f"{e:>10.0e}{a:>14.2e}{b:>14.2e}{a/b:>14,.0f}")
결함 행렬의 기울기       : 0.500   (이론값 0.5)
대각화 가능한 것의 기울기 : 1.000   (이론값 1)

       eps            결함            정상           몇 배
     1e-12      1.00e-06      5.00e-13     1,999,822
     1e-08      1.00e-04      5.00e-09        20,000
     1e-04      1.00e-02      5.00e-05           200

섭동이 작을수록 상대적 피해가 커진다. 보통 기대하는 것과 정반대이다.

크기들 = np.round(np.linspace(-3.5, -0.5, 13), 2)     # log10(eps)
각 = np.linspace(0, 2 * np.pi, 200)
프레임 = []
for lg in 크기들:
    eps = 10.0 ** lg
    위, 아래 = [], []
    for _ in range(500):
        E = rng.normal(size=(2, 2)) * eps
        값 = np.linalg.eigvals(B + E)
        (위 if E[1, 0] > 0 else 아래).extend(값)
    위, 아래 = np.array(위), np.array(아래)
    프레임.append([
        go.Scatter(x=위.real, y=위.imag, mode="markers",
                   marker=dict(color=COLORS["output"], size=4, opacity=0.5),
                   name="E21 > 0 : 실수축으로"),
        go.Scatter(x=아래.real, y=아래.imag, mode="markers",
                   marker=dict(color=COLORS["input"], size=4, opacity=0.5),
                   name="E21 < 0 : 허수축으로"),
        go.Scatter(x=3 + np.sqrt(eps) * np.cos(각), y=np.sqrt(eps) * np.sin(각),
                   mode="lines", line=dict(color="#d62728", width=2, dash="dash"),
                   name="반지름 sqrt(eps)"),
    ])

배치 = layout2d("섭동 크기를 키우면", extent=0.45)
배치["height"] = 560
배치["xaxis"] = dict(title=dict(text="Re"), range=[2.55, 3.45], scaleanchor="y")
배치["yaxis"] = dict(title=dict(text="Im"), range=[-0.45, 0.45])
slider_figure(프레임, 크기들, 배치, prefix="log10(eps) = ", initial=6)
Loading...

6. 컴퓨터에는 보이지 않는다

print(f"{'':>18}{'고유벡터 행렬의 조건수':>26}")
for 이름, X in (("[[2,1],[1,2]]", 정상),
                ("[[5,4],[-1,1]]", A),
                ("[[3,1],[0,3]]", B)):
    _, V = np.linalg.eig(X)
    print(f"{이름:>18}{np.linalg.cond(V):>26.3g}")
print()
print("S^-1 을 쓰는 모든 계산이 이 조건수만큼 오차를 증폭한다.")
                                고유벡터 행렬의 조건수
     [[2,1],[1,2]]                         1
    [[5,4],[-1,1]]                  1.68e+08
     [[3,1],[0,3]]                     3e+15

S^-1 을 쓰는 모든 계산이 이 조건수만큼 오차를 증폭한다.
값, V = np.linalg.eig(B)
print("numpy 가 준 B 의 고유벡터 :\n", np.round(V.real, 6))
print("두 열의 사잇각 :",
      round(float(np.degrees(np.arccos(min(1.0, abs(V[:, 0] @ V[:, 1]).real))))), "도")
print("-> 두 고유벡터가 사실상 같은 방향이다. 독립이 아니다.")
print()
print("배정밀도의 상대오차는 약 2.2e-16 이다.")
print("  그 크기의 섭동만으로 고윳값이", f"{np.sqrt(2.2e-16):.2e}", "만큼 흔들린다.")
print("  = 여덟 자리를 입력 단계에서 잃는다.")
numpy 가 준 B 의 고유벡터 :
 [[ 1. -1.]
 [ 0.  0.]]
두 열의 사잇각 : 0 도
-> 두 고유벡터가 사실상 같은 방향이다. 독립이 아니다.

배정밀도의 상대오차는 약 2.2e-16 이다.
  그 크기의 섭동만으로 고윳값이 1.48e-08 만큼 흔들린다.
  = 여덟 자리를 입력 단계에서 잃는다.

7. 대안 — 슈어 분해

A=QTQHA = QTQ^{\mathsf H} 로 쓴다. QQ 는 유니타리, TT 는 상삼각이다. 모든 행렬에 대해 존재하고 수치적으로 안정하다.

T, Q = schur(A)
print(show_matrix(Q, "Q  (직교)"))
print(show_matrix(T, "T  (상삼각).  대각이 고윳값"))
print("Q^T Q = I 인가 :", np.allclose(Q.T @ Q, np.eye(2)))
print("Q T Q^T = A 인가 :", np.allclose(Q @ T @ Q.T, A))
print("Q 의 조건수 :", np.linalg.cond(Q), "  <- 직교라 언제나 1")
Q  (직교)
[   0.894    0.447 ]
[  -0.447    0.894 ]
T  (상삼각).  대각이 고윳값
[  3   5 ]
[  0   3 ]
Q^T Q = I 인가 : True
Q T Q^T = A 인가 : True
Q 의 조건수 : 1.0   <- 직교라 언제나 1
print("결함 행렬에도 통하는가")
for 이름, X in (("[[3,1],[0,3]]", B), ("[[5,4],[-1,1]]", A),
                ("J_3(2)", J3), ("무작위 5x5", rng.normal(size=(5, 5)))):
    T, Q = schur(X)
    print(f"  {이름:>16} : Q 조건수 {np.linalg.cond(Q):>8.4f},  "
          f"복원 오차 {np.abs(Q @ T @ Q.T - X).max():.2e}")
결함 행렬에도 통하는가
     [[3,1],[0,3]] : Q 조건수   1.0000,  복원 오차 0.00e+00
    [[5,4],[-1,1]] : Q 조건수   1.0000,  복원 오차 8.88e-16
            J_3(2) : Q 조건수   1.0000,  복원 오차 0.00e+00
           무작위 5x5 : Q 조건수   1.0000,  복원 오차 3.55e-15

QQ 의 조건수가 언제나 1이다. 유니타리라 길이를 보존하기 때문이다(L26). 오차가 증폭될 자리가 없다.

마치며...

서술 파트의 내용이 노트북의 코드
유사이면 고윳값 공유무작위 300쌍에서 어긋난 것 0건
대칭성은 안 보존M1SMM^{-1}SM 이 대칭이 아니다
두 중복도랭크로 센다. 3I3I 는 0, BB 는 1
일반화 고유벡터유일하지 않다. 어느 것을 골라도 JJ 는 같다
M1AM=JM^{-1}AM = J손으로 만든 MM 으로 정확히
열별 읽기Av2=3v2+v1A\vv{v}_2 = 3\vv{v}_2 + \vv{v}_1
teλtte^{\lambda t}N2=0N^2 = 0, 급수가 두 항에서 끊긴다
J3J_3 이면 t2t^2N3=0N^3 = 0 까지
결함은 한 점에만ϵ=1012\epsilon = 10^{-12} 에도 사라진다
ϵ\sqrt\epsilon로그-로그 기울기 0.5 대 1
슬라이더섭동을 키우면 십자가 커진다
조건수1 / 108 / 1015
슈어QQ 의 조건수가 언제나 1

더 해 볼 것

  1. 3절의 v2\vv{v}_2 를 아무렇게나 바꿔 보자. v1\vv{v}_1 의 배수를 더하면 JJ 가 그대로인가? v2\vv{v}_2 자체를 상수배하면 어떻게 되는가?

  2. J3(2)J_3(2) 에 섭동을 주면 고윳값이 몇 제곱근만큼 움직이는가? ϵ\sqrt\epsilon 인가 ϵ1/3\epsilon^{1/3} 인가?

  3. L23의 동반행렬 [2110]\begin{bmatrix} -2 & -1 \\ 1 & 0\end{bmatrix} 이 결함 행렬인가? 조르당 형으로 옮겨 보고 tette^{-t} 가 나오는지 확인해 보자.

  4. scipy.linalg.schur 로 얻은 TT 의 대각이 정말 고윳값인가? 결함 행렬에서도 그런가?

다음 강의는 이 교재의 정점이다. 조건이 하나도 없는 분해를 만난다.