共役勾配法（Conjugate Gradient Method）

概要

共役勾配法（CG法）は、対称正定値行列を係数とする線形システムや二次関数の最小化に特化した高効率な反復手法です。最急降下法の欠点であるジグザグ運動を解決し、理論的には $n$ 変数の二次関数に対して $n$ 回の反復で正確解に到達します。

基本概念

共役性（Conjugacy）

行列 $A$ に関して、ベクトル $d_i$ と $d_j$（$i \neq j$）が $A$-共役であるとは：

$$d_i^T A d_j = 0$$

が成り立つことです。これは、$A$ によって定義される内積に関して直交している状態を意味します。

二次関数の最小化問題

$$\min_{x \in \mathbb{R}^n} f(x) = \frac{1}{2}x^T A x - b^T x + c$$

ここで、$A$ は対称正定値行列、$b \in \mathbb{R}^n$、$c \in \mathbb{R}$ です。

勾配は：

$$\nabla f(x) = Ax - b$$

最適解では $\nabla f(x^*) = 0$、すなわち $Ax^* = b$ が成り立ちます。

Fletcher-Reeves共役勾配法

アルゴリズム

Algorithm: Fletcher-Reeves Conjugate Gradient Method
1. 初期点 x₀ ∈ ℝⁿ を選択
2. r₀ = b - Ax₀ (残差ベクトル)
3. d₀ = r₀ (初期探索方向)
4. for k = 0, 1, 2, ... do:
   5. 停止条件: ‖rₖ‖ ≤ ε ならば停止
   6. ステップサイズを計算: αₖ = (rₖᵀrₖ)/(dₖᵀAdₖ)
   7. 解を更新: xₖ₊₁ = xₖ + αₖdₖ
   8. 残差を更新: rₖ₊₁ = rₖ - αₖAdₖ
   9. 共役方向パラメータ: βₖ₊₁ = (rₖ₊₁ᵀrₖ₊₁)/(rₖᵀrₖ)
   10. 新探索方向: dₖ₊₁ = rₖ₊₁ + βₖ₊₁dₖ
   11. k := k + 1

理論的性質

1. 有限終了性: $n$ 次元問題に対して最大 $n$ 回の反復で正確解に到達
2. 共役性の保持: 生成される探索方向 $\{d_0, d_1, \ldots, d_{k-1}\}$ は $A$-共役
3. 残差の直交性: $r_i^T r_j = 0$ （$i \neq j$）

Polak-Ribière共役勾配法

非二次関数に対してより良い性能を示すバリエーションです：

$$\beta_k^{PR} = \frac{r_k^T (r_k - r_{k-1})}{r_{k-1}^T r_{k-1}}$$

実装例

Python実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags

def conjugate_gradient(A, b, x0=None, max_iter=None, tol=1e-6):
    """
    共役勾配法による線形システム Ax = b の解法
    
    Parameters:
    A: 対称正定値行列 (n×n)
    b: 右辺ベクトル (n×1)
    x0: 初期解 (デフォルト: 零ベクトル)
    max_iter: 最大反復回数 (デフォルト: n)
    tol: 許容誤差
    
    Returns:
    x: 解ベクトル
    history: 収束履歴
    """
    n = len(b)
    if x0 is None:
        x0 = np.zeros(n)
    if max_iter is None:
        max_iter = n
    
    x = x0.copy()
    r = b - A @ x  # 残差ベクトル
    d = r.copy()   # 探索方向
    
    history = {
        'residual_norm': [],
        'objective': [],
        'x_values': [x.copy()]
    }
    
    for k in range(max_iter):
        # 残差ノルムの記録
        r_norm = np.linalg.norm(r)
        history['residual_norm'].append(r_norm)
        
        # 目的関数値の記録
        obj_val = 0.5 * x.T @ A @ x - b.T @ x
        history['objective'].append(obj_val)
        
        # 収束判定
        if r_norm < tol:
            print(f"共役勾配法が収束しました (反復回数: {k})")
            break
        
        # ステップサイズの計算
        Ad = A @ d
        alpha = (r.T @ r) / (d.T @ Ad)
        
        # 解の更新
        x = x + alpha * d
        history['x_values'].append(x.copy())
        
        # 残差の更新
        r_new = r - alpha * Ad
        
        # 共役方向パラメータ (Fletcher-Reeves)
        beta = (r_new.T @ r_new) / (r.T @ r)
        
        # 新しい探索方向
        d = r_new + beta * d
        
        # 残差の更新
        r = r_new
    
    return x, history

