Lecture 34. 전체 정리 — 다섯 개의 분해 — 파이썬 실습

import numpy as np
import plotly.graph_objects as go
from scipy.linalg import lu

from linalg_viz import COLORS, arrow, layout3d, line, plane, show_matrix, slider_figure

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

# 이 강의에서 계속 쓸 두 행렬
T = np.array([[2.0, -1.0, 0.0],
              [-1.0, 2.0, -1.0],
              [0.0, -1.0, 2.0]])          # L27 의 삼중대각행렬
A = np.array([[1.0, 1.0, 2.0],
              [1.0, 0.0, 1.0],
              [0.0, 1.0, 1.0]])           # L29 의 C 에 열 하나를 더한 것
print(show_matrix(T, "T   다섯을 전부 통과한다"))
print(show_matrix(A, "A   마지막 그림의 주인공. 3열 = 1열 + 2열"))
print("A 의 앞 두 열이 L29 의 앵커 C 인가 :",
      np.allclose(A[:, :2], np.array([[1.0, 1], [1, 0], [0, 1]])))

def 다섯분해(M, 이름="A", 출력=True):
    """같은 행렬을 다섯 가지로 쪼개고, 각각이 무엇을 알려주는지 한눈에 본다.

    돌려주는 것은 {분해이름: (가능한가, 설명, 조각들)} 형태의 사전이다.
    """
    M = np.asarray(M, dtype=float)
    m, n = M.shape
    r = np.linalg.matrix_rank(M)
    정방 = m == n
    보고 = {}

    # 1) A = LU  —  정방이고 행 교환 없이 소거가 되어야 한다.
    #    선행 주소행렬식으로 재면 안 된다. 그것은 피벗이 n 개 다 잡히는 것까지
    #    요구해서 [[1,2],[2,4]] 같은 특이행렬을 떨어뜨린다. 그 행렬은 교환 없이
    #    소거되고 L @ U 가 정확히 원래 행렬이다. 그래서 직접 소거해 본다.
    if not 정방:
        보고["A = LU"] = (False, "정방이 아니다", ())
    else:
        허용 = 1e-11 * max(1.0, float(np.abs(M).max()))
        Ue, Le, 걸린곳 = M.copy(), np.eye(n), 0
        for k in range(n):
            if abs(Ue[k, k]) <= 허용:
                if np.all(np.abs(Ue[k + 1:, k]) <= 허용):
                    continue                  # 바꿀 행이 없다. 자유열이니 넘어간다
                걸린곳 = k + 1
                break
            곱수 = Ue[k + 1:, k] / Ue[k, k]
            Le[k + 1:, k] = 곱수
            Ue[k + 1:, k:] -= np.outer(곱수, Ue[k, k:])
        피벗 = np.diag(Ue)
        곱 = float(np.prod(피벗))
        곱 = 0.0 if 곱 == 0.0 else 곱          # -0.0 이 찍히는 것을 막는다
        if 걸린곳:
            보고["A = LU"] = (False, f"{걸린곳}번째 피벗 자리가 0이라 행을 바꿔야 한다",
                              lu(M))
        else:
            꼬리 = ("" if float(np.abs(피벗).min()) > 허용
                    else "  (0 인 피벗이 있어 후진 대입은 못 한다)")
            보고["A = LU"] = (True,
                              f"피벗 {np.round(피벗, 4)}, 곱 {곱:.4f}{꼬리}",
                              (np.eye(n), Le, Ue))

    # 2) A = QR  —  열이 독립이어야 한다
    Q, R = np.linalg.qr(M)
    if r == n:
        보고["A = QR"] = (True, f"|det R| = {abs(np.prod(np.diag(R))):.4f}", (Q, R))
    else:
        보고["A = QR"] = (False, f"열이 종속 (랭크 {r} < 열 {n})", (Q, R))

    # 3) S = Q L Q^T  —  대칭이어야 한다
    if 정방 and np.allclose(M, M.T):
        w, V = np.linalg.eigh(M)
        순 = np.argsort(w)[::-1]
        부호 = "전부 양수" if w.min() > 1e-12 else ("전부 음이 아님" if w.min() > -1e-12
                                                  else "부호가 섞임")
        보고["S = Q L Q^T"] = (True, f"lambda {np.round(w[순], 4)}, {부호}",
                               (V[:, 순], w[순]))
    else:
        보고["S = Q L Q^T"] = (False, "대칭이 아니다" if 정방 else "정방이 아니다", ())

    # 4) A = S L S^-1  —  독립인 고유벡터가 n 개여야 한다
    if 정방:
        w2, S = np.linalg.eig(M)
        cond = np.linalg.cond(S)
        복소 = np.max(np.abs(w2.imag)) > 1e-9
        if cond < 1e8:
            꼬리 = " (복소수에서만)" if 복소 else ""
            보고["A = S L S^-1"] = (True, f"cond(S) = {cond:.4g}{꼬리}", (S, w2))
        else:
            보고["A = S L S^-1"] = (False, f"고유벡터가 모자란다 (cond(S) = {cond:.2e})",
                                    (S, w2))
    else:
        보고["A = S L S^-1"] = (False, "정방이 아니다", ())

    # 5) A = U S V^T  —  요구하는 것이 없다
    U_, s, Vt = np.linalg.svd(M)
    조건수 = s[0] / s[-1] if s[-1] > 1e-14 else np.inf
    보고["A = U S V^T"] = (True, f"sigma {np.round(s, 4)}, cond {조건수:.4g}",
                           (U_, s, Vt))

    if 출력:
        print(f"{이름}   ({m}x{n}, 랭크 {r}, "
              f"{'대칭' if 정방 and np.allclose(M, M.T) else '비대칭'})")
        print("-" * 78)
        for 줄, (됨, 설명, _) in 보고.items():
            print(f"  {줄:<14}{'O' if 됨 else 'X'}   {설명}")
        통과 = sum(1 for 됨, _, _ in 보고.values() if 됨)
        print(f"  {'':<14}     -> 다섯 중 {통과}개 통과")
    return 보고

