Lecture 27. 양의 정부호 행렬과 최소값 — 파이썬 실습

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, layout2d, show_matrix, slider_figure

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

def 피벗들(S):
    """행 교환 없이 소거해서 피벗을 뽑는다. 대칭이 깨지지 않는다."""
    U = np.array(S, dtype=float)
    n = U.shape[0]
    for k in range(n - 1):
        if abs(U[k, k]) < 1e-14:
            return None                       # 교환 없이는 진행할 수 없다
        U[k + 1:] -= np.outer(U[k + 1:, k] / U[k, k], U[k])
    return np.diag(U).copy()


def 선행주소행렬식(S):
    """왼쪽 위 모서리부터 잘라 낸 부분행렬들의 행렬식."""
    S = np.asarray(S, dtype=float)
    return np.array([np.linalg.det(S[:k + 1, :k + 1]) for k in range(S.shape[0])])

def 조건1_고윳값(S, 눈감아=1e-10):
    return bool((np.linalg.eigvalsh(S) > 눈감아).all())


def 조건2_주소행렬식(S, 눈감아=1e-10):
    return bool((선행주소행렬식(S) > 눈감아).all())


def 조건3_피벗(S, 눈감아=1e-10):
    p = 피벗들(S)
    return bool(p is not None and (p > 눈감아).all())


def 조건4_에너지(S, rng, 표본=2000):
    """아무 x 나 넣어 봐서 음수가 나오는지 본다. 반례 찾기용이지 증명은 아니다."""
    x = rng.normal(size=(np.shape(S)[0], 표본))
    return bool((np.einsum("ik,ij,jk->k", x, np.asarray(S, float), x) > 0).all())


def 조건5_촐레스키(S):
    """S = R R^T 를 찾을 수 있으면 양정치. 실패하면 그 자리에서 끝난다."""
    try:
        np.linalg.cholesky(S)
        return True
    except np.linalg.LinAlgError:
        return False

S = np.array([[5.0, 4.0],
              [4.0, 5.0]])
print(show_matrix(S, "앵커 S"))
print("고윳값           :", np.linalg.eigvalsh(S), "  (기대 1, 9)")
print("선행 주소행렬식  :", 선행주소행렬식(S), "  (기대 5, 9)")
print("피벗             :", 피벗들(S), "  (기대 5, 1.8)")
print("피벗의 곱        :", np.prod(피벗들(S)), " = det =", np.linalg.det(S))
print()
for 이름, 판정 in (("① 고윳값", 조건1_고윳값), ("② 주소행렬식", 조건2_주소행렬식),
                  ("③ 피벗", 조건3_피벗), ("⑤ 촐레스키", 조건5_촐레스키)):
    print(f"  {이름:>14} : {판정(S)}")
print(f"  {'④ 에너지':>14} : {조건4_에너지(S, rng)}")

T = np.array([[2.0, -1.0, 0.0],
              [-1.0, 2.0, -1.0],
              [0.0, -1.0, 2.0]])
print("고윳값          :", np.linalg.eigvalsh(T))
print("  = 2-sqrt2, 2, 2+sqrt2 인가 :",
      np.allclose(np.linalg.eigvalsh(T), sorted([2-np.sqrt(2), 2, 2+np.sqrt(2)])))
print("선행 주소행렬식 :", 선행주소행렬식(T), "  (기대 2, 3, 4)")
print("피벗            :", 피벗들(T), "  (기대 2, 1.5, 1.3333)")
print("피벗의 곱       :", np.prod(피벗들(T)), " = det =", np.linalg.det(T))

표본 = 600
정확어긋 = []
양, 음 = 0, 0
표본들 = []
for _ in range(표본):
    n = int(rng.integers(2, 6))
    B = rng.normal(size=(n, n))
    M = B + B.T
    if rng.random() < 0.5:
        M = M + n * np.eye(n)                  # 절반쯤은 양정치가 되게 밀어 준다
    표본들.append(M)
    판정 = (조건1_고윳값(M), 조건2_주소행렬식(M), 조건3_피벗(M), 조건5_촐레스키(M))
    if len(set(판정)) != 1:
        정확어긋.append((M, 판정))
    if 판정[0]:
        양 += 1
    else:
        음 += 1

print(f"대칭행렬 {표본}개 : 양정치 {양}개, 아닌 것 {음}개  -> 양쪽이 다 나온다")
print(f"네 판정이 갈린 경우 : {len(정확어긋)}건")

에너지어긋 = []
for M in 표본들:
    if 조건1_고윳값(M) != 조건4_에너지(M, rng):
        에너지어긋.append(M)

