第6章ニュートン法と準ニュートン法

この章の目標

• ニュートン法を 2 次近似の最小化として導き、ヘッセ行列のリプシッツ連続性と正則性のもとで局所 2 次収束を証明できる
• ニュートン法が大域的には発散しうる例を挙げ、減衰ニュートン法で直せることを説明できる
• 非線形最小二乗問題にガウス–ニュートン法とレーベンバーグ–マーカート法を適用し、両者の違いを説明できる
• セカント条件から BFGS 更新を導き、曲率条件のもとで正定値性が保たれることを証明できる
• L-BFGS と、対数バリアと中心パスによる内点法の考え方を説明できる

前提：第5章、01-calculus 第7章（ヘッセ行列・作用素ノルム）。6.7 節では第3章の生産計画の例と第4章の弱双対性を使う。

第5章で見たように、勾配法の反復回数は条件数 $\kappa$ に比例し、悪条件の問題では数万回に達する。ニュートン法はヘッセ行列で各方向の曲がり方を補正するので、変数のスケールの影響を受けず、解の近くでは正しい桁数が毎回ほぼ 2 倍になる。ロジスティック回帰など正準リンクの一般化線形モデルの標準的な計算法（IRLS）もニュートン法である（22 統計学 第6章 定理 6.7）。

代償は二つある。1 回の反復でヘッセ行列を作って $n$ 元の連立一次方程式を解く手間と、解から遠いと発散しうることである。本章では、ニュートン法の局所的な速さを証明したあと、大域的な振る舞いを直す減衰ニュートン法、最小二乗に特化したガウス–ニュートン法とレーベンバーグ–マーカート法、ヘッセ行列を勾配から推定する準ニュートン法（BFGS, L-BFGS）を学ぶ。最後に、制約付き問題にニュートン法を使う内点法の考え方を述べる。

$\lVert A \rVert$ は行列の作用素ノルム $\sup_{\lVert x \rVert \leq 1}\lVert Ax \rVert$ である（01-calculus 第7章。$\lVert AB \rVert \leq \lVert A \rVert\lVert B \rVert$ が成り立つ）。線形収束・2 次収束などの言葉は第5章で定義した。

6.1 ニュートン法

点 $x_k$ のまわりで $f$ を 2 次のテイラー多項式

で近似する。$\nabla^2 f(x_k) \succ 0$ なら $m_k$ は強凸な二次関数で、最小点は $\nabla^2 f(x_k)d = -\nabla f(x_k)$ の解である（第1章 命題 1.21）。これは方程式 $\nabla f(x) = 0$ を $x_k$ で 1 次近似して解くこと（1 変数のニュートン–ラフソン法の多変数版）とも同じである。

定義 6.1（ニュートン法, Newton's method）$\nabla^2 f(x_k)$ が正則のとき、連立一次方程式 $\nabla^2 f(x_k)d_k = -\nabla f(x_k)$ を解いて $x_{k+1} = x_k + d_k$ とする。$d_k$ をニュートン方向 (Newton direction) という。

実装では逆行列を作らず、連立一次方程式を解く（正定値ならコレスキー分解で約 $n^3/3$ 回の浮動小数点演算）。変数を $x = Ay + b$（$A$ は正則）と取り替えてもニュートン法の点列は対応する（問題 6.2）。勾配法と違って、変数のスケールに依存しない。

例 6.2 $f(x) = x - \log x$（$x > 0$、最小点 $x^{\ast} = 1$）では $f'(x) = 1 - 1/x$, $f''(x) = 1/x^2$ なので、$x_{k+1} = x_k - (1 - 1/x_k)x_k^2 = 2x_k - x_k^2$。誤差 $e_k = 1 - x_k$ は $e_{k+1} = e_k^2$ を満たす。$x_0 = 1/2$ なら $x_k = 3/4, 15/16, 255/256, 65535/65536, \dots$ で、$e_5 = 2^{-32} \approx 2.3 \times 10^{-10}$。正しい桁数が毎回 2 倍になる。一方 $x_0 = 3$ なら $x_1 = -3$ で、定義域から出てしまう（問題 6.1）。

6.2 局所 2 次収束

補題 6.3（逆行列の摂動）$A$ を正則な $n$ 次正方行列、$B$ を $n$ 次正方行列とし、$\lVert A^{-1} \rVert\lVert B - A \rVert < 1$ とする。このとき $B$ は正則で

$$\lVert B^{-1} \rVert \leq \frac{\lVert A^{-1} \rVert}{1 - \lVert A^{-1} \rVert\lVert B - A \rVert}$$

