ペナルティ法（Penalty Methods）

概要

ペナルティ法は、制約付き最適化問題を無制約最適化問題の列に変換して解く手法です。制約違反に対してペナルティを課すことで、解が制約条件を満たすよう誘導します。バリア法とは異なり、実行可能領域の外部からも最適解に接近できる外点法です。

問題設定

一般的な制約付き最適化問題：

$$\begin{align} \min_{x \in \mathbb{R}^n} &\quad f(x) \\ \text{subject to} &\quad h_j(x) = 0, \quad j = 1, 2, \ldots, p \\ &\quad g_i(x) \leq 0, \quad i = 1, 2, \ldots, m \end{align}$$

ここで：

• $f: \mathbb{R}^n \to \mathbb{R}$ は目的関数
• $h_j: \mathbb{R}^n \to \mathbb{R}$ は等式制約関数
• $g_i: \mathbb{R}^n \to \mathbb{R}$ は不等式制約関数

ペナルティ関数

基本的なペナルティ関数

制約違反を測定するペナルティ関数：

$$P(x) = \sum_{j=1}^p |h_j(x)|^2 + \sum_{i=1}^m [\max(0, g_i(x))]^2$$

より一般的な形式

$$P(x) = \sum_{j=1}^p \phi_j(h_j(x)) + \sum_{i=1}^m \psi_i(g_i(x))$$

ここで：

• $\phi_j(t) = t^2$ （等式制約用）
• $\psi_i(t) = [\max(0, t)]^2$ （不等式制約用）

外点ペナルティ法

アルゴリズム

Algorithm: Exterior Penalty Method
1. 初期点 x₀ ∈ ℝⁿ を選択 (実行可能性は不要)
2. 初期ペナルティパラメータ r₀ > 0 を設定
3. for k = 0, 1, 2, ... do:
   4. ペナルティ問題を解く:
      xₖ₊₁ = argmin{f(x) + rₖP(x)}
   5. パラメータを更新: rₖ₊₁ = βrₖ (β > 1)
   6. 収束判定: 制約違反が十分小さければ停止

ペナルティ問題

各反復で解くべき無制約最適化問題：

$$\min_{x \in \mathbb{R}^n} F_r(x) = f(x) + r P(x)$$

ここで、$r > 0$ はペナルティパラメータです。

理論的性質

収束性定理

定理（外点ペナルティ法の収束性）: 以下の条件下で：

1. 元問題が実行可能で最適解 $x^*$ を持つ
2. LICQ（線形独立制約限定条件）が成立
3. 二次十分条件が満たされる

$r \to \infty$ のとき、ペナルティ問題の解 $x^*(r)$ は元問題の最適解に収束します。

ラグランジュ乗数との関係

ペナルティ問題の解 $x^*(r)$ において：

$$\begin{align} \lambda_j^*(r) &= r \cdot 2h_j(x^*(r)) \\ \mu_i^*(r) &= r \cdot 2\max(0, g_i(x^*(r))) \end{align}$$

$r \to \infty$ の極限で、これらは元問題のラグランジュ乗数に収束します。

実装例

Python実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

