Lecture 30. 선형변환과 그 행렬 — 파이썬 실습

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

from linalg_viz import COLORS, layout2d, show_matrix, slider_figure

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

def 행렬만들기(T, 입력기저, 출력기저):
    """선형변환 T 의 행렬을 정의대로 만든다.

    입력기저 : 열이 기저벡터인 행렬 (n 개)
    출력기저 : 열이 기저벡터인 행렬 (m 개)
    """
    입력기저 = np.asarray(입력기저, dtype=float)
    출력기저 = np.asarray(출력기저, dtype=float)
    열들 = []
    for j in range(입력기저.shape[1]):
        간곳 = T(입력기저[:, j])                        # j 번째 기저벡터가 가는 곳
        좌표, *_ = np.linalg.lstsq(출력기저, 간곳, rcond=None)   # 출력 기저로 전개
        열들.append(좌표)
    return np.column_stack(열들)


def 깔끔(M, 눈감아=1e-10):
    """눈으로 볼 때 방해가 되는 아주 작은 값을 0 으로 눌러 준다.

    lstsq 는 실수 나눗셈을 하므로 결과가 정확히 0 이 되는 일은 드물다.
    1e-16 같은 값은 '0 인데 반올림 오차가 남은 것' 이라고 읽으면 된다.
    """
    M = np.asarray(M, dtype=float).copy()
    M[np.abs(M) < 눈감아] = 0.0
    return M

A = np.array([[2.0, -1.0], [1.0, 1.0]])
만든것 = 행렬만들기(lambda v: A @ v, np.eye(2), np.eye(2))
print(show_matrix(A, "원래 A"))
print(show_matrix(만든것, "기저벡터에게 물어서 만든 것"))
print("같은가 :", np.allclose(A, 만든것))
print()
print("e1 이 가는 곳 :", A @ np.array([1.0, 0.0]), " = A 의 1열 :", A[:, 0])
print("e2 이 가는 곳 :", A @ np.array([0.0, 1.0]), " = A 의 2열 :", A[:, 1])

def 선형인가(T, n, 횟수=500, 눈감아=1e-9):
    """무작위 v, w, c, d 로 선형성을 두들겨 본다."""
    for _ in range(횟수):
        v, w = rng.normal(size=n), rng.normal(size=n)
        c, d = rng.normal(), rng.normal()
        왼 = T(c * v + d * w)
        오 = c * np.asarray(T(v), float) + d * np.asarray(T(w), float)
        if not np.allclose(왼, 오, atol=눈감아):
            return False
    return True

회전 = np.array([[0.0, -1.0], [1.0, 0.0]])
후보 = (("행렬 곱하기", lambda v: A @ v),
        ("90도 회전", lambda v: 회전 @ v),
        ("y=x 로 투영", lambda v: np.array([v.sum()/2, v.sum()/2])),
        ("평행이동 +(1,2)", lambda v: v + np.array([1.0, 2.0])),
        ("길이 재기", lambda v: np.array([np.linalg.norm(v), 0.0])),
        ("제곱하기", lambda v: v**2))
print(f"{'':>18}{'T(0) = 0 인가':>16}{'선형인가':>12}")
for 이름, T in 후보:
    영 = np.allclose(T(np.zeros(2)), 0.0)
    print(f"{이름:>18}{str(영):>16}{str(선형인가(T, 2)):>12}")

v, w = np.array([3.0, 0.0]), np.array([0.0, 4.0])
평행 = lambda x: x + np.array([1.0, 2.0])
print("T(v+w)    =", 평행(v + w))
print("T(v)+T(w) =", 평행(v) + 평행(w), "  <- 이동벡터를 두 번 더했다")
print()
print("차원을 하나 올려 보자 (동차좌표)")
H = np.array([[1.0, 0.0, 1.0],
              [0.0, 1.0, 2.0],
              [0.0, 0.0, 1.0]])