証明. $\lVert x \rVert = \lVert A^{-1}Ax \rVert \leq \lVert A^{-1} \rVert\lVert Ax \rVert$ より、すべての $x$ で

右辺の係数は正なので $Bx = 0 \Rightarrow x = 0$、すなわち $B$ は単射で正則である。$x = B^{-1}y$ を代入すれば評価を得る。$\square$

定理 6.4（ニュートン法の局所 2 次収束）$U \subset \mathbb{R}^n$ を開集合、$f\colon U \to \mathbb{R}$ を $C^2$ 級とし、$x^{\ast} \in U$ で $\nabla f(x^{\ast}) = 0$ かつ $\nabla^2 f(x^{\ast})$ は正則とする。さらに $r > 0$, $M > 0$ があって、閉球 $\bar{B} = \lbrace x \mid \lVert x - x^{\ast} \rVert \leq r \rbrace \subset U$ 上で

$$\lVert \nabla^2 f(x) - \nabla^2 f(y) \rVert \leq M\lVert x - y \rVert \quad (x, y \in \bar{B})$$

とする（ヘッセ行列のリプシッツ連続性）。$\beta = \lVert \nabla^2 f(x^{\ast})^{-1} \rVert$, $\delta = \min(r, \frac{1}{2\beta M})$ とおく。$\lVert x_0 - x^{\ast} \rVert < \delta$ ならば、ニュートン法の点列はすべて定義されて $\lVert x_k - x^{\ast} \rVert < \delta$ を満たし、

$$\lVert x_{k+1} - x^{\ast} \rVert \leq \beta M\lVert x_k - x^{\ast} \rVert^2 \leq \frac{1}{2}\lVert x_k - x^{\ast} \rVert$$

が成り立つ。特に $x_k \to x^{\ast}$ で、収束は 2 次である。

証明. $\lVert x - x^{\ast} \rVert < \delta$ となる $x$ について、$x^{+} = x - \nabla^2 f(x)^{-1}\nabla f(x)$ が定義されて $\lVert x^{+} - x^{\ast} \rVert \leq \beta M\lVert x - x^{\ast} \rVert^2$ となることを示せばよい（このとき右辺は $\beta M\delta\lVert x - x^{\ast} \rVert \leq \frac{1}{2}\lVert x - x^{\ast} \rVert$ 以下なので、帰納法ですべての主張が従う）。

(i) $\lVert \nabla^2 f(x) - \nabla^2 f(x^{\ast}) \rVert \leq M\lVert x - x^{\ast} \rVert < M\delta \leq \frac{1}{2\beta}$ なので、補題 6.3（$A = \nabla^2 f(x^{\ast})$, $B = \nabla^2 f(x)$）より $\nabla^2 f(x)$ は正則で、$\lVert \nabla^2 f(x)^{-1} \rVert \leq \beta/(1 - \frac{1}{2}) = 2\beta$。

(ii) $h = x^{\ast} - x$、$R = \nabla f(x^{\ast}) - \nabla f(x) - \nabla^2 f(x)h$ とおく。単位ベクトル $w$ について $\varphi(s) = \langle w, \nabla f(x + sh) \rangle$（$0 \leq s \leq 1$。線分は $\bar{B}$ に含まれる）は $C^1$ 級で、連鎖律より $\varphi'(s) = \langle w, \nabla^2 f(x + sh)h \rangle$。微分積分学の基本定理（01-calculus 第5章 定理 5.13）より

$R \neq 0$ なら $w = R/\lVert R \rVert$ として $\lVert R \rVert \leq \frac{M}{2}\lVert x - x^{\ast} \rVert^2$。

(iii) $\nabla f(x^{\ast}) = 0$ より $x^{+} - x^{\ast} = \nabla^2 f(x)^{-1}\bigl(\nabla^2 f(x)(x - x^{\ast}) - \nabla f(x)\bigr) = \nabla^2 f(x)^{-1}R$。(i)(ii) より $\lVert x^{+} - x^{\ast} \rVert \leq 2\beta \cdot \frac{M}{2}\lVert x - x^{\ast} \rVert^2 = \beta M\lVert x - x^{\ast} \rVert^2$。$\square$

$\nabla^2 f(x^{\ast}) \succ 0$ なら $x^{\ast}$ は狭義の局所最小点である（第1章 定理 1.19）。仮定の役割を確かめておく。

