보강 4-3 실습 — 함정을 하나씩 재현한다

import itertools

import numpy as np

np.set_printoptions(precision=4, suppress=True, linewidth=115)


def 같다(번호, 얻은, 적힌, 말="", 허용=1e-9):
    얻은, 적힌 = np.asarray(얻은, dtype=float), np.asarray(적힌, dtype=float)
    차 = np.abs(얻은 - 적힌).max()
    print(f"[{번호}] {'맞다  ' if 차 < 허용 else '틀리다'} "
          f"차이 {차:.2e}   {말}")

S = np.ones((3, 3)) - 3 * np.eye(3)
x = np.array([2.0, 2, 2])
dm = np.array([0.5, 0.25, 0.25])
D = S - np.outer(dm, x)

같다("47", S, S.T, "이론이 말하는 표는 대칭")
같다("47", np.sort(np.linalg.eigvalsh(S)), [-3, -3, 0], "고윳값 0, -3, -3")
같다("47", S @ np.ones(3), np.zeros(3), "S p = 0 (값이 같이 오르면 안 바뀐다)")
print(f"     관측표 D =\n{D}")
print(f"     D 는 대칭인가 : {np.allclose(D, D.T)}   D12={D[0,1]}, D21={D[1,0]}")
print(f"     D 의 행 합 {D.sum(axis=1)}  -> D p 도 0 이 아니다")

# 어긋남이 잡음이라면 규칙적일 리 없다
print()
print(f"{'(i,j)':>8} {'D_ij - D_ji':>12} {'(dm_j - dm_i) x':>18}")
for i, j in ((0, 1), (0, 2), (1, 2)):
    왼 = D[i, j] - D[j, i]
    오 = (dm[j] - dm[i]) * x[i]
    print(f"{f'({i+1},{j+1})':>8} {왼:>12.2f} {오:>18.2f}")
print("  -> 잡음이 아니라 계산 가능한 양이다")

# 표본을 늘려 보면 정말 안 가까워지는지
난수 = np.random.default_rng(40)
print(f"{'표본 수':>10} {'|D - D^T| 최대':>16}")
for n in (50, 500, 5000, 50000):
    잡음 = 난수.normal(scale=1 / np.sqrt(n), size=(3, 3))
    관측 = D + 잡음
    print(f"{n:>10} {np.abs(관측 - 관측.T).max():>16.4f}")
print(f"{'무한':>10} {np.abs(D - D.T).max():>16.4f}   <- 0 으로 안 간다")

def 지렛대(부하):
    A = np.c_[np.ones(len(부하)), np.asarray(부하, dtype=float)]
    return np.diag(A @ np.linalg.inv(A.T @ A) @ A.T)


h = 지렛대([1, 2, 3, 6])
같다("48", h, np.array([30, 18, 14, 50]) / 56, "지렛대값")
같다("48", h.sum(), 2.0, "합 = 랭크 = 2")
print(f"     넷째가 총량의 {h[3]/h.sum():.1%}")

print()
print(f"{'부하':>18} {'h_ii':>28} {'합':>6} {'최대':>7}")
for 부하 in ([1, 2, 3, 6], [1, 2, 3, 4], [1, 2, 3, 4, 6], [1, 2, 3, 4, 5, 6]):
    hh = 지렛대(부하)
    print(f"{str(부하):>18} {str(np.round(hh, 3)):>28} {hh.sum():>6.1f} {hh.max():>7.3f}")
print("  합은 늘 2 다. 나누는 방식만 고를 수 있다.")

부하 = np.array([1.0, 2, 3, 6])
A = np.c_[np.ones(4), 부하]
참 = 10 + 2 * 부하                              # 참 관계
h = 지렛대(부하)

print(f"{'어느 점을 5만큼 흔드나':>24} {'그 점의 잔차':>12} {'h_ii':>8}")
for k in range(4):
    y = 참.copy()
    y[k] += 5
    be = np.linalg.solve(A.T @ A, A.T @ y)
    e = y - A @ be
    print(f"{f'{k+1}번 (부하 {부하[k]:.0f})':>24} {e[k]:>12.3f} {h[k]:>8.3f}")
