Lecture 33. 좌·우 역행렬과 의사역행렬 — 파이썬 실습 - 선형대수

L33 서술 파트의 반전은 하나였다. L16의 최소제곱 공식이 사실은 좌역행렬이었다.

이 노트북에서는 그 재회를 숫자로 확인한다. L16의 앵커를 그대로 가져와 $(A^{\mathsf T}A)^{-1}A^{\mathsf T}$ 와 pinv 가 같은 행렬인지 보고, 해집합을 따라 움직이며 노름이 정말 늘기만 하는지 재고, L16이 미뤄 둔 질문에 답한다.

서술 파트의 내용	여기서 확인하는 방법
네 경우	$A^{+}A$ , $AA^{+}$ 가 $I$ 인지 $P$ 인지
좌역행렬 = 최소제곱	L16 앵커로 대조
$AA^{+}$ = L15의 투영행렬	$P^2=P$ , $P^{\mathsf T}=P$ , 고윳값
우역행렬 = 최소노름	해집합 위를 슬라이더로
$A^{+} = V\Sigma^{+}U^{\mathsf T}$	직접 만들어 numpy 와 대조
$A^{+}A$ , $AA^{+}$ 는 투영	랭크부족 행렬로
$\kappa(A^{\mathsf T}A) = \kappa(A)^2$	정답을 아는 문제로
정규방정식 < QR < SVD	세 방법 오차 겨루기
rcond, ridge	딱 자르기 vs 부드러운 자르기
무어-펜로즈 네 조건	유일성
L30의 적분행렬	$S = D^{+}$

0. 준비¶

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, show_matrix, slider_figure

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

numpy 2.5.2

1. 네 경우¶

$A^{+}A$ 와 $AA^{+}$ 가 각각 $I$ 가 되는지 $P$ 가 되는지가 전부를 가른다.

def 되돌리기(A, 이름=""):
    """A^+ A 와 A A^+ 가 무엇이 되는지 본다."""
    A = np.asarray(A, dtype=float)
    m, n = A.shape
    r = np.linalg.matrix_rank(A)
    P = np.linalg.pinv(A)
    좌, 우 = P @ A, A @ P
    좌I = bool(np.allclose(좌, np.eye(n)))
    우I = bool(np.allclose(우, np.eye(m)))
    이름표 = ("진짜 역행렬" if (좌I and 우I) else
              "좌역행렬만" if 좌I else
              "우역행렬만" if 우I else "의사역행렬만")
    return dict(이름=이름, 크기=(m, n), 랭크=r, 좌I=좌I, 우I=우I,
                종류=이름표, 좌=좌, 우=우)

경우 = (
    ("r = m = n", np.array([[2.0, 1.0], [1.0, 3.0]])),
    ("r = n < m", np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])),
    ("r = m < n", np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 3.0]])),
    ("r < m, n ", np.array([[1.0, 2.0], [2.0, 4.0], [3.0, 6.0]])),
)
print(f"{'':>12}{'크기':>9}{'r':>3}{'A^+A = I':>11}{'AA^+ = I':>11}{'':>4}{'종류':>14}")
for 이름, A in 경우:
    보 = 되돌리기(A, 이름)
    print(f"{이름:>12}{str(보['크기']):>9}{보['랭크']:>3}"
          f"{str(보['좌I']):>11}{str(보['우I']):>11}{'':>4}{보['종류']:>14}")
print()
print("-> 서술 파트의 네 칸이 그대로 나온다.")

                   크기  r   A^+A = I   AA^+ = I                종류
   r = m = n   (2, 2)  2       True       True            진짜 역행렬
   r = n < m   (3, 2)  2       True      False             좌역행렬만
   r = m < n   (2, 3)  2      False       True             우역행렬만
   r < m, n    (3, 2)  1      False      False            의사역행렬만

-> 서술 파트의 네 칸이 그대로 나온다.

