Klein Appendix A.2 — Multivariate Methods for Maximization

1 도입 — 다변수 최적화가 어디에 등장하는가

생존분석 본문에서 다변수 최적화가 등장하는 자리:

위치	모수 차원	형태
Ch.8 — Cox PH (다공변량)	\(\beta \in \mathbb{R}^p\), \(p \geq 2\)	부분우도 Newton-Raphson
Ch.9 — 시간의존 공변량 + 층화	$+ $ 층 모수	동일 알고리즘 확장
Ch.10 Aalen	시점별 OLS — closed form	다변수지만 closed form
Ch.10 Lin-Ying	\(\alpha \in \mathbb{R}^p\) — closed form	다변수지만 closed form
Ch.12 모수 회귀 (Weibull/LL/LN)	\((\mu, \sigma, \gamma) \in \mathbb{R}^{p+2}\)	full 우도 NR
Ch.12 일반화 감마	\((\mu, \sigma, \gamma, \theta) \in \mathbb{R}^{p+3}\)	NR 또는 LM
Ch.13 Frailty MLE (모수적 baseline)	\((\theta, \beta, \mu, \sigma) \in \mathbb{R}^{p+3}\)	다변수 NR
Ch.13 Frailty EM 의 분산 추정	\(D + p + 1\) 차	정보 행렬 역수

§ A.1 (1D) 와 § A.2 (다변수) 의 결정적 차이

1D 영점 찾기: \(f'(x) = 0\) 의 해 — 한 차원의 값.

다변수 영점 찾기: \(\nabla f(\mathbf{x}) = \mathbf{0}\) 의 해 — \(p\) 차원 벡터의 모든 성분이 0 인 점.

왜 다변수가 더 어려운가:

1. 방향 선택: 1D 는 +/- 두 방향만. 다변수는 \(p\) 차원 단위 구면의 모든 방향.
2. Hessian 의 역할: 1D 의 \(f''\) 는 스칼라. 다변수의 \(\mathbf{H}\) 는 \(p \times p\) 행렬 — 역수 비용 \(O(p^3)\).
3. 길쭉한 우도: 한 축의 곡률이 다른 축보다 훨씬 큰 경우 (예: Weibull 의 \(\lambda\) 는 평탄, \(\alpha\) 는 가파름). 단순 그래디언트 방향이 비효율.
4. 수렴 판정: 단일 기준이 부족, 4 가지 (함수값/gradient 노름/모수 변화/반복 한계) 조합.

이 절은 다변수 최적화의 표준 도구 3 가지 (steepest ascent, Newton-Raphson 다변수, Marquardt) 를 알고리즘과 직관으로 다룬다.

2 표기 — Gradient 와 Hessian

2.1 Gradient (1차 도함수 벡터) — 식 (A.3)

\(f: \mathbb{R}^p \to \mathbb{R}\) 의 점 \(\mathbf{x}\) 에서의 그래디언트:

\[ \mathbf{u}(\mathbf{x}) = [u_1(\mathbf{x}), \ldots, u_p(\mathbf{x})]^t , \quad u_j(\mathbf{x}) = \frac{\partial f(\mathbf{x})}{\partial x_j} . \tag{식 A.3} \]

2.2 Hessian (2차 도함수 행렬) — 식 (A.4)

\[ \mathbf{H}(\mathbf{x}) = (H_{ij}(\mathbf{x})) , \quad H_{ij}(\mathbf{x}) = \frac{\partial^2 f(\mathbf{x})}{\partial x_i \partial x_j} . \tag{식 A.4} \]

대칭 행렬 (\(H_{ij} = H_{ji}\), Schwarz 정리). 최대값에서는 음의 정칙 (negative definite, 모든 eigenvalue < 0).

직관 — Gradient 와 Hessian 의 기하적 의미

Gradient \(\mathbf{u}(\mathbf{x})\): 점 \(\mathbf{x}\) 에서 \(f\) 가 가장 빠르게 증가하는 방향을 가리키는 벡터. 길이는 그 증가율. 영점 \(\mathbf{u}(\mathbf{x}^*) = \mathbf{0}\) 인 점이 임계점 (최대/최소/안장점).

Hessian \(\mathbf{H}(\mathbf{x})\): 점 \(\mathbf{x}\) 의 곡률 텐서. \(\mathbf{H}\) 의 eigenvector 들이 우도의 주축 (principal axes), eigenvalue 들이 각 축의 곡률.

최대값의 기하: 모든 방향으로 곡률 음 (\(\mathbf{H}\) 음의 정칙) → 우도가 그 점에서 모든 방향으로 감소 → 진정한 최대.

길쭉한 우도: \(\mathbf{H}\) 의 eigenvalue 들이 크게 다르면 (예: 100 vs 1) 우도 등고선이 길쭉한 타원. 한 축으로는 가파르고 다른 축으로는 평탄.

이 길쭉함이 다변수 최적화의 핵심 도전 — steepest ascent 가 비효율적이 되는 원인.

2.3 Quadratic Form 의 의미

\(f\) 의 점 \(\mathbf{x}_0\) 근처 Taylor 2 차 근사:

\[ f(\mathbf{x}) \approx f(\mathbf{x}_0) + \mathbf{u}(\mathbf{x}_0)^t (\mathbf{x} - \mathbf{x}_0) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_0)^t \mathbf{H}(\mathbf{x}_0) (\mathbf{x} - \mathbf{x}_0) . \]

