직교사영행렬¶

먼저 다음을 계산한다.

그러면

이다. 대각합은 \(\operatorname{tr}(\mathbf{H}) = (5 + 2 + 5)/6 = 12/6 = 2\)로 \(\mathbf{X}\)의 열 개수(\(p = 2\))와 같다. 이는 일반적인 결과 \(\operatorname{tr}(\mathbf{H}) = \operatorname{rank}(\mathbf{H}) = p\)를 확인해 준다.

\(\mathbf{y} = \hat{\mathbf{y}} + \mathbf{e}\)이므로

이다. 교차항은 다음과 같이 사라진다.

\(\mathbf{H}(\mathbf{I} - \mathbf{H}) = \mathbf{H} - \mathbf{H}^2 = \mathbf{H} - \mathbf{H} = \mathbf{0}\)이므로 교차항이 0이 되어 \(\lVert \mathbf{y} \rVert^2 = \lVert \hat{\mathbf{y}} \rVert^2 + \lVert \mathbf{e} \rVert^2\)를 얻는다. \(\square\)

\(\mathbf{z} = \mathbf{X}\boldsymbol{\beta}\)을 \(\text{col}(\mathbf{X})\)의 임의의 벡터라 하자. \(\lVert \mathbf{y} - \hat{\mathbf{y}} \rVert \leq \lVert \mathbf{y} - \mathbf{z} \rVert\)를 보이면 된다.

\(\mathbf{y} - \mathbf{z} = (\mathbf{y} - \hat{\mathbf{y}}) + (\hat{\mathbf{y}} - \mathbf{z})\)로 쓰자. \(\mathbf{y} - \hat{\mathbf{y}} = \mathbf{e} \in \text{col}(\mathbf{X})^\perp\)이고 \(\hat{\mathbf{y}} - \mathbf{z} \in \text{col}(\mathbf{X})\)이므로 이 두 벡터는 직교한다. 피타고라스 정리에 의해

이다. 등호는 \(\mathbf{z} = \hat{\mathbf{y}}\)일 때에 한해 성립하므로 \(\hat{\mathbf{y}}\)이 유일한 최근접점임이 확인된다. \(\square\)

중심화한 두 벡터 \(\tilde{\mathbf{y}} = \mathbf{y} - \bar{y}\mathbf{1}\)과 \(\tilde{\mathbf{y}}_{\!*} = \hat{\mathbf{y}} - \bar{y}\mathbf{1}\)을 잡고 \(\theta\)를 둘 사이의 각이라 하자. \(\tilde{\mathbf{y}} = \tilde{\mathbf{y}}_{\!*} + \mathbf{e}\)이고 \(\mathbf{e}\)가 \(\hat{\mathbf{y}}\)과 \(\mathbf{1}\) 모두에 직교하므로 \(\tilde{\mathbf{y}}_{\!*}^\top\mathbf{e} = 0\)이며, 따라서

이다. 곧 \(R^2 = \text{SSR}/\text{SST} = \cos^2\theta\)로, 평균만 맞히는 모형에서 출발한 반응벡터와 그것을 모형 부분공간 위로 사영한 벡터 사이 각의 코사인의 제곱이다.

\(R^2 = 1\)은 \(\theta = 0\), 곧 \(\mathbf{e} = \mathbf{0}\)이므로 \(\mathbf{y}\)가 정확히 \(\text{col}(\mathbf{X})\) 안에 놓인다는 뜻이다. 자료가 잔차 없이 모형에 완벽히 들어맞는다. \(R^2 = 0\)은 \(\tilde{\mathbf{y}}_{\!*} = \mathbf{0}\), 곧 \(\hat{\mathbf{y}} = \bar{y}\mathbf{1}\)이므로 사영이 절편 방향밖에 쓰지 못한 경우다. 예측변수가 아무것도 설명하지 못한다.