보T = 다섯분해(T, "T = tridiag(-1, 2, -1)")

보A = 다섯분해(A, "A = [[1,1,2],[1,0,1],[0,1,1]]")

P, L, U = 보T["A = LU"][2]
print("T 의 피벗 :", np.diag(U), "  서술의 2, 1.5, 4/3 인가 :",
      np.allclose(np.diag(U), [2.0, 1.5, 4/3]))
print("행 교환이 있었는가 :", not np.allclose(P, np.eye(3)),
      "  L @ U 가 T 인가 :", np.allclose(L @ U, T))
print("피벗의 곱 :", np.prod(np.diag(U)), "  det T :", np.linalg.det(T))
print()
Q, R = 보T["A = QR"][2]
print("|det R| :", abs(np.prod(np.diag(R))), "  <- 직교행렬은 부피를 안 바꾼다")
print()
_, s, _ = 보T["A = U S V^T"][2]
print("T 의 특이값 :", s)
print("2+sqrt2, 2, 2-sqrt2 :", [2 + np.sqrt(2), 2.0, 2 - np.sqrt(2)])
print("cond(T) :", np.linalg.cond(T), "  3 + 2 sqrt2 :", 3 + 2*np.sqrt(2))
print("sigma 의 곱 :", np.prod(s), "  |det T| :", abs(np.linalg.det(T)))

def 부호맞춤(M):
    """열마다 0이 아닌 첫 성분이 양수가 되도록 부호를 고른다."""
    M = np.array(M, dtype=float)
    for j in range(M.shape[1]):
        큰것 = np.flatnonzero(np.abs(M[:, j]) > 1e-8)
        if 큰것.size and M[큰것[0], j] < 0:
            M[:, j] *= -1
    return M