1차 항: 평면 (linear). 2차 항: \(\mathbf{H}\) 가 음의 정칙이면 위로 볼록 (concave) 한 quadratic. 이 quadratic 의 최대점이 Newton 의 한 step 의 도착지.

3 § A.2.1 Steepest Ascent — 식 (A.5)

3.1 알고리즘

각 step 에서 그래디언트 방향으로 직선 탐색 (line search):

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + d_k \, \mathbf{u}(\mathbf{x}_k) , \tag{식 A.5} \]

\(d_k\) 는 그 방향에서 \(f\) 를 최대화하는 step size. 즉:

\[ d_k = \arg\max_{d > 0} f\big(\mathbf{x}_k + d \, \mathbf{u}(\mathbf{x}_k)\big) . \]

이 1D 최대화는 § A.1 의 도구 (bisection, secant, Newton) 로 풀어낸다.

3.1.1 의사코드

INPUT: f, ∇f (gradient), x_0, tol, max_iter
x ← x_0
FOR k = 0, 1, ..., max_iter:
    g ← ∇f(x)                                   # gradient
    IF ||g|| < tol: BREAK                        # 수렴
    # 1D line search: argmax_d f(x + d·g)
    d_star ← line_search(d → f(x + d·g))         # § A.1 도구 사용
    x ← x + d_star · g
RETURN x

직관 — “현재 점에서 가장 가파른 오르막” 만 사용

각 step 에서 gradient 방향만 본다. Hessian 정보는 사용 안 함.

왜 gradient 방향이 가장 가파른가: 단위 벡터 \(\mathbf{v}\) 방향의 \(f\) 증가율은 \(\mathbf{u}^t \mathbf{v}\). 코시-슈바르츠 부등식으로 이 양은 \(\|\mathbf{u}\|\) 에서 최대 (\(\mathbf{v} = \mathbf{u}/\|\mathbf{u}\|\) 일 때). 따라서 gradient 방향이 가장 가파른 오르막.

Step size \(d_k\): 그 방향에서 무한히 멀리 갈 게 아니라 최적 거리 까지만. 1D 최적화로 결정.

비유: 안개 낀 산에서 손으로 더듬어 가장 가파른 방향을 느끼고 (gradient), 그 방향으로 균형 잃지 않을 만큼 (line search) 한 발 디딤. 안개라 멀리 못 봐도 매번 오르막은 올라간다.

3.2 길쭉한 우도와 지그재그 — Steepest Ascent 의 한계

직관 — 왜 길쭉한 우도에서 지그재그가 생기는가

우도 함수가 길쭉한 타원 형태라고 상상하자 (한 축의 곡률이 다른 축의 100 배). 다음 두 사례를 본다.

사례 1: 원형 우도 (Hessian 의 eigenvalue 모두 같음). 모든 방향의 곡률 동일. - 점 \(\mathbf{x}_0\) 의 gradient 는 정확히 최대점 방향을 가리킴. - 한 step 의 line search 로 최대점 도달 (1 step 수렴).

사례 2: 길쭉한 우도 (한 축 곡률 큰 경우). 등고선이 가늘고 긴 타원. - 점 \(\mathbf{x}_0\) 의 gradient 는 등고선에 수직 방향이지만, 이는 최대점 방향이 아니다 (가파른 축 쪽으로 치우침). - Line search 로 그 방향의 최대점에 도달 — 타원의 짧은 축의 다른 끝 근처에 도달. - 새 점에서의 gradient 는 또 등고선에 수직, 이전 방향과 거의 직교. - 이렇게 직교 방향으로 번갈아 지그재그 진행. - 길쭉할수록 지그재그가 심해 수렴이 매우 느림 (\(O(1/k)\) 수렴 — bisection 보다도 느릴 수 있음).

이는 Rosenbrock banana 함수 등에서 극단적으로 보이는 현상 — Steepest ascent 는 이런 함수에서 수백~수천 step 이 필요할 수 있다.

해결: Newton-Raphson 의 \(\mathbf{H}^{-1}\) 정규화 — 길쭉한 우도를 둥글게 펴서 본 후 진행.

3.3 Steepest Ascent 의 사용 시점

조건	적합성
Hessian 계산 불가	O (BFGS 같은 quasi-Newton 도 후보)
우도 단봉 + 거의 원형	O
큰 Hessian (\(p \geq 1000\))	O (\(\mathbf{H}^{-1}\) 비용 \(O(p^3)\) 회피)
길쭉한 우도	X (지그재그)
좋은 시작값 + 빠른 수렴 원할 때	X (Newton 우위)

머신러닝의 SGD 와의 관계

신경망 훈련의 확률적 경사 하강 (SGD) 가 사실상 steepest descent (최소화 버전). \(p\) 가 매우 큼 (수백만~수십억) 이라 Hessian 비현실 → gradient 기반만 사용.

생존분석은 보통 \(p\) 가 작아 (\(p < 100\)) Newton 이 우위. ML 의 high-dim 문제와 정반대 영역. § A.2 의 이론은 둘 다 같지만 실무 default 가 다르다.

4 § A.2.2 Newton-Raphson 다변수 — 식 (A.6)

4.1 알고리즘