print("  h_ii 가 클수록 잔차가 작다. 5만큼 틀렸는데 0.5 밖에 안 남는다.")
같다("48", 1 - h, [abs(5 * (1 - h[k])) / 5 for k in range(4)],
    "잔차 = (1 - h_ii) x 흔든 양")

mS = 3 * np.eye(3) - np.ones((3, 3))
선행 = [np.linalg.det(mS[:k, :k]) for k in (1, 2, 3)]
같다("49", 선행, [2, 3, 0], "선행 주소행렬식")
같다("49", np.sort(np.linalg.eigvalsh(mS)), [0, 3, 3], "고윳값")

모두 = []
for k in (1, 2, 3):
    for idx in itertools.combinations(range(3), k):
        모두.append(round(float(np.linalg.det(mS[np.ix_(idx, idx)])), 9))
print(f"     모든 주소행렬식 {sorted(모두)}  (2^3 - 1 = {2**3-1} 개)")
print(f"     전부 0 이상인가 : {all(v >= -1e-9 for v in 모두)}")

print()
print("반례 : 선행은 통과하는데 준양정부호가 아닌 표")
K = np.array([[0.0, 0], [0, -1]])
print(f"     K = {K.tolist()}")
print(f"     선행 주소행렬식 {[np.linalg.det(K[:k, :k]) for k in (1, 2)]}  -> 음수 없음")
print(f"     고윳값 {np.linalg.eigvalsh(K)}  -> 준양정부호가 아니다")
모두K = [round(float(np.linalg.det(K[np.ix_(i, i)])), 9)
        for k in (1, 2) for i in itertools.combinations(range(2), k)]
print(f"     모든 주소행렬식 {sorted(모두K)}  -> -1 이 있어 걸린다")

print(f"{'n':>4} {'모든 주소행렬식 개수':>20} {'고윳값 개수':>12}")
for n in (3, 5, 10, 20):
    print(f"{n:>4} {2**n - 1:>20} {n:>12}")

def 수율(v, 교차):
    v = np.asarray(v, dtype=float)
    return 80 + 6 * v[1] - 2 * v[0] ** 2 - 2 * v[1] ** 2 + 교차 * v[0] * v[1]


for 교차 in (2.0, 6.0):
    B = np.array([[-2.0, 교차 / 2], [교차 / 2, -2]])
    xs = np.linalg.solve(2 * B, -np.array([0.0, 6]))
    w = np.linalg.eigvalsh(B)
    성격 = "봉우리" if (w < 0).all() else "안장"
    print(f"교차항 {교차:.0f} : B = {B.tolist()},  고윳값 {w},  "
          f"정류점 {np.round(xs, 3)},  y = {수율(xs, 교차):.2f}   -> {성격}")
print("  정류점도 옮겨 가고 성격도 바뀐다. 대각은 둘 다 (-2,-2) 로 그대로다.")

B = np.array([[-2.0, 1], [1, -2]])
xs = np.array([1.0, 2])
같다("50", xs, np.linalg.solve(2 * B, -np.array([0.0, 6])), "정류점 (1,2)")
같다("50", 수율(xs, 2.0), 86.0, "y* = 86")

# 능선 : 방향마다 떨어지는 속도가 lambda 다
w, V = np.linalg.eigh(B)
print(f"{'방향':>12} {'v^T B v':>10} {'lambda |v|^2':>14}")
for k in range(2):
    v = V[:, k]
    print(f"{str(np.round(v, 3)):>12} {v @ B @ v:>10.4f} {w[k] * (v @ v):>14.4f}")
print(f"  느린 방향 {np.round(V[:, int(np.argmax(w))], 3)} 이 능선이다")

print()
print("실험 범위 [-1,1]^2 안에서 가장 좋았던 점과 모형이 가리키는 봉우리 :")
격 = np.linspace(-1, 1, 41)
최고 = max(((수율((a, b), 2.0), (a, b)) for a in 격 for b in 격))
print(f"  범위 안 최고 : {최고[1]} 에서 {최고[0]:.2f}")
print(f"  모형의 봉우리 : (1, 2) 에서 86.00   <- 압력 2 는 시험한 적이 없다")

