Lecture 24. 마코브 행렬 — 파이썬 실습

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(24)
print("numpy", np.__version__)

A = np.array([[0.9, 0.2],
              [0.1, 0.8]])


def 마코브인가(M, 눈감아=1e-12):
    """성분이 음이 아니고 열의 합이 1인지 본다."""
    M = np.asarray(M, dtype=float)
    return bool((M >= -눈감아).all() and np.allclose(M.sum(axis=0), 1.0))


def 확률벡터인가(p, 눈감아=1e-12):
    """성분이 음이 아니고 합이 1인지 본다."""
    p = np.asarray(p, dtype=float)
    return bool((p >= -눈감아).all() and np.isclose(p.sum(), 1.0))

print(show_matrix(A, "A"))
print("열의 합 :", A.sum(axis=0))
print("마코브 행렬인가 :", 마코브인가(A))
print()
p = np.array([0.3, 0.7])
for k in range(4):
    print(f"k={k} :  p = {p},  합 = {p.sum():.12f},  확률벡터 {확률벡터인가(p)}")
    p = A @ p

하나 = np.ones(2)
print("1^T A =", 하나 @ A, "  -> 1^T 인가 :", np.allclose(하나 @ A, 하나))
p = np.array([0.3, 0.7])
print("1^T (A p) =", 하나 @ (A @ p))
print("(1^T A) p =", (하나 @ A) @ p, "  <- 같은 수")

def 무작위마코브(n, rng):
    """성분을 아무렇게나 뽑고 열마다 합으로 나눈다."""
    M = rng.random((n, n))
    return M / M.sum(axis=0)

최악 = 0.0
표본 = 500
for _ in range(표본):
    n = int(rng.integers(2, 9))
    M = 무작위마코브(n, rng)
    값 = np.linalg.eigvals(M)
    최악 = max(최악, float(np.abs(값 - 1.0).min()))     # 1 에 가장 가까운 고윳값과의 거리
print(f"마코브 행렬 {표본}개에서 '1 에 가장 가까운 고윳값' 과 1 의 거리")
print(f"  그중 최댓값 : {최악:.3e}   -> 전부 정확히 1 을 가지고 있다")

M = 무작위마코브(5, rng)
print("A   의 특성다항식 계수 :", np.round(np.poly(M), 10))
print("A^T 의 특성다항식 계수 :", np.round(np.poly(M.T), 10))
print("같은가 :", np.allclose(np.poly(M), np.poly(M.T)))
print()
print("A^T 1 =", M.T @ np.ones(5), "  <- 행의 합이 전부 1")

하나 = np.ones(2)
print("A^T (1,1) =", A.T @ 하나, "  -> (1,1) 그대로.  A^T 의 고유벡터")
print("A   (1,1) =", A @ 하나, "  -> 달라졌다.  A 의 고유벡터가 아니다")
print()
값, V = np.linalg.eig(A)
i = int(np.argmin(np.abs(값 - 1)))
정상 = V[:, i].real
print("np.linalg.eig 가 준 고유벡터 :", 정상, "  <- 길이가 1 로 맞춰져 있다")
정상 = 정상 / 정상.sum()
print("성분의 합으로 나눈 것        :", 정상, "  = (2/3, 1/3)")
print("확률벡터인가 :", 확률벡터인가(정상))

최대크기 = 0.0
최대합 = 0.0
for _ in range(300):
    n = int(rng.integers(2, 8))
    M = 무작위마코브(n, rng)
    값, V = np.linalg.eig(M)
    최대크기 = max(최대크기, float(np.abs(값).max()))
    for j in range(n):
        if abs(값[j] - 1) > 1e-8:                        # lambda != 1 인 것만
            최대합 = max(최대합, float(abs(V[:, j].sum())))
print(f"모든 고윳값의 절댓값 중 최댓값 : {최대크기:.12f}   -> 1 을 넘지 않는다")
print(f"lambda != 1 인 고유벡터의 성분 합의 절댓값 중 최댓값 : {최대합:.3e}   -> 0 이다")

값, V = np.linalg.eig(A)
for j in range(2):
    x = V[:, j].real
    한마디 = "0 이 아니다" if abs(값[j].real - 1) < 1e-9 else "0 이다"
    print(f"lambda = {값[j].real:.4f} :  고유벡터 {np.round(x/abs(x).max(), 4)},"
          f"   성분의 합 {x.sum():>10.6f}   <- {한마디}")

S = np.array([[2/3, 1.0],                              # 열1 = 정규화한 정상상태
              [1/3, -1.0]])
print(f"{'초기 u0':>18}{'c1':>10}{'c2':>10}")
for u0 in ([1.0, 0.0], [0.0, 1.0], [0.5, 0.5], [0.13, 0.87]):
    c = np.linalg.solve(S, np.array(u0))
    print(f"{str(u0):>18}{c[0]:>10.6f}{c[1]:>10.6f}")