Qs, 람 = 보T["S = Q L Q^T"][2]
S2, 람2 = 보T["A = S L S^-1"][2]
순2 = np.argsort(람2.real)[::-1]
Uu, 시그, Vt = 보T["A = U S V^T"][2]

Qs, S2, Uu, Vv = (부호맞춤(Qs), 부호맞춤(S2[:, 순2].real),
                  부호맞춤(Uu), 부호맞춤(Vt.T))
print(show_matrix(Qs, "Q   (eigh)"))
print(show_matrix(S2, "S   (eig)"))
print(show_matrix(Uu, "U   (svd)"))
print(show_matrix(Vv, "V   (svd)"))
print("Q == S :", np.allclose(Qs, S2))
print("Q == U :", np.allclose(Qs, Uu))
print("U == V :", np.allclose(Uu, Vv))
print("lambda == sigma :", np.allclose(np.sort(람)[::-1], 시그))
print()
print("-> 대칭 양정치라서 세 분해가 하나로 무너진다.")

for 줄 in ("A = LU", "A = QR", "S = Q L Q^T", "A = S L S^-1", "A = U S V^T"):
    됨, 설명, _ = 보A[줄]
    print(f"{줄:<14}{'통과' if 됨 else '실패':<6}{설명}")
print()
print("A 의 고윳값 :", np.sort(np.linalg.eigvals(A).real))
print("0, 1-sqrt2, 1+sqrt2 :", sorted([0.0, 1 - np.sqrt(2), 1 + np.sqrt(2)]))
print()
print("-> 랭크가 모자라 QR 은 안 되는데, 고윳값이 서로 달라 대각화는 된다.")
print("   LU 는 통과한다. 교환 없이 소거가 끝나고 마지막 피벗만 0 이다.")

반례 = (
    ("대칭인데 LU 가 없다", np.array([[0.0, 1.0], [1.0, 0.0]]), "S = Q L Q^T", "A = LU"),
    ("LU 는 되는데 대각화가 안 된다", np.array([[3.0, 1.0], [0.0, 3.0]]),
     "A = LU", "A = S L S^-1"),
    ("대각화는 되는데 QR 이 없다", A, "A = S L S^-1", "A = QR"),
)
for 제목, M, 되는것, 안되는것 in 반례:
    보 = 다섯분해(M, 제목, 출력=False)
    print(f"{제목}")
    print(f"   {되는것:<14}{'O' if 보[되는것][0] else 'X'}   {보[되는것][1]}")
    print(f"   {안되는것:<14}{'O' if 보[안되는것][0] else 'X'}   {보[안되는것][1]}")
    print()
print("-> 어느 조건도 다른 조건을 품지 않는다.")

print("대칭 행렬 500개를 만들어 전부 대각화되는지 본다")
깨진것 = 0
for _ in range(500):
    n = int(rng.integers(2, 6))
    B = rng.normal(size=(n, n))
    Sym = B + B.T
    _, V = np.linalg.eigh(Sym)
    if not np.allclose(V @ V.T, np.eye(n)):
        깨진것 += 1
print("  고유벡터가 정규직교가 아닌 경우 :", 깨진것, "/ 500")
print()
print("역은 거짓인가 :  M = [[2,1],[0,3]]")
M = np.array([[2.0, 1.0], [0.0, 3.0]])
w, V = np.linalg.eig(M)
print("  고윳값", w, " cond(V)", f"{np.linalg.cond(V):.4f}", " -> 대각화 O")
print("  대칭인가 :", np.allclose(M, M.T), " -> 대칭 X")

