불량조건 계획행렬¶

OLS 추정량은 행렬 \(\mathbf{X}^\top\mathbf{X}\)의 역행렬을 요구한다. 이 행렬이 특이행렬에 가까우면 역행렬 계산이 수치적으로 불안정해지고, 자료의 작은 교란이 추정 계수의 큰 변동으로 이어진다. 이 불안정성을 이해하면, 통계적인 편향-분산 논증과 별개로 수치적 안정화 도구로서의 정칙화가 보인다.

조건수¶

행렬 \(\mathbf{A}\)의 조건수는 \(\mathbf{A}\mathbf{x} = \mathbf{b}\)의 해가 \(\mathbf{b}\)의 교란에 얼마나 민감한지를 정량화한다. \(\mathbf{X}^\top\mathbf{X}\)처럼 대칭 양반정치 행렬에서는

이며 \(\lambda_{\max}\)와 \(\lambda_{\min}\)은 각각 최대·최소 고윳값이다.

\(\mathbf{X}\)의 특이값이 \(d_1 \geq d_2 \geq \cdots \geq d_p \geq 0\)이면 \(\mathbf{X}^\top\mathbf{X}\)의 고윳값은 \(d_1^2, \ldots, d_p^2\)이므로

이다.

조건수가 1에 가까우면 잘 조건화된 문제다. 조건수가 \(10^k\)이면 연립방정식을 풀 때 정밀도가 대략 \(k\)자리 손실된다. \(\kappa(\mathbf{X}^\top\mathbf{X}) > 10^{10}\)이면 배정밀도 부동소수점 연산으로도 의미 없는 결과가 나올 수 있다.

불량조건의 원인¶

준선형종속. 두 개 이상의 설명변수가 거의 선형종속이면 \(\mathbf{X}^\top\mathbf{X}\)의 최소 고윳값이 0에 가까워지고 조건수가 무한대로 향한다.

고차원. \(p\)가 \(n\)에 가까우면 정확한 선형종속이 없어도 \(\mathbf{X}^\top\mathbf{X}\)가 매우 작은 고윳값을 갖는 경향이 있다. \(p > n\)이면 계수 부족 행렬이 되어 OLS 해가 유일하지 않다.

척도 불일치. 설명변수의 척도가 크게 다르면(하나는 미터, 다른 하나는 킬로미터) \(\mathbf{X}^\top\mathbf{X}\)의 고윳값 폭이 커진다. 회귀 전에 설명변수를 표준화하면 이 원인은 완화된다.

OLS의 교란 민감도¶

OLS 연립방정식 \((\mathbf{X}^\top\mathbf{X})\hat{\boldsymbol{\beta}} = \mathbf{X}^\top\mathbf{y}\)를 생각하자. 우변이 작은 양 \(\delta\)만큼 교란되면 해의 변화는 다음 한계를 만족한다.

즉 \(\hat{\boldsymbol{\beta}}\)의 상대오차가 교란의 상대크기 대비 최대 \(\kappa(\mathbf{X}^\top\mathbf{X})\)배까지 증폭될 수 있다. 조건수가 \(10^6\)이면 자료의 여섯째 소수점에 있는 교란이 계수 추정값의 첫째 자리를 바꿀 수 있다.

안정화 장치로서의 정칙화¶

능형회귀는 역행렬을 취하기 전에 \(\mathbf{X}^\top\mathbf{X}\)에 \(\lambda\mathbf{I}\)를 더한다.

\(\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I}\)의 고윳값은 \(d_j^2 + \lambda\)이므로 조건수는

가 된다. \(\lambda > 0\)이 고윳값 스펙트럼의 바닥을 올리므로 조건수가 줄어든다. 모든 \(\lambda > 0\)에 대해

이다. \(\lambda \gg d_1^2\)인 극단에서는 조건수가 1(단위행렬)에 접근한다. \(d_p^2\)이 0에 가까울 때는 적당한 \(\lambda\)만으로도 조건수를 여러 자릿수 줄일 수 있다.

왼쪽은 \(p = 8\), \(\rho = 0.999\)인 등상관 행렬의 고윳값 여덟 개를 로그 눈금으로 찍은 것이다. 회색 점(\(\lambda = 0\))을 보면 첫 고윳값만 \(7.993\)이고 나머지 일곱 개가 모두 \(0.001\)에 납작하게 깔려 있다. 조건수 \(7.993/0.001 = 7993\)이 바로 이 절벽의 높이다. 여기에 \(\lambda\)를 더한다는 것은 모든 고윳값을 똑같이 위로 밀어 올린다는 뜻인데, 꼭대기에 더하는 \(\lambda\)는 \(7.993\)에 비해 미미하고 바닥에 더하는 \(\lambda\)는 \(0.001\)에 비해 거대하다. 그래서 같은 \(\lambda\)가 절벽 아래만 메운다.

