그람 행렬¶

직접 곱하면

이다.

대칭성은 \((\mathbf{X}^\top\mathbf{X})^\top = \mathbf{X}^\top\mathbf{X}\)로부터 곧바로 따라온다. 양정치성의 경우 선행 주소행렬식이 \(3 > 0\)이고 \(\det = 3 \times 38 - 10^2 = 114 - 100 = 14 > 0\)이므로 \(\mathbf{X}^\top\mathbf{X}\)는 양정치다. 동등하게 \(\mathbf{X}\)의 계수가 2이므로(두 열이 일차독립이므로) \(\mathbf{X}^\top\mathbf{X}\)가 양정치다.

임의의 벡터 \(\mathbf{v} \neq \mathbf{0}\)에 대해

이다. 이것이 0일 필요충분조건은 \(\mathbf{X}\mathbf{v} = \mathbf{0}\), 즉 \(\mathbf{v} \in \ker(\mathbf{X})\)인 것이다. \(\mathbf{X}\)가 완전 열계수를 가지면 \(\ker(\mathbf{X}) = \{\mathbf{0}\}\)이므로 모든 \(\mathbf{v} \neq \mathbf{0}\)에 대해 \(\mathbf{v}^\top(\mathbf{X}^\top\mathbf{X})\mathbf{v} > 0\)이고 이것이 양정치성이다. 역으로 \(\mathbf{X}\)가 완전 열계수를 갖지 않으면 \(\mathbf{X}\mathbf{v} = \mathbf{0}\)인 \(\mathbf{v} \neq \mathbf{0}\)이 존재하여 이차형식이 0이 된다. \(\square\)

\(\mathbf{X}^\top\mathbf{X}\)의 조건수는 \(\kappa = \lambda_{\max}/\lambda_{\min}\)이며, 여기서 \(\lambda_{\max}\)와 \(\lambda_{\min}\)은 최대·최소 고윳값이다.