def 무작위행렬(종류, rng):
    """다섯 족(族)에서 하나씩 뽑는다. (행렬, 대각화 가능한가) 를 돌려준다."""
    n = int(rng.integers(2, 6))
    if 종류 == "일반 정방":
        return rng.normal(size=(n, n)), True
    if 종류 == "대칭":
        B = rng.normal(size=(n, n))
        return B + B.T, True
    if 종류 == "결함":                      # 조르당 블록을 닮음변환으로 감춘다
        J = np.eye(n) * rng.normal()
        for i in range(n - 1):
            J[i, i + 1] = 1.0
        Mv = rng.normal(size=(n, n))
        while abs(np.linalg.det(Mv)) < 1e-2:
            Mv = rng.normal(size=(n, n))
        return Mv @ J @ np.linalg.inv(Mv), False
    if 종류 == "랭크부족":
        r = max(1, n - 1)
        return rng.normal(size=(n, r)) @ rng.normal(size=(r, n)), None
    if 종류 == "직사각":
        m = n + int(rng.integers(1, 4))
        return rng.normal(size=(m, n)), False
    raise ValueError(종류)

종류들 = ["일반 정방", "대칭", "결함", "랭크부족", "직사각"]
줄이름 = ["A = LU", "A = QR", "S = Q L Q^T", "A = S L S^-1", "A = U S V^T"]
N = 400
속은횟수, 조건수들 = 0, []
집계 = {}
rng2 = np.random.default_rng(34)
for 종류 in 종류들:
    합 = np.zeros(5)
    for _ in range(N):
        M, 참 = 무작위행렬(종류, rng2)
        보 = 다섯분해(M, 출력=False)
        표 = [보[k][0] for k in 줄이름]
        if 참 is not None:                       # 대각화는 만들 때의 사실을 쓴다
            표[3] = bool(참) and M.shape[0] == M.shape[1]
        합 += np.array(표, dtype=float)
        if 종류 == "결함":
            _, V = np.linalg.eig(M)
            조건수들.append(np.linalg.cond(V))
            속은횟수 += int(np.linalg.matrix_rank(V) == M.shape[0])
    집계[종류] = 합 / N * 100

print(f"{'':>10}" + "".join(f"{k:>14}" for k in 줄이름))
for 종류 in 종류들:
    print(f"{종류:>10}" + "".join(f"{v:>13.1f}%" for v in 집계[종류]))
전체 = np.mean([집계[k] for k in 종류들], axis=0)
print(f"{'전체':>10}" + "".join(f"{v:>13.1f}%" for v in 전체))
print()
print("서술 파트의 80 / 80 / 20 / 60 / 100 과 맞는가 :", np.round(전체, 1))
print("LU 와 QR 이 똑같이 80% 인데, LU 는 직사각에서 QR 은 랭크부족에서 걸린다.")

print(f"결함 행렬 {N}개에 numpy 를 그대로 믿으면")
print(f"  matrix_rank(V) 가 '꽉 찼다'고 본 비율 : {속은횟수/N*100:.1f}%")
print(f"  고유벡터 행렬 조건수 중앙값 : {np.median(조건수들):.3e}")
print(f"  그중 가장 작은 것          : {min(조건수들):.3e}")
print()
print("-> 조르당 형은 수치적으로 존재하지 않는다는 L28 의 결론 그대로다.")

print(show_matrix(A.T @ A, "A^T A"))
print("눈으로 찾은 고유벡터로 확인")
for v, 배 in (([1.0, 1, 2], 9), ([1.0, -1, 0], 1), ([1.0, 1, -1], 0)):
    v = np.array(v)
    print(f"  A^T A {str(v):<12} = {str(A.T @ A @ v):<16} = {배} x {v}   "
          f"{np.allclose(A.T @ A @ v, 배 * v)}")
print()
print("특이값 :", np.linalg.svd(A, compute_uv=False))
print("sqrt(9), sqrt(1), sqrt(0) :", [3.0, 1.0, 0.0])

v1 = np.array([1.0, 1, 2]) / np.sqrt(6)
v2 = np.array([1.0, -1, 0]) / np.sqrt(2)
v3 = np.array([1.0, 1, -1]) / np.sqrt(3)
시그마 = np.array([3.0, 1.0, 0.0])