def polak_ribiere_cg(A, b, x0=None, max_iter=None, tol=1e-6):
    """Polak-Ribière共役勾配法"""
    n = len(b)
    if x0 is None:
        x0 = np.zeros(n)
    if max_iter is None:
        max_iter = n
    
    x = x0.copy()
    r = b - A @ x
    d = r.copy()
    r_old = r.copy()
    
    history = {'residual_norm': [], 'x_values': [x.copy()]}
    
    for k in range(max_iter):
        r_norm = np.linalg.norm(r)
        history['residual_norm'].append(r_norm)
        
        if r_norm < tol:
            print(f"Polak-Ribière CGが収束しました (反復回数: {k})")
            break
        
        Ad = A @ d
        alpha = (r.T @ r) / (d.T @ Ad)
        
        x = x + alpha * d
        history['x_values'].append(x.copy())
        
        r_new = r - alpha * Ad
        
        # Polak-Ribière公式
        beta = (r_new.T @ (r_new - r)) / (r.T @ r)
        beta = max(0, beta)  # 非負性の保証
        
        d = r_new + beta * d
        r = r_new
    
    return x, history

def cg_for_optimization(f, grad_f, x0, max_iter=1000, tol=1e-6, restart=None):
    """
    一般的な最適化問題に対する共役勾配法
    (非線形共役勾配法)
    """
    if restart is None:
        restart = len(x0)
    
    x = x0.copy()
    g = grad_f(x)
    d = -g.copy()
    
    history = {'f_values': [], 'grad_norms': [], 'x_values': [x.copy()]}
    
    for k in range(max_iter):
        grad_norm = np.linalg.norm(g)
        history['f_values'].append(f(x))
        history['grad_norms'].append(grad_norm)
        
        if grad_norm < tol:
            print(f"非線形CGが収束しました (反復回数: {k})")
            break
        
        # 直線探索 (簡単なバックトラッキング)
        alpha = backtracking_line_search(f, grad_f, x, d)
        
        # 解の更新
        x_new = x + alpha * d
        g_new = grad_f(x_new)
        
        # 再開条件
        if k % restart == 0:
            d = -g_new
        else:
            # Polak-Ribière公式
            beta = max(0, g_new.T @ (g_new - g) / (g.T @ g))
            d = -g_new + beta * d
        
        x = x_new
        g = g_new
        history['x_values'].append(x.copy())
    
    return x, history

def backtracking_line_search(f, grad_f, x, d, alpha_init=1.0, c1=1e-4, rho=0.5):
    """バックトラッキング直線探索"""
    alpha = alpha_init
    f_x = f(x)
    grad_x = grad_f(x)
    
    while f(x + alpha * d) > f_x + c1 * alpha * np.dot(grad_x, d):
        alpha *= rho
    
    return alpha

# 使用例1: 線形システムの解法
def example_linear_system():
    """対称正定値行列の線形システム"""
    n = 100
    # 対角優位行列の生成
    A = diags([1, -2, 1], [-1, 0, 1], shape=(n, n)).toarray()
    A = A.T @ A + np.eye(n)  # 正定値化
    
    # 真の解
    x_true = np.random.randn(n)
    b = A @ x_true
    
    # 初期推定値
    x0 = np.zeros(n)
    
    # 共役勾配法による解法
    x_cg, hist_cg = conjugate_gradient(A, b, x0)
    
    # 結果の比較
    print(f"真の解との誤差: {np.linalg.norm(x_cg - x_true):.2e}")
    print(f"反復回数: {len(hist_cg['residual_norm'])}")
    
    # 収束履歴の可視化
    plt.figure(figsize=(12, 4))
    
    plt.subplot(1, 3, 1)
    plt.semilogy(hist_cg['residual_norm'])
    plt.xlabel('反復回数')
    plt.ylabel('残差ノルム')
    plt.title('収束履歴')
    plt.grid(True)
    
    plt.subplot(1, 3, 2)
    plt.plot(hist_cg['objective'])
    plt.xlabel('反復回数')
    plt.ylabel('目的関数値')
    plt.title('目的関数の減少')
    plt.grid(True)
    
    plt.subplot(1, 3, 3)
    errors = [np.linalg.norm(x - x_true) for x in hist_cg['x_values']]
    plt.semilogy(errors)
    plt.xlabel('反復回数')
    plt.ylabel('真の解との誤差')
    plt.title('解の精度')
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