• 正則性：$f(x) = x^4$ では $f''(0) = 0$ で、ニュートン法は $x_{k+1} = x_k - \frac{4x_k^3}{12x_k^2} = \frac{2}{3}x_k$ となり、線形収束しかしない。
• 初期点の近さ：定理の $\delta$ は保守的なことが多い。例 6.2 では、定理は $\lvert x_0 - 1 \rvert$ が $0.15$ 程度以下なら収束を保証するが（問題 6.1）、実際には $0 < x_0 < 2$ で収束する。それでも、遠い初期点で収束が保証されないことは本質的である（6.3 節）。
• 最小点とは限らない：定理は $\nabla f(x^{\ast}) = 0$ と正則性しか使わないので、極大点や鞍点にも同じ速さで収束する。

6.3 大域的な振る舞いと減衰ニュートン法

例 6.5（狭義凸でも発散する）$f(x) = \sqrt{1 + x^2}$ は $f''(x) = (1 + x^2)^{-3/2} > 0$ で狭義凸、最小点は $0$ だけである（$\lvert x \rvert \to \infty$ で $f''(x) \to 0$ なので強凸ではない）。$f'(x) = x/\sqrt{1 + x^2}$ より $f'(x)/f''(x) = x(1 + x^2)$ で、ニュートン法は

となる。$\lvert x_0 \rvert < 1$ なら非常に速く $0$ に収束するが（$x_0 = 0.5$ から $-0.125, 0.00195, -7.5 \times 10^{-9}$）、$\lvert x_0 \rvert = 1$ なら $\pm 1$ を往復し、$\lvert x_0 \rvert > 1$ なら発散する（$x_0 = 1.1$ から $-1.331, 2.358, -13.11, 2253, \dots$）。遠くでは $f$ がほぼ $\lvert x \rvert$ のように平らで、2 次近似の最小点がはるか遠くに出るからである。

対策は、ニュートン方向を「進む方向」としてだけ使い、進む距離は直線探索で決めることである。

定義 6.6（減衰ニュートン法, damped Newton method）$\nabla^2 f(x_k) \succ 0$ のとき、ニュートン方向 $d_k$ に沿って、$\bar{t} = 1$ から始めるアルミホ条件のバックトラッキング（第5章 定義 5.5）で $t_k$ を選び、$x_{k+1} = x_k + t_kd_k$ とする。

$\nabla^2 f(x_k) \succ 0$ で $\nabla f(x_k) \neq 0$ なら $\nabla f(x_k)^{\top}d_k = -\nabla f(x_k)^{\top}\nabla^2 f(x_k)^{-1}\nabla f(x_k) < 0$ なので、$d_k$ は降下方向で、バックトラッキングは有限回で終わる（第5章 命題 5.6）。$\lambda(x) = \bigl(\nabla f(x)^{\top}\nabla^2 f(x)^{-1}\nabla f(x)\bigr)^{1/2}$ をニュートン減少量 (Newton decrement) という。$x$ での 2 次近似を $m(d)$ とすると $\lambda(x)^2/2 = f(x) - \min_d m(d)$ は 2 次近似で見込まれる減少量なので、停止判定に使われる。

例 6.5 で $x_0 = 3$ から、$c = 1/4$, $\beta = 1/2$ の減衰ニュートン法を行うと、$t_k = 1/8, 1/2, 1, 1, 1$ と選ばれて $x_k = -0.75, -0.164, 0.0044, -8.6 \times 10^{-8}, 6.4 \times 10^{-22}$ となる（多倍長の計算で確かめた。最後の値は倍精度では桁落ちで数 % ずれる）。最初は短い歩幅で近づき、近づくと $t = 1$ が受け入れられて純粋なニュートン法の速さに戻る。この振る舞いは一般に成り立つ。

定理 6.7（減衰ニュートン法の 2 段階の収束、主張）$f\colon \mathbb{R}^n \to \mathbb{R}$ を $C^2$ 級の凸関数とし、初期点の下位集合 $S = \lbrace x \mid f(x) \leq f(x_0) \rbrace$ 上で $\mu I \preceq \nabla^2 f(x) \preceq LI$（$\mu > 0$）を満たし、$\nabla^2 f$ が $S$ 上でリプシッツ連続であるとする。アルミホ条件の定数を $0 < c < 1/2$ とすると、減衰ニュートン法は最小点 $x^{\ast}$ に収束する。さらに、有限回の反復のあとは $t_k = 1$ が常に受け入れられ、収束は 2 次である。$f(x_k) - p^{\ast} \leq \varepsilon$ までの反復回数は $\frac{f(x_0) - p^{\ast}}{\gamma} + \log_2\log_2\frac{\varepsilon_0}{\varepsilon}$ 以下である（$\gamma, \varepsilon_0$ は $\mu, L$、ヘッセ行列のリプシッツ定数と直線探索の定数 $c, \beta$ で決まる正の数。Boyd–Vandenberghe, Convex Optimization 9.5 節）。