열들이 거의 공선적이면 \(\lambda_{\min}\)이 0에 가까워져 \(\kappa\)가 매우 커진다. \(\operatorname{Var}(\hat{\boldsymbol{\beta}}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\)이고 \((\mathbf{X}^\top\mathbf{X})^{-1}\)의 고윳값이 \(1/\lambda_i\)이므로, \(\lambda_{\min}\)이 작으면 그에 대응하는 고유벡터 방향으로 분산 \(\sigma^2/\lambda_{\min}\)이 커진다. 조건수가 크다는 것은 또한 \(\mathbf{y}\)의 작은 섭동이 \(\hat{\boldsymbol{\beta}}\)을 크게 변화시킨다는 뜻이며, 추정값이 수치적으로 불안정해진다.

\(\mathbf{S} = \frac{1}{n-1}\mathbf{X}_c^\top\mathbf{X}_c\)는 그람 행렬의 양의 스칼라배이므로 양반정치성을 물려받는다.

\(\mathbf{S}\)가 양정치일 필요충분조건은 \(\mathbf{X}_c\)가 완전 열계수를 갖는 것이며, 이를 위해서는 \(n - 1 \geq p\)가 필요하다(중심화가 계수를 많아야 1만큼 줄이기 때문이다). 실무적으로는 표본 공분산행렬이 가역이려면 변수보다 관측값이 많아야 한다는(\(n > p\)) 뜻이다. \(p > n\)이면 \(\mathbf{S}\)가 특이행렬이 되어 정칙화나 차원축소 같은 기법이 필요하다. \(\square\)

\([\mathbf{G}]_{ij} = \mathbf{v}_i^\top\mathbf{v}_j\)인데 \(i \neq j\)이면 직교성에 의해 \(\mathbf{v}_i^\top\mathbf{v}_j = 0\)이다. 따라서 비대각 성분이 모두 0이고

이다. 대각행렬의 역행렬은 각 성분의 역수이므로 정규방정식 \(\mathbf{G}\hat{\boldsymbol{\beta}} = \mathbf{X}^\top\mathbf{y}\)가 \(p\)개의 독립적인 한 변수 방정식으로 분리된다.

import numpy as np

X = np.array([[1., 1.], [1., -1.], [1., 1.], [1., -1.]])   # 두 열이 직교
y = np.array([3., 1., 4., 2.])

G = X.T @ X
print("그람 행렬:\n", G)

beta_ls = np.linalg.solve(G, X.T @ y)                       # 정규방정식
beta_sep = np.array([X[:, j] @ y / (X[:, j] @ X[:, j])      # 열마다 따로
                     for j in range(X.shape[1])])
print("정규방정식 해:", beta_ls)
print("열별 계산    :", beta_sep)

출력:

그람 행렬:
 [[4. 0.]
 [0. 4.]]
정규방정식 해: [2.5 1. ]
열별 계산    : [2.5 1. ]

통계적 의미가 크다. 예측변수가 직교하면 한 변수를 모형에 넣거나 빼도 다른 변수의 계수가 전혀 변하지 않는다. 실험계획에서 직교설계를 선호하는 이유이며, 관측자료에서 계수가 모형 설정에 따라 요동치는 이유이기도 하다. \(\square\)

\(\mathbf{X}^\top\mathbf{X}\mathbf{v} = \lambda\mathbf{v}\)이고 \(\lambda \neq 0\)이라 하자. 양변에 왼쪽에서 \(\mathbf{X}\)를 곱하면

이다. \(\lambda \neq 0\)이므로 \(\mathbf{X}\mathbf{v} \neq \mathbf{0}\)이고(만약 \(\mathbf{X}\mathbf{v} = \mathbf{0}\)이면 \(\lambda\mathbf{v} = \mathbf{X}^\top\mathbf{X}\mathbf{v} = \mathbf{0}\)이 되어 모순), 따라서 \(\mathbf{X}\mathbf{v}\)가 \(\mathbf{X}\mathbf{X}^\top\)의 고윳값 \(\lambda\)에 대응하는 고유벡터다. 반대 방향도 \(\mathbf{X}^\top\)를 곱해 같은 방식으로 보인다.

두 행렬의 크기가 다르므로 개수가 많은 쪽에는 0인 고윳값이 채워진다.

import numpy as np

rng = np.random.default_rng(0)
X = rng.normal(size=(5, 3))                 # n=5, p=3

e_small = np.sort(np.linalg.eigvalsh(X.T @ X))[::-1]     # 3개
e_big = np.sort(np.linalg.eigvalsh(X @ X.T))[::-1]       # 5개

print("X^T X 의 고윳값:", e_small.round(6))
print("X X^T 의 고윳값:", e_big.round(6))

출력:

X^T X 의 고윳값: [9.928184 2.921811 0.112128]
X X^T 의 고윳값: [ 9.928184  2.921811  0.112128  0.       -0.      ]

\(\mathbf{X}\mathbf{X}^\top\)의 큰 세 고윳값이 \(\mathbf{X}^\top\mathbf{X}\)의 고윳값과 정확히 같고 나머지 두 개는 0이다.

\(p \gg n\)이면 \(n \times n\)인 \(\mathbf{X}\mathbf{X}^\top\)로 계산해야 한다. 유전체 자료처럼 \(p\)가 수만이고 \(n\)이 수백인 상황에서 \(p \times p\) 행렬은 다루기 어렵지만 \(n \times n\)은 작다. 두 행렬이 같은 정보를 담고 있으므로 손해가 없다. 이것이 커널 트릭의 대수적 근거다. \(\square\)

\(\det(\mathbf{X}^\top\mathbf{X})\)를 그람 행렬식이라 하며 열들이 펼치는 평행육면체 부피의 제곱과 같다.

import numpy as np

# 서로 직교하고 길이가 1, 2인 두 벡터 -> 넓이 2, 그람 행렬식 4
X1 = np.array([[1., 0.], [0., 2.], [0., 0.]])
print("직교:      det =", round(np.linalg.det(X1.T @ X1), 6), " (넓이^2 = 4)")

# 같은 길이지만 방향이 겹치면 넓이가 줄어든다
X2 = np.array([[1., 1.], [0., 1.], [0., 0.]])
print("비스듬함:  det =", round(np.linalg.det(X2.T @ X2), 6))

# 두 열이 일차종속이면 부피가 0
X3 = np.array([[1., 2.], [0., 0.], [0., 0.]])
print("공선:      det =", round(np.linalg.det(X3.T @ X3), 6))

출력:

직교:      det = 4.0  (넓이^2 = 4)
비스듬함:  det = 1.0
공선:      det = 0.0

기하적 의미. \(\det(\mathbf{X}^\top\mathbf{X}) = 0\)은 부피가 0이라는 뜻이고, 이는 열들이 더 낮은 차원의 부분공간에 눌려 있다는 것, 곧 일차종속이라는 뜻이다.

이것이 본문의 양정치성 정리와 정확히 같은 이야기를 기하로 옮긴 것이다. 완전 열계수 \(\iff\) 부피 \(> 0\) \(\iff\) \(\mathbf{X}^\top\mathbf{X} \succ 0\) \(\iff\) 최소제곱해가 유일하다.

부피가 0은 아니지만 매우 작은 경우가 곧 다중공선성이다. 해가 존재하기는 하지만 납작한 평행육면체 위에서 결정되므로 불안정하다. \(\square\)

열 \(j\)를 중심화한 뒤 그 노름으로 나누면 \(\mathbf{z}_j = (\mathbf{v}_j - \bar{v}_j\mathbf{1}) / \lVert \mathbf{v}_j - \bar{v}_j\mathbf{1} \rVert\)이다. 그러면

로 정확히 표본상관계수다. 곧 \(\mathbf{Z}^\top\mathbf{Z}\)가 상관행렬이다.

import numpy as np

rng = np.random.default_rng(2)
X = rng.normal(size=(100, 3)) @ np.array([[1., .8, .2],
                                          [0., .6, .3],
                                          [0., 0., 1.]])
Z = (X - X.mean(0)) / np.linalg.norm(X - X.mean(0), axis=0)

print("Z^T Z =\n", np.round(Z.T @ Z, 4))
print("\ncorrcoef =\n", np.round(np.corrcoef(X, rowvar=False), 4))

출력:

Z^T Z =
 [[1.     0.8133 0.1539]
 [0.8133 1.     0.3154]
 [0.1539 0.3154 1.    ]]

corrcoef =
 [[1.     0.8133 0.1539]
 [0.8133 1.     0.3154]
 [0.1539 0.3154 1.    ]]

세 가지가 같은 대상의 다른 척도임을 알 수 있다.

상관행렬이 언제나 양반정치인 이유도 여기서 나온다. 상관행렬 역시 그람 행렬이기 때문이다. \(\square\)

\(\mathbf{x}_3 = \mathbf{x}_1 + \varepsilon \cdot (\text{잡음})\)으로 두고 \(\varepsilon\)을 줄여 가며 관찰한다.

import numpy as np

rng = np.random.default_rng(0)
n = 50
x1 = rng.normal(size=n)
x2 = rng.normal(size=n)

for eps in (1.0, 0.1, 0.01):
    x3 = x1 + eps * rng.normal(size=n)          # eps가 작을수록 x1과 닮는다
    X = np.column_stack([x1, x2, x3])
    G = X.T @ X
    ev = np.linalg.eigvalsh(G)
    var_factor = np.diag(np.linalg.inv(G)).max()  # Var(beta)/sigma^2 의 최댓값
    print(f"eps={eps:<5}  cond={ev.max()/ev.min():>10.1f}   "
          f"max Var factor={var_factor:>8.3f}")

출력:

eps=1.0    cond=       7.1   max Var factor=   0.048
eps=0.1    cond=     382.7   max Var factor=   2.243
eps=0.01   cond=   27830.5   max Var factor= 162.058

\(\varepsilon\)을 \(1\)에서 \(0.01\)로 줄이면 조건수가 \(7\)에서 \(27{,}830\)으로 약 \(4{,}000\)배 커지고, 계수 분산의 최댓값은 \(0.048\)에서 \(162\)로 약 \(3{,}400\)배 커진다.

대수적 이유. \(\mathbf{x}_3 \approx \mathbf{x}_1\)이면 \(\mathbf{y} = \mathbf{x}_3 - \mathbf{x}_1\) 방향에서 \(\mathbf{X}\mathbf{y} \approx \mathbf{0}\)이므로 \(\mathbf{y}^\top\mathbf{G}\mathbf{y} = \lVert\mathbf{X}\mathbf{y}\rVert^2 \approx 0\)이다. 곧 \(\lambda_{\min} \approx 0\)이고, \(\operatorname{Var}(\hat{\boldsymbol{\beta}}) = \sigma^2\mathbf{G}^{-1}\)의 고윳값 \(\sigma^2/\lambda_{\min}\)이 폭발한다.

자료가 늘어난 것이 아니라 거의 같은 정보를 두 번 넣은 것이므로 추정이 불안정해지는 것이 당연하다. 18장의 능형회귀는 \(\mathbf{G} + \lambda\mathbf{I}\)로 최소고윳값을 끌어올려 이 문제를 완화한다. \(\square\)

\((\subseteq)\) \(\mathbf{X}^\top\mathbf{X}\mathbf{v} = \mathbf{X}^\top(\mathbf{X}\mathbf{v})\)이므로 \(\operatorname{Col}(\mathbf{X}^\top\mathbf{X}) \subseteq \operatorname{Col}(\mathbf{X}^\top)\)은 곧바로 나온다.

\((\supseteq)\) 차원을 센다. 정리 3 에 의해 \(\dim\operatorname{Col}(\mathbf{X}^\top\mathbf{X}) = \operatorname{rank}(\mathbf{X}^\top\mathbf{X}) = \operatorname{rank}(\mathbf{X}) = \operatorname{rank}(\mathbf{X}^\top) = \dim\operatorname{Col}(\mathbf{X}^\top)\)이다. 포함관계가 있고 차원이 같은 두 부분공간은 같으므로 결론이 따라 나온다.

정규방정식은 언제나 무모순이다. \(\mathbf{X}^\top\mathbf{y} \in \operatorname{Col}(\mathbf{X}^\top) = \operatorname{Col}(\mathbf{X}^\top\mathbf{X})\)이므로 \(\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} = \mathbf{X}^\top\mathbf{y}\)를 만족하는 \(\boldsymbol{\beta}\)가 적어도 하나 존재한다. 완전 열계수일 때에만 그것이 유일하다(정리 2). 계수가 \(r < p\)이면 해집합은 \(\ker(\mathbf{X})\) 방향으로 \(p - r\)차원만큼 평행이동한 아핀 집합이 되고, 적합값 \(\hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\beta}\)는 어느 해를 택해도 같다.