# 使用例2: 二次関数の最小化
def example_quadratic_optimization():
    """二次関数の最小化"""
    # 二次関数 f(x) = 1/2 * x^T * A * x - b^T * x
    A = np.array([[4, 1], [1, 3]])
    b = np.array([1, 2])
    
    def f(x):
        return 0.5 * x.T @ A @ x - b.T @ x
    
    def grad_f(x):
        return A @ x - b
    
    x0 = np.array([0.0, 0.0])
    
    # 線形システムとしての解法
    x_exact = np.linalg.solve(A, b)
    
    # 共役勾配法
    x_cg, hist_cg = conjugate_gradient(A, b, x0)
    
    print(f"正確解: {x_exact}")
    print(f"CG解: {x_cg}")
    print(f"誤差: {np.linalg.norm(x_cg - x_exact):.2e}")
    
    # 等高線図での可視化
    plt.figure(figsize=(10, 8))
    
    # 等高線の描画
    x1_range = np.linspace(-1, 1, 100)
    x2_range = np.linspace(-1, 1, 100)
    X1, X2 = np.meshgrid(x1_range, x2_range)
    Z = np.zeros_like(X1)
    
    for i in range(len(x1_range)):
        for j in range(len(x2_range)):
            x = np.array([X1[i,j], X2[i,j]])
            Z[i,j] = f(x)
    
    plt.contour(X1, X2, Z, levels=20, alpha=0.6)
    
    # 最急降下法との比較のため、最急降下法も実行
    x_sd, hist_sd = steepest_descent_simple(f, grad_f, x0, max_iter=50)
    
    # 収束経路の描画
    cg_path = np.array(hist_cg['x_values'])
    sd_path = np.array(hist_sd['x_values'])
    
    plt.plot(cg_path[:, 0], cg_path[:, 1], 'ro-', label='共役勾配法', linewidth=2)
    plt.plot(sd_path[:, 0], sd_path[:, 1], 'bs-', label='最急降下法', alpha=0.7)
    plt.plot(x_exact[0], x_exact[1], 'g*', markersize=15, label='最適解')
    
    plt.xlabel('x₁')
    plt.ylabel('x₂')
    plt.title('収束経路の比較')
    plt.legend()
    plt.grid(True)
    plt.axis('equal')
    plt.show()

def steepest_descent_simple(f, grad_f, x0, max_iter=1000, tol=1e-6):
    """簡単な最急降下法（比較用）"""
    x = x0.copy()
    history = {'x_values': [x.copy()]}
    
    for k in range(max_iter):
        grad = grad_f(x)
        if np.linalg.norm(grad) < tol:
            break
        
        alpha = backtracking_line_search(f, grad_f, x, -grad)
        x = x - alpha * grad
        history['x_values'].append(x.copy())
    
    return x, history

# 実行
if __name__ == "__main__":
    print("=== 線形システムの例 ===")
    example_linear_system()
    
    print("\n=== 二次関数最適化の例 ===")
    example_quadratic_optimization()

MATLAB実装

function [x, history] = conjugate_gradient(A, b, x0, options)
    % 共役勾配法による線形システム Ax = b の解法
    
    if nargin < 3, x0 = zeros(size(b)); end
    if nargin < 4, options = struct(); end
    
    % デフォルトパラメータ
    max_iter = getfield_default(options, 'max_iter', length(b));
    tol = getfield_default(options, 'tol', 1e-6);
    
    x = x0;
    r = b - A * x;  % 残差
    d = r;          % 探索方向
    
    history.residual_norm = [];
    history.x_values = x0;
    
    for k = 1:max_iter
        r_norm = norm(r);
        history.residual_norm(end+1) = r_norm;
        
        % 収束判定
        if r_norm < tol
            fprintf('共役勾配法が収束しました (反復回数: %d)\n', k);
            break;
        end
        
        % ステップサイズ
        Ad = A * d;
        alpha = (r' * r) / (d' * Ad);
        
        % 解の更新
        x = x + alpha * d;
        history.x_values = [history.x_values, x];
        
        % 残差の更新
        r_new = r - alpha * Ad;
        
        % 共役方向パラメータ
        beta = (r_new' * r_new) / (r' * r);
        
        % 新しい探索方向
        d = r_new + beta * d;
        r = r_new;
    end
end