$\log_2\log_2(\varepsilon_0/\varepsilon)$ は $\varepsilon = 10^{-15}\varepsilon_0$ でも $6$ 程度で、2 次収束の段階は実際上数回で終わる。$c < 1/2$ が要るのは、二次関数では $t = 1$ の一歩による減少量がちょうど $\lambda(x)^2/2$ で、アルミホ条件 $f(x + d) \leq f(x) - c\lambda(x)^2$ が $c \leq 1/2$ のときしか成り立たないからである（$c > 1/2$ では解の近くでも $t = 1$ が受け入れられず、収束は線形に落ちる）。前半の段階でも線形収束することは、第5章と同じ議論で示せる（問題 6.7）。

6.4 ガウス–ニュートン法とレーベンバーグ–マーカート法

測定値に非線形なモデルをあてはめる問題（薬の血中濃度の減衰、センサーの較正曲線、需要の飽和曲線など）は、$r\colon \mathbb{R}^n \to \mathbb{R}^m$（残差、$m \geq n$）について

を最小化する非線形最小二乗問題になる。ヤコビ行列を $J(x) = (\partial r_i/\partial x_j)$（$m \times n$）とすると

である。残差が小さいときや $r$ が一次関数に近いときは、第 2 項が小さい。

定義 6.8（ガウス–ニュートン法・レーベンバーグ–マーカート法）ヘッセ行列の第 2 項を無視して $J_k^{\top}J_kd = -J_k^{\top}r_k$（$J_k = J(x_k)$, $r_k = r(x_k)$）を解き、$x_{k+1} = x_k + d$ とする方法をガウス–ニュートン法 (Gauss–Newton method) という。$\lambda > 0$ について $(J_k^{\top}J_k + \lambda I)d = -J_k^{\top}r_k$ を解く方法をレーベンバーグ–マーカート法 (Levenberg–Marquardt method) という。

ガウス–ニュートン方向は、線形化した残差の最小二乗問題 $\min_d\lVert J_kd + r_k \rVert^2$ の解である（02-linear-algebra 第7章 定理 7.19 の正規方程式）。2 階微分が要らず、$J_k^{\top}J_k$ は常に半正定値である。

命題 6.9 $g = J^{\top}r \neq 0$ とする。

1. $J$ の列が一次独立なら、ガウス–ニュートン方向 $d_{\mathrm{GN}} = -(J^{\top}J)^{-1}g$ は降下方向である。
2. $\lambda > 0$ なら、$d(\lambda) = -(J^{\top}J + \lambda I)^{-1}g$ は降下方向であり、$\lVert d \rVert \leq \lVert d(\lambda) \rVert$ の範囲で $\lVert Jd + r \rVert^2$ を最小にする。

証明. 1・2 の前半：対称行列 $P = J^{\top}J$（1 の場合。$v^{\top}Pv = \lVert Jv \rVert^2$ で列が一次独立だから）または $P = J^{\top}J + \lambda I$（2 の場合）は正定値なので $P^{-1}$ も正定値で、$g^{\top}d = -g^{\top}P^{-1}g < 0$。2 の後半：$d(\lambda)$ は強凸な二次関数 $q(d) = \lVert Jd + r \rVert^2 + \lambda\lVert d \rVert^2$ の最小点である（$\nabla q = 2(J^{\top}Jd + J^{\top}r + \lambda d)$）。$\lVert d \rVert \leq \lVert d(\lambda) \rVert$ なら $\lVert Jd + r \rVert^2 = q(d) - \lambda\lVert d \rVert^2 \geq q(d(\lambda)) - \lambda\lVert d \rVert^2 \geq \lVert Jd(\lambda) + r \rVert^2$。$\square$

$\lambda$ が小さければガウス–ニュートン法に近く、$\lambda$ が大きければ $d(\lambda) \approx -g/\lambda$ と短い勾配方向になる（問題 6.4）。命題 6.9 の 2 は、レーベンバーグ–マーカート法が「線形化が信頼できる半径 $\lVert d(\lambda) \rVert$ の中で最良の一歩」をとる信頼領域法 (trust-region method) であることを意味する。実装では、$f$ が減れば一歩を採用して $\lambda$ を小さくし（線形化を信頼する）、減らなければ一歩を捨てて $\lambda$ を大きくする。収束について、ガウス–ニュートン法は解での残差が $0$ で $J$ の列が一次独立なら局所的に 2 次収束し、残差が小さければ速く線形収束するが、残差が大きいと収束しないこともある（主張。Nocedal–Wright 第10章）。