\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{H}(\mathbf{x}_k)^{-1} \mathbf{u}(\mathbf{x}_k) . \tag{식 A.6} \]

4.1.1 의사코드

INPUT: f, ∇f, ∇²f (Hessian), x_0, tol, max_iter
x ← x_0
FOR k = 0, 1, ..., max_iter:
    g ← ∇f(x)
    IF ||g|| < tol: BREAK
    H ← ∇²f(x)
    Δ ← solve(H, g)                               # H · Δ = g
    x_new ← x - Δ
    IF f(x_new) < f(x):                           # step-halving
        Δ ← Δ / 2
        x_new ← x - Δ
    x ← x_new
RETURN x

(solve(H, g) 은 LU 분해 등으로 선형 system 푸는 것 — explicit \(\mathbf{H}^{-1}\) 계산 회피, 더 안정.)

직관 — Hessian 으로 좌표를 정규화한 후 Steepest Ascent

식 (A.6) 의 \(-\mathbf{H}^{-1} \mathbf{u}\) 는 다음과 동치:

1. \(\mathbf{y} = \mathbf{H}^{1/2} (\mathbf{x} - \mathbf{x}^*)\) 로 좌표 변환 (\(\mathbf{H}^{1/2}\) 는 양의 정칙 부분).
2. 이 새 좌표에서 우도가 등방성 (isotropic) — 모든 방향의 곡률 같음.
3. 등방성 좌표에서 steepest ascent 한 step.
4. 원래 좌표로 역변환.

즉 길쭉한 우도를 둥글게 편 다음 steepest ascent. 이것이 Newton 의 본질.

기하적: 우도의 2차 Taylor 근사 \(f(\mathbf{x}_0) + \mathbf{u}^t \Delta + \frac{1}{2} \Delta^t \mathbf{H} \Delta\) 의 최대점이 정확히 \(\Delta = -\mathbf{H}^{-1} \mathbf{u}\). 즉 2차 근사의 최대점으로 한 번에 점프.

우도가 정확히 quadratic 이면 1 step 으로 도달. 일반 우도는 quadratic 근사가 점점 정확해져 quadratic 수렴.

4.2 2차 수렴

Quadratic 수렴의 정량적 의미

오차 \(\mathbf{e}_k = \mathbf{x}_k - \hat{\mathbf{x}}\) 에 대해:

\[ \|\mathbf{e}_{k+1}\| \leq C \|\mathbf{e}_k\|^2 , \]

\(C\) 는 \(\mathbf{H}\) 의 third-order 도함수에 의존하는 상수.

Step	\(\|\mathbf{e}_k\|\)
0	\(10^{-2}\)
1	\(10^{-4}\) (가정 \(C \approx 1\))
2	\(10^{-8}\)
3	\(10^{-16}\) (기계 정밀도)

3-4 step 만에 16 자리 정밀도 — 단변수 Newton 과 동일한 폭발적 가속. 이것이 Newton 의 가장 큰 매력.

4.3 발산 사례 — Newton 의 약점

다변수 Newton 이 실패하는 경우

1. \(\mathbf{H}\) 가 ill-conditioned 또는 singular: Eigenvalue 중 하나가 0 에 매우 가까우면 \(\mathbf{H}^{-1}\) 폭발. 모수 식별 불가능 또는 약한 식별성.

2. \(\mathbf{H}\) 가 양의 정칙이 아님: 안장점 또는 최소값 근처. Newton 은 그런 점에서도 영점으로 인식해 잘못된 방향 진행.

3. Overshoot: \(\mathbf{x}_k\) 가 영점에서 멀고 quadratic 근사가 부정확하면 \(\|\Delta\|\) 가 커 영점 너머로 점프 → 발산.

4. 잘못된 방향 — \(f\) 감소: \(f(\mathbf{x}_{k+1}) < f(\mathbf{x}_k)\) → 최대화에 반대 방향.

대처 — Step-Halving (SAS Cox 의 표준):

Δ ← H⁻¹ u
x_new ← x - Δ
WHILE f(x_new) < f(x) AND ||Δ|| > eps:
    Δ ← Δ / 2
    x_new ← x - Δ

\(f\) 가 증가할 때까지 step 을 절반씩 줄임. 작은 step 이라도 항상 증가 보장.

대안 — Damped Newton 또는 trust region: Step 크기를 사전 제한 (\(\|\Delta\| \leq \delta_k\)). Trust region 은 modern 표준.

5 § A.2.3 Marquardt’s Method — Steepest 와 Newton 의 절충

5.1 알고리즘

\(\mathbf{S}_k\) 를 대각 스케일 행렬 (\(S_{ii} = |H_{ii}(\mathbf{x}_k)|^{-1/2}\)). \(\gamma \geq 0\) 을 blending 상수.

\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{S}_k \big(\mathbf{S}_k \mathbf{H}(\mathbf{x}_k) \mathbf{S}_k + \gamma \mathbf{I}\big)^{-1} \mathbf{S}_k \mathbf{u}(\mathbf{x}_k) . \]

직관 — \(\gamma\) 가 두 방법을 잇는 다이얼

\(\gamma = 0\): \((\mathbf{S}_k \mathbf{H} \mathbf{S}_k)^{-1}\) 만 남음 → 정규화된 Newton-Raphson (스케일 행렬은 단순 normalization).