중심화하지 않은 비 \(\lVert\hat{\mathbf{y}}\rVert^2/\lVert\mathbf{y}\rVert^2\)은 같은 계산을 원점 기준으로 한 것이어서 \(\mathbf{y}\)와 \(\hat{\mathbf{y}}\) 사이 각의 코사인 제곱이다. 절편이 있는 모형에서는 \(\lVert\mathbf{y}\rVert^2\) 안에 평균 수준 \(n\bar{y}^2\)이 들어 있고 그 부분은 절편이 언제나 정확히 맞히므로, 이 비는 중심화한 \(R^2\)보다 크게 나온다. 실제로 \(\lVert\hat{\mathbf{y}}\rVert^2 = \text{SSR} + n\bar{y}^2\), \(\lVert\mathbf{y}\rVert^2 = \text{SST} + n\bar{y}^2\)이므로 분자와 분모에 같은 양을 더한 꼴이며, 두 값이 일치하는 것은 \(\bar{y} = 0\)이거나 \(R^2 = 1\)일 때뿐이다. 연습문제 10에서 같은 자료에 대해 두 값을 나란히 계산한다.

사영 공식 \(\mathbf{P}_{\mathbf{X}} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\)에서 \(\mathbf{X} = \mathbf{a}\)(\(n \times 1\))로 두면 \(\mathbf{X}^\top\mathbf{X} = \mathbf{a}^\top\mathbf{a}\)가 스칼라이므로 역행렬이 곧 역수다.

멱등성은 \(\mathbf{P}^2 = \dfrac{\mathbf{a}(\mathbf{a}^\top\mathbf{a})\mathbf{a}^\top}{(\mathbf{a}^\top\mathbf{a})^2} = \mathbf{P}\), 대칭성은 \((\mathbf{a}\mathbf{a}^\top)^\top = \mathbf{a}\mathbf{a}^\top\)에서 나온다. 계수가 1이므로 \(\operatorname{tr}(\mathbf{P}) = 1\)이다.

\(\mathbf{P}\mathbf{y} = \mathbf{a}\dfrac{\mathbf{a}^\top\mathbf{y}}{\mathbf{a}^\top\mathbf{a}}\)로 쓰면 계수가 바로 읽힌다.

이것이 절편 없는 단순회귀 \(y_i = \beta a_i + \varepsilon_i\)의 최소제곱추정량이다.

import numpy as np

a = np.array([1., 2., 2.])
y = np.array([3., 1., 4.])

P = np.outer(a, a) / (a @ a)
print("P =\n", P.round(4))
print("멱등:", np.allclose(P @ P, P), " 대칭:", np.allclose(P, P.T),
      " tr:", round(np.trace(P), 6))
print("beta_hat =", round(a @ y / (a @ a), 6))
print("P y =", (P @ y).round(6), " = beta_hat * a =", (a @ y / (a @ a) * a).round(6))

출력:

P =
 [[0.1111 0.2222 0.2222]
 [0.2222 0.4444 0.4444]
 [0.2222 0.4444 0.4444]]
멱등: True  대칭: True  tr: 1.0
beta_hat = 1.444444
P y = [1.444444 2.888889 2.888889]  = beta_hat * a = [1.444444 2.888889 2.888889]

벡터 하나 위로의 사영은 앞으로 계속 쓰인다. 중심화행렬의 여집합 \(\frac{1}{n}\mathbf{J}\)가 \(\mathbf{1}\) 위로의 사영이고, 그람–슈미트의 각 단계도 이 형태다. \(\square\)

\(\mathbf{H}\)가 대칭 멱등이므로 \(\mathbf{H} = \mathbf{H}^2 = \mathbf{H}^\top\mathbf{H}\)이고, 따라서

이다. 오른쪽의 합이 음이 아니므로 \(h_{ii} \ge h_{ii}^2\), 곧 \(h_{ii}(1 - h_{ii}) \ge 0\)이고 \(0 \le h_{ii} \le 1\)이다. 합은 \(\sum_i h_{ii} = \operatorname{tr}(\mathbf{H}) = p\)다.

\(h_{ii} = 1\)인 경우. 위 식에서 \(j \neq i\)인 모든 \(h_{ji} = 0\)이어야 한다. 그러면 \(\hat{y}_i = \sum_j h_{ij}y_j = y_i\)이므로 그 관측값은 정확히 맞춰진다. 잔차가 항상 0이고, 그 점은 회귀선을 자기 쪽으로 완전히 끌어당긴다.

import numpy as np

rng = np.random.default_rng(0)
n = 40
X = np.column_stack([np.ones(n), rng.normal(size=n), rng.normal(size=n)])
H = X @ np.linalg.inv(X.T @ X) @ X.T
h = np.diag(H)