次のコードは、$y = 10e^{-0.5t}$ に雑音を加えた 10 点のデータに $y \approx ae^{-bt}$ をあてはめる。

import numpy as np

rng = np.random.default_rng(0)
t = np.arange(10.0)                                    # 計測時刻
y = 10.0 * np.exp(-0.5 * t) + 0.2 * rng.standard_normal(10)
r = lambda p: p[0] * np.exp(-p[1] * t) - y            # 残差、p = (a, b)
J = lambda p: np.column_stack([np.exp(-p[1] * t), -p[0] * t * np.exp(-p[1] * t)])
F = lambda p: 0.5 * np.sum(r(p) ** 2)

def gauss_newton(p, iters=20):
    for _ in range(iters):
        p = p + np.linalg.lstsq(J(p), -r(p), rcond=None)[0]
    return p

def levenberg_marquardt(p, iters=20, lam=1.0):
    for _ in range(iters):
        Jp = J(p)
        d = np.linalg.solve(Jp.T @ Jp + lam * np.eye(2), -Jp.T @ r(p))
        if F(p + d) < F(p):
            p, lam = p + d, lam / 10      # 成功：ガウス–ニュートン法に近づける
        else:
            lam = lam * 10                # 失敗：短い勾配方向に近づける
    return p

for p0 in [(1.0, 0.0), (1.0, 1.0)]:
    p_gn, p_lm = gauss_newton(np.array(p0)), levenberg_marquardt(np.array(p0))
    print(f"start {p0}:  GN {np.round(p_gn, 4)} F={F(p_gn):.4f}"
          f"  LM {np.round(p_lm, 4)} F={F(p_lm):.4f}")

start (1.0, 0.0):  GN [10.0022  0.4933] F=0.1066  LM [10.0022  0.4933] F=0.1066
start (1.0, 1.0):  GN [-0.     -7.1409] F=79.8532  LM [10.0022  0.4933] F=0.1066

$(a, b) = (1, 0)$ からはどちらも最小点 $(10.0022, 0.4933)$ に着く。$(1, 1)$ からのガウス–ニュートン法は、最初の一歩で $b \approx -7.14$（減衰ではなく急増するモデル）に飛んで $F \approx 3.2 \times 10^{57}$ となり、その後 $a \approx 0$（$-0.$ は $-1.7 \times 10^{-29}$ を丸めた表示）で動けなくなる。このときモデル $ae^{7.14t}$ は、最後の時刻 $t = 9$ の測定値（雑音で負の $-0.14$ になっている）だけに合い、ほかの時刻ではほぼ $0$ である。$J$ の第 1 列 $e^{7.14t}$ は $1$ から約 $8 \times 10^{27}$ まで変わり、$J$ の 2 つの特異値は約 $8 \times 10^{27}$ と $1 \times 10^{-4}$ になる。lstsq は小さいほうの特異値を数値的に $0$ とみなして捨てるので、一歩は事実上 $0$ になる（勾配は $0$ でなく、停留点で止まったのではない）。レーベンバーグ–マーカート法は値が増える一歩を捨てるので、同じ初期点から最小点に着く。

6.5 準ニュートン法と BFGS 更新

ヘッセ行列が計算できない、または計算が高すぎるとき、勾配の変化からヘッセ行列を推定しながら進むのが準ニュートン法 (quasi-Newton method) である。$s_k = x_{k+1} - x_k$, $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$ とおく。

定義 6.10（セカント条件, secant condition）対称行列 $B_{k+1}$ が $B_{k+1}s_k = y_k$ を満たすとき、$B_{k+1}$ はセカント条件を満たすという。

微分積分学の基本定理を成分ごとに使うと $y_k = \bigl(\int_0^1 \nabla^2 f(x_k + \tau s_k)\ d\tau\bigr)s_k$ なので、セカント条件は「$B_{k+1}$ が線分上の平均のヘッセ行列と $s_k$ 方向で一致する」ことを要求する。1 変数なら $B_{k+1} = (f'(x_{k+1}) - f'(x_k))/(x_{k+1} - x_k)$ で、割線（セカント）の傾きである。