u1 = A @ v1 / 시그마[0]              # u 는 반드시 A v / sigma 로 얻는다
u2 = A @ v2 / 시그마[1]
u3 = np.array([-1.0, 1, 1]) / np.sqrt(3)

V손 = np.column_stack([v1, v2, v3])
U손 = np.column_stack([u1, u2, u3])
print("u1 =", u1, "  (2,1,1)/sqrt6 =", np.array([2.0, 1, 1])/np.sqrt(6))
print("u2 =", u2, "  (0,1,-1)/sqrt2 =", np.array([0.0, 1, -1])/np.sqrt(2))
print()
print("V 가 직교인가 :", np.allclose(V손.T @ V손, np.eye(3)))
print("U 가 직교인가 :", np.allclose(U손.T @ U손, np.eye(3)))
print("U S V^T 로 A 가 복원되는가 :",
      np.allclose(U손 @ np.diag(시그마) @ V손.T, A),
      "  오차", f"{np.linalg.norm(U손 @ np.diag(시그마) @ V손.T - A):.2e}")
print()
print("A v3     =", np.round(A @ v3, 12), " (영공간)")
print("A^T u3   =", np.round(A.T @ u3, 12), " (좌영공간)")

w, U따로 = np.linalg.eigh(A @ A.T)
순 = np.argsort(w)[::-1]
U따로 = U따로[:, 순]
print(show_matrix(U따로, "A A^T 의 고유벡터를 따로 뽑은 것"))
print(show_matrix(U손, "A v / sigma 로 얻은 것"))
print()
복원 = U따로 @ np.diag(시그마) @ V손.T
print("따로 뽑은 U 로 복원하면 A 가 되는가 :", np.allclose(복원, A))
print("복원 오차 :", f"{np.linalg.norm(복원 - A):.4f}")
print("어느 열의 부호가 뒤집혔는가 :",
      [j for j in range(3) if not np.allclose(U따로[:, j], U손[:, j])])

A플러스 = np.linalg.pinv(A)
print(show_matrix(A플러스 * 9, "9 A^+   <- 서술의 [[1,5,-4],[1,-4,5],[2,1,1]] 인가"))
print("정수 행렬과 같은가 :",
      np.allclose(A플러스 * 9, [[1, 5, -4], [1, -4, 5], [2, 1, 1]]))
print()
print(f"{'':>6}{'A v = sigma u':>34}{'A^+ u = v / sigma':>34}")
for i, (v, u, s) in enumerate(zip([v1, v2, v3], [u1, u2, u3], 시그마), 1):
    앞 = np.allclose(A @ v, s * u)
    뒤 = np.allclose(A플러스 @ u, v / s) if s > 0 else np.allclose(A플러스 @ u, 0)
    print(f"  i={i} {str(np.round(A @ v, 4)):>28}{str(np.round(A플러스 @ u, 4)):>34}"
          f"   {앞 and 뒤}")
print()
print("i=3 은 양쪽 다 0 이다. 그것이 편도라는 뜻이다.")

왼투영, 오른투영 = A플러스 @ A, A @ A플러스
print(show_matrix(왼투영, "A^+ A   행공간으로의 투영"))
print(show_matrix(오른투영, "A A^+   열공간으로의 투영"))
print("I - v3 v3^T 와 같은가 :", np.allclose(왼투영, np.eye(3) - np.outer(v3, v3)))
print("I - u3 u3^T 와 같은가 :", np.allclose(오른투영, np.eye(3) - np.outer(u3, u3)))
print()
for 이름, P in (("A^+A", 왼투영), ("AA^+", 오른투영)):
    print(f"  {이름} : P^2=P {np.allclose(P@P, P)}, P^T=P {np.allclose(P.T, P)}, "
          f"고윳값 {np.round(np.linalg.eigvalsh(P), 10)}")
print()
print("-> 살렸다, 살렸다, 죽였다.")