print(show_matrix(H, "H"))
print("(3,0,1) -> ", H @ np.array([3.0, 0.0, 1.0]), "  = (4,2) 로 옮겨졌다")
print("R^3 에서는 선형인가 :", 선형인가(lambda x: H @ x, 3))
print()
print("-> 평행이동을 선형으로 만든 것이 아니라,")
print("   한 차원 위에서 선형인 것을 z=1 평면에서 잘라 본 것이다.")

B = np.array([[0.0, 1.0], [1.0, 0.0]])          # x, y 맞바꾸기
합성 = 행렬만들기(lambda v: A @ (B @ v), np.eye(2), np.eye(2))
print(show_matrix(합성, "합성을 정의대로 만든 것"))
print(show_matrix(A @ B, "A @ B"))
print("같은가 :", np.allclose(합성, A @ B))
print()
print("순서를 바꾸면 :", np.allclose(행렬만들기(lambda v: B @ (A @ v),
                                          np.eye(2), np.eye(2)), B @ A))
print("A@B 와 B@A 가 같은가 :", np.allclose(A @ B, B @ A),
      "  <- 합성은 순서를 탄다")

def 미분행렬(n):
    """P_(n-1) -> P_(n-2) 의 미분 행렬. 기저는 1, x, ..., x^(n-1)."""
    D = np.zeros((n - 1, n))
    for j in range(1, n):
        D[j - 1, j] = j                       # d/dx x^j = j x^(j-1)
    return D


def 적분행렬(n):
    """P_(n-2) -> P_(n-1) 의 적분 행렬. 적분상수는 붙이지 않는다."""
    S = np.zeros((n, n - 1))
    for j in range(n - 1):
        S[j + 1, j] = 1.0 / (j + 1)           # int x^j dx = x^(j+1)/(j+1)
    return S

D, S = 미분행렬(4), 적분행렬(4)
print(show_matrix(D, "D  (3x4)  미분"))
print(show_matrix(S * 6, "6 S  (4x3)  적분 (6배 해서 정수로 보기)"))
print()
p = np.array([5.0, 3.0, 0.0, 2.0])            # 5 + 3x + 2x^3
print("p  =", p, " -> 5 + 3x + 2x^3")
print("Dp =", D @ p, " -> 3 + 6x^2      (p' = 3 + 6x^2 이 맞는가)")

def 다항식미분(계수):
    """계수 벡터를 미분한다. (numpy 의 polyder 를 쓰지 않고 직접)"""
    계수 = np.asarray(계수, dtype=float)
    return np.array([j * 계수[j] for j in range(1, len(계수))])

정의대로 = 행렬만들기(다항식미분, np.eye(4), np.eye(3))
print(show_matrix(정의대로, "기저 1, x, x^2, x^3 에게 물어 만든 행렬"))
print("D 와 같은가 :", np.allclose(정의대로, D))
print()
print("각 기저벡터가 가는 곳을 하나씩 :")
for j, 이름 in enumerate(["1", "x", "x^2", "x^3"]):
    기저 = np.eye(4)[:, j]
    print(f"  d/dx {이름:>4}  ->  계수 {다항식미분(기저)}   = D 의 {j+1}열 {D[:, j]}")

영 = null_space(D)
print("D 의 영공간 차원 :", 영.shape[1])
print("기저 :", np.round(영.ravel(), 6), " -> 상수함수")
print("정말 죽는가 :", D @ 영.ravel())
print()
print(f"{'공간':>8}{'D 의 크기':>12}{'랭크':>6}{'영공간 차원':>12}")
for n in (4, 6, 8, 12):
    Dn = 미분행렬(n)
    r = np.linalg.matrix_rank(Dn)
    print(f"{'P_' + str(n-1):>8}{str(Dn.shape):>12}{r:>6}{n - r:>12}")
print()
print("-> 차원을 아무리 올려도 죽는 것은 상수 하나뿐이다.")
print("   그래서 미분 행렬은 어느 차원에서도 정방이 될 수 없다.")