문제는 \(\lambda\)의 크기다. 파란 사각형(\(\lambda = 0.01\))은 바닥을 \(0.001\)에서 \(0.011\)로 열 배 올려 조건수를 \(7993\)에서 \(728\)로 줄인다. 주황 마름모(\(\lambda = 1\))는 바닥을 \(1.001\)까지 올려 조건수를 \(8.98\)로 만든다. 무려 890배 감소다. 꼭대기 고윳값 \(7.993\)은 \(\lambda = 1\)을 더해도 \(8.993\)으로 거의 그대로라는 점을 눈으로 확인하라. 정칙화가 잘 결정된 방향을 거의 건드리지 않는다는 말의 정확한 의미가 이것이다.

오른쪽은 \(\lambda\)를 연속적으로 바꾸며 조건수를 추적한 것이다. 두 곡선 모두 왼쪽에서는 평평하다. \(\rho = 0.999\) 곡선은 \(\lambda < 0.001\) 구간에서 \(7993\)에 머물고, \(\rho = 0.99\) 곡선은 \(\lambda < 0.01\) 구간에서 \(793\)에 머문다. 세로 점선으로 표시한 각자의 \(\lambda_{\min}\)을 넘어서야 비로소 내려가기 시작한다. \(\lambda\)가 \(\lambda_{\min}\)보다 작으면 편향만 더할 뿐 안정화는 하나도 못 한다. 이 장의 연습문제 3에서 \(\lambda = 0.01\)이 아무 도움이 되지 않았던 이유가 여기 있다.

두 곡선이 \(\lambda \approx 0.3\) 부근에서 만나 이후 겹친다는 점도 읽을 만하다. \(\lambda \gg \lambda_{\min}\)이면 조건수는 \((\lambda_{\max} + \lambda)/\lambda \approx \lambda_{\max}/\lambda\)로, 원래 행렬이 얼마나 나빴는지와 무관해진다. 충분히 큰 \(\lambda\) 앞에서는 모든 불량조건이 똑같아진다. 그러므로 \(\lambda\)를 고를 때 기준으로 삼아야 할 눈금은 \(\lambda_{\min} = d_p^2\)이지 \(1.0\) 같은 기본값이 아니다.

고윳값 스펙트럼과 유효계수¶

\(\mathbf{X}^\top\mathbf{X}\)의 고윳값 스펙트럼이 불량조건의 심각도를 드러낸다. 문턱 \(\tau\)에서의 유효계수를

로 정의한다. 유효계수가 \(p\)보다 훨씬 작으면 계수공간의 많은 방향이 자료로 잘 결정되지 않는다. 능형회귀는 SVD 축소인자 \(d_j^2/(d_j^2 + \lambda)\)가 말하는 대로 이렇게 잘 결정되지 않는 방향을 가장 강하게 축소한다.

실무에서 불량조건 탐지하기¶

분산팽창인자(VIF)는 다중공선성에 관한 다음 절에서 자세히 다룬다.

연습문제¶

import numpy as np
p = 8
for rho in (0.99, 0.999):
    S = rho*np.ones((p, p)) + (1-rho)*np.eye(p)
    ev = np.sort(np.linalg.eigvalsh(S))[::-1]
    for lam in (0.001, 0.01, 0.1, 1.0):
        print(rho, lam, (ev[0]+lam)/(ev[-1]+lam))

0.99 0.001 721.0000000001435
0.99 0.01 397.0000000000435
0.99 0.1 73.00000000000145
0.99 1.0 8.841584158415863
0.999 0.001 3997.0000000013506
0.999 0.01 727.5454545454992
0.999 0.1 80.12871287128766
0.999 1.0 8.984015984015992

import numpy as np
from numpy.linalg import cond, svd
rng = np.random.default_rng(11)
n, p, rho = 40, 8, 0.999
S = rho*np.ones((p, p)) + (1-rho)*np.eye(p)
X = rng.normal(size=(n, p)) @ np.linalg.cholesky(S).T
y = X @ np.ones(p) + rng.normal(0, 1, n)
A = X.T @ X; g = X.T @ y; b0 = np.linalg.solve(A, g)
print(cond(A))          # 16601.7

16601.655795266535

정리하며¶

조건수 \(\kappa(\mathbf{X}^\top\mathbf{X})\)는 OLS 계수가 자료의 작은 교란에 얼마나 민감하게 반응하는지를 잰다. 준선형종속, 고차원, 척도 불일치가 모두 조건수를 부풀려 OLS를 수치적으로 불안정하게 만든다. 능형 정칙화는 \(\mathbf{X}^\top\mathbf{X}\)에 \(\lambda\mathbf{I}\)를 더해 모든 고윳값을 \(\lambda\)만큼 올리고 조건수를 줄인다. 이 안정화 덕분에, 통계적 편향-분산 논증만으로는 정당화되지 않는 상황에서도 정칙화가 값지다.