# I 가 아닌 것들은 정말 투영행렬인가
print("I 가 아닌 자리에 앉은 행렬이 투영인지 확인")
for 이름, A in 경우:
    보 = 되돌리기(A, 이름)
    for 쪽, M, 맞음 in (("A^+A", 보["좌"], 보["좌I"]), ("AA^+", 보["우"], 보["우I"])):
        if 맞음:
            continue
        print(f"  {이름} 의 {쪽} : P^2=P {np.allclose(M@M, M)}, "
              f"P^T=P {np.allclose(M.T, M)}, 랭크 {np.linalg.matrix_rank(M)}, "
              f"고윳값 {np.round(np.linalg.eigvalsh(M), 6)}")

I 가 아닌 자리에 앉은 행렬이 투영인지 확인
  r = n < m 의 AA^+ : P^2=P True, P^T=P True, 랭크 2, 고윳값 [-0.  1.  1.]
  r = m < n 의 A^+A : P^2=P True, P^T=P True, 랭크 2, 고윳값 [-0.  1.  1.]
  r < m, n  의 A^+A : P^2=P True, P^T=P True, 랭크 1, 고윳값 [-0.  1.]
  r < m, n  의 AA^+ : P^2=P True, P^T=P True, 랭크 1, 고윳값 [-0.  0.  1.]

2. 좌역행렬은 L16의 최소제곱이었다¶

L16의 앵커를 그대로 가져온다.

A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
b = np.array([1.0, 3.0, 2.0])
print(show_matrix(A, "A  (L16 의 그 행렬)"))
print("b =", b)
print()
왼 = np.linalg.inv(A.T @ A) @ A.T
print(show_matrix(왼, "(A^T A)^-1 A^T   손으로 만든 좌역행렬"))
print(show_matrix(np.linalg.pinv(A), "np.linalg.pinv(A)"))
print("같은 행렬인가 :", np.allclose(왼, np.linalg.pinv(A)))

A  (L16 의 그 행렬)
[  1   1 ]
[  1   2 ]
[  1   3 ]
b = [1. 3. 2.]

(A^T A)^-1 A^T   손으로 만든 좌역행렬
[    1.33    0.333   -0.667 ]
[    -0.5        0      0.5 ]
np.linalg.pinv(A)
[       1.33       0.333      -0.667 ]
[       -0.5   -3.02e-17         0.5 ]
같은 행렬인가 : True

x햇 = 왼 @ b
p = A @ x햇
e = b - p
print("x_hat =", x햇, "   L16 의 (1, 0.5) 인가 :", np.allclose(x햇, [1.0, 0.5]))
print("p     =", p, "   L16 의 (1.5, 2, 2.5) 인가 :", np.allclose(p, [1.5, 2.0, 2.5]))
print("e     =", e, "   |e|^2 =", float(e @ e))
print()
print("잔차가 열공간에 수직인가 : A^T e =", np.round(A.T @ e, 12))
print("lstsq 와 같은가 :", np.allclose(x햇, np.linalg.lstsq(A, b, rcond=None)[0]))

x_hat = [1.  0.5]    L16 의 (1, 0.5) 인가 : True
p     = [1.5 2.  2.5]    L16 의 (1.5, 2, 2.5) 인가 : True
e     = [-0.5  1.  -0.5]    |e|^2 = 1.5

잔차가 열공간에 수직인가 : A^T e = [0. 0.]
lstsq 와 같은가 : True

$AA^{+}$ 는 L15의 투영행렬이다¶

P = A @ 왼
print(show_matrix(P, "A A^+"))
print("L15 의 A(A^T A)^-1 A^T 와 같은가 :",
      np.allclose(P, A @ np.linalg.inv(A.T @ A) @ A.T))
print()
print("P^2 = P    :", np.allclose(P @ P, P))
print("P^T = P    :", np.allclose(P.T, P))
print("랭크       :", np.linalg.matrix_rank(P), " (열공간의 차원)")
print("고윳값     :", np.round(np.linalg.eigvalsh(P), 12))
print()
print("A^+ A =", np.round(왼 @ A, 12).tolist(), " <- 이쪽은 I 다")

A A^+
[   0.833    0.333   -0.167 ]
[   0.333    0.333    0.333 ]
[  -0.167    0.333    0.833 ]
L15 의 A(A^T A)^-1 A^T 와 같은가 : True

P^2 = P    : True
P^T = P    : True
랭크       : 2  (열공간의 차원)
고윳값     : [0. 1. 1.]