def exterior_penalty_method(f, grad_f, h=None, grad_h=None, g=None, grad_g=None,
                           x0=None, r0=1.0, beta=10.0, max_outer=20, 
                           max_inner=100, tol=1e-6):
    """
    外点ペナルティ法の実装
    
    Parameters:
    f: 目的関数
    grad_f: 目的関数の勾配
    h: 等式制約関数のリスト [h1, h2, ..., hp]
    grad_h: 等式制約関数の勾配のリスト
    g: 不等式制約関数のリスト [g1, g2, ..., gm]
    grad_g: 不等式制約関数の勾配のリスト
    x0: 初期点
    r0: 初期ペナルティパラメータ
    beta: rの増加率 (beta > 1)
    """
    
    if h is None:
        h = []
        grad_h = []
    if g is None:
        g = []
        grad_g = []
    
    def penalty_function(x):
        """ペナルティ関数"""
        penalty = 0.0
        
        # 等式制約のペナルティ
        for hj in h:
            penalty += hj(x)**2
        
        # 不等式制約のペナルティ
        for gi in g:
            penalty += max(0, gi(x))**2
        
        return penalty
    
    def penalty_gradient(x):
        """ペナルティ関数の勾配"""
        grad = np.zeros_like(x)
        
        # 等式制約の寄与
        for hj, grad_hj in zip(h, grad_h):
            grad += 2 * hj(x) * grad_hj(x)
        
        # 不等式制約の寄与
        for gi, grad_gi in zip(g, grad_g):
            if gi(x) > 0:
                grad += 2 * gi(x) * grad_gi(x)
        
        return grad
    
    def penalty_objective(x, r):
        """ペナルティ目的関数"""
        return f(x) + r * penalty_function(x)
    
    def penalty_objective_grad(x, r):
        """ペナルティ目的関数の勾配"""
        return grad_f(x) + r * penalty_gradient(x)
    
    # 初期化
    if x0 is None:
        x0 = np.zeros(2)  # デフォルトの次元
    
    x = x0.copy()
    r = r0
    
    history = {
        'x_values': [x.copy()],
        'f_values': [f(x)],
        'r_values': [r],
        'penalty_values': [penalty_function(x)],
        'constraint_violations': []
    }
    
    def constraint_violation(x):
        """制約違反の測定"""
        violation = 0.0
        for hj in h:
            violation += abs(hj(x))
        for gi in g:
            violation += max(0, gi(x))
        return violation
    
    print(f"初期点: {x}, f(x0) = {f(x):.6f}")
    print(f"初期制約違反: {constraint_violation(x):.6f}")
    
    for outer_iter in range(max_outer):
        print(f"\n外側反復 {outer_iter}: r = {r:.2f}")
        
        # ペナルティ問題を解く
        def objective(x_var):
            return penalty_objective(x_var, r)
        
        def objective_grad(x_var):
            return penalty_objective_grad(x_var, r)
        
        try:
            # 無制約最適化
            result = minimize(objective, x, method='BFGS', jac=objective_grad,
                            options={'maxiter': max_inner, 'gtol': tol})
            
            x = result.x
            
            # 履歴の更新
            history['x_values'].append(x.copy())
            history['f_values'].append(f(x))
            history['r_values'].append(r)
            history['penalty_values'].append(penalty_function(x))
            
            cv = constraint_violation(x)
            history['constraint_violations'].append(cv)
            
            print(f"ペナルティ問題の解: {x}")
            print(f"目的関数値: {f(x):.6f}")
            print(f"制約違反: {cv:.6f}")
            
            if len(h) > 0:
                print(f"等式制約値: {[hj(x) for hj in h]}")
            if len(g) > 0:
                print(f"不等式制約値: {[gi(x) for gi in g]}")
            
            # 収束判定
            if cv < tol:
                print(f"\n収束しました (制約違反 = {cv:.2e})")
                break
            
            # rの更新
            r *= beta
            
        except Exception as e:
            print(f"最適化でエラーが発生: {e}")
            break
    
    return x, history

def quadratic_penalty_method(f, grad_f, h, grad_h, g, grad_g, x0, r0=1.0, 
                            beta=10.0, max_outer=20, tol=1e-6):
    """二次ペナルティ法（より一般的な実装）"""
    
    def penalty_function(x):
        penalty = 0.0
        
        # 等式制約: Σ h_j(x)²
        for hj in h:
            penalty += hj(x)**2
        
        # 不等式制約: Σ [max(0, g_i(x))]²
        for gi in g:
            penalty += max(0, gi(x))**2
        
        return penalty
    
    def penalty_gradient(x):
        grad = np.zeros_like(x)
        
        # 等式制約の勾配
        for hj, grad_hj in zip(h, grad_h):
            grad += 2 * hj(x) * grad_hj(x)
        
        # 不等式制約の勾配
        for gi, grad_gi in zip(g, grad_g):
            if gi(x) > 0:
                grad += 2 * gi(x) * grad_gi(x)
        
        return grad
    
    # 以下の実装は exterior_penalty_method と同様
    # ... (省略)

