선형대수 표기와 관례¶

(1) 해석적으로. 네 성질이 모두 \((\mathbf{X}^\top\mathbf{X})^{-1}\)의 존재 하나에서 나온다.

대칭성. \(\mathbf{X}^\top\mathbf{X}\)가 대칭이므로 그 역행렬도 대칭이고, \((\mathbf{A}\mathbf{B})^\top = \mathbf{B}^\top\mathbf{A}^\top\)를 두 번 쓰면

이다.

멱등성. 가운데에서 \((\mathbf{X}^\top\mathbf{X})^{-1}(\mathbf{X}^\top\mathbf{X}) = \mathbf{I}_p\)가 지워진다.

대각합. 순환 성질 \(\operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{B}\mathbf{A})\)를 \(\mathbf{A} = \mathbf{X}\), \(\mathbf{B} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\)로 쓰면

이다. \(n \times n\) 행렬의 대각합이 \(n\)에 무관하게 \(p\)로 정해지는 것이 요점이다.

고윳값. \(\mathbf{H}\mathbf{v} = \lambda\mathbf{v}\)이고 \(\mathbf{v} \ne \mathbf{0}\)이면 멱등성에서 \(\lambda\mathbf{v} = \mathbf{H}\mathbf{v} = \mathbf{H}^2\mathbf{v} = \lambda^2\mathbf{v}\)이므로 \(\lambda^2 = \lambda\), 곧 \(\lambda \in \{0, 1\}\)이다. \(\mathbf{H}\)가 대칭이므로 스펙트럼 정리가 실수 고윳값 \(n\)개를 보장하고, 대각합이 고윳값의 합이라는 사실에서 \(1\)의 개수가 \(\operatorname{tr}(\mathbf{H}) = p\)개, 나머지 \(n - p\)개가 \(0\)이다. 이 \(n - p\)가 잔차 자유도 \(d\)다.

(2) 수치적으로. \(\sigma = 0.5\)인 잡음을 얹은 자료로 네 성질을 확인한다.

import numpy as np

rng = np.random.default_rng(42)
n, p = 50, 3
X = np.column_stack([np.ones(n), rng.standard_normal((n, p - 1))])
beta_true = np.array([2.0, -1.0, 0.5])
sigma = 0.5
y = X @ beta_true + rng.standard_normal(n) * sigma

# === 정규방정식으로 최소제곱 풀기 ===
beta_hat = np.linalg.solve(X.T @ X, X.T @ y)
print("참값   beta:", beta_true)
print("최소제곱 beta:", beta_hat.round(4))

# 추정값이 참값에서 얼마나 벗어났는지를 표준오차 단위로 재어 본다.
se = sigma * np.sqrt(np.diag(np.linalg.inv(X.T @ X)))
print("표준오차      :", se.round(4))
print("(추정-참)/SE  :", ((beta_hat - beta_true) / se).round(2))

# === 모자 행렬의 성질 ===
H = X @ np.linalg.inv(X.T @ X) @ X.T
print(f"대칭:   {np.allclose(H, H.T)}")
print(f"멱등:   {np.allclose(H @ H, H)}")
print(f"tr(H) = {np.trace(H):.1f}  (p = {p} 와 같아야 한다)")

# 고윳값은 1 이 p 개, 0 이 n-p 개여야 한다.
eig_H = np.linalg.eigvalsh(H)
print(f"H 의 고윳값 중 1 에 가까운 것: {np.sum(np.isclose(eig_H, 1))} 개")
print(f"H 의 고윳값 중 0 에 가까운 것: {np.sum(np.isclose(eig_H, 0))} 개")

# === X'X 의 스펙트럼 분해 (대칭이고 여기서는 양정치다) ===
eigvals, eigvecs = np.linalg.eigh(X.T @ X)
print(f"고윳값: {eigvals.round(2)}")
print(f"고윳값의 합 {eigvals.sum():.2f} = tr(X'X) {np.trace(X.T @ X):.2f}")
print(f"조건수: {eigvals.max() / eigvals.min():.2f}")
reconstruction = eigvecs @ np.diag(eigvals) @ eigvecs.T
print(f"스펙트럼 복원이 X'X 와 일치: {np.allclose(reconstruction, X.T @ X)}")

출력:

참값   beta: [ 2.  -1.   0.5]
최소제곱 beta: [ 1.9711 -0.95    0.4159]
표준오차      : [0.071  0.0882 0.1041]
(추정-참)/SE  : [-0.41  0.57 -0.81]
대칭:   True
멱등:   True
tr(H) = 3.0  (p = 3 와 같아야 한다)
H 의 고윳값 중 1 에 가까운 것: 3 개
H 의 고윳값 중 0 에 가까운 것: 47 개
고윳값: [20.34 38.64 51.  ]
고윳값의 합 109.97 = tr(X'X) 109.97
조건수: 2.51
스펙트럼 복원이 X'X 와 일치: True

네 성질이 모두 맞는다. \(\operatorname{tr}(\mathbf{H}) = 3.0\)이 \(p = 3\)과 같고, \(50 \times 50\) 행렬인 \(\mathbf{H}\)의 고윳값은 \(1\)이 \(3\)개, \(0\)이 \(47 = n - p\)개다. 유도한 대로 \(\mathbf{H}\)는 \(\mathbb{R}^{50}\)을 3차원 열공간과 47차원 잔차공간으로 쪼개며, \(n\)이 아무리 커져도 열공간의 차원은 \(p\)에 묶여 있다.

추정값 \((1.9711, -0.9500, 0.4159)\)는 참값 \((2, -1, 0.5)\)과 다르다. 유도가 틀린 것이 아니라 자료에 잡음이 있어서다. \(\operatorname{Cov}(\hat{\boldsymbol\beta}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\)이 주는 표준오차가 \((0.0710, 0.0882, 0.1041)\)인데, 실제 벗어난 폭을 그것으로 나누면 \((-0.41, 0.57, -0.81)\)로 셋 다 1 표준오차 안이다. 이만큼 벗어나는 것이 정상이라는 뜻이다.

\(\mathbf{X}^\top\mathbf{X}\) 쪽에서는 고윳값의 합 \(109.97\)이 대각합과 정확히 같다는 것, 그리고 조건수가 \(2.51\)로 작다는 것을 볼 수 있다. 조건수가 작으니 열들이 서로 충분히 갈라져 있어 \((\mathbf{X}^\top\mathbf{X})^{-1}\)을 안심하고 쓸 수 있다. 앞에서 본 오른쪽 그림처럼 한 열이 다른 열의 배수에 가까워지면 가장 작은 고윳값이 \(0\)으로 내려가고 조건수가 터져, 같은 공식이 수치적으로 못 쓰게 된다.

(a) \(\mathbf{X}^\top \mathbf{X} = \begin{pmatrix} 3 & 12 \\ 12 & 56 \end{pmatrix}\), \(\mathbf{X}^\top \mathbf{y} = \begin{pmatrix} 27 \\ 124 \end{pmatrix}\).

(b) \(\det(\mathbf{X}^\top \mathbf{X}) = 168 - 144 = 24\)이므로

(c) \(\hat{\mathbf{y}} = (5, 9, 13)^\top = \mathbf{y}\)이고 잔차는 모두 0이다. 세 점이 한 직선 위에 있으므로 적합이 정확하다.

직접 계산한다.

순서를 바꾸는 데에는 합이 유한하다는 사실만 쓰였다. \(\square\)

대칭성: \((\mathbf{X}^\top \mathbf{X})^{-1}\)은 대칭이므로(대칭행렬의 역행렬)

멱등성:

대각합: 순환 성질에 의해

\(\square\)

순환 성질과 \(\mathbf{Q}^\top \mathbf{Q} = \mathbf{I}\)를 쓰면

이다. 행렬식의 경우, \(\det\)는 곱셈적이고 직교행렬 \(\mathbf{Q}\)에 대해 \(\det(\mathbf{Q}) = \pm 1\)이므로

이다. \(\square\)

\(\mathbb{E}[\mathbf{Z}] = 0\)이고 \(\operatorname{Cov}(\mathbf{Z}) = \boldsymbol{\Sigma}\)인 \(\mathbf{Z}\)를 써서 \(\mathbf{X} = \boldsymbol{\mu} + \mathbf{Z}\)로 놓자. 전개하면(\(\mathbf{A}\)의 대칭성을 이용한다)

이다. 기댓값을 취하면 \(\mathbb{E}[\mathbf{Z}] = 0\)이므로 가운데 항이 사라진다. 마지막 항의 경우 \(\mathbf{Z}^\top \mathbf{A} \mathbf{Z}\)가 스칼라이므로 자기 자신의 대각합과 같고,

이다. 여기에 결정론적인 항을 더하면 결과를 얻는다. \(\square\)

(a) 대수적으로: 특이성은 \(\det(\mathbf{X}^\top \mathbf{X}) = 0\)을 뜻하므로 \((\mathbf{X}^\top \mathbf{X})^{-1}\)이 존재하지 않는다. 정규방정식 \(\mathbf{X}^\top \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y}\)은 (우변이 \(\mathbf{X}^\top \mathbf{X}\)의 열공간 안에 있으므로) 해를 갖지만 무한히 많다. \(\boldsymbol{\beta}^*\)가 해이고 \(\mathbf{v} \in \ker(\mathbf{X})\)이면 \(\mathbf{X}\mathbf{v} = \mathbf{0}\)이므로 \(\boldsymbol{\beta}^* + \mathbf{v}\)도 이 방정식을 만족한다.

(b) 기하적으로: \(\hat{\mathbf{y}} = \mathbf{H}\mathbf{y}\)는 여전히 \(\mathbf{y}\)를 \(\operatorname{col}(\mathbf{X})\) 위로 직교사영한 유일한 벡터다. 그러나 \(\mathbf{X}\)의 열들이 선형종속이면 그 사영을 열들의 선형결합으로 나타내는 방법이 무한히 많고, 각각이 타당한 \(\hat{\boldsymbol{\beta}}\)을 준다. 적합값은 식별되지만 계수는 식별되지 않는다. 해결책은 종속인 열을 제거하거나, 능형회귀(\(\mathbf{X}^\top \mathbf{X}\)에 \(\lambda \mathbf{I}\)를 더해 가역성을 회복한다), 또는 유사역행렬(최소 노름 해를 준다)을 쓰는 것이다. \(\square\)

셋째 열이 첫째와 둘째 열의 합이므로 \(\operatorname{rank}(\mathbf{X}) = 2\)이다.

import numpy as np

X = np.array([[1., 2., 3.],
              [1., 3., 4.],
              [1., 4., 5.],
              [1., 5., 6.]])          # 3열 = 1열 + 2열

r = np.linalg.matrix_rank(X)
U, s, Vt = np.linalg.svd(X)
print(f"rank = {r},  특잇값 = {s.round(6)}")

null_basis = Vt[r:]                    # 0 특잇값에 딸린 오른쪽 특이벡터
row_basis = Vt[:r]
print(f"영공간 기저: {null_basis.round(4)}")
print(f"X v = {(X @ null_basis.T).ravel().round(12)}")
print(f"rank + nullity = {r} + {null_basis.shape[0]} = {r + null_basis.shape[0]} = 열 개수")
print(f"행공간 ⟂ 영공간: {np.allclose(row_basis @ null_basis.T, 0)}")

출력:

rank = 2,  특잇값 = [11.982576  0.646436  0.      ]
영공간 기저: [[-0.5774 -0.5774  0.5774]]
X v = [0. 0. 0. 0.]
rank + nullity = 2 + 1 = 3 = 열 개수
행공간 ⟂ 영공간: True

세 번째 특잇값이 정확히 \(0\)이라는 것이 계수 결손의 신호다. 영공간 기저 \(\propto (1, 1, -1)\)은 "첫째 열 더하기 둘째 열 빼기 셋째 열 = 0"이라는 종속관계를 그대로 읽어 준다.

회귀를 돌리면. \(\mathbf{X}^\top\mathbf{X}\)가 특이행렬이라 \(\hat{\boldsymbol{\beta}}\)이 유일하지 않다(연습문제 6). \(\boldsymbol{\beta}^*\)가 해이면 임의의 \(c\)에 대해 \(\boldsymbol{\beta}^* + c(1,1,-1)^\top\)도 해이며, 넷 모두 똑같은 적합값을 준다. np.linalg.lstsq는 오류를 내지 않고 최소 노름 해 하나를 조용히 돌려주므로 더 위험하다. 계수의 부호나 크기를 해석하려 들면 아무 의미 없는 수를 해석하게 된다.

실무에서 이 상황은 완전한 다중공선성으로 나타난다. 가장 흔한 원인은 범주형 변수를 원-핫 인코딩하면서 기준 범주를 빼지 않은 것이다(더미변수 함정). \(\square\)

\(\boldsymbol{\mu} = \mathbb{E}[\mathbf{X}]\)라 하면 \(\mathbb{E}[\mathbf{A}\mathbf{X}] = \mathbf{A}\boldsymbol{\mu}\)이므로

이다. \(\mathbf{A}\)가 상수라 기댓값 밖으로 빠져나온 것이 전부다.

양반정치성. \(\mathbf{A}\) 자리에 행벡터 \(\mathbf{c}^\top\)를 넣으면 좌변이 스칼라 확률변수의 분산이므로

이고, 이것이 모든 \(\mathbf{c}\)에 대해 성립한다. 즉 공분산행렬의 양반정치성은 "분산은 음수가 될 수 없다"를 행렬로 옮겨 적은 것일 뿐이다. 등호는 \(\mathbf{c}^\top\mathbf{X}\)가 상수일 때, 곧 변수들 사이에 완전한 선형관계가 있을 때만 성립한다. \(\square\)

import numpy as np

Sigma = np.array([[0.04,  0.006,  0.000],
                  [0.006, 0.09,  -0.012],
                  [0.000, -0.012, 0.16]])

w = np.full(3, 1 / 3)
print(f"개별 표준편차:        {np.sqrt(np.diag(Sigma)).round(3)}")
print(f"개별 표준편차의 평균: {np.sqrt(np.diag(Sigma)).mean():.4f}")
print(f"동일가중 포트폴리오:   분산 {w @ Sigma @ w:.6f}, 표준편차 {np.sqrt(w @ Sigma @ w):.4f}")

Si = np.linalg.inv(Sigma)
one = np.ones(3)
w_mv = Si @ one / (one @ Si @ one)          # 최소분산 포트폴리오
print(f"\n최소분산 가중치: {w_mv.round(4)}")
print(f"최소분산 표준편차: {np.sqrt(w_mv @ Sigma @ w_mv):.4f}")

출력:

개별 표준편차:        [0.2 0.3 0.4]
개별 표준편차의 평균: 0.3000
동일가중 포트폴리오:   분산 0.030889, 표준편차 0.1758

최소분산 가중치: [0.5721 0.2561 0.1718]
최소분산 표준편차: 0.1563

동일가중 포트폴리오의 표준편차 \(0.1758\)은 개별 표준편차의 평균 \(0.3\)보다 훨씬 작다. 이것이 분산투자다. 비대각 성분이 작거나 음수라서 \(\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w}\)가 대각 성분의 가중평균보다 작아지는 것이다.