print()
print("c1 은 사실 초기벡터의 성분 합이다. 확률벡터가 아닌 것을 넣어 보자.")
for u0 in ([2.0, 1.0], [5.0, 0.0], [0.4, 0.4]):
    c1 = np.linalg.solve(S, np.array(u0))[0]
    print(f"  u0 = {str(u0):>14} :  성분 합 {sum(u0):>5.1f},  c1 = {c1:>7.4f}")

정상 = np.array([2/3, 1/3])
u = np.array([1.0, 0.0])
print(f"{'k':>4}{'u_k':>24}{'공식':>24}{'합':>6}")
for k in range(6):
    공식 = 정상 + (1/3) * (0.7 ** k) * np.array([1.0, -1.0])
    print(f"{k:>4}{str(np.round(u, 4)):>24}{str(np.round(공식, 4)):>24}{u.sum():>6.1f}")
    u = A @ u
print()
print("A^50 =")
print(show_matrix(np.linalg.matrix_power(A, 50), "  랭크 1, 두 열이 모두 정상상태"))

P = np.array([[0.8, 0.1, 0.1],
              [0.1, 0.7, 0.3],
              [0.1, 0.2, 0.6]])
print("열의 합 :", P.sum(axis=0), "  마코브인가 :", 마코브인가(P))
값 = np.linalg.eigvals(P)
print("고윳값 :", np.round(np.sort(값.real)[::-1], 6))
print("두 번째 크기 |lambda2| :", np.sort(np.abs(값))[::-1][1])

_, V = np.linalg.eig(P)
s = V[:, int(np.argmin(np.abs(값 - 1)))].real
s = s / s.sum()
print("정상상태 :", np.round(s, 6), "  = (1/3, 7/18, 5/18) 인가 :",
      np.allclose(s, [1/3, 7/18, 5/18]))

꼭짓 = np.array([[0.0, 0.0], [1.0, 0.0], [0.5, np.sqrt(3)/2]])
구름 = rng.dirichlet(np.ones(3), size=200)

테두리 = np.vstack([꼭짓, 꼭짓[:1]])
끝점 = s @ 꼭짓

프레임 = []
걸음들 = list(range(0, 16))
점 = 구름.copy()
for k in 걸음들:
    좌표 = 점 @ 꼭짓
    프레임.append([
        go.Scatter(x=테두리[:, 0], y=테두리[:, 1], mode="lines",
                   line=dict(color="#bbbbbb", width=2), showlegend=False),
        go.Scatter(x=좌표[:, 0], y=좌표[:, 1], mode="markers",
                   marker=dict(color=COLORS["input"], size=6, opacity=0.65),
                   name="분포 200개"),
        go.Scatter(x=[끝점[0]], y=[끝점[1]], mode="markers",
                   marker=dict(color="#d62728", size=14), name="정상상태"),
    ])
    점 = 점 @ P.T                                      # 행마다 P p 를 한꺼번에

배치 = layout2d("걸음을 옮기면", extent=1.0)
배치["height"] = 560
배치["xaxis"] = dict(visible=False, range=[-0.1, 1.1], scaleanchor="y")
배치["yaxis"] = dict(visible=False, range=[-0.1, 0.98])
slider_figure(프레임, 걸음들, 배치, prefix="k = ")

점 = 구름.copy()
print(f"{'k':>4}{'구름의 지름':>14}{'0.7^k':>10}")
초기지름 = None
for k in range(0, 16, 3):
    좌표 = np.linalg.matrix_power(P, k) @ 점.T
    지름 = float(np.abs(좌표 - (s[:, None])).max())
    if 초기지름 is None:
        초기지름 = 지름
    print(f"{k:>4}{지름:>14.6f}{초기지름 * 0.7 ** k:>10.6f}")

B = np.array([[0.0, 1.0], [1.0, 0.0]])              # 반드시 번갈아 오간다
print("B 는 마코브인가 :", 마코브인가(B), "  고윳값 :", np.linalg.eigvals(B).real)
u = np.array([1.0, 0.0])
궤적 = [np.linalg.matrix_power(B, k) @ u for k in range(7)]
print("궤적 :", " -> ".join(f"({v[0]:.0f},{v[1]:.0f})" for v in 궤적))
print("-> 정상상태 (0.5, 0.5) 는 있지만 (1,0) 에서 출발하면 닿지 못한다")
print("   B (0.5,0.5) =", B @ np.array([0.5, 0.5]), " <- 정상상태 자체는 맞다")
print()
I = np.eye(3)
print("I 는 마코브인가 :", 마코브인가(I), "  고윳값 :", np.linalg.eigvals(I).real)
print("-> 모든 확률벡터가 정상상태. 유일하지 않다")
print("   아무 것이나 :", I @ np.array([0.2, 0.3, 0.5]))

최대둘째 = 0.0
for _ in range(300):
    n = int(rng.integers(2, 8))
    M = 무작위마코브(n, rng)                          # 성분이 전부 양수
    크기 = np.sort(np.abs(np.linalg.eigvals(M)))[::-1]
    최대둘째 = max(최대둘째, float(크기[1]))