def example_equality_constrained():
    """等式制約付き最適化の例"""
    
    # 問題: min f(x) = (x1-1)² + (x2-2)²
    # subject to: h(x) = x1² + x2² - 1 = 0
    
    def f(x):
        return (x[0] - 1)**2 + (x[1] - 2)**2
    
    def grad_f(x):
        return np.array([2*(x[0] - 1), 2*(x[1] - 2)])
    
    def h1(x):
        return x[0]**2 + x[1]**2 - 1
    
    def grad_h1(x):
        return np.array([2*x[0], 2*x[1]])
    
    h = [h1]
    grad_h = [grad_h1]
    
    # 実行不可能な初期点から開始
    x0 = np.array([0.5, 0.5])
    
    # ペナルティ法を実行
    x_opt, history = exterior_penalty_method(f, grad_f, h=h, grad_h=grad_h, x0=x0)
    
    print(f"\n最適解: {x_opt}")
    print(f"最適値: {f(x_opt):.6f}")
    print(f"制約値: {h1(x_opt):.6f}")
    
    # 解析的最適解（ラグランジュ乗数法）
    # 制約円上での最小点
    direction = np.array([1, 2]) - np.array([0, 0])
    direction = direction / np.linalg.norm(direction)
    x_analytical = direction  # 単位円上の点
    
    print(f"解析解: {x_analytical}")
    print(f"解析解での目的関数値: {f(x_analytical):.6f}")
    print(f"誤差: {np.linalg.norm(x_opt - x_analytical):.6f}")
    
    # 結果の可視化
    plt.figure(figsize=(15, 5))
    
    # 収束経路
    plt.subplot(1, 3, 1)
    x_history = np.array(history['x_values'])
    
    # 等高線の描画
    x1_range = np.linspace(-1.5, 1.5, 100)
    x2_range = np.linspace(-1.5, 2.5, 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)):
            Z[i,j] = f([X1[i,j], X2[i,j]])
    
    plt.contour(X1, X2, Z, levels=20, alpha=0.6)
    
    # 制約円の描画
    theta = np.linspace(0, 2*np.pi, 100)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    plt.plot(circle_x, circle_y, 'r-', linewidth=2, label='制約: x₁² + x₂² = 1')
    
    # 収束経路
    plt.plot(x_history[:, 0], x_history[:, 1], 'bo-', linewidth=2, markersize=6)
    plt.plot(x_history[0, 0], x_history[0, 1], 'go', markersize=10, label='開始点')
    plt.plot(x_history[-1, 0], x_history[-1, 1], 'ro', markersize=10, label='最適解')
    plt.plot(x_analytical[0], x_analytical[1], 'k*', markersize=15, label='解析解')
    
    plt.xlabel('x₁')
    plt.ylabel('x₂')
    plt.title('ペナルティ法の収束経路')
    plt.legend()
    plt.grid(True)
    plt.axis('equal')
    
    # 目的関数値と制約違反
    plt.subplot(1, 3, 2)
    plt.plot(history['f_values'], 'b-o', label='目的関数値')
    plt.axhline(y=f(x_analytical), color='k', linestyle='--', label='解析解の値')
    plt.xlabel('外側反復回数')
    plt.ylabel('目的関数値')
    plt.title('目的関数値の収束')
    plt.legend()
    plt.grid(True)
    
    plt.subplot(1, 3, 3)
    plt.semilogy(history['constraint_violations'][1:], 'r-o', label='制約違反')
    plt.semilogy(history['r_values'][1:], 'b-s', label='ペナルティパラメータ')
    plt.xlabel('外側反復回数')
    plt.ylabel('値 (対数スケール)')
    plt.title('制約違反とペナルティパラメータ')
    plt.legend()
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    return x_opt, history