print(f"지렛대 최소 {h.min():.4f}, 최대 {h.max():.4f}")
print(f"합 = {h.sum():.6f}  (= p = {X.shape[1]})")
print(f"평균 = {h.mean():.4f}  (= p/n = {X.shape[1]/n:.4f})")

출력:

지렛대 최소 0.0285, 최대 0.2371
합 = 3.000000  (= p = 3)
평균 = 0.0750  (= p/n = 0.0750)

지렛대의 평균이 \(p/n\)이므로 \(2p/n\)을 넘는 점을 주의해서 보라는 실무 규칙이 여기서 나온다. 2.3절 이상치와 지렛대점의 논의와 이어진다. \(\square\)

임의의 \(\mathbf{x}\)에 대해 \(\mathbf{H}_1\mathbf{x} \in \operatorname{col}(\mathbf{X}_1) \subseteq \operatorname{col}(\mathbf{X}_2)\)이다. 사영은 자기 치역의 벡터를 그대로 두므로 \(\mathbf{H}_2(\mathbf{H}_1\mathbf{x}) = \mathbf{H}_1\mathbf{x}\)이고, 따라서 \(\mathbf{H}_2\mathbf{H}_1 = \mathbf{H}_1\)이다.

두 행렬 모두 대칭이므로 전치를 취하면

이다.

F 검정의 근거. 이 성질에서 \(\mathbf{H}_2 - \mathbf{H}_1\)도 사영임이 따라 나온다.

게다가 \((\mathbf{H}_2 - \mathbf{H}_1)(\mathbf{I} - \mathbf{H}_2) = \mathbf{O}\)이므로 두 사영이 직교한다. 그래서 큰 모형과 작은 모형의 제곱합 차이 \(\lVert(\mathbf{H}_2 - \mathbf{H}_1)\mathbf{y}\rVert^2\)과 잔차제곱합 \(\lVert(\mathbf{I}-\mathbf{H}_2)\mathbf{y}\rVert^2\)이 (정규성 아래에서) 독립인 카이제곱이 되고, 그 비가 \(F\) 분포를 따른다. 내포모형 F 검정의 기하가 바로 이것이다. \(\square\)

\(\mathbf{P}_1\)과 \(\mathbf{P}_2\)가 모두 \(\mathcal{V}\) 위로의 직교사영이라 하자. 연습문제 7의 논법을 양방향으로 쓴다. 두 치역이 같으므로 \(\mathbf{P}_2\mathbf{P}_1 = \mathbf{P}_1\)이고 \(\mathbf{P}_1\mathbf{P}_2 = \mathbf{P}_2\)다.

대칭성에서

이므로 둘은 같다.

왜 중요한가. 사영 공식 \(\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\)은 겉보기에 \(\mathbf{X}\)에 의존하지만, 실제로는 \(\operatorname{col}(\mathbf{X})\)에만 의존한다. 같은 열공간을 주는 다른 계획행렬(예: 예측변수를 재척도화하거나 선형결합한 것)을 써도 모자 행렬은 똑같다. 회귀에서 적합값 \(\hat{\mathbf{y}}\)과 \(R^2\)이 모수화 방식에 영향받지 않는 이유가 이것이다. 반면 계수 \(\hat{\boldsymbol{\beta}}\)은 모수화에 따라 달라진다.

빗각 사영에서는 이 유일성이 성립하지 않는다. 치역이 같아도 눌러 없애는 방향이 다르면 다른 사영이다. \(\square\)

이것이 Frisch–Waugh–Lovell 정리다. 다중회귀의 한 계수는 "다른 변수들의 영향을 걷어낸 뒤"의 단순회귀 계수와 같다.

import numpy as np

rng = np.random.default_rng(0)
n = 40
x1 = np.ones(n)
x2 = rng.normal(size=n)
x3 = 0.5 * x2 + rng.normal(size=n)          # x2 와 상관된 변수
y = 1 + 2 * x2 - 1.5 * x3 + rng.normal(size=n)

# (1) 전체 다중회귀
X = np.column_stack([x1, x2, x3])
beta_full = np.linalg.lstsq(X, y, rcond=None)[0]