\(\gamma \to \infty\): 분모의 \(\gamma \mathbf{I}\) 가 지배 → \(\mathbf{S}_k \cdot \frac{\mathbf{S}_k \mathbf{u}}{\gamma}\) → steepest ascent (Hessian 영향 사라짐).

중간 \(\gamma\): 둘의 보간. \(\mathbf{H}\) 가 양의 정칙이 아닐 때도 \(\mathbf{S}_k \mathbf{H} \mathbf{S}_k + \gamma \mathbf{I}\) 가 양의 정칙이 되어 안정.

핵심 통찰 (Marquardt 1963): “\(\mathbf{H}\) 가 신뢰할 만하면 Newton 쪽 (\(\gamma\) 작게), 못 믿겠으면 steepest 쪽 (\(\gamma\) 크게).”

5.2 적응 절차 (Adaptive)

INPUT: f, ∇f, ∇²f, x_0, γ_0 = 0.001, tol
x ← x_0
γ ← γ_0
FOR k = 0, 1, ...:
    g ← ∇f(x); H ← ∇²f(x); S ← diag(|H_ii|^(-1/2))
    M ← S·H·S + γ·I
    Δ ← S · solve(M, S·g)
    x_new ← x - Δ
    IF f(x_new) > f(x):                        # 진전
        x ← x_new
        γ ← γ / 10                             # 다음에 더 Newton 쪽
    ELSE:                                       # 후퇴
        γ ← γ × 10                             # 더 steepest 쪽으로
        # x 는 그대로, γ 만 늘려 재시도
    IF ||g|| < tol AND γ small: BREAK
    # 마지막 step 은 γ = 0 (순수 Newton) 으로 정확 마무리

비유 — 운전 경험에 따라 액셀 강도 조절

• 직선 도로 (Newton 신뢰): \(\gamma\) 작게 → 액셀 풀로 (빠름).
• 굽은 도로 (Newton 의심): \(\gamma\) 크게 → 안전 운전 (steepest, 느리지만 안전).

이 적응이 단일 method 보다 거의 항상 빠름 — 처음에는 안전, 끝에서는 빠르게.

LM 의 다른 응용: 비선형 최소제곱 문제 (nls() in R, scipy.optimize.least_squares(method='lm')) 의 표준. 회귀 잔차 제곱합 최소화에서는 Hessian 이 이론상 양의 정칙이지만 작은 표본 / 모형 오설정에서 문제 발생 → LM 의 견고성이 결정적.

6 수렴 판정 — 4 가지 기준

다변수 최적화에서는 단일 기준이 부족, 여러 기준 조합 사용:

수렴 판정 4 가지

1. 함수값 변화 (\(\epsilon = 10^{-6}\) 권장):

\[ |f(\mathbf{x}_{k+1}) - f(\mathbf{x}_k)| < \epsilon \quad \text{또는} \quad \frac{|f(\mathbf{x}_{k+1}) - f(\mathbf{x}_k)|}{|f(\mathbf{x}_k)|} < \epsilon . \]

2. Gradient 노름:

\[ \sum_j u_j(\mathbf{x}_{k+1})^2 < \epsilon \quad \text{또는} \quad \max_j |u_j(\mathbf{x}_{k+1})| < \epsilon . \]

3. 모수 변화:

\[ \sum_j (x_{k+1, j} - x_{k, j})^2 < \epsilon \quad \text{또는} \quad \max_j |x_{k+1, j} - x_{k, j}| < \epsilon . \]

4. 반복 한계 (안전장치):

\[ k \geq k_{\max} \quad (\text{보통 } 100 \text{ 또는 } 500) . \]

실무: (1) + (2) 둘 다 점검. (3) 은 보조 — 함수가 평탄하면 모수가 변해도 진정한 최적이 아닐 수 있음. (4) 는 발산 방지.

기준 충돌 사례: \(\|\mathbf{u}\|\) 작은데 \(\|\Delta \mathbf{x}\|\) 큼 → 평탄 영역에서 조심해야 함. 또는 \(|\Delta f|\) 작은데 \(\|\mathbf{u}\|\) 큼 → quadratic 근사 부정확.

7 Example A.2 — Weibull (λ, α) MLE

7.1 데이터와 모형

§ A.1 의 Example A.1 과 같은 10 obs:

t = (2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37)

\(n = 10\), \(\sum t_j \approx 9.76\).

이번엔 두 모수 \(S(t) = \exp(-\lambda t^\alpha)\) 모두 추정.

7.1.1 Log-likelihood, Score, Hessian

\[ \ell(\lambda, \alpha) = n \ln \lambda + n \ln \alpha + (\alpha - 1) \sum \ln t_j - \lambda \sum t_j^\alpha . \]

\[ \mathbf{u}(\lambda, \alpha) = \begin{pmatrix} u_\lambda \\ u_\alpha \end{pmatrix} = \begin{pmatrix} n/\lambda - \sum t_j^\alpha \\ n/\alpha + \sum \ln t_j - \lambda \sum t_j^\alpha \ln t_j \end{pmatrix} . \]

\[ \mathbf{H}(\lambda, \alpha) = \begin{pmatrix} -n/\lambda^2 & -\sum t_j^\alpha \ln t_j \\ -\sum t_j^\alpha \ln t_j & -n/\alpha^2 - \lambda \sum t_j^\alpha (\ln t_j)^2 \end{pmatrix} . \]