traces = [plane(spans=[v1, v2], extent=2.4, color=COLORS["colspace"],
                opacity=0.30, name="행공간 (v1, v2 의 평면)"),
          line(v3, extent=2.8, color=COLORS["nullspace"], name="영공간 (v3 방향)")]
traces += arrow(np.zeros(3), 1.4 * v1, color=COLORS["input"], name="v1")
traces += arrow(np.zeros(3), 1.4 * v2, color=COLORS["second"], name="v2")
traces += arrow(np.zeros(3), 1.4 * v3, color=COLORS["nullspace"], name="v3")
go.Figure(traces, layout=layout3d("정의역 R^3 — 행공간(평면) + 영공간(직선)",
                                  extent=2.6, axis_names=("x1", "x2", "x3")))

traces = [plane(spans=[u1, u2], extent=2.4, color=COLORS["output"],
                opacity=0.25, name="열공간 (u1, u2 의 평면)"),
          line(u3, extent=2.8, color=COLORS["nullspace"], name="좌영공간 (u3 방향)")]
traces += arrow(np.zeros(3), 3.0 * u1, color=COLORS["output"], name="A v1 = 3 u1")
traces += arrow(np.zeros(3), 1.0 * u2, color=COLORS["second"], name="A v2 = 1 u2")
traces += arrow(np.zeros(3), 1.4 * u3, color=COLORS["nullspace"], name="u3 (아무것도 안 온다)")
go.Figure(traces, layout=layout3d("치역 R^3 — 열공간(평면) + 좌영공간(직선)",
                                  extent=2.6, axis_names=("y1", "y2", "y3")))

표본 = rng.normal(size=(600, 3))
표본 /= np.linalg.norm(표본, axis=1, keepdims=True)
프레임, 이름표 = [], []
for t in np.linspace(0.0, 1.0, 11):
    Y = (1 - t) * 표본 + t * (표본 @ A.T)
    프레임.append([
        plane(spans=[u1, u2], extent=2.6, color=COLORS["output"], opacity=0.22,
              grid=6, name="열공간"),          # 평평한 면이라 격자는 성기게
        go.Scatter3d(x=Y[:, 0], y=Y[:, 1], z=Y[:, 2], mode="markers",
                     marker=dict(size=2.4, color=COLORS["input"]),
                     name=f"점 600개 (t = {t:.1f})"),
    ])
    이름표.append(f"{t:.1f}")
배치 = layout3d("구면이 열공간의 평면으로 눌린다", extent=3.2,
                axis_names=("y1", "y2", "y3"))
배치["height"] = 600
slider_figure(프레임, 이름표, 배치, prefix="t = ", initial=0)

Y = 표본 @ A.T
거리 = np.abs(Y @ u3)                       # 열공간까지의 거리
print("t = 1 에서 모든 점이 열공간 위에 있는가")
print("  u3 방향 성분의 최댓값 :", f"{거리.max():.3e}")
print("  (0 이면 평면 위에 정확히 앉았다는 뜻이다)")
print()
print("부피는 얼마나 남았는가 : |det A| =", f"{abs(np.linalg.det(A)):.3e}")
print("  0 이다. 3차원 덩어리가 2차원 평면이 되었으니 부피가 남을 리 없다.")

import math
print(f"{'n':>4}{'크래머 (n+1)n!':>18}{'여인수 n!':>16}{'소거 n^3/3':>14}{'배수':>12}")
for n in (5, 10, 15, 20, 25):
    크 = float(n + 1) * math.factorial(n)
    여 = float(math.factorial(n))
    소 = n ** 3 / 3
    print(f"{n:>4}{크:>18.3e}{여:>16.3e}{소:>14.1f}{크/소:>12.2e}")
print()
print("-> n = 20 에서 5.1e19 대 2667 이다. 아름다움에는 값이 붙는다.")