def example_mixed_constraints():
    """等式・不等式混合制約の例"""
    
    # 問題: min f(x) = x1² + x2²
    # subject to: h(x) = x1 + x2 - 1 = 0
    #            g(x) = -x1 ≤ 0 (i.e., x1 ≥ 0)
    
    def f(x):
        return x[0]**2 + x[1]**2
    
    def grad_f(x):
        return np.array([2*x[0], 2*x[1]])
    
    def h1(x):
        return x[0] + x[1] - 1
    
    def grad_h1(x):
        return np.array([1, 1])
    
    def g1(x):
        return -x[0]  # x1 ≥ 0
    
    def grad_g1(x):
        return np.array([-1, 0])
    
    h = [h1]
    grad_h = [grad_h1]
    g = [g1]
    grad_g = [grad_g1]
    
    # 制約を違反する初期点
    x0 = np.array([-0.5, 0.8])
    
    x_opt, history = exterior_penalty_method(f, grad_f, h=h, grad_h=grad_h, 
                                           g=g, grad_g=grad_g, x0=x0)
    
    print(f"\n混合制約問題の最適解: {x_opt}")
    print(f"最適値: {f(x_opt):.6f}")
    print(f"等式制約値: {h1(x_opt):.6f}")
    print(f"不等式制約値: {g1(x_opt):.6f}")
    
    # 解析解: x1 = 0.5, x2 = 0.5 (制約 x1 + x2 = 1, x1 ≥ 0 下での最小値)
    x_analytical = np.array([0.5, 0.5])
    print(f"解析解: {x_analytical}")
    print(f"誤差: {np.linalg.norm(x_opt - x_analytical):.6f}")
    
    return x_opt, history

def l1_penalty_method(f, grad_f, h, grad_h, g, grad_g, x0, r0=1.0, beta=10.0):
    """L1ペナルティ法（非微分可能）"""
    
    def l1_penalty_function(x):
        penalty = 0.0
        
        # 等式制約: Σ |h_j(x)|
        for hj in h:
            penalty += abs(hj(x))
        
        # 不等式制約: Σ max(0, g_i(x))
        for gi in g:
            penalty += max(0, gi(x))
        
        return penalty
    
    def l1_penalty_objective(x, r):
        return f(x) + r * l1_penalty_function(x)
    
    # L1ペナルティは非微分可能なので、微分不要な最適化手法を使用
    x = x0.copy()
    r = r0
    
    history = {'x_values': [x.copy()], 'f_values': [f(x)]}
    
    for outer_iter in range(20):
        # Nelder-Mead法など微分不要な手法を使用
        def objective(x_var):
            return l1_penalty_objective(x_var, r)
        
        result = minimize(objective, x, method='Nelder-Mead')
        x = result.x
        
        history['x_values'].append(x.copy())
        history['f_values'].append(f(x))
        
        print(f"L1反復 {outer_iter}: r = {r:.2f}, x = {x}, f(x) = {f(x):.6f}")
        
        # 制約違反のチェック
        violation = sum(abs(hj(x)) for hj in h) + sum(max(0, gi(x)) for gi in g)
        if violation < 1e-6:
            break
        
        r *= beta
    
    return x, history

# 実行例
if __name__ == "__main__":
    print("=== 等式制約付き最適化の例 ===")
    x_opt1, history1 = example_equality_constrained()
    
    print("\n=== 混合制約問題の例 ===")
    x_opt2, history2 = example_mixed_constraints()

MATLAB実装