7.1.2 시작값

Method-of-moments 추정: \(\alpha_0 = 1\) (지수분포 가정), \(\lambda_0 = 10/\sum t_j \approx 1.024\). 초기 log-likelihood \(\ell(\lambda_0, \alpha_0) = -9.757\).

정지 기준: \(\max(|u_\lambda|, |u_\alpha|) < 0.1\).

7.2 Steepest Ascent — 5 step Trace

Step \(k\)	\(\lambda_k\)	\(\alpha_k\)	\(\ell\)	\(u_\lambda\)	\(u_\alpha\)	\(d_k\)
0	1.024	1.000	\(-9.757\)	0.001	7.035	0.098
1	1.025	1.693	\(-7.491\)	0.001	\(-1.80\)	0.089
2	0.865	1.694	\(-7.339\)	0.661	0.001	0.126
3	0.865	1.777	\(-7.311\)	0.000	\(-0.363\)	0.073
4	0.839	1.777	\(-7.307\)	0.121	0.000	0.128
5	0.839	1.792	\(-7.306\)	0.000	\(-0.072\)	0.007

5 step 후 \((\hat{\lambda}, \hat{\alpha}) = (0.839, 1.792)\).

Trace 의 패턴 — 지그재그 명확

각 step 의 변화:

• Step 0 → 1: \(\alpha\) 만 거의 변화 (1.000 → 1.693). \(\lambda\) 거의 그대로 (1.024 → 1.025). \(u_\alpha = 7.035\) 가 \(u_\lambda = 0.001\) 의 압도적으로 크므로 gradient 가 거의 \(\alpha\) 방향.
• Step 1 → 2: \(\lambda\) 만 거의 변화 (1.025 → 0.865). \(\alpha\) 거의 그대로 (1.693 → 1.694). 새 gradient 가 거의 \(\lambda\) 방향.
• Step 2 → 3: \(\alpha\) 변화 (1.694 → 1.777). \(\lambda\) 그대로.
• Step 3 → 4: \(\lambda\) 변화. \(\alpha\) 그대로.
• Step 4 → 5: \(\alpha\) 변화. \(\lambda\) 그대로.

즉 각 step 이 한 모수씩 번갈아 큰 변화 — 이것이 길쭉한 우도 함수의 지그재그 패턴이다.

수렴은 결국 도달하지만 step 수가 많고 (5 step) 진행이 비효율적. Newton 의 Hessian 정규화가 이 비효율을 해결한다 (다음 절).

7.3 Newton-Raphson 다변수 — 3 step Trace

같은 시작값 + 정지 기준.

Step \(k\)	\(\lambda_k\)	\(\alpha_k\)	\(u_\lambda\)	\(u_\alpha\)	\(H_{\lambda\lambda}\)	\(H_{\alpha\alpha}\)	\(H_{\lambda\alpha}\)
0	1.024	1.000	0.001	7.035	\(-9.537\)	\(-13.449\)	\(-1.270\)
1	0.954	1.530	\(-0.471\)	1.684	\(-10.987\)	\(-8.657\)	\(-3.34\)
2	0.838	1.769	0.035	0.181	\(-14.223\)	\(-7.783\)	\(-1.220\)
3	0.832	1.796	\(-0.001\)	0.001	\(-14.431\)	\(-7.750\)	\(-4.539\)

3 step 후 \((\hat{\lambda}, \hat{\alpha}) = (0.832, 1.796)\).

Newton 의 첫 step 분석

Step 0 → 1: \(\Delta = -\mathbf{H}^{-1} \mathbf{u}\). Hessian:

\[ \mathbf{H} = \begin{pmatrix} -9.537 & -1.270 \\ -1.270 & -13.449 \end{pmatrix}, \quad \det(\mathbf{H}) = 9.537 \times 13.449 - 1.270^2 = 126.65 . \]

\(\mathbf{H}^{-1} = \frac{1}{126.65} \begin{pmatrix} -13.449 & 1.270 \\ 1.270 & -9.537 \end{pmatrix}\).

\(\Delta = -\mathbf{H}^{-1} \mathbf{u} = -\frac{1}{126.65} \begin{pmatrix} -13.449 \cdot 0.001 + 1.270 \cdot 7.035 \\ 1.270 \cdot 0.001 + (-9.537) \cdot 7.035 \end{pmatrix} \approx \begin{pmatrix} -0.07 \\ 0.530 \end{pmatrix}\).

따라서 새 점: \((1.024 - 0.07, 1.000 + 0.530) = (0.954, 1.530)\) — 표 첫 행과 일치.

Step 1 의 \(\alpha\) 변화 \(1.000 → 1.530\) 은 steepest 의 \(1.000 → 1.693\) 보다 적다 (overshoot 회피). 그러나 동시에 \(\lambda\) 도 적절히 줄임 (\(1.024 → 0.954\)) — 두 모수를 동시에 합리적으로 update, steepest 의 한 모수씩 변화 패턴 회피.

이것이 Hessian 의 cross-term (\(H_{\lambda\alpha}\)) 이 만드는 차이 — 두 모수의 상관을 고려해 동시 update.

7.4 Marquardt — 5 step Trace (\(\gamma = 0.5\))