# (2) x1, x2 의 영향을 걷어낸 뒤 단순회귀
X1 = np.column_stack([x1, x2])
M1 = np.eye(n) - X1 @ np.linalg.inv(X1.T @ X1) @ X1.T
y_res, x3_res = M1 @ y, M1 @ x3
beta_fwl = (x3_res @ y_res) / (x3_res @ x3_res)

print("다중회귀의 x3 계수 :", round(beta_full[2], 8))
print("FWL 로 구한 계수   :", round(beta_fwl, 8))

출력:

다중회귀의 x3 계수 : -1.44939995
FWL 로 구한 계수   : -1.44939995

두 값이 소수점 여덟째 자리까지 같다.

해석적 의미가 크다. 다중회귀의 계수 \(\hat\beta_3\)은 "\(x_3\)이 한 단위 늘 때 \(y\)의 변화"가 아니라 "다른 변수로 설명되지 않는 부분의 \(x_3\)이 한 단위 늘 때, 다른 변수로 설명되지 않는 부분의 \(y\)의 변화"다. 이것이 다중회귀 계수를 "다른 변수를 통제했을 때의 효과"로 읽는 근거이며, 13.5절의 부분회귀 그림(추가변수 그림)의 바탕이기도 하다. 그 그림은 정확히 \(\mathbf{M}_1\mathbf{y}\)를 \(\mathbf{M}_1\mathbf{x}_2\)에 대해 찍은 산점도다. \(\square\)

import numpy as np

rng = np.random.default_rng(1)
n, p = 30, 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)

H = X @ np.linalg.inv(X.T @ X) @ X.T
y_hat = H @ y
e = y - y_hat

print("대칭   :", np.allclose(H, H.T))
print("멱등   :", np.allclose(H @ H, H))
print("tr(H)  :", round(np.trace(H), 6), " (= p =", p, ")")
print("y_hat . e :", round(float(y_hat @ e), 12), " (직교)")

print("\n||y||^2          =", round(y @ y, 6))
print("||y_hat||^2 + ||e||^2 =", round(y_hat @ y_hat + e @ e, 6))
print("원점 기준 cos^2(theta) =", round((y_hat @ y_hat) / (y @ y), 6))

# 중심화: 1 위로의 사영을 걷어낸 뒤 같은 분해를 반복한다
yt = y - y.mean()
yt_hat = y_hat - y.mean()
sst, ssr, sse = yt @ yt, yt_hat @ yt_hat, e @ e
print("\nSST =", round(sst, 6), " SSR =", round(ssr, 6), " SSE =", round(sse, 6))
print("SSR + SSE =", round(ssr + sse, 6))
print("R^2 = SSR/SST     =", round(ssr / sst, 6))
print("R^2 = 1 - SSE/SST =", round(1 - sse / sst, 6))

출력:

대칭   : True
멱등   : True
tr(H)  : 3.0  (= p = 3 )
y_hat . e : -0.0  (직교)

||y||^2          = 158.337856
||y_hat||^2 + ||e||^2 = 158.337856
원점 기준 cos^2(theta) = 0.885583

SST = 121.938838  SSR = 103.822325  SSE = 18.116513
SSR + SSE = 121.938838
R^2 = SSR/SST     = 0.85143
R^2 = 1 - SSE/SST = 0.85143

네 성질이 모두 확인된다. 잔차와 적합값의 내적은 반올림 오차 수준에서 0이고, 피타고라스 분해는 원점 기준(\(158.337856\))과 평균 기준(\(121.938838\)) 양쪽에서 모두 성립한다.

두 \(R^2\)이 다르다. 원점 기준 값은 \(0.885583\)이고 보고되는 통상의 \(R^2\)은 \(0.851430\)이다. 차이는 \(\bar{y} = 1.1015\) 때문이다. \(\lVert\mathbf{y}\rVert^2 = \text{SST} + n\bar{y}^2 = 121.938838 + 30 \times 1.1015^2 = 158.34\)에서 보듯, 원점 기준 비는 절편이 공짜로 맞히는 평균 수준까지 설명력으로 세는 셈이다. \(R^2\)을 보고할 때는 \(\mathbf{1}\) 위로의 사영을 먼저 걷어낸 중심화된 쪽을 쓴다. \(\square\)