function [x, history] = nonlinear_cg(f, grad_f, x0, options)
    % 非線形最適化のための共役勾配法
    
    if nargin < 4, options = struct(); end
    
    max_iter = getfield_default(options, 'max_iter', 1000);
    tol = getfield_default(options, 'tol', 1e-6);
    restart = getfield_default(options, 'restart', length(x0));
    
    x = x0;
    g = grad_f(x);
    d = -g;
    
    history.f_values = [];
    history.grad_norms = [];
    history.x_values = x0;
    
    for k = 1:max_iter
        grad_norm = norm(g);
        history.f_values(end+1) = f(x);
        history.grad_norms(end+1) = grad_norm;
        
        if grad_norm < tol
            fprintf('非線形CGが収束しました (反復回数: %d)\n', k);
            break;
        end
        
        % 直線探索
        alpha = line_search(f, grad_f, x, d);
        
        % 解の更新
        x_new = x + alpha * d;
        g_new = grad_f(x_new);
        
        % 再開条件
        if mod(k, restart) == 0
            d = -g_new;
        else
            % Polak-Ribière公式
            beta = max(0, g_new' * (g_new - g) / (g' * g));
            d = -g_new + beta * d;
        end
        
        x = x_new;
        g = g_new;
        history.x_values = [history.x_values, x];
    end
end