import numpy as np

X = np.array([[1., 2., 3.], [1., 0., 1.], [1., 1., 2.], [1., 3., 4.]])  # 3열 = 1열 + 2열
y = np.array([6., 2., 4., 8.])
G, b = X.T @ X, X.T @ y

print("rank(X) =", np.linalg.matrix_rank(X), "  rank(X^T X) =", np.linalg.matrix_rank(G))
print("X^T X 의 고윳값:", np.linalg.eigvalsh(G).round(6))
print("X^T y 가 Col(X^T X) 안에 있는가:",
      np.linalg.matrix_rank(G) == np.linalg.matrix_rank(np.column_stack([G, b])))

beta_min = np.linalg.lstsq(G, b, rcond=None)[0]           # 최소노름 해
beta_alt = beta_min + 5.0 * np.array([1., 1., -1.])        # ker(X) 방향으로 밀어낸 다른 해
print("\n해 1:", beta_min.round(6), " 잔차", round(np.linalg.norm(G @ beta_min - b), 10))
print("해 2:", beta_alt.round(6), " 잔차", round(np.linalg.norm(G @ beta_alt - b), 10))
print("두 해의 적합값이 같은가:", np.allclose(X @ beta_min, X @ beta_alt))

출력:

rank(X) = 2   rank(X^T X) = 2
X^T X 의 고윳값: [-0.        1.284367 46.715633]
X^T y 가 Col(X^T X) 안에 있는가: True