print(show_matrix(D @ S, "D S"))
print("I_3 인가 :", np.allclose(D @ S, np.eye(3)))
print()
print(show_matrix(S @ D, "S D"))
print("I_4 인가 :", np.allclose(S @ D, np.eye(4)))
print("어디가 다른가 :", np.argwhere(~np.isclose(S @ D, np.eye(4))).tolist(),
      " <- 왼쪽 위 한 칸")

print("p 를 미분했다가 적분해 보자")
print("  p      =", p)
print("  D p    =", D @ p)
print("  S D p  =", S @ (D @ p), "  <- 상수 5 가 사라졌다")
print()
print("적분했다가 미분하면")
q = np.array([3.0, 0.0, 6.0])                 # 3 + 6x^2
print("  q      =", q)
print("  S q    =", S @ q)
print("  D S q  =", D @ (S @ q), "  <- 그대로 돌아왔다")

P = S @ D
print("P^2 = P 인가 :", np.allclose(P @ P, P))
print("P^T = P 인가 :", np.allclose(P.T, P))
print("고윳값 :", np.round(np.linalg.eigvalsh(P), 12), " <- 0 과 1 뿐 (L21)")
print("랭크   :", np.linalg.matrix_rank(P), " = D 의 랭크 :",
      np.linalg.matrix_rank(D))
print()
print("무엇으로의 투영인가 : 상수항이 없는 다항식들")
print("  (0,3,0,2) 를 넣으면 :", P @ np.array([0.0, 3.0, 0.0, 2.0]), " 그대로")
print("  (5,0,0,0) 을 넣으면 :", P @ np.array([5.0, 0.0, 0.0, 0.0]), " 통째로 죽는다")

print(show_matrix(np.linalg.pinv(D) * 6, "6 * pinv(D)"))
print(show_matrix(S * 6, "6 * S  (우리가 만든 적분 행렬)"))
print("같은가 :", np.allclose(np.linalg.pinv(D), S))
print()
print(f"{'n':>4}{'pinv(D) == S 인가':>20}")
for n in (4, 6, 8, 12):
    print(f"{n:>4}{str(np.allclose(np.linalg.pinv(미분행렬(n)), 적분행렬(n))):>20}")
print()
print("L33 에서 이 행렬에 이름을 붙인다.")

투영 = lambda v: np.array([v.sum()/2, v.sum()/2])
기저들 = (("표준기저 e1, e2", np.eye(2)),
          ("직선에 맞춘 (1,1), (1,-1)", np.array([[1.0, 1.0], [1.0, -1.0]])),
          ("순서를 바꾼 (1,-1), (1,1)", np.array([[1.0, 1.0], [-1.0, 1.0]])),
          ("길이를 바꾼 (2,2), (3,-3)", np.array([[2.0, 3.0], [2.0, -3.0]])),
          ("비스듬한 (1,0), (1,1)", np.array([[1.0, 1.0], [0.0, 1.0]])))
for 이름, M in 기저들:
    표현 = 행렬만들기(투영, M, M)
    print(show_matrix(깔끔(표현), 이름))

print(f"{'':>28}{'대각합':>10}{'행렬식':>10}{'랭크':>8}{'고윳값':>16}")
for 이름, M in 기저들:
    표현 = 행렬만들기(투영, M, M)
    고 = np.sort(np.linalg.eigvals(표현).real)
    print(f"{이름:>28}{np.trace(표현):>10.4f}{np.linalg.det(표현):>10.4f}"
          f"{np.linalg.matrix_rank(표현):>8}{str(np.round(고, 4)):>16}")
print()
print("-> 대각합, 행렬식, 랭크, 고윳값은 좌표계를 타지 않는다.")
print("   이것들은 변환의 성질이지 기저의 성질이 아니다.")

print("대칭성은 어떤가")
for 이름, M in 기저들:
    표현 = 행렬만들기(투영, M, M)
    print(f"  {이름:>28} : 대칭인가 {np.allclose(표현, 표현.T)}")
print()
print("-> 비스듬한 기저에서는 대칭성이 깨진다.")
print("   '대칭' 은 기저가 직교라는 약속 위에서만 뜻이 있다.")