正定値な $B_{k+1}$ がセカント条件を満たすには、$s_k^{\top}y_k = s_k^{\top}B_{k+1}s_k > 0$ でなければならない。この曲率条件 (curvature condition) $s_k^{\top}y_k > 0$ は、$f$ が $\mu$-強凸なら $s_k^{\top}y_k \geq \mu\lVert s_k \rVert^2$（第2章 定理 2.22 の 3）で自動的に成り立ち、一般の $f$ では次の直線探索で保証される。

定義 6.11（ウルフ条件, Wolfe conditions）$0 < c_1 < c_2 < 1$ とし、$d$ を降下方向とする。ステップ幅 $t$ がアルミホ条件 $f(x + td) \leq f(x) + c_1t\nabla f(x)^{\top}d$ と $\nabla f(x + td)^{\top}d \geq c_2\nabla f(x)^{\top}d$ を満たすとき、ウルフ条件を満たすという。

$t$ がウルフ条件を満たせば、$s = td$, $y = \nabla f(x + td) - \nabla f(x)$ について $s^{\top}y = t(\nabla f(x + td) - \nabla f(x))^{\top}d \geq t(c_2 - 1)\nabla f(x)^{\top}d > 0$ である。$f$ が $C^1$ 級でその直線上で下に有界なら、ウルフ条件を満たす $t$ は存在する（主張。Nocedal–Wright 第3章）。準ニュートン法では $c_2 = 0.9$ がよく使われる。

セカント条件だけでは $B_{k+1}$ は決まらない（$n(n + 1)/2$ 個の未知数に $n$ 個の式）。そこで $B_k$ を少しだけ修正する。

BFGS 更新の導出. $B_k$ を対称正定値とし、$B_{k+1}$ を対称な 2 階の修正 $B_{k+1} = B_k + a\ uu^{\top} + b\ vv^{\top}$（$a, b \in \mathbb{R}$, $u, v \in \mathbb{R}^n$）の形で探す。セカント条件は $B_ks_k + a(u^{\top}s_k)u + b(v^{\top}s_k)v = y_k$ となる。左辺の $B_ks_k$ を打ち消して $y_k$ を作るために $u = y_k$, $v = B_ks_k$ とし、$a(y_k^{\top}s_k) = 1$, $b(s_k^{\top}B_ks_k) = -1$ とすれば満たされる。こうして

を得る。第 2 項は $B_k$ がもっていた $s_k$ 方向の曲率の情報を取り除き、第 3 項は観測した曲率 $y_k$ を入れる。これを BFGS 更新（ブロイデン–フレッチャー–ゴールドファーブ–シャノ）という。同じ考え方を逆行列 $H_k = B_k^{-1}$ について行うと DFP 更新が得られる。BFGS 更新の逆行列 $H_{k+1}$（定理 6.12 の 3）は、$Hy_k = s_k$（逆行列の形のセカント条件）を満たす対称行列 $H$ のうち、（ある重みつきのフロベニウスノルムで）$H_k$ に最も近いものとしても特徴づけられる（主張。Nocedal–Wright 第6章）。

定理 6.12（BFGS 更新の性質）$B_k$ を対称正定値、$s_k^{\top}y_k > 0$ とし、$B_{k+1}$ を (6.1) で定める。$\rho_k = 1/(y_k^{\top}s_k)$, $H_k = B_k^{-1}$ とする。

1. $B_{k+1}$ は対称で、セカント条件 $B_{k+1}s_k = y_k$ を満たす。
2. $B_{k+1}$ は正定値である。
3. $B_{k+1}$ の逆行列は $H_{k+1} = (I - \rho_ks_ky_k^{\top})H_k(I - \rho_ky_ks_k^{\top}) + \rho_ks_ks_k^{\top}$ である。

証明. 添字 $k$ を省く。1：$B_{k+1}s = Bs - Bs\frac{s^{\top}Bs}{s^{\top}Bs} + y\frac{y^{\top}s}{y^{\top}s} = y$。

2：$z \neq 0$ について

内積 $\langle u, v \rangle_B = u^{\top}Bv$ についてのコーシー–シュワルツの不等式 $(s^{\top}Bz)^2 \leq (s^{\top}Bs)(z^{\top}Bz)$ より括弧内は $0$ 以上で、等号は $z$ が $s$ の定数倍のときに限る。そのとき $z = \alpha s$（$\alpha \neq 0$）で、最後の項は $\alpha^2y^{\top}s > 0$。いずれにしても $z^{\top}B_{k+1}z > 0$。

