バリア法（Barrier Methods）

概要

バリア法は、不等式制約付き最適化問題を無制約最適化問題の列に変換して解く手法です。制約条件の境界近くで急激に増加するバリア関数を目的関数に加えることで、解が実行可能領域の内部を通って最適解に近づくよう誘導します。

問題設定

不等式制約付き最適化問題：

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

ここで、$f: \mathbb{R}^n \to \mathbb{R}$ は目的関数、$g_i: \mathbb{R}^n \to \mathbb{R}$ は制約関数です。

実行可能領域：

$$\mathcal{F} = \{x \in \mathbb{R}^n : g_i(x) \leq 0, \, i = 1, \ldots, m\}$$

バリア関数

対数バリア関数

最も一般的なバリア関数は対数バリア関数です：

$$\phi(x) = -\sum_{i=1}^m \log(-g_i(x))$$

この関数は以下の性質を持ちます：

• $g_i(x) < 0$ （実行可能領域の内部）で有限値
• $g_i(x) \to 0$ のとき $\phi(x) \to +\infty$
• 制約境界への接近を防ぐ

その他のバリア関数

逆数バリア関数

$$\phi(x) = \sum_{i=1}^m \frac{1}{-g_i(x)}$$

二次バリア関数

$$\phi(x) = \sum_{i=1}^m \frac{1}{(-g_i(x))^2}$$

バリア法のアルゴリズム

基本アルゴリズム

Algorithm: Barrier Method
1. 初期点 x₀ ∈ int(F) を選択 (実行可能領域の内部)
2. 初期パラメータ μ₀ > 0 を設定
3. for k = 0, 1, 2, ... do:
   4. バリア問題を解く:
      xₖ₊₁ = argmin{f(x) + μₖφ(x)}
   5. パラメータを更新: μₖ₊₁ = σμₖ (0 < σ < 1)
   6. 収束判定: μₖ₊₁ < ε ならば停止

バリア問題

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

$$\min_{x \in \text{int}(\mathcal{F})} B_\mu(x) = f(x) + \mu \phi(x)$$

ここで、$\mu > 0$ はバリアパラメータです。

理論的性質

中心経路（Central Path）

バリアパラメータ $\mu > 0$ に対する各バリア問題の解 $x^*(\mu)$ は、中心経路を形成します：

$$\{x^*(\mu) : \mu > 0\}$$

KKT条件との関係

$x^*(\mu)$ は以下の条件を満たします：

$$\begin{align} \nabla f(x^*(\mu)) + \sum_{i=1}^m \frac{\mu}{-g_i(x^*(\mu))} \nabla g_i(x^*(\mu)) &= 0 \\ g_i(x^*(\mu)) &< 0, \quad i = 1, \ldots, m \end{align}$$

$\mu \to 0$ の極限で、これは元問題のKKT条件に収束します：

$$\lambda_i^* = \frac{\mu}{-g_i(x^*(\mu))}$$

収束性

定理（バリア法の収束性）:

• 元問題が実行可能で最適解を持つ
• Slater条件が成立する
• 適切な仮定の下で

$\mu \to 0$ のとき、$x^*(\mu)$ は元問題の最適解に収束します。

実装例

Python実装

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

