ニュートン法と修正ニュートン法


## 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実装

```python
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}$ について、以下の条件は同値です：

1. $\mathbf{A}$ は正定値である（positive definite、$\mathbf{A} \succ 0$）
2. すべての非ゼロベクトル $\mathbf{v}$ に対して $\mathbf{v}^T\mathbf{A}\mathbf{v} > 0$
3. $\mathbf{A}$ のすべての固有値（eigenvalues）が正である
4. $\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）

コレスキー分解を利用した修正法では、以下のステップで行列を正定値化します：

1. ヘッセ行列 $\mathbf{H}$ のコレスキー分解を試みます
2. 分解に失敗した場合（正定値でない）、$\mathbf{H} + \mu\mathbf{I}$ の形で修正します
3. コレスキー分解が成功するまで $\mu$ を調整します

この方法のメリットは、最小限の修正でヘッセ行列を正定値化できることです。$\mu$ が小さければ元の二次収束性を近似的に保持できます。

```python
@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$）になります。これは元の行列の構造を尊重しつつ、必要最小限の修正を加える方法です。

```python
@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）の実装

実装が最も簡単なバックトラッキングラインサーチを示します：

```python
@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ラインサーチを組み合わせています：

```python
@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 大規模問題への対応

大規模問題では、以下のアプローチが考えられます：

1. 準ニュートン法（Quasi-Newton Methods、BFGS, L-BFGS） - ヘッセ行列を低コストで近似
2. Hessian-Free 法（Hessian-Free Methods） - ヘッセ行列を明示的に計算せずにベクトル積だけを利用
3. 部分空間法（Subspace Methods） - 低次元の部分空間で最適化を行う