보강 2. 수치선형대수 입문 — 파이썬 실습

import math

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, show_matrix, slider_figure

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(36)
eps = np.finfo(float).eps
print("numpy", np.__version__, "  eps =", f"{eps:.3e}")

def 삼중(n):
    """tridiag(-1, 2, -1). 고윳값이 닫힌 꼴로 알려져 있다."""
    return (np.diag(2.0 * np.ones(n)) + np.diag(-np.ones(n - 1), 1)
            + np.diag(-np.ones(n - 1), -1))


def 참고윳값(n):
    return np.array([4 * np.sin(k * np.pi / (2 * (n + 1))) ** 2
                     for k in range(1, n + 1)])

for n in (3, 20):
    T = 삼중(n)
    print(f"n={n:>3} : 수치와 공식이 일치 "
          f"{np.allclose(np.sort(np.linalg.eigvalsh(T)), np.sort(참고윳값(n)))}"
          f"   cond = {np.linalg.cond(T):.4f}")
print(show_matrix(삼중(4), "T  (n=4)"))

print(f"{'n':>4}{'소거 n^3/3':>14}{'여인수 n!':>24}{'크래머 (n+1)n!':>28}")
for n in (5, 10, 15, 20, 25):
    print(f"{n:>4}{n**3//3:>14,}{math.factorial(n):>24,}"
          f"{(n+1)*math.factorial(n):>28,}")
print()
초 = 26 * math.factorial(25) / 1e9
print(f"n=25 에서 크래머를 1 GFLOP/s 로 돌리면 {초/3.15e7:.2e} 년")
print(f"우주의 나이는 약 1.4e10 년이다.")

def 후방오차(A, b, x햇):
    """잔차만으로 후방오차를 계산한다. 참해를 몰라도 된다."""
    r = b - A @ x햇
    ΔA = np.outer(r, x햇) / (x햇 @ x햇)
    맞는가 = np.allclose((A + ΔA) @ x햇, b)
    return np.linalg.norm(ΔA, 2) / np.linalg.norm(A, 2), 맞는가

A = 삼중(6)
b = rng.normal(size=6)
x햇 = np.linalg.solve(A, b)
값, 맞는가 = 후방오차(A, b, x햇)
print("(A + dA) x_hat = b 인가 :", 맞는가)
print(f"후방오차 = {값:.4e}   eps = {eps:.4e}   비 = {값/eps:.2f}")
print()
r = b - A @ x햇
print("잔차만으로 계산한 값과 같은가 :",
      np.isclose(값, np.linalg.norm(r)/(np.linalg.norm(A,2)*np.linalg.norm(x햇))))
print()
# 일부러 나쁜 답을 넣어 본다
나쁜x = x햇 * 1.01
값2, _ = 후방오차(A, b, 나쁜x)
print(f"1% 틀린 답의 후방오차 : {값2:.4e}   <- 훨씬 크다")

def 정수문제(n, 폭=5, seed=0):
    """A 와 x 를 정수로 잡아 b = Ax 가 반올림 없이 정확하게 한다."""
    r = np.random.default_rng(seed)
    while True:
        A = r.integers(-폭, 폭 + 1, size=(n, n)).astype(float)
        if abs(np.linalg.det(A)) > 0.5:
            break
    x = r.integers(-9, 10, size=n).astype(float)
    return A, A @ x, x

print("정수 문제 400개에서 부등식이 깨지는가")
깨짐, 최대비, 셈 = 0, 0.0, 0
for s in range(400):
    n = int(np.random.default_rng(s).integers(3, 8))
    A0, b0, x참 = 정수문제(n, seed=s)
    κ = np.linalg.cond(A0)
    if κ > 1e12:
        continue
    셈 += 1
    x = np.linalg.solve(A0, b0)
    후, _ = 후방오차(A0, b0, x)
    전 = np.linalg.norm(x - x참) / np.linalg.norm(x)
    한계 = κ * 후
    if 한계 > 0:
        깨짐 += 전 > 한계 * (1 + 1e-9)
        최대비 = max(최대비, 전 / 한계)
    elif 전 > 0:
        깨짐 += 1
        최대비 = np.inf