def log_barrier_method(f, grad_f, g, grad_g, x0, mu0=1.0, sigma=0.1, 
                      max_outer=50, max_inner=100, tol=1e-6):
    """
    対数バリア法の実装
    
    Parameters:
    f: 目的関数
    grad_f: 目的関数の勾配
    g: 制約関数のリスト [g1, g2, ..., gm]
    grad_g: 制約関数の勾配のリスト [grad_g1, grad_g2, ..., grad_gm]
    x0: 初期点 (実行可能領域の内部)
    mu0: 初期バリアパラメータ
    sigma: μの減少率 (0 < sigma < 1)
    """
    
    def barrier_function(x):
        """対数バリア関数"""
        barrier = 0.0
        for gi in g:
            gi_val = gi(x)
            if gi_val >= 0:
                return np.inf  # 実行可能領域外
            barrier -= np.log(-gi_val)
        return barrier
    
    def barrier_gradient(x):
        """バリア関数の勾配"""
        grad = np.zeros_like(x)
        for i, (gi, grad_gi) in enumerate(zip(g, grad_g)):
            gi_val = gi(x)
            if gi_val >= 0:
                return np.full_like(x, np.inf)
            grad += grad_gi(x) / (-gi_val)
        return grad
    
    def barrier_objective(x, mu):
        """バリア目的関数"""
        return f(x) + mu * barrier_function(x)
    
    def barrier_objective_grad(x, mu):
        """バリア目的関数の勾配"""
        return grad_f(x) + mu * barrier_gradient(x)
    
    # 初期化
    x = x0.copy()
    mu = mu0
    history = {
        'x_values': [x.copy()],
        'f_values': [f(x)],
        'mu_values': [mu],
        'barrier_values': []
    }
    
    print(f"初期点: {x}, f(x0) = {f(x):.6f}")
    
    for outer_iter in range(max_outer):
        print(f"\n外側反復 {outer_iter}: μ = {mu:.6f}")
        
        # バリア問題を解く
        def objective(x_var):
            return barrier_objective(x_var, mu)
        
        def objective_grad(x_var):
            return barrier_objective_grad(x_var, mu)
        
        # 制約チェック
        feasible = all(gi(x) < 0 for gi in g)
        if not feasible:
            print("警告: 現在の点が実行可能領域外です")
            break
        
        try:
            # 無制約最適化
            result = minimize(objective, x, method='BFGS', jac=objective_grad,
                            options={'maxiter': max_inner, 'gtol': tol * mu})
            
            x = result.x
            
            # 履歴の更新
            history['x_values'].append(x.copy())
            history['f_values'].append(f(x))
            history['mu_values'].append(mu)
            history['barrier_values'].append(barrier_function(x))
            
            print(f"バリア問題の解: {x}")
            print(f"目的関数値: {f(x):.6f}")
            print(f"制約値: {[gi(x) for gi in g]}")
            
            # 収束判定
            if mu < tol:
                print(f"\n収束しました (μ = {mu:.2e})")
                break
            
            # μの更新
            mu *= sigma
            
        except Exception as e:
            print(f"最適化でエラーが発生: {e}")
            break
    
    return x, history

def inverse_barrier_method(f, grad_f, g, grad_g, x0, mu0=1.0, sigma=0.1,
                          max_outer=50, max_inner=100, tol=1e-6):
    """逆数バリア法の実装"""
    
    def barrier_function(x):
        """逆数バリア関数"""
        barrier = 0.0
        for gi in g:
            gi_val = gi(x)
            if gi_val >= 0:
                return np.inf
            barrier += 1.0 / (-gi_val)
        return barrier
    
    def barrier_gradient(x):
        """逆数バリア関数の勾配"""
        grad = np.zeros_like(x)
        for gi, grad_gi in zip(g, grad_g):
            gi_val = gi(x)
            if gi_val >= 0:
                return np.full_like(x, np.inf)
            grad += grad_gi(x) / (gi_val**2)
        return grad
    
    # 以下の実装は log_barrier_method と同様
    # ... (省略)