Pl = np.array([[0.5, 0.5], [0.5, 0.5]])       # 표준기저에서의 투영
M = np.array([[1.0, 1.0], [1.0, -1.0]])       # 새 기저를 열에 세운다

print("새 좌표 (3,1) 은 어떤 벡터인가")
print("  3*(1,1) + 1*(1,-1) =", 3*M[:, 0] + 1*M[:, 1])
print("  M @ (3,1)          =", M @ np.array([3.0, 1.0]), "  같다")
print()
print("표준 좌표 (4,2) 를 새 좌표로 :", np.linalg.inv(M) @ np.array([4.0, 2.0]))
print()
print(show_matrix(np.linalg.inv(M) @ Pl @ M, "M^-1 P M"))
print(show_matrix(행렬만들기(투영, M, M), "정의대로 만든 것"))
print("같은가 :", np.allclose(np.linalg.inv(M) @ Pl @ M, 행렬만들기(투영, M, M)))

c = np.array([3.0, 1.0])                       # 새 좌표
단계 = (("① 새 좌표", c),
        ("② M 을 곱해 표준 좌표로", M @ c),
        ("③ P 로 변환 (표준 좌표에서)", Pl @ (M @ c)),
        ("④ M^-1 로 다시 새 좌표로", np.linalg.inv(M) @ (Pl @ (M @ c))))
for 이름, 값 in 단계:
    print(f"  {이름:>28} : {np.round(값, 6)}")
print()
print("한 번에 :", np.round((np.linalg.inv(M) @ Pl @ M) @ c, 6))

Ssym = np.array([[3.0, 1.0], [1.0, 3.0]])
값, 벡 = np.linalg.eigh(Ssym)
print("고윳값 :", 값)
print(show_matrix(벡, "고유벡터를 열에 세운 것"))
print()
print(show_matrix(행렬만들기(lambda v: Ssym @ v, np.eye(2), np.eye(2)),
                  "표준기저에서"))
print(show_matrix(깔끔(행렬만들기(lambda v: Ssym @ v, 벡, 벡)), "고유기저에서"))
print("대각인가 :", np.allclose(행렬만들기(lambda v: Ssym @ v, 벡, 벡),
                              np.diag(값)))

# 3x3 로도. L25 의 앵커를 다시 쓴다.
S3 = np.array([[6.0, 0.0, -3.0], [0.0, 6.0, -3.0], [-3.0, -3.0, 9.0]])
값3, 벡3 = np.linalg.eigh(S3)
print("고윳값 :", 값3)
print(show_matrix(깔끔(행렬만들기(lambda v: S3 @ v, 벡3, 벡3)),
                  "고유기저에서의 S3"))
print()
# 평면 x+y+z=0 으로의 투영도 같은 방법으로
b1, b2, n = (np.array([1.0, -1.0, 0.0]), np.array([1.0, 1.0, -2.0]),
             np.array([1.0, 1.0, 1.0]))
Pm = np.eye(3) - np.outer(n, n) / 3.0
좋은기저 = np.column_stack([b1, b2, n])
print(show_matrix(Pm * 3, "표준기저에서의 P  (3배)"))
print(show_matrix(깔끔(행렬만들기(lambda v: Pm @ v, 좋은기저, 좋은기저)),
                  "b1, b2, n 기저에서의 P"))
print("diag(1,1,0) 인가 :",
      np.allclose(행렬만들기(lambda v: Pm @ v, 좋은기저, 좋은기저),
                  np.diag([1.0, 1.0, 0.0])))
print("고윳값 :", np.sort(np.linalg.eigvalsh(Pm)), "  <- 살렸다, 살렸다, 죽였다")
print()
민것 = 행렬만들기(lambda v: Pm @ v, 좋은기저, 좋은기저)
print("깔끔() 을 씌우기 전의 (3,3) 성분 :", 민것[2, 2])
print("-> 0 이어야 하는데 1e-16 쯤 남는다. lstsq 가 실수 나눗셈을 하기 때문이다.")
print("   L28 에서 본 그대로다. 컴퓨터의 0 은 '충분히 작은 것' 이다.")