최소분산 포트폴리오는 여기서 더 나아가 \(0.1563\)까지 낮추는데, 변동성이 가장 낮은 첫 자산에 가장 큰 비중(\(0.57\))을 준다. 상관이 완전하다면(\(\rho = 1\)) 비대각 성분이 커져 이런 이득이 사라진다. 위기 때 상관이 \(1\)로 몰리면서 분산투자 효과가 증발하는 현상이 바로 이 계산에 들어 있다.

기울기 규칙. \(\mathbf{a}^\top\mathbf{x} = \sum_j a_j x_j\)이므로 \(\partial/\partial x_k = a_k\), 곧 \(\nabla = \mathbf{a}\)이다. 이차형식은 \(\mathbf{x}^\top\mathbf{A}\mathbf{x} = \sum_{i,j}a_{ij}x_ix_j\)이므로

이고, \(\mathbf{A}\)가 대칭이면 \(2\mathbf{A}\mathbf{x}\)가 된다.

유도. 목적함수를 전개하면

이다(\(\mathbf{y}^\top\mathbf{X}\boldsymbol{\beta}\)가 스칼라라 전치와 같으므로 두 교차항이 합쳐졌다). 위의 두 규칙을 쓰면

최소점 확인. 헤세행렬은 \(\nabla^2 S = 2\mathbf{X}^\top\mathbf{X}\)이고, 임의의 \(\mathbf{v} \ne \mathbf{0}\)에 대해