def example_constrained_optimization():
    """制約付き最適化の例"""
    
    # 問題: min f(x) = (x1-2)² + (x2-1)²
    # subject to: g1(x) = x1² + x2² - 4 ≤ 0
    #            g2(x) = -x1 ≤ 0
    #            g3(x) = -x2 ≤ 0
    
    def f(x):
        return (x[0] - 2)**2 + (x[1] - 1)**2
    
    def grad_f(x):
        return np.array([2*(x[0] - 2), 2*(x[1] - 1)])
    
    def g1(x):
        return x[0]**2 + x[1]**2 - 4
    
    def grad_g1(x):
        return np.array([2*x[0], 2*x[1]])
    
    def g2(x):
        return -x[0]
    
    def grad_g2(x):
        return np.array([-1, 0])
    
    def g3(x):
        return -x[1]
    
    def grad_g3(x):
        return np.array([0, -1])
    
    g = [g1, g2, g3]
    grad_g = [grad_g1, grad_g2, grad_g3]
    
    # 実行可能な初期点
    x0 = np.array([1.0, 0.5])
    
    # バリア法を実行
    x_opt, history = log_barrier_method(f, grad_f, g, grad_g, x0)
    
    print(f"\n最適解: {x_opt}")
    print(f"最適値: {f(x_opt):.6f}")
    print(f"制約値: {[gi(x_opt) for gi in g]}")
    
    # 結果の可視化
    plt.figure(figsize=(15, 5))
    
    # 収束経路
    plt.subplot(1, 3, 1)
    x_history = np.array(history['x_values'])
    
    # 等高線の描画
    x1_range = np.linspace(-0.5, 2.5, 100)
    x2_range = np.linspace(-0.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 = 2 * np.cos(theta)
    circle_y = 2 * np.sin(theta)
    plt.plot(circle_x, circle_y, 'r-', linewidth=2, label='x₁² + x₂² = 4')
    plt.axhline(y=0, color='r', linewidth=2, label='x₂ = 0')
    plt.axvline(x=0, color='r', linewidth=2, label='x₁ = 0')
    
    # 収束経路
    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.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')
    plt.xlabel('外側反復回数')
    plt.ylabel('目的関数値')
    plt.title('目的関数値の収束')
    plt.grid(True)
    
    # バリアパラメータとバリア値
    plt.subplot(1, 3, 3)
    plt.semilogy(history['mu_values'][1:], 'r-o', label='μ値')
    plt.semilogy(history['barrier_values'], 'b-s', label='バリア関数値')
    plt.xlabel('外側反復回数')
    plt.ylabel('値 (対数スケール)')
    plt.title('バリアパラメータの変化')
    plt.legend()
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    return x_opt, history

def compare_barrier_functions():
    """異なるバリア関数の比較"""
    
    # 単純な1次元問題: min x² subject to x ≥ 1
    def f(x):
        return x[0]**2
    
    def g(x):
        return 1 - x[0]  # x ≥ 1 → 1 - x ≤ 0
    
    x_vals = np.linspace(1.01, 3, 100)
    
    # 各バリア関数の計算
    log_barrier = [-np.log(-g([x])) for x in x_vals]
    inv_barrier = [1/(-g([x])) for x in x_vals]
    quad_barrier = [1/(-g([x]))**2 for x in x_vals]
    
    plt.figure(figsize=(12, 4))
    
    plt.subplot(1, 3, 1)
    plt.plot(x_vals, log_barrier, label='対数バリア')
    plt.xlabel('x')
    plt.ylabel('φ(x)')
    plt.title('対数バリア関数')
    plt.grid(True)
    plt.ylim(0, 10)
    
    plt.subplot(1, 3, 2)
    plt.plot(x_vals, inv_barrier, label='逆数バリア')
    plt.xlabel('x')
    plt.ylabel('φ(x)')
    plt.title('逆数バリア関数')
    plt.grid(True)
    plt.ylim(0, 10)
    
    plt.subplot(1, 3, 3)
    plt.plot(x_vals, quad_barrier, label='二次バリア')
    plt.xlabel('x')
    plt.ylabel('φ(x)')
    plt.title('二次バリア関数')
    plt.grid(True)
    plt.ylim(0, 50)
    
    plt.tight_layout()
    plt.show()

# 実行例
if __name__ == "__main__":
    print("=== 制約付き最適化の例 ===")
    x_opt, history = example_constrained_optimization()
    
    print("\n=== バリア関数の比較 ===")
    compare_barrier_functions()

MATLAB実装

function [x, history] = log_barrier_method(f, grad_f, g, grad_g, x0, options)
    % 対数バリア法の実装
    
    if nargin < 6, options = struct(); end
    
    % デフォルトパラメータ
    mu0 = getfield_default(options, 'mu0', 1.0);
    sigma = getfield_default(options, 'sigma', 0.1);
    max_outer = getfield_default(options, 'max_outer', 50);
    max_inner = getfield_default(options, 'max_inner', 100);
    tol = getfield_default(options, 'tol', 1e-6);
    
    x = x0;
    mu = mu0;
    
    history.x_values = x0;
    history.f_values = f(x0);
    history.mu_values = mu0;
    
    fprintf('初期点: [%.4f, %.4f], f(x0) = %.6f\n', x(1), x(2), f(x));
    
    for outer_iter = 1:max_outer
        fprintf('\n外側反復 %d: μ = %.6f\n', outer_iter, mu);
        
        % バリア問題の定義
        barrier_obj = @(x_var) barrier_objective(x_var, f, g, mu);
        barrier_grad = @(x_var) barrier_gradient(x_var, grad_f, g, grad_g, mu);
        
        % 制約チェック
        feasible = true;
        for i = 1:length(g)
            if g{i}(x) >= 0
                feasible = false;
                break;
            end
        end
        
        if ~feasible
            warning('現在の点が実行可能領域外です');
            break;
        end
        
        % 無制約最適化 (簡単な勾配法)
        x_prev = x;
        for inner_iter = 1:max_inner
            grad = barrier_grad(x);
            
            if norm(grad) < tol * mu
                break;
            end
            
            % 簡単な直線探索
            alpha = backtrack_line_search(barrier_obj, barrier_grad, x, -grad);
            x = x - alpha * grad;
        end
        
        % 履歴の更新
        history.x_values = [history.x_values, x];
        history.f_values = [history.f_values, f(x)];
        history.mu_values = [history.mu_values, mu];
        
        fprintf('バリア問題の解: [%.6f, %.6f]\n', x(1), x(2));
        fprintf('目的関数値: %.6f\n', f(x));
        
        % 収束判定
        if mu < tol
            fprintf('\n収束しました (μ = %.2e)\n', mu);
            break;
        end
        
        % μの更新
        mu = sigma * mu;
    end
end

function obj = barrier_objective(x, f, g, mu)
    obj = f(x);
    for i = 1:length(g)
        gi_val = g{i}(x);
        if gi_val >= 0
            obj = inf;
            return;
        end
        obj = obj - mu * log(-gi_val);
    end
end

function grad = barrier_gradient(x, grad_f, g, grad_g, mu)
    grad = grad_f(x);
    for i = 1:length(g)
        gi_val = g{i}(x);
        if gi_val >= 0
            grad = inf(size(x));
            return;
        end
        grad = grad + mu * grad_g{i}(x) / (-gi_val);
    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

実用的な考慮事項

初期点の選択

1. 実行可能領域の内部: $g_i(x_0) < 0$ を満たす点が必要
2. 制約境界からの距離: 境界に近すぎると数値的問題が発生
3. Phase I問題: 実行可能な初期点を見つけるための補助問題

パラメータ調整

バリアパラメータの減少率

• $\sigma$ が小さい場合: 収束が速いが、数値的不安定性
• $\sigma$ が大きい場合: 安定だが収束が遅い
• 推奨値: $\sigma = 0.1 \sim 0.5$

内側問題の収束判定

各バリア問題の解の精度：

$$\|\nabla B_\mu(x)\| \leq \tau \mu$$

ここで、$\tau$ は相対的な許容誤差（例：$\tau = 0.1$）。

数値的安定性

1. 対数の引数: $g_i(x) \geq -\epsilon$ の場合の処理
2. オーバーフロー防止: 大きなバリア値の取り扱い
3. 勾配の計算: 数値微分 vs 解析的勾配

バリア法の変形

自己一致バリア関数（Self-Concordant Barriers）

関数 $\phi(x)$ が $\nu$-自己一致バリア関数であるとは：

1. $\phi$ は凸関数
2. $\phi$ は3回連続微分可能
3. 自己一致性条件を満たす

この性質により、ニュートン法の収束が保証されます。

原始-双対内点法

バリア法と双対理論を組み合わせた手法：

Algorithm: Primal-Dual Interior Point Method
1. 初期点 (x₀, λ₀, s₀) を選択
2. for k = 0, 1, 2, ... do:
   3. ニュートン方向 (Δx, Δλ, Δs) を計算
   4. ステップサイズ α を決定
   5. (xₖ₊₁, λₖ₊₁, sₖ₊₁) = (xₖ, λₖ, sₖ) + α(Δx, Δλ, Δs)
   6. 収束判定

線形計画法への応用

線形計画問題：

$$\begin{align} \min &\quad c^T x \\ \text{s.t.} &\quad Ax = b \\ &\quad x \geq 0 \end{align}$$

対数バリア関数：

$$\phi(x) = -\sum_{i=1}^n \log x_i$$

バリア問題：

$$\min_{Ax=b} c^T x - \mu \sum_{i=1}^n \log x_i$$

まとめ

バリア法の主な特徴：

1. 内点アプローチ: 実行可能領域の内部を通って最適解に接近
2. 理論的基盤: 中心経路とKKT条件との明確な関係
3. 数値的効率: 大規模問題への適用可能性
4. 汎用性: 様々な制約タイプへの適用
5. 発展性: 現代的な内点法の基礎

利点

• 不等式制約の自然な取り扱い
• 多項式時間での収束（適切な条件下）
• 大規模問題への適用可能性

欠点

• 実行可能な初期点が必要
• 等式制約の直接的な扱いが困難
• 数値的な不安定性の可能性

実装のポイント

バリア法を実装する際は、制約境界への接近を適切に制御し、数値的安定性を保つためのパラメータ調整が重要です。特に、バリアパラメータの減少スケジュールと内側問題の収束判定に注意してください。

参考文献

1. Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. Springer.
2. Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press.
3. Nesterov, Y., & Nemirovskii, A. (1994). Interior-Point Polynomial Algorithms in Convex Programming. SIAM.
4. Wright, S. J. (1997). Primal-Dual Interior-Point Methods. SIAM.