표 = {"주 표": [[.5, .1, .3], [.2, .5, .1], [.1, .2, .4]],
      "잘못 입력한 표": [[.5, .1, .3], [.2, .5, .1], [.6, .6, .4]],
      "열 합이 1 을 넘는 표": [[.5, .1, .3], [.2, .5, .1], [.4, .5, .4]]}
print(f"{'표':>22} {'열 합':>18} {'선행 주소행렬식':>26} {'rho':>8} {'성립':>6}")
for 이름, M in 표.items():
    A = np.array(M)
    IA = np.eye(3) - A
    선행 = [np.linalg.det(IA[:k, :k]) for k in (1, 2, 3)]
    rho = max(abs(np.linalg.eigvals(A)))
    성립 = all(v > 0 for v in 선행)
    print(f"{이름:>22} {str(np.round(A.sum(axis=0), 2)):>18} "
          f"{str(np.round(선행, 4)):>26} {rho:>8.4f} {str(성립):>6}")

print()
print("열 합이 rho 를 가두는 범위 (성분이 음이 아닐 때) :")
for 이름, M in 표.items():
    A = np.array(M)
    열 = A.sum(axis=0)
    rho = max(abs(np.linalg.eigvals(A)))
    print(f"  {이름:>22} : {열.min():.2f} <= {rho:.4f} <= {열.max():.2f}  "
          f"{'참' if 열.min() - 1e-9 <= rho <= 열.max() + 1e-9 else '거짓'}")
print("  -> 최소가 1 을 넘어야 rho > 1 이 확정된다. 일부만 넘는 것으로는 모른다.")

e0 = np.array([2.0, 0])
print("이산 : e_(t+1) = M e_t")
for 이름, M in [("갑", [[.4, .3], [.3, .4]]), ("을", [[.4, .8], [.8, .4]])]:
    M = np.array(M)
    w = np.linalg.eigvalsh(M)
    e5 = np.linalg.matrix_power(M, 5) @ e0
    print(f"  {이름} : 대각 {np.diag(M)}  대각합 {np.trace(M):.1f}  det {np.linalg.det(M):+.2f}  "
          f"고윳값 {w}  5년 뒤 {np.round(e5, 4)}")
같다("52", np.linalg.matrix_power(np.array([[.4, .3], [.3, .4]]), 5) @ e0,
    [0.7**5, 0.7**5], "갑은 0.7^5 로 줄어든다", 허용=1e-4)
print("  대각합이 둘 다 0.8 이다. 가른 것은 행렬식이다.")

print("연속 : dp/dt = J p,  세 도시 모두 대각 (-2,-2), 대각합 -4")
for 이름, J in [("(가)", [[-2., 1], [1, -2]]), ("(나)", [[-2., -3], [3, -2]]),
               ("(다)", [[-2., 9], [1, -2]])]:
    J = np.array(J)
    w = np.linalg.eigvals(J)
    판별 = np.trace(J) ** 2 - 4 * np.linalg.det(J)
    거동 = ("발산" if max(np.real(w)) > 0
           else ("나선 수렴" if abs(np.imag(w)).max() > 1e-9 else "단조 수렴"))
    print(f"  {이름} det {np.linalg.det(J):+6.1f}  판별식 {판별:+7.1f}  "
          f"고윳값 {np.round(w, 3)}  -> {거동}")
print("  대각합이 같으니 남은 자유도는 행렬식 하나뿐이다.")

from scipy.linalg import expm
print()
print(f"{'t':>5} {'(가) 거리':>11} {'(나) 거리':>11} {'(다) 거리':>13}")
p0 = np.array([1.0, 0])
for t in (0.5, 1, 2, 4):
    줄 = [np.linalg.norm(expm(np.array(J) * t) @ p0)
          for J in ([[-2., 1], [1, -2]], [[-2., -3], [3, -2]], [[-2., 9], [1, -2]])]
    print(f"{t:>5.1f} {줄[0]:>11.5f} {줄[1]:>11.5f} {줄[2]:>13.4f}")