function alpha = line_search(f, grad_f, x, d)
    % 簡単なバックトラッキング直線探索
    alpha = 1.0;
    c1 = 1e-4;
    rho = 0.5;
    
    f_x = f(x);
    grad_x = grad_f(x);
    
    while f(x + alpha * d) > f_x + c1 * alpha * (grad_x' * d)
        alpha = rho * alpha;
    end
end

function val = getfield_default(s, field, default)
    if isfield(s, field)
        val = s.(field);
    else
        val = default;
    end
end

理論的性質と証明のスケッチ

有限終了性

定理: $n$ 次元の二次関数に対して、共役勾配法は最大 $n$ 回の反復で正確解に到達する。

証明のスケッチ:

1. 生成される探索方向 $\{d_0, d_1, \ldots, d_{k-1}\}$ は $A$-共役
2. これらの方向は線形独立であり、最大 $n$ 個まで生成可能
3. $k$ 回目の反復後、解 $x_k$ は初期点 $x_0$ から $\text{span}\{d_0, \ldots, d_{k-1}\}$ 内での最適点
4. $k = n$ のとき、全空間での最適化となり、正確解に到達

共役性の保持

各反復で生成される探索方向は、過去のすべての方向と $A$-共役関係を保持します：

$$d_i^T A d_j = 0, \quad \text{for } i \neq j$$

残差の直交性

残差ベクトル $r_k = b - Ax_k$ は相互に直交します：

$$r_i^T r_j = 0, \quad \text{for } i \neq j$$

条件数と収束速度

収束率の評価

対称正定値行列 $A$ の条件数を $\kappa(A) = \lambda_{\max}/\lambda_{\min}$ とすると、共役勾配法の収束率は：

$$\|x_k - x^*\|_A \leq 2\left(\frac{\sqrt{\kappa(A)} - 1}{\sqrt{\kappa(A)} + 1}\right)^k \|x_0 - x^*\|_A$$

ここで、$\|x\|_A = \sqrt{x^T A x}$ は $A$-ノルムです。

実用的な収束予測

条件数が小さい（$\kappa(A) \ll n$）場合、共役勾配法は $O(\sqrt{\kappa(A)})$ 回の反復で収束し、これは理論上限の $n$ 回よりも大幅に少なくなります。

前処理付き共役勾配法

前処理の概念

前処理行列 $M$ を用いて、元の問題：

$$Ax = b$$

を等価な問題：

$$M^{-1}Ax = M^{-1}b$$

に変換します。$M \approx A$ かつ $M^{-1}$ が計算しやすい行列を選択します。

PCG アルゴリズム

Algorithm: Preconditioned Conjugate Gradient (PCG)
1. 初期点 x₀, r₀ = b - Ax₀
2. z₀ = M⁻¹r₀, d₀ = z₀
3. for k = 0, 1, 2, ... do:
   4. αₖ = (rₖᵀzₖ)/(dₖᵀAdₖ)
   5. xₖ₊₁ = xₖ + αₖdₖ
   6. rₖ₊₁ = rₖ - αₖAdₖ
   7. zₖ₊₁ = M⁻¹rₖ₊₁
   8. βₖ₊₁ = (rₖ₊₁ᵀzₖ₊₁)/(rₖᵀzₖ)
   9. dₖ₊₁ = zₖ₊₁ + βₖ₊₁dₖ

非線形問題への拡張

Fletcher-Reeves型

$$\beta_k^{FR} = \frac{\|\nabla f(x_{k+1})\|^2}{\|\nabla f(x_k)\|^2}$$

Polak-Ribière型

$$\beta_k^{PR} = \frac{\nabla f(x_{k+1})^T (\nabla f(x_{k+1}) - \nabla f(x_k))}{\|\nabla f(x_k)\|^2}$$

Hestenes-Stiefel型

$$\beta_k^{HS} = \frac{\nabla f(x_{k+1})^T (\nabla f(x_{k+1}) - \nabla f(x_k))}{d_k^T (\nabla f(x_{k+1}) - \nabla f(x_k))}$$

実用的な考慮事項

メモリ効率

共役勾配法は以下の利点があります：

• 低メモリ使用量: 行列 $A$ を明示的に保存する必要がない（行列ベクトル積のみ）
• スパース行列対応: 大規模スパース問題に適用可能

数値的安定性

1. 再開戦略: 蓄積誤差を避けるため、定期的に $d_k = -\nabla f(x_k)$ でリセット
2. 直交性の監視: 残差の直交性 $r_i^T r_j$ をチェック
3. 条件数の改善: 前処理による条件数の改善

実装のコツ

# 数値的安定性のための改良版実装
def robust_cg(A, b, x0=None, M=None, max_iter=None, tol=1e-6, restart=None):
    """
    数値的に安定な共役勾配法
    
    Parameters:
    M: 前処理行列 (デフォルト: 恒等行列)
    restart: 再開間隔 (デフォルト: 50回)
    """
    n = len(b)
    if x0 is None:
        x0 = np.zeros(n)
    if M is None:
        M = np.eye(n)
    if max_iter is None:
        max_iter = n
    if restart is None:
        restart = min(50, n)
    
    x = x0.copy()
    r = b - A @ x
    
    for restart_iter in range(max_iter // restart + 1):
        # 再開時の初期化
        if restart_iter > 0:
            r = b - A @ x
        
        z = np.linalg.solve(M, r)
        d = z.copy()
        
        for k in range(min(restart, max_iter - restart_iter * restart)):
            if np.linalg.norm(r) < tol:
                return x, True
            
            Ad = A @ d
            alpha = (r.T @ z) / (d.T @ Ad)
            
            x = x + alpha * d
            r = r - alpha * Ad
            
            z_new = np.linalg.solve(M, r)
            beta = (r.T @ z_new) / (r.T @ z)
            
            d = z_new + beta * d
            z = z_new
    
    return x, False

実際の応用例

有限要素法での応用

def fem_example():
    """有限要素法における剛性行列の解法"""
    # 1次元熱伝導方程式の離散化
    n = 1000
    h = 1.0 / (n + 1)
    
    # 剛性行列 (三重対角行列)
    main_diag = 2 / h * np.ones(n)
    off_diag = -1 / h * np.ones(n - 1)
    A = diags([off_diag, main_diag, off_diag], [-1, 0, 1]).toarray()
    
    # 右辺 (熱源項)
    x_grid = np.linspace(h, 1 - h, n)
    b = h * np.sin(np.pi * x_grid)
    
    # 共役勾配法による解法
    import time
    start_time = time.time()
    u_cg, hist = conjugate_gradient(A, b)
    cg_time = time.time() - start_time
    
    # 直接法との比較
    start_time = time.time()
    u_direct = np.linalg.solve(A, b)
    direct_time = time.time() - start_time
    
    print(f"CG解法時間: {cg_time:.4f}秒")
    print(f"直接解法時間: {direct_time:.4f}秒")
    print(f"解の誤差: {np.linalg.norm(u_cg - u_direct):.2e}")
    
    # 解析解との比較
    u_exact = np.sin(np.pi * x_grid) / (np.pi**2)
    print(f"解析解との誤差: {np.linalg.norm(u_cg - u_exact):.2e}")

まとめ

共役勾配法の重要な特徴：

1. 理論的保証: 二次関数に対して有限回での収束
2. 計算効率: $O(n^2)$ の計算量（直接法の $O(n^3)$ より効率的）
3. メモリ効率: $O(n)$ のメモリ使用量
4. スケーラビリティ: 大規模問題への適用可能性
5. 前処理対応: 条件数改善による収束加速

実装のポイント

共役勾配法を実装する際は、数値的安定性を保つために定期的な再開、前処理の活用、および適切な収束判定が重要です。特に、残差の直交性が保たれているかを監視してください。

参考文献

1. Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. Springer.
2. Shewchuk, J. R. (1994). An introduction to the conjugate gradient method without the agonizing pain. Carnegie Mellon University.
3. Hestenes, M. R., & Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards.
4. Fletcher, R., & Reeves, C. M. (1964). Function minimization by conjugate gradients. The Computer Journal.