뉴턴-랩슨과 반복재가중최소제곱¶

반복 알고리즘이 필요한 이유¶

선형회귀에서 최소제곱 추정량은 닫힌 형태의 해 \(\hat{\boldsymbol{\theta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\)를 갖는다. 로지스틱 회귀에는 그런 호사가 없다. 로그가능도가 모수에 대해 비선형이므로 최대가능도 추정치는 반복적 수치 최적화로 찾아야 한다. 이 절에서는 가장 널리 쓰이는 두 알고리즘, 뉴턴-랩슨과 반복재가중최소제곱(IRLS)을 전개한다.

로지스틱 회귀의 로그가능도¶

가능도 절에서 보았듯이, \(y_i \in \{0,1\}\)인 독립 관측치 \((y_i, \mathbf{x}_i)\) \(n\)개에 대한 로그가능도는

이며, \(\boldsymbol{\theta} \in \mathbb{R}^{p}\)는 절편과 기울기 계수를 모은 것이다. 관측치 \(i\)의 예측확률은 \(\hat{p}_i = \sigma(\mathbf{x}_i^\top\boldsymbol{\theta})\)이고 \(\sigma\)는 시그모이드 함수다.

점수방정식¶

점수(로그가능도의 기울기)는

이다. 여기서 \(\mathbf{X}\)는 \(n \times p\) 계획행렬이고 \(\hat{\mathbf{p}} = (\hat{p}_1, \ldots, \hat{p}_n)^\top\)이다. 점수를 0으로 놓으면 \(p\)개의 비선형 방정식계를 얻는다. \(\hat{\mathbf{p}}\)가 시그모이드를 통해 \(\boldsymbol{\theta}\)에 의존하므로 닫힌 형태의 해가 없다.

헤세행렬¶

로그가능도의 이차도함수가 헤세행렬이다.

여기서 \(\mathbf{W} = \operatorname{diag}\bigl(\hat{p}_1(1-\hat{p}_1),\ldots, \hat{p}_n(1-\hat{p}_n)\bigr)\)는 \(n \times n\) 대각 가중행렬이다. \(0 < \hat{p}_i < 1\)인 한 \(\mathbf{W}\)의 대각원소는 모두 양수이므로, \(\mathbf{X}\)가 완전열계수이면 헤세행렬은 음정치다. 그러면 로그가능도가 강오목이므로 국소최대점은 유일한 전역최대점이 된다.

뉴턴-랩슨 알고리즘¶

뉴턴-랩슨은 다음 반복으로 점수방정식의 근을 찾는다.

위에서 유도한 점수와 헤세행렬을 대입하면

이며, \(\mathbf{W}^{(t)}\)와 \(\hat{\mathbf{p}}^{(t)}\)는 현재 반복값 \(\boldsymbol{\theta}^{(t)}\)에서 평가한 것이다.

수렴¶

로그가능도가 강오목이므로 뉴턴-랩슨은 어떤 출발점에서든 유일한 MLE로 수렴한다. 최적점 근처에서 수렴은 이차적이다. 즉 반복마다 정확한 자릿수가 대략 두 배가 된다. 실무에서는 3회에서 8회면 충분한 경우가 많다.

피셔 점수법과 IRLS¶

피셔 점수법은 관측 헤세행렬을 그 기댓값으로 대체한다. 로지스틱 회귀에서는 두 양이 일치한다.

여기서 \(\mathcal{I}(\boldsymbol{\theta})\)는 피셔 정보행렬이다. 따라서 로지스틱 회귀에서는 피셔 점수법과 뉴턴-랩슨이 같은 반복열을 만든다.

IRLS 해석¶

뉴턴-랩슨 갱신식을 다시 쓰면

이고, 작업 반응변수는

이다. 각 반복은 가중치 \(\mathbf{W}^{(t)}\)와 반응변수 \(\mathbf{z}^{(t)}\)를 갖는 가중최소제곱 문제를 푼다. 가중치가 단계마다 바뀌므로 이 절차를 반복재가중최소제곱이라 부른다.

알고리즘 요약¶

분리와 비수렴¶

두 범주가 완전히 분리되어 있으면 — 특성공간의 어떤 초평면이 모든 훈련점을 정확히 분류하면 — MLE가 존재하지 않는다. \(\lVert\boldsymbol{\theta}\rVert \to \infty\)일 때 로그가능도가 상한에 접근할 뿐이고, 뉴턴-랩슨은 발산한다. 계수가 한없이 커지고 표준오차가 매우 커지는 것이 증상이다.

경사하강과의 관계¶

경사하강은 고정된 보폭 \(\eta\)로 \(\boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + \eta\,\mathbf{s}(\boldsymbol{\theta}^{(t)})\) 를 갱신한다. 뉴턴-랩슨은 헤세행렬의 역이 정해 주는 적응적 보폭을 쓰는 경사하강으로 볼 수 있다. 이 곡률 정보 덕분에 빠르게 수렴하지만, 매 단계 \(p \times p\) 헤세행렬을 만들고 역을 구해야 하는 비용이 든다. \(p\)가 큰 고차원 문제에서는 L-BFGS 같은 준뉴턴법이 헤세행렬을 근사하여 계산 비용을 줄인다.

두 알고리즘을 나란히 보기¶

"이차 수렴"과 "발산"은 말로만 들으면 정도의 차이처럼 들린다. 실제로 재 보면 차원이 다른 이야기다. 아래 왼쪽은 19.1절의 공부 시간 자료 300건에서 최적값까지 남은 음의 로그가능도를 반복 횟수에 대해 그린 것이고(두 축 모두 로그 눈금), 오른쪽은 이 절 앞의 표에 나온 네 점짜리 분리 자료에서 반복이 만들어 내는 확률 곡선이다.

왼쪽에서 뉴턴-랩슨의 남은 거리는 \(84.5 \to 7.91 \to 0.576 \to 0.0059 \to 7.8\times10^{-7} \to 1.4\times10^{-14}\)로 줄어든다. 정확한 자릿수가 반복마다 대략 두 배씩 늘어나는 이차 수렴이 눈에 그대로 보인다. 여섯 걸음이면 배정도 부동소수점의 한계에 닿는다. 같은 자료에서 학습률 \(\alpha=0.001\)의 경사하강은 같은 정확도에 이르는 데 약 2,900회, \(\alpha=0.0001\)이면 2만 회를 넘어도 \(10^{-9}\) 근처에 머문다. 곡률 정보를 쓰느냐 마느냐가 세 자릿수의 차이를 만든다.

한 걸음의 비용은 다르다. 뉴턴 단계는 \(p\times p\) 행렬을 만들고 풀어야 하므로 \(O(np^2 + p^3)\)이고, 경사하강 한 걸음은 \(O(np)\)다. 설명변수가 두어 개인 문제에서는 뉴턴이 압도적이지만, \(p\)가 수천인 문제에서는 \(p^3\)이 감당되지 않아 L-BFGS나 확률적 경사하강으로 갈아타게 된다. 그 선택은 수렴 속도가 아니라 한 걸음의 가격이 정한다.

오른쪽은 같은 알고리즘이 멈추지 못하는 경우다. 자료가 \(x=2.5\)를 기준으로 완전히 갈리므로 반복이 진행될수록 곡선이 점점 가팔라진다. \(t=1\)에서 \(\theta_1=-1.60\)이던 기울기가 \(t=5\)에서 \(-8.63\), \(t=20\)에서 \(-38.67\)이 되고, 곡선은 계단함수로 수렴해 간다. 로그가능도는 계속 올라가지만 상한 0에 도달하지는 못하므로 최대점이 존재하지 않는다. 알고리즘이 고장 난 것이 아니라, 찾으라고 시킨 것이 없는 것이다. 수렴 판정이 계수의 변화량에 걸려 있다면 이 반복은 영원히 멈추지 않는다.

연습문제¶

import numpy as np

x = np.array([1., 2., 3., 4., 5., 6.])
y = np.array([0., 0., 1., 0., 1., 1.])
X = np.column_stack([np.ones(6), x])

theta = np.zeros(2)
history = [theta.copy()]
for t in range(10):
    p_hat = 1 / (1 + np.exp(-X @ theta))
    W = p_hat * (1 - p_hat)
    H = X.T @ (X * W[:, None])          # X^T W X
    s = X.T @ (y - p_hat)               # score
    theta = theta + np.linalg.solve(H, s)
    history.append(theta.copy())

theta_hat = theta
for t, th in enumerate(history):
    print(f"{t}  {th}  err = {np.linalg.norm(th - theta_hat):.3e}")

0  [0. 0.]  err = 4.419e+00
1  [-2.8  0.8]  err = 1.507e+00
2  [-3.88421522  1.10977578]  err = 3.795e-01
3  [-4.22143618  1.20612462]  err = 2.877e-02
4  [-4.24892948  1.21397985]  err = 1.738e-04
5  [-4.24909654  1.21402758]  err = 6.365e-09
6  [-4.24909655  1.21402759]  err = 9.155e-16
7  [-4.24909655  1.21402759]  err = 2.747e-15
8  [-4.24909655  1.21402759]  err = 0.000e+00
9  [-4.24909655  1.21402759]  err = 2.747e-15
10  [-4.24909655  1.21402759]  err = 0.000e+00

정리하며¶

반복 알고리즘으로 해를 찾는다.

• 뉴턴–랩슨이 표준이다. 헤세행렬을 쓰므로 이차 수렴하며, 오목한 로그가능도에서 빠르고 안정적이다.
• IRLS 가 그 다른 이름이다. 뉴턴 갱신을 다시 쓰면 가중최소제곱을 반복하는 형태가 되며, 가중치가 \(p_i(1-p_i)\) 다. 13장의 WLS 와 직접 이어진다.
• 기울기 하강은 느리지만 확장성이 좋다. \(p\) 가 아주 크면 헤세행렬 계산이 부담이므로 확률적 기울기 하강 계열을 쓴다.
• 수렴하지 않으면 분리를 의심한다. 계수가 계속 커지고 표준오차가 폭발하면 완전분리나 준분리이며, 정칙화나 퍼스 보정이 처방이다.
• 초기값은 대개 \(\mathbf 0\) 으로 충분하다. 오목하므로 어디서 출발해도 같은 곳에 도달한다.

다음 절 검정으로 넘어간다.