R = np.array([[0.0, -1.0], [1.0, 0.0]])       # 90도 회전
print("고윳값 :", np.linalg.eigvals(R), "  <- 실수가 아니다")
print()
print("실수 고유벡터가 정말 하나도 없는지 두들겨 보자.")
최대 = 0.0
for _ in range(200000):
    x = rng.normal(size=2)
    x /= np.linalg.norm(x)
    최대 = max(최대, abs(float(x @ (R @ x))))   # 같은 방향이면 |cos| = 1
print(f"단위벡터 20만개 중 |x . Rx| 의 최댓값 : {최대:.3e}")
print("-> 정확히 0 이다. 90도 회전은 모든 벡터를 자기 자신과 수직으로 보낸다.")
print("   자기 방향을 지키는 벡터가 하나도 없으니 대각이 될 수가 없다.")

print("무작위 기저 2000개로 회전을 적어 봐도 대각이 되는 일은 없다.")
가장대각 = np.inf
for _ in range(2000):
    Mr = rng.normal(size=(2, 2))
    if abs(np.linalg.det(Mr)) < 1e-6:
        continue
    표현 = np.linalg.inv(Mr) @ R @ Mr
    비대각 = abs(표현[0, 1]) + abs(표현[1, 0])
    가장대각 = min(가장대각, 비대각)
print(f"비대각 성분 크기의 최솟값 : {가장대각:.4f}   (0 이면 대각화 성공)")
print()
결함 = np.array([[3.0, 1.0], [0.0, 3.0]])
print("결함 행렬 [[3,1],[0,3]] 도 마찬가지다.")
print("  고윳값 :", np.linalg.eigvals(결함))
print("  독립인 고유벡터 개수 :",
      np.linalg.matrix_rank(np.linalg.eig(결함)[1], tol=1e-8), " (2 여야 대각화 가능)")

각도들 = np.arange(0, 91, 7.5)
프레임, 이름표 = [], []
for 각 in 각도들:
    t = np.radians(각)
    b1v = np.array([np.cos(t), np.sin(t)])
    b2v = np.array([-np.sin(t), np.cos(t)])
    Mb = np.column_stack([b1v, b2v])
    표현 = np.linalg.inv(Mb) @ Pl @ Mb
    s = np.linspace(-2.0, 2.0, 2)
    프레임.append([
        go.Scatter(x=s, y=s, mode="lines", name="y = x",
                   line=dict(color="#bbbbbb", width=2)),
        go.Scatter(x=[0, b1v[0]], y=[0, b1v[1]], mode="lines+markers",
                   name="b1", line=dict(color=COLORS["input"], width=5)),
        go.Scatter(x=[0, b2v[0]], y=[0, b2v[1]], mode="lines+markers",
                   name="b2", line=dict(color=COLORS["output"], width=5)),
        go.Scatter(x=[0], y=[-1.72], mode="text",
                   text=[f"[[{표현[0,0]:+.2f}, {표현[0,1]:+.2f}], "
                         f"[{표현[1,0]:+.2f}, {표현[1,1]:+.2f}]]"],
                   textfont=dict(size=18), showlegend=False),
    ])
    이름표.append(f"{각:.0f}")
배치 = layout2d("같은 투영, 기저만 돌린다", extent=2.05)
배치["height"] = 600
slider_figure(프레임, 이름표, 배치, prefix="기저 각도 ")

print(f"{'각도':>8}{'대각합':>10}{'행렬식':>10}{'대칭인가':>10}")
for 각 in (0, 15, 30, 45, 60, 90):
    t = np.radians(각)
    Mb = np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])
    표현 = np.linalg.inv(Mb) @ Pl @ Mb
    print(f"{각:>8}{np.trace(표현):>10.4f}{np.linalg.det(표현):>10.4f}"
          f"{str(np.allclose(표현, 표현.T)):>10}")
print()
print("직교 기저만 쓰면 대칭성도 지켜진다. 5절에서 깨졌던 것은 비스듬한 기저였다.")