print(f"  {셈}개 중 깨진 횟수 : {깨짐}")
print(f"  전방오차 / 한계 의 최댓값 : {최대비:.4f}")
print()
print("-> 깨진다. 그런데 유도에는 근사가 없었다. 무엇이 잘못되었는가.")

from fractions import Fraction


def 정확한잔차(A, b, x햇):
    """A, b 가 정수일 때 r = b - A x_hat 을 유리수로 정확히 계산한다."""
    F = [Fraction(v) for v in x햇]
    out = []
    for i in range(len(b)):
        s = Fraction(int(round(b[i])))
        for j in range(len(F)):
            s -= Fraction(int(round(A[i, j]))) * F[j]
        out.append(float(s))
    return np.array(out)

print(f"{'잔차를 어떻게 재는가':>22}{'깨진 횟수':>12}{'전방/한계 최댓값':>18}")
for 이름, 잔차재기 in (("보통 (float)", lambda A, b, x: b - A @ x),
                      ("정확 (Fraction)", 정확한잔차)):
    깨짐, 최대비, 셈 = 0, 0.0, 0
    for s in range(400):
        n = int(np.random.default_rng(s).integers(3, 8))
        A0, b0, x참 = 정수문제(n, seed=s)
        κ = np.linalg.cond(A0)
        if κ > 1e12:
            continue
        셈 += 1
        x = np.linalg.solve(A0, b0)
        r = 잔차재기(A0, b0, x)
        후 = np.linalg.norm(r) / (np.linalg.norm(A0, 2) * np.linalg.norm(x))
        전 = np.linalg.norm(x - x참) / np.linalg.norm(x)
        한계 = κ * 후
        if 한계 > 0:
            깨짐 += 전 > 한계 * (1 + 1e-9)
            최대비 = max(최대비, 전 / 한계)
        elif 전 > 0:
            깨짐 += 1
            최대비 = np.inf
    print(f"{이름:>22}{깨짐:>12}{최대비:>18.4f}")
print()
print("-> 정확히 재면 400개 중 0번 깨진다. 부등식은 옳았다.")
print("   그리고 최댓값이 0.99 다. 부등식이 옳을 뿐 아니라 **꽉 차 있다.**")
print()
print("   이것이 반복 개선(iterative refinement)이 잔차를 더 높은 정밀도로")
print("   계산하는 이유다. 잔차를 같은 정밀도로 재면 개선할 것이 안 보인다.")

print("한 문제를 자세히 보자")
A0, b0, x참 = 정수문제(5, seed=3)
x = np.linalg.solve(A0, b0)
r보 = b0 - A0 @ x
r정 = 정확한잔차(A0, b0, x)
print("  보통 잔차 :", " ".join(f"{v:+.4e}" for v in r보))
print("  정확 잔차 :", " ".join(f"{v:+.4e}" for v in r정))
다른곳 = int(np.sum(r보 != r정))
비 = np.linalg.norm(r보) / np.linalg.norm(r정)
print(f"\n  다섯 성분 중 {다른곳}개가 다르다.")
print(f"  노름 : 보통 {np.linalg.norm(r보):.6e},  정확 {np.linalg.norm(r정):.6e}")
print(f"  보통 쪽이 {비:.4f} 배, 곧 {(비-1)*100:+.1f}% 만큼 어긋났다.")
print()
print("  잔차 자체가 1e-15 수준이라 성분 하나의 마지막 자리가 통째로 흔들린다.")
print("  그 흔들림이 한계를 몇 % 움직이고, 그것이 부등식을 깨 보이게 한다.")

print("힐베르트 행렬 : 조건수가 아주 나쁜 대표 선수")
print(f"{'n':>4}{'조건수':>12}{'후방오차':>12}{'전방오차':>12}{'조건수x후방':>14}")
for n in (4, 6, 8, 10, 12):
    H = np.array([[1.0/(i+j+1) for j in range(n)] for i in range(n)])
    x참 = np.ones(n)
    bb = H @ x참
    x = np.linalg.solve(H, bb)
    후, _ = 후방오차(H, bb, x)
    전 = np.linalg.norm(x - x참) / np.linalg.norm(x참)
    κ = np.linalg.cond(H)
    print(f"{n:>4}{κ:>12.2e}{후:>12.2e}{전:>12.2e}{κ*후:>14.2e}")