Step \(k\)	\(\lambda_k\)	\(\alpha_k\)	\(u_\lambda\)	\(u_\alpha\)	\(H_{\lambda\lambda}\)	\(H_{\alpha\alpha}\)	\(H_{\lambda\alpha}\)
0	1.024	1.000	0.001	7.035	\(-9.537\)	\(-13.449\)	\(-1.270\)
1	0.993	1.351	\(-0.357\)	3.189	\(-10.136\)	\(-9.534\)	\(-2.565\)
2	0.930	1.585	\(-0.394\)	1.295	\(-11.557\)	\(-8.424\)	\(-3.599\)
3	0.883	1.701	\(-0.275\)	0.523	\(-12.813\)	\(-8.049\)	\(-4.176\)
4	0.858	1.753	\(-0.162\)	0.218	\(-13.591\)	\(-7.891\)	\(-4.453\)
5	0.845	1.777	\(-0.087\)	0.094	\(-14.013\)	\(-7.817\)	\(-4.581\)

5 step 후 \((\hat{\lambda}, \hat{\alpha}) = (0.845, 1.777)\).

Marquardt 의 보수적 진행

\(\gamma = 0.5\) 의 고정값 (적응 안 함). 효과:

• Newton 의 첫 step (1.024, 1.000) → (0.954, 1.530) 보다 더 보수적 — Marquardt 는 (0.993, 1.351) 로 절반쯤만 진행.
• 매 step 의 변화도 작음 — 안전하게 천천히 도달.
• 결국 5 step 필요 (Newton 의 3 step 보다 많음).

\(\gamma\) 적응이 없을 때의 Marquardt 가 보여주는 안전 운전 패턴. 적응 (\(\gamma\) 자동 조절) 을 사용하면 끝 step 에서 \(\gamma \to 0\) 으로 Newton 과 동등해져 더 빠르다.

7.5 세 방법 종합 비교 — Example A.2

방법	Step	\((\hat{\lambda}, \hat{\alpha})\)	마지막 \(\max(|u_\lambda|, |u_\alpha|)\)
Steepest Ascent	5	\((0.839, 1.792)\)	0.072
Newton-Raphson	3	\((0.832, 1.796)\)	0.001
Marquardt (\(\gamma = 0.5\))	5	\((0.845, 1.777)\)	0.094

통찰 — 세 방법의 trade-off 가 명확히 보임

Newton 의 우월성: - Step 수: 3 (Marquardt/Steepest 의 5 보다 적음). - 마지막 정확도: \(|u| = 0.001\) — Marquardt/Steepest 의 0.07-0.09 보다 100 배 정확. - 점추정: \((0.832, 1.796)\) — 다른 방법과 약간 다름 (정지 기준 안에서의 우연 차이).

Steepest 의 비효율: - 한 모수씩 번갈아 변화 (지그재그). - 같은 step 수로 Marquardt 와 비슷.

Marquardt 의 안전: - 보수적 진행, 중간 step 의 안정성 좋음. - 적응 (\(\gamma\) 조절) 추가하면 Newton 과 동등 속도까지 가능.

실무 권고: - 단일 우도 함수, 좋은 시작값 → Newton. - 복잡한 모형, 시작값 불안정 → Marquardt 적응. - 매우 큰 \(p\), Hessian 비싼 경우 → BFGS / L-BFGS (다음 절).

8 Quasi-Newton (BFGS) — Klein 미언급의 표준

현대 표준 — Quasi-Newton

Klein 부록 A 가 작성된 시점 (2003) 부터 표준 도구는 거의 변하지 않았지만, 한 가지 빠진 부분이 Quasi-Newton.

아이디어: Hessian 명시적 계산 대신 gradient 변화 로 Hessian 의 근사 update.

BFGS 갱신 공식:

\[ \mathbf{B}_{k+1} = \mathbf{B}_k + \frac{\mathbf{y}_k \mathbf{y}_k^t}{\mathbf{y}_k^t \mathbf{s}_k} - \frac{\mathbf{B}_k \mathbf{s}_k \mathbf{s}_k^t \mathbf{B}_k}{\mathbf{s}_k^t \mathbf{B}_k \mathbf{s}_k} , \]

여기서 \(\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k\), \(\mathbf{y}_k = \mathbf{u}(\mathbf{x}_{k+1}) - \mathbf{u}(\mathbf{x}_k)\).

장점: - Hessian 명시적 계산 불필요 → \(f''\) 비싼 우도에 적합. - Super-linear 수렴 (Newton 의 quadratic 보다 약간 느림). - \(\mathbf{B}_k\) 가 자동으로 양의 정칙 유지 → 안전.

단점: - \(p \times p\) 행렬 저장 — \(p\) 매우 클 때 (예: 신경망의 수백만 모수) 부담. - 해결: L-BFGS (limited-memory BFGS) — 최근 \(m\) step 만 저장. 메모리 \(O(mp)\).

구현: - R: optim(method = "BFGS") 또는 optim(method = "L-BFGS-B") (with bounds). - Python: scipy.optimize.minimize(method='BFGS') 또는 'L-BFGS-B'.

실무 권고: \(p \leq 100\) 이면 BFGS, \(p > 100\) 이면 L-BFGS. 생존분석은 거의 항상 \(p \leq 50\) 이라 BFGS 가 default.

9 코드 예시

9.1 Step 1 — Python: scipy 의 통합 API

import numpy as np
from scipy.optimize import minimize

# Example A.2 데이터
t = np.array([2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37])
n = len(t)