A^+ A = [[1.0, -0.0], [0.0, 1.0]]  <- 이쪽은 I 다

3. 우역행렬은 가장 짧은 해를 고른다¶

A2 = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 3.0]])
b2 = np.array([6.0, 14.0])
오른 = A2.T @ np.linalg.inv(A2 @ A2.T)
print("A^T (A A^T)^-1 == pinv 인가 :", np.allclose(오른, np.linalg.pinv(A2)))
x별 = 오른 @ b2
print()
print("최소노름 해 x+ =", x별, "  노름", f"{np.linalg.norm(x별):.6f}")
print("정말 해인가 : A x+ =", A2 @ x별, " b =", b2,
      " 같은가", np.allclose(A2 @ x별, b2))
print()
영 = np.array([1.0, -2.0, 1.0])
print("영공간 방향 n =", 영, "   A n =", A2 @ 영)
print("x+ 가 영공간에 수직인가 : x+ . n =", float(x별 @ 영))

A^T (A A^T)^-1 == pinv 인가 : True

최소노름 해 x+ = [1. 2. 3.]   노름 3.741657
정말 해인가 : A x+ = [ 6. 14.]  b = [ 6. 14.]  같은가 True

영공간 방향 n = [ 1. -2.  1.]    A n = [0. 0.]
x+ 가 영공간에 수직인가 : x+ . n = 0.0

print("해집합을 따라 움직여 본다  (x = x+ + t n)")
print(f"{'t':>7}{'해':>24}{'A x':>16}{'노름':>11}")
for t in (-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0):
    x = x별 + t * 영
    표시 = "  <- 최소" if t == 0.0 else ""
    print(f"{t:>7.1f}{str(np.round(x, 3)):>24}{str(A2 @ x):>16}"
          f"{np.linalg.norm(x):>11.6f}{표시}")
print()
print("피타고라스로 확인 : |x+ + tn|^2 = |x+|^2 + t^2|n|^2 인가")
for t in (0.7, -1.3, 2.5):
    왼쪽 = np.linalg.norm(x별 + t*영)**2
    오른쪽 = np.linalg.norm(x별)**2 + t**2 * np.linalg.norm(영)**2
    print(f"  t={t:>5} : {왼쪽:.9f} vs {오른쪽:.9f}  "
          f"같은가 {np.isclose(왼쪽, 오른쪽)}")

해집합을 따라 움직여 본다  (x = x+ + t n)
      t                       해             A x         노름
   -2.0           [-1.  6.  1.]       [ 6. 14.]   6.164414
   -1.0           [-0.  4.  2.]       [ 6. 14.]   4.472136
   -0.5           [0.5 3.  2.5]       [ 6. 14.]   3.937004
    0.0              [1. 2. 3.]       [ 6. 14.]   3.741657  <- 최소
    0.5           [1.5 1.  3.5]       [ 6. 14.]   3.937004
    1.0           [ 2. -0.  4.]       [ 6. 14.]   4.472136
    2.0           [ 3. -2.  5.]       [ 6. 14.]   6.164414

피타고라스로 확인 : |x+ + tn|^2 = |x+|^2 + t^2|n|^2 인가
  t=  0.7 : 16.940000000 vs 16.940000000  같은가 True
  t= -1.3 : 24.140000000 vs 24.140000000  같은가 True
  t=  2.5 : 51.500000000 vs 51.500000000  같은가 True

ts = np.linspace(-3.0, 3.0, 41)
프레임, 이름표 = [], []
for t in ts[::4]:
    x = x별 + t * 영
    선 = x별[:, None] + 영[:, None] * np.linspace(-3.2, 3.2, 40)
    프레임.append([
        go.Scatter3d(x=선[0], y=선[1], z=선[2], mode="lines",
                     line=dict(color=COLORS["output"], width=6),
                     name="해집합"),
        go.Scatter3d(x=[0, x[0]], y=[0, x[1]], z=[0, x[2]], mode="lines+markers",
                     line=dict(color="#d62728", width=7),
                     marker=dict(size=5), name=f"|x| = {np.linalg.norm(x):.3f}"),
        go.Scatter3d(x=[0], y=[0], z=[0], mode="markers",
                     marker=dict(size=5, color="black"), name="원점"),
    ])
    이름표.append(f"{t:+.1f}")