print()
print("-> 후방오차는 늘 eps 언저리다. 알고리즘은 제 몫을 했다.")
print("   그런데 전방오차가 n=12 에서 1 을 넘는다. 문제가 나쁜 것이다.")

def 쪼개기(A):
    """A = D + L + U 로 나눈다."""
    D = np.diag(np.diag(A))
    return D, np.tril(A, -1), np.triu(A, 1)


def 반복(A, b, M, N, 횟수):
    x = np.zeros(len(b))
    자취 = [1.0]
    for _ in range(횟수):
        x = np.linalg.solve(M, N @ x + b)
        자취.append(np.linalg.norm(A @ x - b) / np.linalg.norm(b))
    return x, np.array(자취)


def 켤레기울기(A, b, 횟수):
    """대칭 양의 정부호에서만 쓴다."""
    x = np.zeros(len(b))
    r = b - A @ x
    p = r.copy()
    자취 = [1.0]
    for _ in range(횟수):
        Ap = A @ p
        a = (r @ r) / (p @ Ap)
        x = x + a * p
        r새 = r - a * Ap
        beta = (r새 @ r새) / (r @ r)
        p = r새 + beta * p
        r = r새
        자취.append(np.linalg.norm(A @ x - b) / np.linalg.norm(b))
    return x, np.array(자취)

n = 20
T = 삼중(n)
bb = rng.normal(size=n)
D, L_, U_ = 쪼개기(T)

print("스펙트럼 반지름이 1 보다 작은가")
for 이름, M, N in (("야코비", D, -(L_ + U_)),
                   ("가우스-자이델", D + L_, -U_)):
    ρ = max(abs(np.linalg.eigvals(np.linalg.solve(M, N))))
    print(f"  {이름:>12} : rho = {ρ:.6f}   수렴 {ρ < 1}")
야ρ = max(abs(np.linalg.eigvals(np.linalg.solve(D, -(L_+U_)))))
가ρ = max(abs(np.linalg.eigvals(np.linalg.solve(D+L_, -U_))))
print()
print(f"야코비의 rho 가 cos(pi/(n+1)) 인가 : {야ρ:.9f} vs "
      f"{np.cos(np.pi/(n+1)):.9f}", np.isclose(야ρ, np.cos(np.pi/(n+1))))
print(f"가우스-자이델이 그 제곱인가       : {가ρ:.9f} vs {야ρ**2:.9f}",
      np.isclose(가ρ, 야ρ**2))

_, 야 = 반복(T, bb, D, -(L_ + U_), 250)
_, 가 = 반복(T, bb, D + L_, -U_, 250)
xCG, 씨 = 켤레기울기(T, bb, 30)
x참 = np.linalg.solve(T, bb)

print(f"{'':>14}{'50회 뒤':>12}{'250회 뒤':>12}{'1e-10 도달':>12}")
for 이름, 자취 in (("야코비", 야), ("가우스-자이델", 가)):
    도달 = int(np.argmax(자취 < 1e-10)) if (자취 < 1e-10).any() else -1
    print(f"{이름:>14}{자취[50]:>12.2e}{자취[250]:>12.2e}"
          f"{(도달 if 도달 > 0 else '못 감'):>12}")
도달 = int(np.argmax(씨 < 1e-10))
print(f"{'켤레기울기':>14}{'':>12}{씨[-1]:>12.2e}{도달:>12}")
print()
print(f"CG 가 n={n} 번째에서 : 잔차 {씨[n]:.3e},  "
      f"해 오차 {np.linalg.norm(xCG-x참)/np.linalg.norm(x참):.3e}")
print("-> 정확히 n 걸음에서 끝난다. 이론이 말한 그대로다.")