z1 = np.array([3.0, 4, 0, -4, -3])
z2 = np.array([4.0, 3, 0, -3, -4])
Z = np.c_[z1, z2]
G = Z.T @ Z
같다("53", G, [[50, 48], [48, 50]], "Z^T Z")
같다("53", np.sort(np.linalg.eigvalsh(G)), [2, 98], "고윳값 2, 98")
s = np.linalg.svd(Z)[1]
같다("53", s, [7 * np.sqrt(2), np.sqrt(2)], "특이값")
같다("53", np.linalg.cond(Z), 7.0, "kappa(Z) = 7")
같다("53", np.linalg.cond(G), 49.0, "kappa(Z^T Z) = 49")
같다("53", s ** 2, np.linalg.svd(G)[1], "sigma(Z^T Z) = sigma(Z)^2")

print()
print(f"{'우변':>16} {'합':>6} {'차':>6} {'계수':>16}")
for rhs in ([98., 98], [100., 96], [99., 97], [102., 94]):
    be = np.linalg.solve(G, rhs)
    print(f"{str(rhs):>16} {rhs[0]+rhs[1]:>6.0f} {rhs[0]-rhs[1]:>6.0f} "
          f"{str(np.round(be, 3)):>16}")
print("  합은 거의 안 움직이는데 차만 벌어진다")
print(f"  합은 98 로 나뉘고 차는 2 로 나뉜다 -> {98/2:.0f}배 예민하다")

x = np.array([3.0, 1, -1, -3])
v = np.array([1.0, -1, 1, -1])          # 숨은 요인 (현실에서는 안 보인다)
z = np.array([1.0, 1, -1, -1])          # 도구
y = 2 * x + v

같다("56", y, [7, 1, -1, -7], "관측된 사용량")
be = x @ y / (x @ x)
같다("56", be, 2.2, "통상 계산 (참값 2 보다 위)")
같다("56", be, 2 + (x @ v) / (x @ x), "편의 = x.v / x.x")
print(f"     x.v = {x @ v:.0f}  <- 0 이 아닌 것이 편의의 정체")

e = y - be * x
같다("56", x @ e, 0.0, "잔차가 x 와 직교")
print("     -> 그런데 이 직교는 be 를 그렇게 골랐기 때문이다")

# 어떤 자료를 넣어도 직교가 나온다는 것을 확인한다
난 = np.random.default_rng(56)
어긋남 = 0
for _ in range(20000):
    xx = 난.normal(size=6)
    yy = 난.normal(size=6)                # 아무 관계도 없는 자료
    b = xx @ yy / (xx @ xx)
    if abs(xx @ (yy - b * xx)) > 1e-9:
        어긋남 += 1
print(f"무작위 자료 2만 쌍 중 잔차가 직교하지 않은 경우 : {어긋남}개")
print("  -> 직교는 모형의 증거가 아니다. 나눗셈을 제대로 했다는 뜻뿐이다.")

# 도구변수가 참값을 되찾는다
같다("56", z @ y / (z @ x), 2.0, "도구로 계산")
같다("56", z @ v, 0.0, "도구가 숨은 요인과 직교")
print(f"     z.x = {z @ x:.0f}  <- 0 이 아니어야 나눗셈이 된다")

xh = z * (z @ x) / (z @ z)
같다("56", xh, [2, 2, -2, -2], "1단계 : x 를 z 위로 사영")
같다("56", xh @ y / (xh @ xh), 2.0, "2단계도 같은 답")

zp = np.array([1.0, -1, -1, 1])
print(f"\n     잘못 고른 도구 z' : z'.x = {zp @ x:.0f}  -> 나눗셈 불가")

x3 = np.array([420.0, 320, 260])
Z3 = np.array([[210.0, 32, 78], [84, 160, 26], [42, 64, 104]])
A3 = Z3 @ np.diag(1 / x3)
B3 = np.diag(1 / x3) @ Z3