이다. \(\mathbf{X}\)가 완전 열계수이면 \(\mathbf{X}\mathbf{v} \ne \mathbf{0}\)이라 엄격한 양정치가 되어 \(S\)가 순볼록이고, 정류점은 유일한 전역 최소점이다. \(\square\)

import numpy as np

rng = np.random.default_rng(1)
n, p = 40, 3
X = np.column_stack([np.ones(n), rng.normal(size=(n, p - 1))])
y = X @ np.array([1., 2., -1.]) + rng.normal(size=n)

beta = np.array([0.3, -0.7, 1.1])          # 아무 점에서나
S = lambda b: ((y - X @ b) ** 2).sum()

grad_analytic = -2 * X.T @ y + 2 * X.T @ X @ beta
h = 1e-6
grad_numeric = np.array([(S(beta + h * e) - S(beta - h * e)) / (2 * h)
                         for e in np.eye(p)])
print("해석적 기울기:", grad_analytic.round(6))
print("수치적 기울기:", grad_numeric.round(6))

hess = 2 * X.T @ X
print(f"\n헤세행렬 고윳값: {np.linalg.eigvalsh(hess).round(3)}  (모두 > 0)")
print(f"정류점 = OLS 해: "
      f"{np.allclose(np.linalg.solve(X.T @ X, X.T @ y), np.linalg.lstsq(X, y, rcond=None)[0])}")