print("감소 인자를 n 에 따라")
print(f"{'n':>5}{'cond':>12}{'야코비':>12}{'최급강하':>12}{'CG':>12}{'같은가':>8}")
for nn in (5, 10, 20, 50, 100):
    Tn = 삼중(nn)
    κ = np.linalg.cond(Tn)
    야 = np.cos(np.pi/(nn+1))
    최 = (κ-1)/(κ+1)
    씨 = (np.sqrt(κ)-1)/(np.sqrt(κ)+1)
    print(f"{nn:>5}{κ:>12.1f}{야:>12.6f}{최:>12.6f}{씨:>12.6f}"
          f"{str(np.isclose(야, 최)):>8}")
print()
print("-> 야코비와 최급강하가 이 행렬에서는 정확히 같은 수다. 항등식이다.")
print("   lam_max + lam_min = 4 이고 lam_max - lam_min = 4cos(pi/(n+1)) 이라 그렇다.")
print()
κ = np.linalg.cond(삼중(20))
print(f"n=20 에서 자릿수 하나(10배)를 더 얻는 데 드는 걸음 수")
for 이름, 인자 in (("야코비", np.cos(np.pi/21)), ("CG", (np.sqrt(κ)-1)/(np.sqrt(κ)+1))):
    print(f"  {이름:>8} : {np.log(0.1)/np.log(인자):.1f} 걸음")

def 거듭제곱법(A, 횟수=40, 씨=0):
    """A 를 자꾸 곱하고 길이를 1 로 맞춘다. L22 의 그 식이다."""
    r = np.random.default_rng(씨)
    v = r.normal(size=len(A))
    v = v / np.linalg.norm(v)
    자취 = []
    for _ in range(횟수):
        v = A @ v
        v = v / np.linalg.norm(v)
        자취.append((v.copy(), float(v @ A @ v)))
    return 자취

A4 = np.array([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 2.0]])
값, 벡 = np.linalg.eigh(A4)
순 = np.argsort(값)[::-1]
값, 벡 = 값[순], 벡[:, 순]
q1 = 벡[:, 0]
비 = abs(값[1] / 값[0])
print("고윳값 :", 값, "   |lam2/lam1| =", f"{비:.6f}")
print()
자취 = 거듭제곱법(A4, 40)
벡오차 = np.array([min(np.linalg.norm(v-q1), np.linalg.norm(v+q1)) for v, _ in 자취])
값오차 = np.array([abs(l - 값[0]) for _, l in 자취])
기 = np.polyfit(np.arange(20, 38), np.log(벡오차[20:38]), 1)[0]
기2 = np.polyfit(np.arange(15, 30), np.log(값오차[15:30]), 1)[0]
print(f"고유벡터 수렴 인자 실측 {np.exp(기):.6f}   이론 |lam2/lam1| = {비:.6f}")
print(f"고윳값   수렴 인자 실측 {np.exp(기2):.6f}   이론 (lam2/lam1)^2 = {비**2:.6f}")
print()
print(f"{'k':>4}{'고유벡터 오차':>16}{'레일리 몫 오차':>18}")
for k in (5, 10, 20, 30):
    print(f"{k:>4}{벡오차[k]:>16.3e}{값오차[k]:>18.3e}")
print()
print("-> 고유벡터 오차가 delta 일 때 고윳값 오차가 delta^2 이다. 두 배 빠르다.")
print(f"   30번째에서 레일리 몫 = {자취[29][1]:.12f},  참값 = {값[0]:.12f}")

def QR알고리즘(A, 시프트=False, 횟수=60):
    """A_(k+1) = R_k Q_k. 시프트를 넣으면 훨씬 빠르다."""
    Ak = A.astype(float).copy()
    n = len(A)
    자취 = [abs(Ak[-1, -2])]
    사진 = [Ak.copy()]
    for _ in range(횟수):
        mu = Ak[-1, -1] if 시프트 else 0.0
        Q, R = np.linalg.qr(Ak - mu * np.eye(n))
        Ak = R @ Q + mu * np.eye(n)
        자취.append(abs(Ak[-1, -2]))
        사진.append(Ak.copy())
    return Ak, np.array(자취), 사진