function [x, history] = exterior_penalty_method(f, grad_f, h, grad_h, g, grad_g, x0, options)
    % 外点ペナルティ法の実装
    
    if nargin < 8, options = struct(); end
    
    % デフォルトパラメータ
    r0 = getfield_default(options, 'r0', 1.0);
    beta = getfield_default(options, 'beta', 10.0);
    max_outer = getfield_default(options, 'max_outer', 20);
    max_inner = getfield_default(options, 'max_inner', 100);
    tol = getfield_default(options, 'tol', 1e-6);
    
    if isempty(h), h = {}; grad_h = {}; end
    if isempty(g), g = {}; grad_g = {}; end
    
    x = x0;
    r = r0;
    
    history.x_values = x0;
    history.f_values = f(x0);
    history.r_values = r0;
    
    fprintf('初期点: [%.4f, %.4f], f(x0) = %.6f\n', x(1), x(2), f(x));
    
    for outer_iter = 1:max_outer
        fprintf('\n外側反復 %d: r = %.2f\n', outer_iter, r);
        
        % ペナルティ問題の定義
        penalty_obj = @(x_var) penalty_objective(x_var, f, h, g, r);
        penalty_grad = @(x_var) penalty_gradient(x_var, grad_f, h, grad_h, g, grad_g, r);
        
        % 無制約最適化
        x_prev = x;
        for inner_iter = 1:max_inner
            grad = penalty_grad(x);
            
            if norm(grad) < tol
                break;
            end
            
            % 簡単な勾配法
            alpha = backtrack_line_search(penalty_obj, penalty_grad, x, -grad);
            x = x - alpha * grad;
        end
        
        % 履歴の更新
        history.x_values = [history.x_values, x];
        history.f_values = [history.f_values, f(x)];
        history.r_values = [history.r_values, r];
        
        fprintf('ペナルティ問題の解: [%.6f, %.6f]\n', x(1), x(2));
        fprintf('目的関数値: %.6f\n', f(x));
        
        % 制約違反の計算
        violation = constraint_violation(x, h, g);
        fprintf('制約違反: %.6f\n', violation);
        
        % 収束判定
        if violation < tol
            fprintf('\n収束しました (制約違反 = %.2e)\n', violation);
            break;
        end
        
        % rの更新
        r = beta * r;
    end
end

function obj = penalty_objective(x, f, h, g, r)
    obj = f(x);
    
    % 等式制約のペナルティ
    for i = 1:length(h)
        obj = obj + r * h{i}(x)^2;
    end
    
    % 不等式制約のペナルティ
    for i = 1:length(g)
        obj = obj + r * max(0, g{i}(x))^2;
    end
end

function grad = penalty_gradient(x, grad_f, h, grad_h, g, grad_g, r)
    grad = grad_f(x);
    
    % 等式制約の寄与
    for i = 1:length(h)
        grad = grad + r * 2 * h{i}(x) * grad_h{i}(x);
    end
    
    % 不等式制約の寄与
    for i = 1:length(g)
        if g{i}(x) > 0
            grad = grad + r * 2 * g{i}(x) * grad_g{i}(x);
        end
    end
end

function violation = constraint_violation(x, h, g)
    violation = 0;
    
    for i = 1:length(h)
        violation = violation + abs(h{i}(x));
    end
    
    for i = 1:length(g)
        violation = violation + max(0, g{i}(x));
    end
end

function alpha = backtrack_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;
        if alpha < 1e-10
            break;
        end
    end
end

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

ペナルティ法の変形

L1ペナルティ法

等式制約に対して絶対値ペナルティを使用：

$$P(x) = \sum_{j=1}^p |h_j(x)| + \sum_{i=1}^m \max(0, g_i(x))$$

利点: 有限のペナルティパラメータで正確解が得られる場合がある 欠点: 非微分可能なため、特殊な最適化手法が必要

拡張ラグランジュ法（Augmented Lagrangian Method）

ペナルティ項とラグランジュ項を組み合わせ：

$$L_r(x, \lambda, \mu) = f(x) + \sum_{j=1}^p \left(\lambda_j h_j(x) + \frac{r}{2} h_j(x)^2\right) + \sum_{i=1}^m \left(\mu_i g_i(x) + \frac{r}{2} [\max(0, g_i(x))]^2\right)$$

適応的ペナルティ法

ペナルティパラメータを制約違反に応じて動的に調整：