print(f"성분이 전부 양수인 마코브 행렬 300개의 |lambda2| 중 최댓값 : {최대둘째:.6f}")
print("전부 1 보다 작은가 :", 최대둘째 < 1.0)

링크 = {0: [1, 2], 1: [2], 2: [0], 3: [0, 1],
        4: [],                       # 막다른 페이지
        5: [4, 6], 6: [5], 7: [0, 5]}


def 링크행렬(링크):
    """A[i, j] = 1/d_j  (j 가 i 로 링크했을 때). 막다른 열은 0 으로 둔다."""
    n = len(링크)
    A = np.zeros((n, n))
    for j, 나가는곳 in 링크.items():
        for i in 나가는곳:
            A[i, j] = 1.0 / len(나가는곳)
    return A

A0 = 링크행렬(링크)
print("열의 합 :", A0.sum(axis=0))
print("마코브인가 :", 마코브인가(A0), "  <- 4번 열이 0 이라 아니다")

A1 = A0.copy()
막다른 = A1.sum(axis=0) == 0
A1[:, 막다른] = 1.0 / len(링크)
print("고친 뒤 열의 합 :", A1.sum(axis=0), "  마코브인가 :", 마코브인가(A1))
값 = np.linalg.eigvals(A1)
print("고윳값의 절댓값 :", np.round(np.sort(np.abs(값))[::-1], 4))
print("|lambda2| =", np.sort(np.abs(값))[::-1][1], " <- 1 에 가까우면 느리다")

def 구글행렬(A, 알파=0.85):
    """G = alpha A + (1-alpha) J/n. 모든 성분이 양수가 된다."""
    n = A.shape[0]
    return 알파 * A + (1 - 알파) * np.ones((n, n)) / n

알파 = 0.85
G = 구글행렬(A1, 알파)
print("마코브인가 :", 마코브인가(G), "   모든 성분이 양수인가 :", (G > 0).all())
크기 = np.sort(np.abs(np.linalg.eigvals(G)))[::-1]
print("고윳값의 절댓값 :", np.round(크기, 6))
print(f"|lambda2| = {크기[1]:.6f}  <=  alpha = {알파} 인가 :", 크기[1] <= 알파 + 1e-9)

값, V = np.linalg.eig(G)
순위 = V[:, int(np.argmin(np.abs(값 - 1)))].real
순위 = 순위 / 순위.sum()

# 고윳값 분해 대신 거듭 곱하기만 해도 같은 답이 나온다 (L22 의 지배 고윳값)
p = np.ones(len(링크)) / len(링크)
for _ in range(60):
    p = G @ p
print("고유벡터로 구한 순위 :", np.round(순위, 5))
print("거듭 곱해 구한 순위 :", np.round(p, 5))
print("최대 차이 :", np.abs(순위 - p).max())
print()
for i in np.argsort(순위)[::-1]:
    막대 = "#" * int(round(순위[i] * 200))
    print(f"  page {i} : {순위[i]:.4f}  {막대}")

print(f"{'alpha':>8}{'|lambda2|':>12}{'60걸음 오차':>14}")
for 알파 in (0.5, 0.7, 0.85, 0.95, 0.99):
    G2 = 구글행렬(A1, 알파)
    크기 = np.sort(np.abs(np.linalg.eigvals(G2)))[::-1]
    값2, V2 = np.linalg.eig(G2)
    참 = V2[:, int(np.argmin(np.abs(값2 - 1)))].real
    참 = 참 / 참.sum()
    p = np.ones(len(링크)) / len(링크)
    for _ in range(60):
        p = G2 @ p
    print(f"{알파:>8}{크기[1]:>12.6f}{np.abs(p - 참).max():>14.2e}")

Ac = np.array([[-1.0, 2.0], [1.0, -2.0]])           # L23 의 앵커. 열의 합 0
print("Ac 의 열의 합 :", Ac.sum(axis=0))
print(f"{'h':>6}{'열의 합':>16}{'성분>=0':>10}{'고윳값':>24}{'기대 1+h*lambda':>20}")
for h in (0.05, 0.1, 0.25, 0.5, 0.6):
    Ad = np.eye(2) + h * Ac
    값 = np.sort(np.linalg.eigvals(Ad).real)[::-1]
    print(f"{h:>6}{str(Ad.sum(axis=0)):>16}{str((Ad>=0).all()):>10}"
          f"{str(np.round(값, 4)):>24}{str([1.0, round(1-3*h, 4)]):>20}")
print()
print("h = 0.1 이 이번 강의의 앵커인가 :", np.allclose(np.eye(2) + 0.1*Ac, A))

for h in (0.1, 0.3, 0.5):
    Ad = np.eye(2) + h * Ac
    for x, 이름 in ((np.array([2.0, 1.0]), "(2,1)"), (np.array([1.0, -1.0]), "(1,-1)")):
        결과 = Ad @ x
        배율 = 결과[0] / x[0]
        print(f"h={h}:  Ad {이름:>6} = {np.round(결과,4)} = {배율:>6.2f} x "
              f"  (1 + h*lambda = {1 + h*(0 if 이름=='(2,1)' else -3):.2f})")