print("조르당 형 : eps 를 흔들면 고유벡터 행렬이 얼마나 나빠지는가")
print(f"{'eps':>10}{'고윳값':>28}{'cond(S)':>14}{'1/sqrt(eps)':>14}")
for e in (0.0, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4):
    J = np.array([[3.0, 1.0], [e, 3.0]])
    w, V = np.linalg.eig(J)
    이론 = np.inf if e == 0 else 1 / np.sqrt(e)
    갈림 = f"3 +- {abs(w[0] - w[1]) / 2:.3e}" if abs(w[0] - w[1]) > 0 else "3, 3 (구별 못 함)"
    print(f"{e:>10.0e}{갈림:>28}{np.linalg.cond(V):>14.3e}{이론:>14.3e}")
print()
print("-> cond(S) 가 1/sqrt(eps) 를 따라간다. L28 의 그 법칙이다.")

print("정규방정식 : 조건수가 정말 제곱되는가")
print(f"{'cond(A)':>12}{'cond(A^T A)':>16}{'cond(A)^2':>16}{'비':>10}")
for k in (1e2, 1e4, 1e6, 1e8):
    m, n = 40, 6
    Uq, _ = np.linalg.qr(rng.normal(size=(m, m)))
    Vq, _ = np.linalg.qr(rng.normal(size=(n, n)))
    X = Uq[:, :n] @ np.diag(np.logspace(0, -np.log10(k), n)) @ Vq.T
    k1, k2 = np.linalg.cond(X), np.linalg.cond(X.T @ X)
    print(f"{k1:>12.4e}{k2:>16.4e}{k1**2:>16.4e}{k2/k1**2:>10.4f}")
print()
print("-> 마지막 칸이 1 근처다. inv(A.T @ A) @ A.T 를 손으로 쓰지 마라.")

print("같은 일을 어느 기저에서 하느냐")
print(f"{'N':>10}{'직접 DFT N^2':>16}{'FFT N log2 N':>16}{'배수':>12}")
for N in (2**10, 2**12, 2**16, 2**20):
    직 = float(N) ** 2
    빠 = N * np.log2(N)
    print(f"{N:>10}{직:>16.3e}{빠:>16.3e}{직/빠:>12.1f}")
print()
print("대각화 : A^100 을 어떻게 구하는가")
B = np.array([[0.9, 0.2], [0.1, 0.8]])          # 마코브 행렬 (L24)
w, S = np.linalg.eig(B)
직접 = np.linalg.matrix_power(B, 100)
기저 = (S @ np.diag(w ** 100) @ np.linalg.inv(S)).real
print("  행렬곱 99번 :", np.round(직접, 8).tolist())
print("  스칼라 100제곱 2번 :", np.round(기저, 8).tolist())
print("  같은가 :", np.allclose(직접, 기저))
print("  고윳값 :", w, " -> |lambda|<1 인 쪽이 사라지고 1 인 쪽만 남는다")

개수, 복소인것 = 400, 0
시그마모음 = []
rng34 = np.random.default_rng(34)      # 그림 4와 같은 표본을 쓴다
for _ in range(개수):
    B = rng34.normal(size=(4, 4))
    w = np.linalg.eigvals(B)
    시그마모음.extend(np.linalg.svd(B, compute_uv=False))
    복소인것 += int(np.max(np.abs(w.imag)) > 1e-9)
print(f"무작위 4x4 실행렬 {개수}개 중 복소 고윳값을 가진 것 : "
      f"{복소인것} ({복소인것/개수*100:.0f}%)")
print("특이값 중 음수이거나 복소인 것 :",
      int(np.sum(np.array(시그마모음) < 0)), "개")
print()
X = rng34.normal(size=(4, 3))
print("4x3 행렬의 고윳값을 구하려 하면")
try:
    np.linalg.eigvals(X)
except np.linalg.LinAlgError as 오류:
    print("  LinAlgError :", 오류)
print("  특이값은 :", np.linalg.svd(X, compute_uv=False))