배치 = dict(title=dict(text="해집합 위를 움직이면 원점까지의 거리가 어떻게 되는가"),
           scene=dict(xaxis=dict(range=[-4, 5], title="x1"),
                      yaxis=dict(range=[-4, 7], title="x2"),
                      zaxis=dict(range=[-1, 7], title="x3"),
                      aspectmode="cube"),
           height=620, margin=dict(l=10, r=10, t=60, b=10))
slider_figure(프레임, 이름표, 배치, prefix="t = ", initial=len(이름표)//2)

4. 일반형 — 정의대로 만들어 본다¶

def 손으로pinv(A, 눈감아=None):
    """A^+ = V Sigma^+ U^T 를 정의대로 만든다."""
    A = np.asarray(A, dtype=float)
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    if 눈감아 is None:
        눈감아 = max(A.shape) * np.finfo(float).eps * (s[0] if s.size else 0.0)
    s플러스 = np.where(s > 눈감아, 1.0 / np.where(s > 눈감아, s, 1.0), 0.0)
    return Vt.T @ np.diag(s플러스) @ U.T

시험 = (("과결정 3x2", A), ("부족결정 2x3", A2),
        ("랭크부족 3x2", np.array([[1.0, 2.0], [2.0, 4.0], [3.0, 6.0]])),
        ("정방 가역", np.array([[2.0, 1.0], [1.0, 3.0]])),
        ("영행렬 2x3", np.zeros((2, 3))),
        ("무작위 5x7", rng.normal(size=(5, 7))))
print(f"{'':>16}{'numpy 와 일치':>14}{'A^+A 랭크':>11}{'AA^+ 랭크':>11}{'r':>4}")
for 이름, X in 시험:
    P = 손으로pinv(X)
    print(f"{이름:>16}{str(np.allclose(P, np.linalg.pinv(X))):>14}"
          f"{np.linalg.matrix_rank(P @ X):>11}{np.linalg.matrix_rank(X @ P):>11}"
          f"{np.linalg.matrix_rank(X):>4}")
print()
print("-> A^+A 와 AA^+ 의 랭크가 언제나 r 이다. 둘 다 r 차원으로의 투영이다.")

                    numpy 와 일치    A^+A 랭크    AA^+ 랭크   r
         과결정 3x2          True          2          2   2
        부족결정 2x3          True          2          2   2
        랭크부족 3x2          True          1          1   1
           정방 가역          True          2          2   2
         영행렬 2x3          True          0          0   0
         무작위 5x7          True          5          5   5

-> A^+A 와 AA^+ 의 랭크가 언제나 r 이다. 둘 다 r 차원으로의 투영이다.

랭크가 모자라면 양쪽 다 $I$ 가 아니다¶

X = np.array([[1.0, 2.0], [2.0, 4.0], [3.0, 6.0]])      # 랭크 1
P = np.linalg.pinv(X)
print(show_matrix(X, "X  (3x2, 랭크 1)"))
print(show_matrix(P @ X, "X^+ X   <- 행공간으로의 투영"))
print(show_matrix(X @ P, "X X^+   <- 열공간으로의 투영"))
print()
행 = np.array([1.0, 2.0]) / np.sqrt(5)
영방향 = np.array([2.0, -1.0]) / np.sqrt(5)
print("행공간 벡터에 X^+X 를 곱하면 :", np.round((P @ X) @ 행, 6), " (그대로)")
print("영공간 벡터에 X^+X 를 곱하면 :", np.round((P @ X) @ 영방향, 12), " (0 으로)")
print()
print("-> X 가 영공간에서 죽인 것은 X^+ 가 되살리지 못한다.")

X  (3x2, 랭크 1)
[  1   2 ]
[  2   4 ]
[  3   6 ]
X^+ X   <- 행공간으로의 투영
[  0.2   0.4 ]
[  0.4   0.8 ]
X X^+   <- 열공간으로의 투영
[  0.0714    0.143    0.214 ]
[   0.143    0.286    0.429 ]
[   0.214    0.429    0.643 ]

행공간 벡터에 X^+X 를 곱하면 : [0.4472 0.8944]  (그대로)
영공간 벡터에 X^+X 를 곱하면 : [0. 0.]  (0 으로)

-> X 가 영공간에서 죽인 것은 X^+ 가 되살리지 못한다.

5. L16이 미뤄 둔 질문의 답¶

“우리는 $A^{\mathsf T}A$ 를 뒤집었다. 이 계산은 안전한가.”

print("조건수가 정말 제곱되는가")
print(f"{'':>28}{'cond(A)':>12}{'cond(A^T A)':>14}{'cond(A)^2':>13}{'비':>8}")
후보 = (("L16 의 A", A),
        ("[[1,1],[1,1.01],[1,1]]", np.array([[1.,1.],[1.,1.01],[1.,1.]])),
        ("[[1,1],[1,1+1e-5],[1,1]]", np.array([[1.,1.],[1.,1+1e-5],[1.,1.]])),
        ("힐베르트 6x4",
         np.array([[1/(i+j+1) for j in range(4)] for i in range(6)])),
        ("무작위 20x5", rng.normal(size=(20, 5))))
for 이름, X in 후보:
    k = np.linalg.cond(X)
    k2 = np.linalg.cond(X.T @ X)
    print(f"{이름:>28}{k:>12.4e}{k2:>14.4e}{k**2:>13.4e}{k2/k**2:>8.4f}")
print()
print("-> 마지막 칸이 전부 1.0000 근처다. 조건수가 정확히 제곱된다.")
print("   (조건수가 아주 크면 cond(A^T A) 자체를 재는 일이 이미 부정확해")
print("    소수점 아래가 조금 흔들린다.)")

조건수가 정말 제곱되는가
                                 cond(A)   cond(A^T A)    cond(A)^2       비
                     L16 의 A  6.7930e+00    4.6145e+01   4.6145e+01  1.0000
      [[1,1],[1,1.01],[1,1]]  4.2568e+02    1.8121e+05   1.8121e+05  1.0000
    [[1,1],[1,1+1e-5],[1,1]]  4.2427e+05    1.8000e+11   1.8000e+11  1.0000
                    힐베르트 6x4  6.4854e+03    4.2060e+07   4.2060e+07  1.0000
                    무작위 20x5  2.1857e+00    4.7775e+00   4.7775e+00  1.0000

-> 마지막 칸이 전부 1.0000 근처다. 조건수가 정확히 제곱된다.
   (조건수가 아주 크면 cond(A^T A) 자체를 재는 일이 이미 부정확해
    소수점 아래가 조금 흔들린다.)

def 정답아는문제(kappa, m=40, n=6, seed=7):
    """최소제곱의 정답이 x_true 인 문제를 만든다.

    b 에 좌영공간 성분만 섞으면 최소제곱해가 x_true 그대로다.
    """
    r = np.random.default_rng(seed)
    U, _ = np.linalg.qr(r.normal(size=(m, m)))
    V, _ = np.linalg.qr(r.normal(size=(n, n)))
    s = np.logspace(0, -np.log10(kappa), n)
    X = U[:, :n] @ np.diag(s) @ V.T
    x참 = r.normal(size=n)
    잡 = U[:, n:] @ r.normal(size=m - n)
    return X, X @ x참 + 잡 / np.linalg.norm(잡) * 0.3, x참

print("같은 문제를 세 가지로 풀어 정답과 견준다")
print("(뽑기 하나는 흔들리므로 씨앗 21개의 중앙값을 쓴다)")
print(f"{'cond(A)':>10}{'정규방정식':>14}{'QR (lstsq)':>14}{'SVD (pinv)':>14}{'배수':>8}")
for k in (1e2, 1e4, 1e6, 1e7, 1e8, 1e9):
    정규들, qr들, sv들 = [], [], []
    for 씨 in range(21):
        X, y, x참 = 정답아는문제(k, seed=씨)
        기준 = np.linalg.norm(x참)
        try:
            정규들.append(np.linalg.norm(np.linalg.solve(X.T@X, X.T@y) - x참) / 기준)
        except np.linalg.LinAlgError:
            정규들.append(np.inf)
        qr들.append(np.linalg.norm(np.linalg.lstsq(X, y, rcond=None)[0] - x참) / 기준)
        sv들.append(np.linalg.norm(np.linalg.pinv(X) @ y - x참) / 기준)
    정규, qr, sv = np.median(정규들), np.median(qr들), np.median(sv들)
    print(f"{k:>10.0e}{정규:>14.3e}{qr:>14.3e}{sv:>14.3e}{정규/qr:>8.1f}")
print()
print("-> 정규방정식이 늘 위에 있고, 조건수가 커질수록 격차가 벌어진다.")
print("   cond(A) 가 1e8 을 넘으면 cond(A^T A) 가 1e16 = 1/eps 를 넘는다.")
print("   그 지점부터는 자릿수가 남지 않는다. inv(A.T@A)@A.T 를 손으로 쓰지 마라.")

같은 문제를 세 가지로 풀어 정답과 견준다
(뽑기 하나는 흔들리므로 씨앗 21개의 중앙값을 쓴다)
   cond(A)         정규방정식    QR (lstsq)    SVD (pinv)      배수
     1e+02     1.481e-13     6.841e-15     9.026e-15    21.7
     1e+04     1.053e-09     3.162e-11     3.168e-11    33.3
     1e+06     1.573e-05     3.703e-07     3.703e-07    42.5
     1e+07     2.269e-03     2.968e-05     2.968e-05    76.5
     1e+08     1.245e-01     6.314e-03     6.314e-03    19.7
     1e+09     2.553e+00     4.763e-01     4.763e-01     5.4

-> 정규방정식이 늘 위에 있고, 조건수가 커질수록 격차가 벌어진다.
   cond(A) 가 1e8 을 넘으면 cond(A^T A) 가 1e16 = 1/eps 를 넘는다.
   그 지점부터는 자릿수가 남지 않는다. inv(A.T@A)@A.T 를 손으로 쓰지 마라.

6. 작은 특이값을 어떻게 할 것인가¶

U6, _ = np.linalg.qr(rng.normal(size=(6, 6)))
V6, _ = np.linalg.qr(rng.normal(size=(4, 4)))
sig = np.array([1.0, 0.5, 1e-3, 1e-6])
Ar = U6[:, :4] @ np.diag(sig) @ V6.T
br = rng.normal(size=6)
print("특이값 :", " ".join(f"{v:.0e}" for v in sig))
print()
print("[rcond] 딱 자르기")
print(f"{'rcond':>10}{'살아남은 sigma':>16}{'해의 노름':>14}{'잔차':>12}")
for rc in (1e-8, 1e-5, 1e-4, 1e-2, 1e-1):
    x = np.linalg.pinv(Ar, rcond=rc) @ br
    print(f"{rc:>10.0e}{int((sig > rc*sig[0]).sum()):>16}"
          f"{np.linalg.norm(x):>14.4e}{np.linalg.norm(Ar@x-br):>12.6f}")
print()
print("-> 문턱을 낮추면 잔차는 줄지만 해가 폭발한다. 사람이 정해야 하는 판단이다.")

특이값 : 1e+00 5e-01 1e-03 1e-06

[rcond] 딱 자르기
     rcond      살아남은 sigma         해의 노름          잔차
     1e-08               4    2.6065e+05    0.570913
     1e-05               3    5.9623e+01    0.627599
     1e-04               3    5.9623e+01    0.627599
     1e-02               2    3.4586e+00    0.630415
     1e-01               2    3.4586e+00    0.630415

-> 문턱을 낮추면 잔차는 줄지만 해가 폭발한다. 사람이 정해야 하는 판단이다.

print("[ridge] 부드러운 자르기 — 필터 인자가 sigma/(sigma^2+lam) 인가")
lam = 1e-6
직접 = np.linalg.solve(Ar.T @ Ar + lam*np.eye(4), Ar.T @ br)
필터 = V6 @ np.diag(sig/(sig**2 + lam)) @ U6[:, :4].T @ br
print("  정규방정식으로 푼 것과 필터 공식이 같은가 :", np.allclose(직접, 필터))
print()
print(f"  lambda = {lam:.0e} 이면 문턱 sqrt(lambda) = {np.sqrt(lam):.0e}")
print(f"{'sigma':>10}{'1/sigma':>14}{'sigma/(s^2+lam)':>18}{'비율':>9}")
for s in sig:
    print(f"{s:>10.0e}{1/s:>14.4e}{s/(s**2+lam):>18.4e}{(s/(s**2+lam))*s:>9.4f}")
print()
print("-> sigma 가 문턱보다 크면 1/sigma 그대로(비율 1), 작으면 0 쪽으로 눌린다.")

[ridge] 부드러운 자르기 — 필터 인자가 sigma/(sigma^2+lam) 인가
  정규방정식으로 푼 것과 필터 공식이 같은가 : True

  lambda = 1e-06 이면 문턱 sqrt(lambda) = 1e-03
     sigma       1/sigma   sigma/(s^2+lam)       비율
     1e+00    1.0000e+00        1.0000e+00   1.0000
     5e-01    2.0000e+00        2.0000e+00   1.0000
     1e-03    1.0000e+03        5.0000e+02   0.5000
     1e-06    1.0000e+06        1.0000e+00   0.0000

-> sigma 가 문턱보다 크면 1/sigma 그대로(비율 1), 작으면 0 쪽으로 눌린다.

7. 무어-펜로즈 네 조건¶

def 네조건(A, X):
    """AXA=A, XAX=X, (AX)^T=AX, (XA)^T=XA."""
    return (bool(np.allclose(A @ X @ A, A)),
            bool(np.allclose(X @ A @ X, X)),
            bool(np.allclose((A @ X).T, A @ X)),
            bool(np.allclose((X @ A).T, X @ A)))

print(f"{'':>16}{'pinv 가 네 조건을 만족하는가':>30}")
for 이름, X in 시험:
    if X.size == 0:
        continue
    print(f"{이름:>16}{str(네조건(X, np.linalg.pinv(X))):>30}")
print()
가짜 = np.linalg.pinv(A) + 0.01 * rng.normal(size=(2, 3))
print("살짝 흔든 것 :", 네조건(A, 가짜), "  <- 넷이 한꺼번에 무너진다")
print()
print("가역 행렬이면 A^+ = A^-1 인가")
G = np.array([[2.0, 1.0], [1.0, 3.0]])
print("  pinv(G) == inv(G) :", np.allclose(np.linalg.pinv(G), np.linalg.inv(G)))
print("  inv(G) 가 네 조건 :", 네조건(G, np.linalg.inv(G)))

                            pinv 가 네 조건을 만족하는가
         과결정 3x2      (True, True, True, True)
        부족결정 2x3      (True, True, True, True)
        랭크부족 3x2      (True, True, True, True)
           정방 가역      (True, True, True, True)
         영행렬 2x3      (True, True, True, True)
         무작위 5x7      (True, True, True, True)

살짝 흔든 것 : (False, False, False, False)   <- 넷이 한꺼번에 무너진다

가역 행렬이면 A^+ = A^-1 인가
  pinv(G) == inv(G) : True
  inv(G) 가 네 조건 : (True, True, True, True)

8. L30에서 이름을 안 붙이고 넘어간 그것¶

def 미분행렬(n):
    """1, x, ..., x^(n-1) 기저에서의 미분. (n-1) x n."""
    D = np.zeros((n - 1, n))
    for j in range(1, n):
        D[j - 1, j] = j
    return D


def 적분행렬(n):
    """상수항을 0 으로 두는 적분. n x (n-1)."""
    S = np.zeros((n, n - 1))
    for j in range(n - 1):
        S[j + 1, j] = 1.0 / (j + 1)
    return S

for n in (4, 6, 8):
    D, S = 미분행렬(n), 적분행렬(n)
    print(f"n={n} : S == pinv(D) 인가 {np.allclose(S, np.linalg.pinv(D))},"
          f"  네 조건 {네조건(D, S)}")
print()
D, S = 미분행렬(4), 적분행렬(4)
print(show_matrix(D, "D  (미분)"))
print(show_matrix(S, "S  (적분)"))
print(show_matrix(D @ S, "D S   <- I_3 다. 미분하고 적분하면 돌아온다"))
print(show_matrix(S @ D, "S D   <- I 가 아니다. 첫 칸이 죽었다"))
print()
print("첫 칸은 상수항이다. 미분이 죽였으니 적분이 되살릴 수 없다.")
print("S D 가 무엇인가 : P^2 = P 인가", np.allclose((S@D)@(S@D), S@D),
      ", 랭크", np.linalg.matrix_rank(S@D), "  -> 행공간으로의 투영")

n=4 : S == pinv(D) 인가 True,  네 조건 (True, True, True, True)
n=6 : S == pinv(D) 인가 True,  네 조건 (True, True, True, True)
n=8 : S == pinv(D) 인가 True,  네 조건 (True, True, True, True)

D  (미분)
[  0   1   0   0 ]
[  0   0   2   0 ]
[  0   0   0   3 ]
S  (적분)
[      0       0       0 ]
[      1       0       0 ]
[      0     0.5       0 ]
[      0       0   0.333 ]
D S   <- I_3 다. 미분하고 적분하면 돌아온다
[  1   0   0 ]
[  0   1   0 ]
[  0   0   1 ]
S D   <- I 가 아니다. 첫 칸이 죽었다
[  0   0   0   0 ]
[  0   1   0   0 ]
[  0   0   1   0 ]
[  0   0   0   1 ]

첫 칸은 상수항이다. 미분이 죽였으니 적분이 되살릴 수 없다.
S D 가 무엇인가 : P^2 = P 인가 True , 랭크 3   -> 행공간으로의 투영

마치며...¶

서술 파트의 내용	이 노트북의 코드
네 경우	$A^{+}A$ , $AA^{+}$ 로 정확히 갈린다
좌역행렬 = 최소제곱	손으로 만든 것과 `pinv` 가 같은 행렬
L16 앵커	$\hat{x}=(1,0.5)$ , $p$ , $e$ 가 그대로
$AA^{+}$ = L15의 $P$	$P^2=P$ , $P^{\mathsf T}=P$ , 고윳값 $0,1,1$
우역행렬 = 최소노름	해집합을 따라 움직이면 늘기만
피타고라스	$\lVert x^{+}+tn\rVert^2 = \lVert x^{+}\rVert^2 + t^2\lVert n\rVert^2$
$V\Sigma^{+}U^{\mathsf T}$	여섯 가지 모양에서 numpy 와 일치
둘 다 랭크 $r$ 투영	랭크부족 행렬로 확인
$\kappa(A^{\mathsf T}A)=\kappa(A)^2$	상대오차 10^-6 이내로 일치
정규방정식이 무너지는 곳	$\kappa(A)=10^{8}$
rcond / ridge	딱 자르기 vs $\sqrt\lambda$ 문턱
무어-펜로즈	살짝 흔들면 넷이 한꺼번에 무너진다
L30의 $S$	$S = D^{+}$ , $SD$ 는 행공간 투영

더 해 볼 것¶

3절의 앵커에서 $\vv{b}$ 를 바꿔 보자. 최소노름 해가 여전히 행공간 안에 있는가? $\vv{b}$ 를 어떻게 잡아도 그런가?
5절의 실험에서 $m$ 과 $n$ 을 바꿔 보자. 정규방정식이 무너지는 조건수가 달라지는가? 왜 달라지지 않는가?
pinv 의 rcond 를 아주 크게(예: 0.5) 주면 어떻게 되는가? 랭크가 몇으로 보이는가? 그때 해는 무엇을 뜻하는가?
8절의 미분행렬에서 적분 상수를 0이 아닌 다른 값으로 두는 행렬 $S'$ 을 만들어 보자. $DS' = I$ 는 여전히 성립하는가? 네 조건 중 어느 것이 깨지는가?
무어-펜로즈 조건 중 ③과 ④만 빼고 만족하는 행렬을 찾아보자. 그런 것을 "일반화 역행렬"이라 부르는데, 왜 유일하지 않은가?

다음 강의는 마지막이다. 서른세 편이 어떻게 다섯 개의 분해로 정리되는지 살펴보자.

Lecture 33. 좌·우 역행렬과 의사역행렬 — 파이썬 실습