print(f"표본으로 판정한 에너지가 실제와 어긋난 경우 : {len(에너지어긋)}건 / {표본}")
print()
for M in 에너지어긋[:3]:
    값 = np.linalg.eigvalsh(M)
    print(f"  n={len(M)}, 고윳값 {np.round(값, 4)}")
    print(f"     음수인 고윳값의 크기 / 가장 큰 고윳값 = "
          f"{abs(값).min()/abs(값).max():.1e}   <- 아주 얇은 방향 하나뿐")

def 조건4_정확(S, 눈감아=1e-10):
    """최소 고유벡터 방향만 보면 된다. 거기가 에너지가 가장 작은 곳이다."""
    값, V = np.linalg.eigh(S)
    v = V[:, 0]
    return bool(v @ np.asarray(S, float) @ v > 눈감아)


어긋 = sum(1 for M in 표본들 if 조건1_고윳값(M) != 조건4_정확(M))
print(f"최소 고유벡터로 판정하면 어긋난 경우 : {어긋}건 / {표본}")
print()
for M in 에너지어긋[:3]:
    값, V = np.linalg.eigh(M)
    v = V[:, 0]
    print(f"  최소 고유벡터 방향의 에너지 : {v @ M @ v:>11.4e}"
          f"   = 최소 고윳값 {값[0]:>11.4e}")

A = np.array([[3.0, 7.0],
              [-1.0, 2.0]])                      # 대칭이 아니다
대칭부 = (A + A.T) / 2
비대칭부 = (A - A.T) / 2
print(show_matrix(A, "A  (비대칭)"))
print(show_matrix(대칭부, "대칭 부분"))
print(show_matrix(비대칭부, "비대칭 부분"))

x = rng.normal(size=(2, 5))
왼 = np.array([xx @ A @ xx for xx in x.T])
오 = np.array([xx @ 대칭부 @ xx for xx in x.T])
비 = np.array([xx @ 비대칭부 @ xx for xx in x.T])
print()
print("x^T A x        :", np.round(왼, 6))
print("x^T (대칭부) x :", np.round(오, 6), "  <- 같다")
print("x^T (비대칭부) x:", np.round(비, 12), "  <- 언제나 0")

print(f"{'t':>6}{'고윳값':>20}{'det':>10}{'판정':>28}")
for t in (0.0, 1.0, 1.9, 2.0, 2.5, 4.0):
    M = np.array([[2.0, t], [t, 2.0]])
    값 = np.linalg.eigvalsh(M)
    if (값 > 1e-9).all():
        판정 = "양의 정부호 — 그릇"
    elif (값 > -1e-9).all():
        판정 = "준정부호 — 골짜기"
    else:
        판정 = "부정부호 — 안장"
    print(f"{t:>6}{str(np.round(값,3)):>20}{np.linalg.det(M):>10.2f}{판정:>28}")

격자 = np.linspace(-2.2, 2.2, 160)
X, Y = np.meshgrid(격자, 격자)
t들 = np.round(np.linspace(0.0, 4.0, 21), 2)

프레임 = []
for t in t들:
    M = np.array([[2.0, t], [t, 2.0]])
    Z = 2 * X**2 + 2 * t * X * Y + 2 * Y**2
    값 = np.linalg.eigvalsh(M)
    if (값 > 1e-9).all():
        말, 색 = "양의 정부호 — 그릇", "#2ca02c"
    elif (값 > -1e-9).all():
        말, 색 = "준정부호 — 골짜기", "#7f7f7f"
    else:
        말, 색 = "부정부호 — 안장", "#d62728"
    프레임.append([
        go.Contour(x=격자, y=격자, z=Z, contours=dict(coloring="lines",
                   start=-16, end=32, size=3), line=dict(width=2),
                   colorscale="RdBu", reversescale=True, showscale=False),
        go.Scatter(x=[0], y=[-2.75], mode="text",
                   text=[f"고윳값 {np.round(값,2)}   {말}"],
                   textfont=dict(size=16, color=색), showlegend=False),
    ])

배치 = layout2d("S(t) = [[2, t], [t, 2]] 의 등고선", extent=3.0)
배치["height"] = 600
slider_figure(프레임, t들, 배치, prefix="t = ", initial=5)

x = rng.normal(size=(2, 6))
왼 = np.array([xx @ S @ xx for xx in x.T])
오 = 5 * (x[0] + 0.8 * x[1]) ** 2 + 1.8 * x[1] ** 2
print("x^T S x      :", np.round(왼, 6))
print("제곱완성 결과 :", np.round(오, 6))
print("최대 차이 :", np.abs(왼 - 오).max())