print("랭크가 떨어지는 순간 무엇이 무너지는가")
print(f"{'t':>10}{'세 번째 피벗':>16}{'|R33|':>14}{'sigma_3':>14}{'sigma_3/t':>12}")
for t in (1e-1, 1e-3, 1e-5, 1e-7, 1e-9, 0.0):
    M = A.copy()
    M[2, 2] += t
    _, _, Um = lu(M)
    _, Rm = np.linalg.qr(M)
    s3 = np.linalg.svd(M, compute_uv=False)[-1]
    비 = s3 / t if t else float("nan")
    print(f"{t:>10.0e}{abs(Um[2,2]):>16.3e}{abs(Rm[2,2]):>14.3e}{s3:>14.3e}"
          f"{비:>12.4f}")
print()
print("-> 셋 다 0 으로 간다. 다만 LU 와 QR 은 '나눌 수 없다'로 끝나고,")
print("   sigma_3 은 '얼마나 가까운가'를 숫자로 남긴다.")

동물원 = [
    ("L2  치환 [[0,1],[1,0]]", np.array([[0.0, 1.0], [1.0, 0.0]])),
    ("L11 랭크 1 [[1,2],[2,4]]", np.array([[1.0, 2.0], [2.0, 4.0]])),
    ("L16 최소제곱 3x2", np.array([[1.0, 1], [1, 2], [1, 3]])),
    ("L20 [[3,0],[4,5]]", np.array([[3.0, 0.0], [4.0, 5.0]])),
    ("L21 90도 회전", np.array([[0.0, -1.0], [1.0, 0.0]])),
    ("L24 마코브 [[.9,.2],[.1,.8]]", np.array([[0.9, 0.2], [0.1, 0.8]])),
    ("L27 S = [[5,4],[4,5]]", np.array([[5.0, 4.0], [4.0, 5.0]])),
    ("L27 T = tridiag", T),
    ("L28 결함 [[3,1],[0,3]]", np.array([[3.0, 1.0], [0.0, 3.0]])),
    ("L29 C = 3x2", np.array([[1.0, 1], [1, 0], [0, 1]])),
    ("L34 A = 3x3 랭크 2", A),
]
줄이름 = ["A = LU", "A = QR", "S = Q L Q^T", "A = S L S^-1", "A = U S V^T"]
머리 = ["LU", "QR", "QLQ^T", "SLS^-1", "USV^T"]
print(f"{'':>30}" + "".join(f"{h:>9}" for h in 머리) + f"{'통과':>7}")
합계 = np.zeros(5)
for 이름, M in 동물원:
    보 = 다섯분해(M, 출력=False)
    표 = np.array([보[k][0] for k in 줄이름], dtype=float)
    합계 += 표
    print(f"{이름:>30}" + "".join(f"{'O' if v else 'X':>9}" for v in 표)
          + f"{int(표.sum()):>7}")
print(f"{'합계':>30}" + "".join(f"{int(v):>9}" for v in 합계)
      + f"{'/ ' + str(len(동물원)):>7}")
print()
print("-> 마지막 열만 꽉 찼다. 이것이 34강의 결론이다.")
print("   (SLS^-1 열은 복소수까지 허용해 센 것이다. 90도 회전이 그런 경우다.)")

print("마지막으로, 이 책의 행렬 전부에 대해 네 부분공간의 차원을 세어 본다")
print(f"{'':>30}{'m':>4}{'n':>4}{'r':>4}{'행공간':>8}{'영공간':>8}"
      f"{'열공간':>8}{'좌영공간':>10}")
for 이름, M in 동물원:
    m, n = M.shape
    r = np.linalg.matrix_rank(M)
    print(f"{이름:>30}{m:>4}{n:>4}{r:>4}{r:>8}{n-r:>8}{r:>8}{m-r:>10}")
print()
print("네 칸의 차원이 전부 (m, n, r) 하나로 정해진다. L10 의 차원의 기본정리다.")