1. ニュートン法（Newton's Method）
1.1 基本原理と導出
最適化問題の目標は、目的関数 f(x) の最小値を見つけることです：
$$\min_{\mathbf{x}} f(\mathbf{x})$$
最適解では勾配（gradient）がゼロになるという条件を利用して、解を探索します。ニュートン法では、現在の点 x_k から次の点を効率的に求めるために、目的関数の二次近似を使います。
二次近似をテイラー展開（Taylor expansion）で表すと：
$$f(\mathbf{x}_k + \Delta \mathbf{x}) \approx f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^T \Delta \mathbf{x} + \frac{1}{2}\Delta \mathbf{x}^T \mathbf{H}_k \Delta \mathbf{x}$$
ここで $\mathbf{H}_k$ はヘッセ行列（Hessian matrix）です。
この近似関数を最小化する Δx を求めるために、勾配をゼロとおきます：
$$\nabla f(\mathbf{x}_k) + \mathbf{H}_k \Delta \mathbf{x} = \mathbf{0}$$
これを解くと最適な更新量が得られます：
$$\Delta \mathbf{x} = -\mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k)$$
よって、更新式は次のようになります：
$$\mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k)$$
この方法が強力なのは、凸関数（convex function）の場合、一度の反復で最適解に到達できる可能性があり、一般的に収束速度が速いためです。
1.2 JAX実装
import jax
import jax.numpy as jnp
from jax import grad, hessian, jit
@jit
def newton_step(f, x):
g = grad(f)(x)
H = hessian(f)(x)
d = -jnp.linalg.solve(H, g) # 数値的に安定
return d
@jit
def newton_method(f, x0, tol=1e-6, max_iter=100):
x = x0
for i in range(max_iter):
g = grad(f)(x)
if jnp.linalg.norm(g) < tol:
return x, i, True
d = newton_step(f, x)
x = x + d
return x, max_iter, False
2. ニュートン法の問題点と修正ニュートン法（Modified Newton Method）
2.1 問題点と必要性
標準的なニュートン法には重大な問題があります。ヘッセ行列が正定値でない場合（$\mathbf{H}_k \not\succ 0$）、更新方向が降下方向（descent direction）とならない可能性があるのです。
これがなぜ問題なのかを見てみましょう。二次近似の勾配から：
$$\nabla f(\mathbf{x}_k + \Delta \mathbf{x}) \approx \nabla f(\mathbf{x}_k) + \mathbf{H}_k \Delta \mathbf{x}$$
ニュートン方向 $\Delta \mathbf{x} = -\mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k)$ を代入すると：
$$\nabla f(\mathbf{x}_k + \Delta \mathbf{x}) \approx \nabla f(\mathbf{x}_k) - \mathbf{H}_k \mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k) = \mathbf{0}$$
これは二次近似の範囲内では、一度の更新で勾配がゼロになることを示しています。しかし、この更新が実際に目的関数を減少させるには、更新方向が降下方向である必要があります。
更新方向 $\mathbf{d} = -\mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k)$ が降下方向であるためには：
$$\mathbf{d}^T \nabla f(\mathbf{x}_k) < 0$$
これを展開すると：
$$(-\mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k))^T \nabla f(\mathbf{x}_k) = -\nabla f(\mathbf{x}_k)^T \mathbf{H}_k^{-1} \nabla f(\mathbf{x}_k) < 0$$
ヘッセ行列が正定値（positive definite、$\mathbf{H}_k \succ 0$）の場合、その逆行列も正定値となり、上の不等式が常に成立するため降下方向が保証されます。しかし、ヘッセ行列が正定値でない場合、この不等式は成立せず、更新方向が上り方向になってしまう可能性があります。
非凸関数（non-convex function）や鞍点（saddle point）付近では、ヘッセ行列は正定値でないことがあり、この問題が顕著になります。
2.2 行列の正定値性（Positive Definiteness）
2.2.1 正定値性の定義と性質
実対称行列 $\mathbf{A}$ について、以下の条件は同値です：
- $\mathbf{A}$ は正定値である（positive definite、$\mathbf{A} \succ 0$）
- すべての非ゼロベクトル $\mathbf{v}$ に対して $\mathbf{v}^T\mathbf{A}\mathbf{v} > 0$
- $\mathbf{A}$ のすべての固有値（eigenvalues）が正である
- $\mathbf{A} = \mathbf{L}\mathbf{L}^T$ となるコレスキー分解（Cholesky decomposition）が存在する
同様に、半正定値行列（positive semidefinite、$\mathbf{A} \succeq 0$）は、これらの条件の不等号を $\geq$ に置き換えたものと同値です。
2.2.2 コレスキー分解と正定値性の関係
コレスキー分解が存在することと行列が正定値であることの同値性は、修正ニュートン法の基礎となる重要な性質です。
証明:
実対称行列 $\mathbf{A}$ について：
$\mathbf{A} \succ 0 \Rightarrow$ コレスキー分解が存在：
- $\mathbf{A}$ が正定値の場合、ガウス消去法（Gaussian elimination）で $\mathbf{A} = \mathbf{L}\mathbf{D}\mathbf{L}^T$ と分解できます
- ここで $\mathbf{D}$ は正の対角成分のみを持つ対角行列です
- $\tilde{\mathbf{L}} = \mathbf{L}\mathbf{D}^{1/2}$ とすると $\mathbf{A} = \tilde{\mathbf{L}}\tilde{\mathbf{L}}^T$ となります
コレスキー分解が存在 $\Rightarrow \mathbf{A} \succ 0$：
- $\mathbf{A} = \mathbf{L}\mathbf{L}^T$ と表せる場合
- 任意の非ゼロベクトル $\mathbf{v}$ について $\mathbf{v}^T\mathbf{A}\mathbf{v} = \mathbf{v}^T\mathbf{L}\mathbf{L}^T\mathbf{v} = |\mathbf{L}^T\mathbf{v}|^2$
- $\mathbf{L}$ が正則なら $\mathbf{L}^T\mathbf{v} \neq \mathbf{0}$ なので $|\mathbf{L}^T\mathbf{v}|^2 > 0$
- よって $\mathbf{A}$ は正定値
2.3 修正手法
修正ニュートン法の目的は、ヘッセ行列が正定値でない場合でも降下方向を保証することです。主なアプローチを見ていきましょう。
2.3.1 修正コレスキー分解法（Modified Cholesky Factorization）
コレスキー分解を利用した修正法では、以下のステップで行列を正定値化します：
- ヘッセ行列 $\mathbf{H}$ のコレスキー分解を試みます
- 分解に失敗した場合（正定値でない）、$\mathbf{H} + \mu\mathbf{I}$ の形で修正します
- コレスキー分解が成功するまで $\mu$ を調整します
この方法のメリットは、最小限の修正でヘッセ行列を正定値化できることです。$\mu$ が小さければ元の二次収束性を近似的に保持できます。
@jit
def modified_cholesky(H, tau=1e-6, max_iter=10):
n = H.shape[0]
mu = 0.0
identity = jnp.eye(n)
for i in range(max_iter):
H_mod = H + mu * identity
try:
L = jnp.linalg.cholesky(H_mod)
return L, mu
except:
mu = max(2.0 * mu, tau) if mu == 0.0 else 2.0 * mu
H_mod = H + mu * identity
L = jnp.linalg.cholesky(H_mod)
return L, mu
@jit
def newton_direction_cholesky(f, x, tau=1e-6):
g = grad(f)(x)
H = hessian(f)(x)
L, mu = modified_cholesky(H, tau)
# LLᵀ y = -g を解く（前進・後退代入）
y = jnp.linalg.solve(L, -g)
d = jnp.linalg.solve(L.T, y)
return d, g, mu
2.3.2 固有値修正法（Eigenvalue Modification）
ヘッセ行列の固有値分解（eigenvalue decomposition）を利用する方法も効果的です：
$$\mathbf{H} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T$$
ここで $\mathbf{Q}$ は直交行列（orthogonal matrix）、$\mathbf{\Lambda}$ は固有値を対角成分とする対角行列です。
負または小さな固有値を持つ場合、次のように修正します：
$$\mathbf{H}' = \mathbf{Q}\mathbf{\Lambda}'\mathbf{Q}^T$$
ここで $\mathbf{\Lambda}' = \text{diag}(\max(\lambda_1, \epsilon), \ldots, \max(\lambda_n, \epsilon))$ であり、$\epsilon > 0$ は閾値です。
この修正により、元の固有ベクトルの方向は保持しながら、$\mathbf{H}'$ が正定値（$\mathbf{H}' \succ 0$）になります。これは元の行列の構造を尊重しつつ、必要最小限の修正を加える方法です。
@jit
def newton_direction_eigenval(f, x, epsilon=1e-6):
g = grad(f)(x)
H = hessian(f)(x)
eigvals, eigvecs = jnp.linalg.eigh(H)
eigvals_mod = jnp.maximum(eigvals, epsilon)
H_inv = eigvecs @ jnp.diag(1.0 / eigvals_mod) @ eigvecs.T
d = -H_inv @ g
return d, g
3. ラインサーチ（Line Search）
3.1 必要性と原理
修正ニュートン法でヘッセ行列を正定値化したとしても、二次近似の有効範囲は限られています。大きすぎるステップは発散の原因になり得ます。そこでラインサーチが必要になります。
ラインサーチでは、更新方向 $\mathbf{d}_k$ を固定し、ステップサイズ $\alpha_k$ を決定します：
$$\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{d}_k$$
適切なステップサイズを選ぶために、以下の条件を使います。
3.2 Armijo条件（Armijo Condition）
Armijo条件は十分な減少（sufficient decrease）を保証します：
$$f(\mathbf{x}_k + \alpha_k \mathbf{d}_k) \leq f(\mathbf{x}_k) + c_1 \alpha_k \nabla f(\mathbf{x}_k)^T \mathbf{d}_k$$
この条件は、実際の関数の減少量が、線形予測による減少量の少なくとも一部（$c_1$ の割合）を達成することを要求します。これにより、各ステップで十分な進歩があることを保証します。
3.3 Wolfe条件（Wolfe Conditions）
Wolfe条件はステップサイズが十分に大きいことを保証します：
$$\nabla f(\mathbf{x}_k + \alpha_k \mathbf{d}_k)^T \mathbf{d}_k \geq c_2 \nabla f(\mathbf{x}_k)^T \mathbf{d}_k$$
ここで $c_2 \in (c_1, 1)$ です。この条件は、ステップの終点での方向微分（directional derivative）が、初期点での方向微分よりもフラットになることを要求します。これにより、ステップサイズが小さすぎるのを防ぎます。
3.4 バックトラッキングラインサーチ（Backtracking Line Search）の実装
実装が最も簡単なバックトラッキングラインサーチを示します：
@jit
def armijo_line_search(f, x, d, g, alpha=1.0, c1=1e-4, rho=0.5):
f_x = f(x)
slope = jnp.dot(g, d) # 方向微分
def cond_fun(state):
alpha, _ = state
f_new = f(x + alpha * d)
return f_new > f_x + c1 * alpha * slope
def body_fun(state):
alpha, _ = state
return alpha * rho, None # αを縮小
alpha, _ = jax.lax.while_loop(cond_fun, body_fun, (alpha, None))
return alpha
このアルゴリズムでは、ステップサイズを最初は1として、Armijo条件を満たすまで繰り返し縮小します。
4. 完全な修正ニュートン法実装
以上の要素を組み合わせた実装です。固有値修正によるヘッセ行列の正定値化とArmijoラインサーチを組み合わせています：
@jit
def modified_newton_step(f, x, reg=1e-4):
g = grad(f)(x)
H = hessian(f)(x)
# 固有値分解と修正
eigvals, eigvecs = jnp.linalg.eigh(H)
modified_eigvals = jnp.maximum(eigvals, reg)
H_inv = eigvecs @ jnp.diag(1.0 / modified_eigvals) @ eigvecs.T
d = -H_inv @ g
return d, g
@jit
def modified_newton_method(f, x0, tol=1e-6, max_iter=100, reg=1e-4):
def cond_fun(state):
_, _, iter_num, g_norm, converged = state
return (iter_num < max_iter) & (~converged)
def body_fun(state):
x, _, iter_num, _, _ = state
# 探索方向計算
d, g = modified_newton_step(f, x, reg)
# ラインサーチでステップサイズ決定
alpha = armijo_line_search(f, x, d, g)
# 更新
x_new = x + alpha * d
# 収束判定
g_new = grad(f)(x_new)
g_norm = jnp.linalg.norm(g_new)
converged = g_norm < tol
return x_new, x, iter_num + 1, g_norm, converged
# 初期状態
g0 = grad(f)(x0)
g0_norm = jnp.linalg.norm(g0)
init_state = (x0, x0, 0, g0_norm, g0_norm < tol)
# 反復
final_state = jax.lax.while_loop(cond_fun, body_fun, init_state)
x_opt, _, iter_num, _, converged = final_state
return x_opt, iter_num, converged
5. 収束判定基準（Convergence Criteria）
収束判定には様々な基準があり、それぞれにメリットとデメリットがあります：
5.1 勾配ノルムによる判定（Gradient Norm Criterion）
最適性の必要条件に基づく理論的に妥当な基準です：
$$|\nabla f(\mathbf{x}_k)| < \epsilon_g$$
5.2 変位による判定（Step Size Criterion）
連続する反復間の変化が小さくなったことで収束を判断します：
$$|\mathbf{x}_{k+1} - \mathbf{x}_k| < \epsilon_x$$
5.3 関数値変化による判定（Function Value Criterion）
目的関数の相対的な変化が小さくなったかを判断します：
$$\frac{|f(\mathbf{x}_{k+1}) - f(\mathbf{x}_k)|}{\max(|f(\mathbf{x}_k)|, 1)} < \epsilon_f$$
実際には、これらの基準を組み合わせて使うことが望ましいでしょう。特に勾配ノルムは理論的に最も正当な基準であり、最適解の必要条件を直接チェックします。
6. 理論と実践
6.1 収束速度（Convergence Rate）
修正ニュートン法は、条件が整うと局所的な2次収束性（quadratic convergence）を持ちます：
$$|\mathbf{x}_{k+1} - \mathbf{x}^| \leq C |\mathbf{x}_k - \mathbf{x}^|^2, \exists C > 0$$
この2次収束性は、最適解に十分近く、ヘッセ行列が最適解の近傍で正定値かつリプシッツ連続（Lipschitz continuous）である場合に保証されます。
6.2 計算コスト（Computational Cost）
ニュートン法の主な計算コストは：
- ヘッセ行列計算: $O(n^2)$
- 行列分解/逆行列: $O(n^3)$
次元 $n$ が大きい場合、このコストは非常に高くなります。
6.3 大規模問題への対応
大規模問題では、以下のアプローチが考えられます：
- 準ニュートン法（Quasi-Newton Methods、BFGS, L-BFGS） - ヘッセ行列を低コストで近似
- Hessian-Free 法（Hessian-Free Methods） - ヘッセ行列を明示的に計算せずにベクトル積だけを利用
- 部分空間法（Subspace Methods） - 低次元の部分空間で最適化を行う