def LDL(S):
    """행 교환 없이 소거하며 배수를 L 에, 피벗을 D 에 모은다."""
    U = np.array(S, dtype=float)
    n = U.shape[0]
    L = np.eye(n)
    for k in range(n - 1):
        for i in range(k + 1, n):
            배수 = U[i, k] / U[k, k]
            L[i, k] = 배수
            U[i] -= 배수 * U[k]
    return L, np.diag(np.diag(U))

L, D = LDL(S)
print(show_matrix(L, "L  (배수)"))
print(show_matrix(D, "D  (피벗)"))
print("L D L^T = S 인가 :", np.allclose(L @ D @ L.T, S))
print()
L3, D3 = LDL(T)
print(show_matrix(L3, "3x3 의 L"))
print(show_matrix(D3, "3x3 의 D"))
print("복원되는가 :", np.allclose(L3 @ D3 @ L3.T, T))
print("D 의 대각이 피벗과 같은가 :", np.allclose(np.diag(D3), 피벗들(T)))

x3 = rng.normal(size=(3, 5))
y3 = L3.T @ x3
왼 = np.array([xx @ T @ xx for xx in x3.T])
오 = (np.diag(D3)[:, None] * y3 ** 2).sum(axis=0)
print("x^T T x    :", np.round(왼, 6))
print("sum d_i y_i^2:", np.round(오, 6))
print("최대 차이 :", np.abs(왼 - 오).max())

값, Q = np.linalg.eigh(S)
print("고윳값 :", 값)
print("반지름 1/sqrt(lambda) :", 1 / np.sqrt(값), "  (기대 1, 1/3)")
print()
for j in range(2):
    끝 = Q[:, j] / np.sqrt(값[j])                # 축의 끝점
    print(f"  lambda={값[j]:.0f} 축 끝 {np.round(끝,4)},  "
          f"거기서 x^T S x = {끝 @ S @ 끝:.6f}  <- 1 이어야 한다")

t = np.linspace(0, 2 * np.pi, 400)
타원 = Q @ np.diag(1 / np.sqrt(값)) @ np.vstack([np.cos(t), np.sin(t)])
자료 = [go.Scatter(x=타원[0], y=타원[1], mode="lines", name="x^T S x = 1",
                  line=dict(color=COLORS["output"], width=4))]
for j, 색 in enumerate(("#2ca02c", "#d62728")):
    끝 = Q[:, j] / np.sqrt(값[j])
    자료.append(go.Scatter(x=[-끝[0], 끝[0]], y=[-끝[1], 끝[1]], mode="lines",
                          line=dict(color=색, width=3),
                          name=f"q{j+1} 축,  lambda={값[j]:.0f},  "
                               f"반지름 {1/np.sqrt(값[j]):.3f}"))
go.Figure(data=자료, layout=layout2d("축은 고유벡터, 반지름은 1/sqrt(lambda)",
                                    extent=1.3))

print(f"{'A 의 크기':>12}{'열 독립':>10}{'최소 고윳값':>16}{'양정치':>10}{'가역':>8}")
for m, n in ((5, 3), (4, 4), (3, 5), (6, 2)):
    A = rng.normal(size=(m, n))
    G = A.T @ A
    독립 = np.linalg.matrix_rank(A) == n
    최소 = np.linalg.eigvalsh(G).min()
    print(f"{f'{m}x{n}':>12}{str(독립):>10}{최소:>16.3e}"
          f"{str(조건1_고윳값(G)):>10}{str(abs(np.linalg.det(G)) > 1e-12):>8}")

# 에너지가 정말 |Ax|^2 인가
A = rng.normal(size=(5, 3))
x = rng.normal(size=(3, 4))
왼 = np.array([xx @ (A.T @ A) @ xx for xx in x.T])
오 = (np.linalg.norm(A @ x, axis=0)) ** 2
print("x^T A^T A x :", np.round(왼, 8))
print("|A x|^2     :", np.round(오, 8))
print("최대 차이 :", np.abs(왼 - 오).max())
print()
print("A^T A 가 대칭인가 :", np.allclose(A.T @ A, (A.T @ A).T))

A = rng.normal(size=(5, 3))
A[:, 2] = A[:, 0] + 2 * A[:, 1]                  # 3열을 종속으로
G = A.T @ A
print("A 의 랭크 :", np.linalg.matrix_rank(A))
print("A^T A 의 고윳값 :", np.linalg.eigvalsh(G))
print("0 인 고윳값의 개수 :", int((np.abs(np.linalg.eigvalsh(G)) < 1e-9).sum()))
print("  = 3 - 랭크 =", 3 - np.linalg.matrix_rank(A))
print()
_, V = np.linalg.eigh(G)
평평 = V[:, 0]
print("평평한 방향 :", np.round(평평 / abs(평평).max(), 4), "  (1, 2, -1 의 배수)")
print("  그 방향의 에너지 :", 평평 @ G @ 평평)
print("  A 를 곱하면 :", np.round(A @ 평평, 12))