3：$H_{k+1}y = (I - \rho sy^{\top})H(y - \rho y(s^{\top}y)) + \rho s(s^{\top}y) = s$ なので $B_{k+1}H_{k+1}y = B_{k+1}s = y$。$s^{\top}z = 0$ となる $z$ では $H_{k+1}z = Hz - \rho(y^{\top}Hz)s$ で、$B_{k+1}Hz = z - Bs\frac{s^{\top}z}{s^{\top}Bs} + \rho y(y^{\top}Hz) = z + \rho(y^{\top}Hz)y$ と $B_{k+1}s = y$ より $B_{k+1}H_{k+1}z = z$。$s^{\top}y > 0$ より $y$ は $s$ の直交補空間に含まれないので、$y$ とこの補空間で $\mathbb{R}^n$ が張られ、$B_{k+1}H_{k+1} = I$。$\square$

逆に $s_k^{\top}y_k \leq 0$ なら、セカント条件を満たす正定値行列は存在しない（問題 6.5）。準ニュートン法でウルフ条件を使うのはこのためである。

BFGS 法は、$H_0 \succ 0$（たとえば $I$）から始め、$d_k = -H_k\nabla f(x_k)$ に沿ってウルフ条件を満たす $t_k$ で進み、定理 6.12 の 3 で $H_k$ を更新する。$d_k$ は常に降下方向で、1 回の反復の手間は行列とベクトルの積と 2 階の更新の $O(n^2)$ で済み、連立一次方程式を解く必要もない。$f$ が $C^2$ 級で、初期点の下位集合 $\lbrace x \mid f(x) \leq f(x_0) \rbrace$ が凸で、その上で $\mu I \preceq \nabla^2 f \preceq LI$（$\mu > 0$）なら、ウルフ条件の BFGS 法は任意の $H_0 \succ 0$ から最小点に収束する。さらにヘッセ行列が最小点の近くでリプシッツ連続で、直線探索が $t = 1$ を最初に試し $c_1 < 1/2$ なら、収束は超 1 次である（主張。Nocedal–Wright 第6章）。

次のコードは、列のスケールが大きく異なる特徴量をもつリッジ正則化つきロジスティック回帰（第1章 例 1.6・例 1.13）で、勾配のノルムが初期値の $10^{-6}$ 倍以下になるまでの反復回数を比べる。直線探索はアルミホ条件のバックトラッキングだけを使っているが、この目的関数は強凸なので曲率条件は自動的に成り立つ。

import numpy as np

rng = np.random.default_rng(0)
m, n = 500, 10
scale = np.logspace(0, 1.5, n)                      # 列のスケールが 1〜約 32 倍
A = rng.standard_normal((m, n)) * scale
y = (rng.random(m) < 1 / (1 + np.exp(-A @ (rng.standard_normal(n) / scale)))).astype(float)
sig = lambda s: 1 / (1 + np.exp(-s))
f = lambda x: np.sum(np.logaddexp(0, A @ x) - y * (A @ x)) + x @ x / 2   # リッジ正則化つき
grad = lambda x: A.T @ (sig(A @ x) - y) + x
hess = lambda x: A.T @ (A * (sig(A @ x) * (1 - sig(A @ x)))[:, None]) + np.eye(n)

def run(method):
    x, H, I, k = np.zeros(n), np.eye(n), np.eye(n), 0
    g = grad(x); tol = 1e-6 * np.linalg.norm(g)
    while np.linalg.norm(g) > tol:
        if method == "gradient": d = -g
        elif method == "BFGS": d = -H @ g
        else: d = -np.linalg.solve(hess(x), g)
        t = 1.0                                            # アルミホ条件のバックトラッキング
        while f(x + t * d) > f(x) + 1e-4 * t * (g @ d):
            t /= 2
        s = t * d; x = x + s
        g_new = grad(x); yk = g_new - g; g, k = g_new, k + 1
        if method == "BFGS":                               # 強凸なので s^T y > 0
            rho = 1 / (yk @ s)
            H = (I - rho * np.outer(s, yk)) @ H @ (I - rho * np.outer(yk, s)) + rho * np.outer(s, s)
    return k

for method in ["gradient", "BFGS", "Newton"]:
    print(f"{method:8s} {run(method)}")

gradient 9037
BFGS     21
Newton   6