같다("57", A3, [[.5, .1, .3], [.2, .5, .1], [.1, .2, .4]], "레온티예프 A")
같다("57", B3, np.diag(1 / x3) @ A3 @ np.diag(x3), "B = xhat^-1 A xhat (유사변환)")
같다("57", np.sort(abs(np.linalg.eigvals(A3))), np.sort(abs(np.linalg.eigvals(B3))),
    "고윳값이 같다")

print(f"\n{'':>16} {'대각합':>10} {'행렬식':>12} {'랭크':>6}")
for 이름, M in (("A", A3), ("B", B3)):
    print(f"{이름:>16} {np.trace(M):>10.4f} {np.linalg.det(M):>12.6f} "
          f"{np.linalg.matrix_rank(M):>6}")
print("  -> 유사변환이 지키는 것들")

print(f"\n{'':>16} {'열 합':>26} {'행 합':>26}")
for 이름, M in (("A", A3), ("B", B3)):
    print(f"{이름:>16} {str(np.round(M.sum(axis=0), 4)):>26} "
          f"{str(np.round(M.sum(axis=1), 4)):>26}")
print("  -> A 의 열 합은 균일한데 B 의 행 합은 아니다. 지키지 않는 것들")
print(f"  A12 = {A3[0,1]:.4f},  B12 = {B3[0,1]:.4f}  <- 성분도 다르다")

def 부분행렬식(M):
    n = len(M)
    값 = set()
    for k in range(1, n + 1):
        for r in itertools.combinations(range(n), k):
            for c in itertools.combinations(range(n), k):
                값.add(round(float(np.linalg.det(M[np.ix_(r, c)]))))
    return sorted(값)


M3 = np.array([[1.0, 0, 1], [1, 1, 0], [0, 1, 1]])        # 삼각형 (홀수 사이클)
M4 = np.array([[1.0, 0, 0, 1], [1, 1, 0, 0],
               [0, 1, 1, 0], [0, 0, 1, 1]])                # 4-사이클 (짝수)
같다("59", np.linalg.det(M3), 2.0, "삼각형의 행렬식")
같다("59", np.linalg.solve(M3, np.ones(3)), [.5, .5, .5], "답이 0.5 회")
print(f"     3-사이클 부분행렬식 {부분행렬식(M3)}  <- 2 가 있다")
print(f"     4-사이클 부분행렬식 {부분행렬식(M4)}  <- 0, +-1 뿐 (완전단모듈)")

print()
print(f"{'사이클 길이':>12} {'det':>6} {'부분행렬식':>16} {'해':>26}")
for n in range(3, 8):
    M = np.eye(n) + np.eye(n, k=-1)
    M[0, -1] = 1
    if abs(np.linalg.det(M)) < 1e-9:
        print(f"{n:>12} {0:>6} {'특이':>16} {'유일하지 않다':>26}")
        continue
    해 = np.linalg.solve(M, np.ones(n))
    print(f"{n:>12} {np.linalg.det(M):>6.0f} {str(부분행렬식(M)):>16} "
          f"{str(np.round(해, 3)):>26}")
print("  홀수 사이클에서만 0.5 가 나온다")
print()
print("짝수 사이클은 특이행렬이라 해가 유일하지 않다. 그런데 정수해는 있다 :")
M4b = np.array([[1.0, 0, 0, 1], [1, 1, 0, 0], [0, 1, 1, 0], [0, 0, 1, 1]])
후보 = np.array([1.0, 0, 1, 0])
print(f"  M4 x = {M4b @ 후보}  일 때  x = {후보.astype(int)}")
print("  완전단모듈이 보장하는 것은 '유일한 해'가 아니라 '해가 있으면 정수'다.")

A = np.array([[-0.5, 1.0], [1.0, -0.5]])
d = np.array([1.0, 1])
w = np.linalg.eigvalsh(A)
같다("60", np.sort(w), [-1.5, 0.5], "고윳값")