B = np.array([[4.0, 1.0, 2.0], [1.0, 3.0, 1.0], [2.0, 1.0, 5.0]])
참 = np.sort(np.linalg.eigvalsh(B))[::-1]
끝, 민자취, 사진 = QR알고리즘(B, False, 60)
print("60회 뒤 A_k :")
print(show_matrix(끝, ""))
print("대각선 :", np.sort(np.diag(끝))[::-1])
print("참 고윳값 :", 참)
print("일치 :", np.allclose(np.sort(np.diag(끝)), np.sort(참)))
print()
print("매 단계가 유사변환인가 (고윳값이 안 움직이는가)")
for k in (0, 1, 5, 20, 60):
    print(f"  k={k:>3} : 고윳값 {np.sort(np.linalg.eigvals(사진[k]).real)}")

비들 = [abs(참[i+1]/참[i]) for i in range(len(참)-1)]
print("이웃 고윳값의 비 :", np.round(비들, 6), "  최댓값 :", f"{max(비들):.6f}")
기 = np.polyfit(np.arange(20, 50), np.log(민자취[20:50]), 1)[0]
print(f"실측 감소 인자 {np.exp(기):.6f}")
print()
print("-> |lam2/lam1| = 0.375 가 아니라 max_i |lam_(i+1)/lam_i| 가 정한다.")
print("   가까운 고윳값이 하나라도 있으면 전체가 느려진다.")

_, 시자취, _ = QR알고리즘(B, True, 12)
print("시프트를 넣으면")
print(f"{'k':>4}{'시프트 없음':>16}{'시프트 있음':>16}")
for k in range(8):
    없 = f"{민자취[k]:.3e}"
    있 = f"{시자취[k]:.3e}" if k < len(시자취) else ""
    print(f"{k:>4}{없:>16}{있:>16}")
print()
없도달 = int(np.argmax(민자취 < 1e-12)) if (민자취 < 1e-12).any() else -1
있도달 = int(np.argmax(시자취 < 1e-12))
print(f"1e-12 도달 : 시프트 없음 {없도달 if 없도달>0 else '60회 내 못함'}, "
      f"시프트 있음 {있도달}회")
print("-> 자릿수가 매번 세 배로 늘어난다. 세제곱 수렴이다.")

프레임, 이름표 = [], []
for k in (0, 1, 2, 3, 5, 8, 12, 20, 30, 45, 60):
    프레임.append([go.Heatmap(z=np.abs(사진[k])[::-1], colorscale="Oranges",
                             zmin=0, zmax=7, showscale=False,
                             text=np.round(사진[k], 3)[::-1],
                             texttemplate="%{text}")])
    이름표.append(str(k))
배치 = dict(title=dict(text="A_k 가 대각으로 굳어 간다"),
           xaxis=dict(visible=False, scaleanchor="y"),
           yaxis=dict(visible=False),
           height=460, margin=dict(l=60, r=60, t=60, b=40))
slider_figure(프레임, 이름표, 배치, prefix="k = ", initial=0)

근 = np.arange(1.0, 21.0)
계수 = np.poly(근)
print(f"근이 1..20 인 다항식.  x^19 의 계수 = {계수[1]:.0f}")
print()
print(f"{'상대섭동':>12}{'복소수 근':>12}{'최대 이동':>14}")
원근 = np.sort_complex(np.roots(계수))
for 크기 in (1e-12, 1e-10, 1e-8, 1e-6):
    c = 계수.copy()
    c[1] *= (1 + 크기)
    r = np.sort_complex(np.roots(c))
    print(f"{크기:>12.0e}{int(np.sum(np.abs(r.imag) > 1e-8)):>12}"
          f"{np.max(np.abs(r - 원근)):>14.4f}")
print()
print("같은 크기로 행렬 쪽을 흔들면")
D20 = np.diag(근)
print(f"{'상대섭동':>12}{'복소수 고윳값':>16}{'최대 이동':>14}")
for 크기 in (1e-12, 1e-10, 1e-8, 1e-6):
    E = D20 + 크기 * 20.0 * rng.normal(size=(20, 20))
    w = np.linalg.eigvals(E)
    print(f"{크기:>12.0e}{int(np.sum(np.abs(w.imag) > 1e-8)):>16}"
          f"{np.max(np.abs(np.sort(w.real) - 근)):>14.4e}")
print()
print("-> 같은 크기의 섭동인데 한쪽은 근이 복소평면으로 날아가고")
print("   다른 쪽은 제자리를 지킨다. 행렬을 다항식으로 바꾸면 문제가 나빠진다.")