출력:

해석적 기울기: [ -69.943347 -130.877197   77.167931]
수치적 기울기: [ -69.943347 -130.877197   77.167931]

헤세행렬 고윳값: [44.717 64.702 84.593]  (모두 > 0)
정류점 = OLS 해: True

수치미분이 해석적 기울기와 소수점 여섯째 자리까지 맞는다. 새 손실함수를 유도했을 때 이런 대조는 언제나 해 볼 만한 값싼 점검이다.

import numpy as np

rng = np.random.default_rng(7)
n, p = 12, 3
X = np.column_stack([np.ones(n), rng.normal(size=(n, p - 1))])
y = rng.normal(size=n)

XtXi = np.linalg.inv(X.T @ X)
H = X @ XtXi @ X.T
e = y - X @ (XtXi @ X.T @ y)
h = np.diag(H)

loo = np.empty(n)                       # 무식하게: n 번 다시 적합
for i in range(n):
    keep = np.ones(n, bool)
    keep[i] = False
    loo[i] = y[i] - X[i] @ np.linalg.lstsq(X[keep], y[keep], rcond=None)[0]

print(f"공식과 실제 재적합의 최대 차이: {np.abs(loo - e / (1 - h)).max():.2e}")
print(f"지렛값 합 = {h.sum():.4f}  (= p = {p})")
print(f"\nSSE   = {(e ** 2).sum():.4f}   (같은 자료로 평가)")
print(f"PRESS = {(loo ** 2).sum():.4f}   (빼놓은 자료로 평가)")

출력:

공식과 실제 재적합의 최대 차이: 4.55e-15
지렛값 합 = 3.0000  (= p = 3)

SSE   = 7.3015   (같은 자료로 평가)
PRESS = 13.1653   (빼놓은 자료로 평가)

왜 중요한가. \(n\)번 다시 적합하는 대신 한 번의 적합에서 나온 \(e_i\)와 \(h_{ii}\)만으로 하나 빼기 교차검증을 끝낼 수 있다. 회귀에서 LOOCV가 사실상 공짜인 이유다.

\(h_{ii} \to 1\)일 때. 분모가 \(0\)으로 가면서 예측오차가 폭발한다. 뜻은 이렇다. 관측 \(i\)가 계획공간에서 워낙 외따로 떨어져 있어 그 점을 빼면 그 위치를 짚어 줄 다른 관측이 없다. 실제로 \(h_{ii} = 1\)이면 회귀선이 그 점을 정확히 통과해 \(e_i = 0\)이 되고, 모형은 그 점에 대해 아무것도 검증하지 못한다.

이 때문에 \(h_{ii} > 2p/n\)인 관측을 높은 지렛점으로 보아 따로 살핀다. 쿡의 거리도 \(e_i\)와 \(h_{ii}\)를 같은 방식으로 결합해 만든다.

위 출력에서 PRESS가 SSE의 두 배쯤 된다는 점도 눈여겨보라. 같은 자료로 적합하고 같은 자료로 평가하면 성능이 낙관적으로 나온다는 사실이 한 줄로 드러난다. \(\square\)