고정 = np.linalg.solve(np.eye(2) - A, d)
같다("60", 고정, [2, 2], "고정점")
print(f"     이산 조건 max|lambda| = {max(abs(w)):.1f} < 1 ?  "
      f"{'수렴' if max(abs(w)) < 1 else '발산'}")
print(f"     연속 조건 max Re(lambda) = {max(w):.1f} < 1 ?  "
      f"{'수렴' if max(w) < 1 else '발산'}")

from scipy.linalg import expm
x0 = np.array([2.0, 0])
xi = x0.copy()
print(f"\n{'t':>4} {'이산 거리':>12} {'연속 거리':>12}")
for t in range(1, 9):
    xi = A @ xi + d
    xc = expm((A - np.eye(2)) * t) @ (x0 - 고정) + 고정
    if t in (1, 2, 4, 6, 8):
        print(f"{t:>4} {np.abs(xi - 고정).max():>12.4f} "
              f"{np.abs(xc - 고정).max():>12.6f}")
print("  같은 표인데 결론이 반대다. 갱신 방식이 다르기 때문이다.")

# 두 조건이 복소평면에서 어떤 영역인지
print("이산은 원 안(|lambda| < 1), 연속은 수직선 왼쪽(Re lambda < 1)")
print(f"{'lambda':>16} {'|lambda|':>10} {'Re':>7} {'이산':>6} {'연속':>6}")
for lam in (0.5, -1.5, 0.9, -0.9, 1.5, complex(0, 2), complex(-0.5, 0.5)):
    이산 = "수렴" if abs(lam) < 1 else "발산"
    연속 = "수렴" if np.real(lam) < 1 else "발산"
    print(f"{str(lam):>16} {abs(lam):>10.3f} {np.real(lam):>7.2f} "
          f"{이산:>6} {연속:>6}")
print("  원은 수직선 왼쪽에 들어간다 -> 이산이 수렴하면 연속도 수렴한다")
print("  반대는 성립하지 않는다. 그 틈이 이 문제다.")

import collections
import pathlib
import re

뿌리 = pathlib.Path.cwd()
while not (뿌리 / "lectures").exists() and 뿌리 != 뿌리.parent:
    뿌리 = 뿌리.parent

셈 = collections.Counter()
문제수 = 0
for 편 in ("L38", "L39", "L40"):
    p = 뿌리 / "lectures" / 편 / f"{편}_theory.md"
    본문 = p.read_text(encoding="utf-8")
    조각 = 본문.split("## 이 편의 회수표")
    for 줄 in 조각[1].splitlines():
        if not 줄.startswith("|") or 줄.count("|") < 5:
            continue
        칸 = [c.strip() for c in 줄.strip().strip("|").split("|")]
        if not re.match(r"^\[\d+\]\(#prob-", 칸[0]):
            continue                      # 전체 색인 표가 섞이지 않게
        문제수 += 1
        for n in re.findall(r"L(\d{2})", 칸[-1]):
            셈[int(n)] += 1

print(f"문제 {문제수}개,  회수 {sum(셈.values())}건,  닿은 강의 {len(셈)}/37")
print(f"한 번도 안 나온 강의 : {[k for k in range(1, 38) if k not in 셈]}")
print()
print("많이 불린 것 열 개")
for 강의, n in sorted(셈.items(), key=lambda t: (-t[1], t[0]))[:10]:
    print(f"  L{강의:02d} {'#' * n} {n}")

막 = [(1, 4, "1막 계산"), (5, 13, "2막 공간"), (14, 17, "3막 각도"),
      (18, 20, "4막 행렬식"), (21, 25, "5막 고윳값"), (26, 29, "6막 조건 버리기"),
      (30, 34, "7막 변환"), (35, 37, "8막 보강")]
print(f"{'막':>18} {'강의 수':>8} {'회수':>6} {'강의당':>8}")
for 시, 끝, 이름 in 막:
    n = sum(셈.get(k, 0) for k in range(시, 끝 + 1))
    개 = 끝 - 시 + 1
    print(f"{이름:>18} {개:>8} {n:>6} {n/개:>8.1f}")