print("그래서 numpy 는 반대로 간다 — 동반행렬을 만들어 고윳값을 구한다")
p = np.array([1.0, -6.0, 11.0, -6.0])       # (x-1)(x-2)(x-3)
동반 = np.zeros((3, 3))
동반[0] = -p[1:] / p[0]
동반[1, 0] = 동반[2, 1] = 1.0
print(show_matrix(동반, "동반행렬"))
print("동반행렬의 고윳값 :", np.sort(np.linalg.eigvals(동반).real))
print("np.roots 의 답    :", np.sort(np.roots(p)))
print("같은가 :", np.allclose(np.sort(np.linalg.eigvals(동반).real),
                             np.sort(np.roots(p))))

def 고윳값조건수(A):
    """1 / |y^T x|. 좌우 고유벡터의 내적이 정한다."""
    w, V = np.linalg.eig(A)
    wL, W = np.linalg.eig(A.T)
    순, 순L = np.argsort(-w.real), np.argsort(-wL.real)
    w, V, W = w[순], V[:, 순], W[:, 순L]
    out = []
    for i in range(len(w)):
        x = V[:, i] / np.linalg.norm(V[:, i])
        y = W[:, i] / np.linalg.norm(W[:, i])
        s = abs(np.vdot(y, x))
        out.append((w[i].real, np.inf if s < 1e-14 else 1.0 / s))
    return out

print(f"{'':>28}{'lam_1':>12}{'조건수':>14}")
경우 = (("대칭 [[4,1],[1,3]]", np.array([[4.,1.],[1.,3.]])),
        ("비대칭 [[4,1],[0,3]]", np.array([[4.,1.],[0.,3.]])),
        ("거의 결함 [[3,1],[1e-6,3]]", np.array([[3.,1.],[1e-6,3.]])),
        ("거의 결함 [[3,1],[1e-12,3]]", np.array([[3.,1.],[1e-12,3.]])))
for 이름, M in 경우:
    l, κ = 고윳값조건수(M)[0]
    print(f"{이름:>28}{l:>12.6f}{κ:>14.4e}")
print()
print("-> 대칭이면 정확히 1 이다. 좌우 고유벡터가 같기 때문이다.")
print("   스펙트럼 정리(L25)가 주는 또 하나의 선물이다.")

print("결함 행렬은 delta^(1/m) 로 튄다  (m = 조르당 블록 크기)")
J = np.array([[3.0, 1.0], [0.0, 3.0]])
print(f"{'섭동 delta':>14}{'고윳값 간격':>16}{'sqrt(delta)':>14}{'비':>8}")
for d in (1e-12, 1e-10, 1e-8, 1e-6, 1e-4):
    Jd = J.copy(); Jd[1, 0] = d
    w = np.linalg.eigvals(Jd)
    간격 = float(abs(w[0] - w[1]))
    print(f"{d:>14.0e}{간격:>16.6e}{2*np.sqrt(d):>14.6e}"
          f"{간격/(2*np.sqrt(d)):>8.4f}")
print()
print("-> 간격이 정확히 2 sqrt(delta) 다. 로그-로그 기울기가 0.5 라는 뜻이고,")
print("   delta = eps = 2.2e-16 이면 고윳값이 1.5e-8 만큼 튄다.")
print("   이것이 L28 에서 관찰만 하고 이름을 못 붙였던 그 법칙이다.")

from scipy.linalg import schur
print("대안 : 슈어 분해 A = Q T Q^H")
Tq, Q = schur(np.array([[3.0, 1.0], [1e-13, 3.0]]))
print(show_matrix(Q, "Q  (직교)"))
print("Q 가 직교인가 :", np.allclose(Q.T @ Q, np.eye(2)),
      "  cond(Q) =", f"{np.linalg.cond(Q):.6f}")
print(show_matrix(Tq, "T  (위삼각, 대각선이 고윳값)"))
print()
print("-> Q 의 조건수가 정확히 1 이다. 고유벡터 행렬처럼 터지지 않는다.")
print("   블록 구조는 못 얻지만, 얻는 것은 믿을 수 있다.")