# 음의 log-likelihood (minimize 기준)
def neg_loglik(params):
    lam, alpha = params
    if lam <= 0 or alpha <= 0:
        return np.inf
    return -(n*np.log(lam) + n*np.log(alpha)
             + (alpha-1)*np.sum(np.log(t)) - lam*np.sum(t**alpha))

# Gradient (식 A.3)
def neg_grad(params):
    lam, alpha = params
    return -np.array([
        n/lam - np.sum(t**alpha),
        n/alpha + np.sum(np.log(t)) - lam*np.sum(t**alpha * np.log(t))
    ])

# Hessian (식 A.4)
def neg_hess(params):
    lam, alpha = params
    H = np.array([
        [-n/lam**2, -np.sum(t**alpha * np.log(t))],
        [-np.sum(t**alpha * np.log(t)),
         -n/alpha**2 - lam*np.sum(t**alpha * np.log(t)**2)]
    ])
    return -H  # neg_loglik 의 Hessian

# 1. Newton-Raphson (Newton-CG)
res_nr = minimize(neg_loglik, x0=[1.024, 1.0], jac=neg_grad, hess=neg_hess,
                  method='Newton-CG', options={'xtol': 1e-6})
print(f"Newton-CG: λ̂={res_nr.x[0]:.4f}, α̂={res_nr.x[1]:.4f}, iter={res_nr.nit}")

# 2. BFGS (Hessian 자동 근사)
res_bfgs = minimize(neg_loglik, x0=[1.024, 1.0], jac=neg_grad,
                    method='BFGS', options={'gtol': 1e-6})
print(f"BFGS: λ̂={res_bfgs.x[0]:.4f}, α̂={res_bfgs.x[1]:.4f}, iter={res_bfgs.nit}")

# 3. L-BFGS-B (with bounds)
res_lb = minimize(neg_loglik, x0=[1.024, 1.0], jac=neg_grad,
                  method='L-BFGS-B', bounds=[(0.01, None), (0.01, None)])
print(f"L-BFGS-B: λ̂={res_lb.x[0]:.4f}, α̂={res_lb.x[1]:.4f}, iter={res_lb.nit}")

# 4. Steepest descent (직접 구현)
def steepest_ascent(f_obj, grad, x0, max_iter=50, tol=0.1):
    x = np.array(x0, dtype=float)
    history = [x.copy()]
    for k in range(max_iter):
        g = -grad(x)  # f_obj 는 neg, gradient 도 neg
        if np.max(np.abs(g)) < tol:
            break
        # 1D line search
        from scipy.optimize import minimize_scalar
        res = minimize_scalar(lambda d: f_obj(x + d*g), method='brent',
                               bracket=(0, 0.1, 0.5))
        x = x + res.x * g
        history.append(x.copy())
    return x, k, history

x_sa, k_sa, _ = steepest_ascent(neg_loglik, neg_grad, [1.024, 1.0])
print(f"Steepest: λ̂={x_sa[0]:.4f}, α̂={x_sa[1]:.4f}, iter={k_sa}")

9.2 Step 2 — R: optim() 의 다양한 method

t <- c(2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37)
n <- length(t)

# 음의 log-likelihood
neg_loglik <- function(params) {
  lam <- params[1]; alpha <- params[2]
  if (lam <= 0 || alpha <= 0) return(Inf)
  -(n*log(lam) + n*log(alpha) + (alpha-1)*sum(log(t)) - lam*sum(t^alpha))
}

# Gradient
neg_grad <- function(params) {
  lam <- params[1]; alpha <- params[2]
  c(-(n/lam - sum(t^alpha)),
    -(n/alpha + sum(log(t)) - lam*sum(t^alpha * log(t))))
}

# 1. BFGS (default for smooth)
fit_bfgs <- optim(c(1.024, 1.0), neg_loglik, gr = neg_grad,
                  method = "BFGS", control = list(reltol = 1e-8))
cat("BFGS:", fit_bfgs$par, "iters:", fit_bfgs$counts, "\n")

# 2. L-BFGS-B (with bounds)
fit_lbfgs <- optim(c(1.024, 1.0), neg_loglik, gr = neg_grad,
                   method = "L-BFGS-B", lower = c(0.01, 0.01))
cat("L-BFGS-B:", fit_lbfgs$par, "\n")

# 3. Nelder-Mead (gradient-free, robust)
fit_nm <- optim(c(1.024, 1.0), neg_loglik, method = "Nelder-Mead")
cat("Nelder-Mead:", fit_nm$par, "\n")

# 4. CG (Conjugate Gradient — large-scale 대안)
fit_cg <- optim(c(1.024, 1.0), neg_loglik, gr = neg_grad, method = "CG")
cat("CG:", fit_cg$par, "\n")

# 5. nlminb (Quasi-Newton + box constraint, 기본 default)
fit_nlm <- nlminb(c(1.024, 1.0), neg_loglik, gradient = neg_grad,
                  lower = c(0.01, 0.01))
cat("nlminb:", fit_nlm$par, "iter:", fit_nlm$iterations, "\n")

9.3 Step 3 — Marquardt 직접 구현