def adaptive_penalty_update(r, constraint_violation, prev_violation, 
                          beta_increase=10.0, beta_decrease=0.5, 
                          tau_increase=0.9, tau_decrease=0.1):
    """適応的ペナルティパラメータ更新"""
    
    if constraint_violation > tau_increase * prev_violation:
        # 制約違反が十分減少していない
        return r * beta_increase
    elif constraint_violation < tau_decrease * prev_violation:
        # 制約違反が大幅に減少
        return max(r * beta_decrease, 1e-6)
    else:
        # 現状維持
        return r

実用的な考慮事項

数値的問題

1. 条件数の悪化: $r$ が大きくなると、ヘッセ行列の条件数が悪化
2. 丸め誤差: 大きなペナルティパラメータでの計算精度の低下
3. 収束の遅化: 高精度解に到達するのに多くの反復が必要

パラメータ選択

初期ペナルティパラメータ $r_0$

• 小さすぎる場合: 制約違反が大きいまま
• 大きすぎる場合: 数値的困難
• 推奨範囲: $r_0 = 1 \sim 100$

増加率 $\beta$

• 推奨範囲: $\beta = 5 \sim 20$
• 適応的調整: 制約違反の減少率に応じて動的変更

収束判定

複数の基準を組み合わせ：

1. 制約違反: $\|h(x)\| + \|\max(0, g(x))\| < \epsilon_c$
2. 目的関数: $|f(x_{k+1}) - f(x_k)| < \epsilon_f$
3. 勾配ノルム: $\|\nabla L_r(x)\| < \epsilon_g$
4. ステップサイズ: $\|x_{k+1} - x_k\| < \epsilon_x$

他手法との比較

ペナルティ法 vs バリア法

特徴	ペナルティ法	バリア法
初期点	任意	実行可能領域内
接近経路	外部から	内部から
パラメータ変化	$r \to \infty$	$\mu \to 0$
実装の容易さ	容易	やや困難
数値的安定性	やや不安定	比較的安定

ペナルティ法 vs SQP法

• ペナルティ法: 単純だが収束が遅い
• SQP法: 複雑だが高速収束

大規模問題への応用

並列化

def parallel_penalty_subproblems(f_list, constraint_list, x0_list, r):
    """複数のペナルティ問題を並列処理"""
    from multiprocessing import Pool
    
    def solve_subproblem(args):
        f, constraints, x0 = args
        # 各サブ問題を独立に解く
        return penalty_solver(f, constraints, x0, r)
    
    with Pool() as pool:
        results = pool.map(solve_subproblem, 
                          zip(f_list, constraint_list, x0_list))
    
    return results

スパース性の活用

大規模問題では制約の構造（ブロック対角など）を活用：

def block_penalty_method(f, constraints_blocks, x0):
    """ブロック構造を持つペナルティ法"""
    
    # 各ブロックを独立に最適化
    x_blocks = []
    for i, (f_block, constraints_block, x0_block) in enumerate(constraints_blocks):
        x_block = penalty_solver(f_block, constraints_block, x0_block)
        x_blocks.append(x_block)
    
    # 全体解の構成
    return np.concatenate(x_blocks)

まとめ

ペナルティ法の主な特徴：

1. 外点アプローチ: 実行可能領域外からの接近が可能
2. 実装の簡単さ: 無制約最適化の繰り返し
3. 汎用性: 様々な制約タイプに適用可能
4. 理論的基盤: ラグランジュ乗数法との明確な関係
5. 数値的課題: 大きなペナルティパラメータでの困難

利点

• 初期点選択の自由度
• アルゴリズムの単純さ
• 既存の無制約最適化手法の利用

欠点

• 収束の遅さ
• 数値的安定性の問題
• 高精度解の困難

実装のポイント

ペナルティ法を実装する際は、ペナルティパラメータの適切な増加スケジュールと数値的安定性の確保が重要です。また、制約違反と目的関数値の両方を監視して収束を判定してください。

参考文献

1. Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. Springer.
2. Bertsekas, D. P. (2016). Nonlinear Programming. Athena Scientific.
3. Fletcher, R. (2013). Practical Methods of Optimization. John Wiley & Sons.
4. Gill, P. E., Murray, W., & Wright, M. H. (2019). Practical Optimization. SIAM.