def 경사하강(M, 시작=(1.1, 0.35), 눈감아=1e-6, 최대=200_000):
    """f = x^T M x 를 최적 고정 걸음으로 내려간다. 걸음 수와 경로를 돌려준다."""
    값 = np.linalg.eigvalsh(M)
    보폭 = 1.0 / (값.max() + 값.min())            # 표준 최적 걸음
    v = np.array(시작, dtype=float)
    경로 = [v.copy()]
    while np.linalg.norm(v) > 눈감아 and len(경로) <= 최대:
        v = v - 보폭 * (2 * M @ v)                # 기울기는 2 M v
        경로.append(v.copy())
    return len(경로) - 1, np.array(경로)

print(f"{'':>22}{'조건수':>10}{'(k-1)/(k+1)':>14}{'걸음':>8}{'이론':>8}")
for 이름, M in (("[[1,0],[0,1.2]]", np.array([[1.0, 0.0], [0.0, 1.2]])),
                ("앵커 [[5,4],[4,5]]", S),
                ("[[1,0],[0,30]]", np.array([[1.0, 0.0], [0.0, 30.0]])),
                ("[[1,0],[0,300]]", np.array([[1.0, 0.0], [0.0, 300.0]]))):
    값 = np.linalg.eigvalsh(M)
    k = 값.max() / 값.min()
    r = (k - 1) / (k + 1)
    걸음, _ = 경사하강(M)
    이론 = np.log(1e-6 / np.linalg.norm([1.1, 0.35])) / np.log(r)
    print(f"{이름:>22}{k:>10.1f}{r:>14.4f}{걸음:>8}{이론:>8.0f}")

자료 = []
for 이름, M, 색 in (("조건수 1.2", np.array([[1.0, 0.0], [0.0, 1.2]]), "#2ca02c"),
                   ("조건수 9 (앵커)", S, COLORS["input"]),
                   ("조건수 30", np.array([[1.0, 0.0], [0.0, 30.0]]), "#d62728")):
    _, 경로 = 경사하강(M)
    보일 = 경로[:40]
    자료.append(go.Scatter(x=보일[:, 0], y=보일[:, 1], mode="lines+markers",
                          name=이름, line=dict(color=색, width=2),
                          marker=dict(size=5)))
go.Figure(data=자료, layout=layout2d("같은 자리에서 출발한 세 궤적 (앞 40걸음)",
                                    extent=1.3))

R = np.linalg.cholesky(T)
print(show_matrix(R, "R  (하삼각)"))
print("R R^T = T 인가 :", np.allclose(R @ R.T, T))
print("R 의 대각 :", np.round(np.diag(R), 6))
print("그 제곱   :", np.round(np.diag(R) ** 2, 6), " = 피벗", np.round(피벗들(T), 6))
print()
print("A = R^T 로 두면 A^T A = T 인가 :", np.allclose(R @ R.T, T),
      " <- 조건 ⑤ 가 구성적으로 확인된다")

나쁨 = np.array([[1.0, 10.0], [10.0, 1.0]])
print(show_matrix(나쁨, "[[1,10],[10,1]]  — 대각은 양수인데"))
print("고윳값          :", np.linalg.eigvalsh(나쁨), "  <- 하나가 음수")
print("선행 주소행렬식 :", 선행주소행렬식(나쁨), "  <- 두 번째에서 걸린다")
print("피벗            :", 피벗들(나쁨))
print("촐레스키        :", 조건5_촐레스키(나쁨))
print()
x = np.array([1.0, -1.0])
print("x = (1,-1) 의 에너지 :", x @ 나쁨 @ x, "  <- 음수. 반례를 손에 쥐었다")

print(f"{'c':>6}{'[[1,c],[c,1]] 의 고윳값':>26}{'양정치':>10}")
for c in (0.0, 0.5, 0.9, 1.0, 1.5, 10.0):
    M = np.array([[1.0, c], [c, 1.0]])
    print(f"{c:>6}{str(np.round(np.linalg.eigvalsh(M), 3)):>26}"
          f"{str(조건1_고윳값(M)):>10}")
print()
print("고윳값이 1 ± c 이므로 |c| < 1 일 때만 양정치이다.")