def marquardt(f_obj, grad, hess, x0, gamma_0=0.001, max_iter=50, tol=1e-6):
    """
    Marquardt's algorithm with adaptive γ.
    f_obj: function to MAXIMIZE (positive log-likelihood).
    """
    x = np.array(x0, dtype=float)
    gamma = gamma_0
    f_curr = f_obj(x)

    for k in range(max_iter):
        g = grad(x)                                # gradient
        if np.max(np.abs(g)) < tol:
            break
        H = hess(x)                                # Hessian
        S = np.diag(np.abs(np.diag(H))**(-0.5))    # 스케일 행렬
        M = S @ H @ S + gamma * np.eye(len(x))
        delta = S @ np.linalg.solve(M, S @ g)
        x_new = x - delta

        f_new = f_obj(x_new)
        if f_new > f_curr:
            x = x_new
            f_curr = f_new
            gamma /= 10                            # Newton 쪽으로
        else:
            gamma *= 10                            # Steepest 쪽으로
            if gamma > 1e10:
                break                              # 발산 방지

    return x, k, gamma

10 핵심 요약

§ A.2 한 줄 요약

다변수 최적화는 gradient \(\mathbf{u}\) (식 A.3) 와 Hessian \(\mathbf{H}\) (식 A.4) 를 도구로 사용한다. Steepest Ascent (식 A.5) 는 gradient 방향 + 1D line search — Hessian 불필요하지만 길쭉한 우도에서 지그재그로 수렴 매우 느림. Newton-Raphson 다변수 (식 A.6, \(-\mathbf{H}^{-1} \mathbf{u}\)) 는 Hessian 으로 좌표를 정규화한 후 한 step 으로 점프 — 2차 수렴이지만 \(\mathbf{H}\) 비싸거나 ill-conditioned 시 위험. Marquardt (\(\gamma\) blending) 는 둘의 절충 — \(\gamma = 0\) 이면 Newton, \(\gamma \to \infty\) 면 steepest. Klein Example A.2 의 Weibull (λ, α) MLE 에서 NR 3 step, Steepest 5 step, Marquardt 5 step. 현대 표준은 BFGS / L-BFGS (Hessian 을 gradient 변화로 자동 근사) — optim / scipy.minimize 의 default.

방법	식	정보 요구	수렴 차수	견고성
Steepest Ascent	식 (A.5)	\(\mathbf{u}\) (+ 1D line search)	선형	매우 강
Newton-Raphson	식 (A.6)	\(\mathbf{u}, \mathbf{H}\)	2 (이차)	약
Marquardt	\(\gamma\) blending	\(\mathbf{u}, \mathbf{H}, \gamma\)	적응	중간-강
BFGS (Quasi-Newton)	rank-1 update	\(\mathbf{u}\) + Hessian 근사	super-linear	강

5 가지 실무 권고

1. Default 는 BFGS / L-BFGS: \(f''\) 명시적 계산 부담 회피, super-linear 수렴, 양의 정칙 자동 유지. optim / scipy.minimize default.

2. 좋은 시작값이 핵심: Method-of-moments, OLS, profile likelihood 등으로 사전 추정. 시작값에서 최적까지 거리가 멀수록 발산 위험.

3. Step-Halving 안전장치: Newton 의 단순 구현 위험. \(f\) 값 모니터링 + step 절반 fallback.

4. 비선형 회귀에는 Marquardt (LM): nls(), minpack.lm::nlsLM(), scipy.optimize.least_squares(method='lm'). 작은 표본 / 모형 오설정에서 가장 견고.

5. 수렴 후 검증: \(\mathbf{H}(\hat{\mathbf{x}})\) 가 음의 정칙 (eigenvalue < 0) 인지 확인. 안 그러면 안장점 또는 local min — 다른 시작값으로 재시도.

11 관련 주제

Klein Appendix A 시리즈

• Appendix A — Numerical Techniques for Maximization — overview (단변수 + 다변수 종합)
• § A.1 — Univariate Methods — Bisection, Secant, Newton-Raphson 의 알고리즘과 수렴 차수

본문 응용

• Ch.8 — Cox 비례위험 모형 — 다변수 부분우도 Newton-Raphson + step-halving
• Ch.12 — 모수적 회귀 모형 — Weibull/LL/LN 다변수 MLE
• § 13.3-13.4 — Gamma Frailty EM & Marginal Model — EM M-step 의 다변수 NR

관련 개념

• BFGS / L-BFGS (Wikipedia)
• Levenberg-Marquardt 알고리즘 (Wikipedia)
• Trust region method (Wikipedia) — Marquardt 의 현대 일반화

12 참고 문헌

• Klein, J. P., & Moeschberger, M. L. (2003). Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.). Springer. Appendix A.2.
• Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. J. SIAM, 11(2), 431-441. (Marquardt 의 원논문)
• Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. Quart. Appl. Math., 2(2), 164-168. (LM 의 한쪽 역사)
• Broyden, C. G., Fletcher, R., Goldfarb, D., & Shanno, D. F. (1970년대). BFGS 의 4 명 독립 발견. (Wikipedia 참조)
• Nocedal, J., & Wright, S. J. (2006). Numerical Optimization (2nd ed.). Springer. (현대 표준 — line search, trust region, BFGS, L-BFGS, conjugate gradient 종합)
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes (3rd ed.). Cambridge. (Marquardt, BFGS 의 실용 구현)
• Thisted, R. A. (1988). Elements of Statistical Computing: Numerical Computation. Chapman and Hall. (Klein 인용 표준 교재)