勾配法（バックトラッキングつき）は 9037 回かかる。BFGS 法は最初の 10 回ほどは短い歩幅で曲率を学習し、その後は $t = 1$ が受け入れられて勾配のノルムが急速に小さくなる。ニュートン法は 6 回で、最後の 3 回で勾配のノルムが $29 \to 2.1 \to 0.014 \to 5.8 \times 10^{-7}$ と 2 次収束の様子を示す。1 回の反復の手間は、勾配法が $O(mn)$、BFGS 法が $O(mn + n^2)$、ニュートン法がヘッセ行列を作る $O(mn^2)$ と解く $O(n^3)$ である。

6.6 L-BFGS（紹介）

$n = 10^6$ では $n \times n$ の $H_k$ を記憶することさえできない（$10^{12}$ 個の成分）。L-BFGS（記憶制限つき BFGS, limited-memory BFGS）は、直近の $m$ 組（数個から 20 個程度）の $(s_i, y_i)$ だけを記憶し、$H_k^0 = \gamma_kI$（$\gamma_k = s_{k-1}^{\top}y_{k-1}/y_{k-1}^{\top}y_{k-1}$ がよく使われる）から始めて、記憶した組で定理 6.12 の 3 の更新を古い順に $m$ 回施した行列 $H_k$ を、行列を作らずに $\nabla f(x_k)$ に掛ける。この積は 2 重のループによる再帰で $O(mn)$ の手間で計算できる（Nocedal–Wright 第7章）。記憶容量は $O(mn)$ で、勾配法とほぼ同じ手間で曲率の情報を使えるため、変数の多い滑らかな問題（大規模なロジスティック回帰、物理シミュレーションのパラメータ推定など）の標準的な解法の一つになっている。

6.7 内点法の考え方：対数バリアと中心パス

制約付き問題にニュートン法を使うために、制約を目的関数の「壁」に置き換える。不等式形の線形計画問題

を考え、$a_i^{\top}x < b_i$（すべての $i$）を満たす点があり、実行可能領域は有界であるとする。対数バリア関数 (logarithmic barrier) を

と定める。$\phi$ は凸で、境界に近づくと $+\infty$ に発散する。$t > 0$ について、$tc^{\top}x + \phi(x)$ を実行可能領域の内部で最小化する点を $x^{\ast}(t)$ とする（境界で $+\infty$ に発散し領域が有界なので最小点が存在する。また、$a_i^{\top}d = 0$（すべての $i$）となる $d \neq 0$ があれば直線 $x + sd$ が領域に含まれて有界性に反するので、$a_i$ たちは $\mathbb{R}^n$ を張り、$\nabla^2\phi = \sum_i a_ia_i^{\top}/(b_i - a_i^{\top}x)^2 \succ 0$ より最小点は一つである）。曲線 $t \mapsto x^{\ast}(t)$ を中心パス (central path) という。$t$ が小さいと壁の効果が強く $x^{\ast}(t)$ は領域の「中心」にあり、$t$ が大きいと目的関数が優勢になって最適解に近づく。

命題 6.13（中心パス上の双対ギャップ）$\lambda_i(t) = \dfrac{1}{t(b_i - a_i^{\top}x^{\ast}(t))} > 0$ とすると $A^{\top}\lambda(t) + c = 0$（$A$ は $a_i^{\top}$ を行とする行列）であり、最適値 $p^{\ast}$ について

$$c^{\top}x^{\ast}(t) - p^{\ast} \leq c^{\top}x^{\ast}(t) + b^{\top}\lambda(t) = \frac{m}{t}$$

証明. $x^{\ast}(t)$ で勾配が $0$ なので $tc + \sum_i\frac{a_i}{b_i - a_i^{\top}x^{\ast}(t)} = 0$ で、$t$ で割れば $c + A^{\top}\lambda = 0$。実行可能な任意の $x$ について、$\lambda \geq 0$ と $Ax \leq b$ より $c^{\top}x = -\lambda^{\top}Ax \geq -\lambda^{\top}b$ なので $p^{\ast} \geq -b^{\top}\lambda$（第4章 定理 4.11 の弱双対性の特別な場合）。また $c^{\top}x^{\ast}(t) + b^{\top}\lambda = \lambda^{\top}(b - Ax^{\ast}(t)) = \sum_i\frac{1}{t} = \frac{m}{t}$。$\square$

$\lambda(t)$ は双対問題の実行可能解で、KKT 条件（第4章 定義 4.2）の相補性 $\lambda_i(b_i - a_i^{\top}x) = 0$ を $\lambda_i(b_i - a_i^{\top}x) = 1/t$ にゆるめたものを満たしている。$t \to \infty$ で KKT 条件に近づく。