해 1: [0.666667 0.666667 1.333333]  잔차 0.0
해 2: [ 5.666667  5.666667 -3.666667]  잔차 0.0
두 해의 적합값이 같은가: True

계수가 3 에서 2 로 떨어져 고윳값 하나가 0 이 되었지만, 정규방정식은 여전히 무모순이다. 서로 다른 두 해가 같은 적합값을 주는 것이 핵심이다. 추정 불가능한 것은 \(\boldsymbol{\beta}\) 자체이고 \(\mathbf{X}\boldsymbol{\beta}\)는 여전히 유일하게 추정된다. 회귀에서 완전 공선성이 있을 때 "개별 계수는 해석할 수 없지만 예측은 할 수 있다"고 말하는 근거가 이것이다.

복소행렬에서는 다르다. 정리 3 의 증명이 실행렬임을 썼으므로 이 모든 논의도 실행렬에 한한다. 열벡터 \(\mathbf{X} = (1, i)^\top\)를 잡으면 \(\mathbf{X}^\top\mathbf{X} = 1 + i^2 = 0\)이라 계수가 \(1\)에서 \(0\)으로 떨어지고, 켤레전치를 쓴 \(\mathbf{X}^*\mathbf{X} = 1 + \lvert i\rvert^2 = 2\)만이 계수를 보존한다. \(\square\)