確率的最適化アルゴリズムの実装と比較：遺伝的アルゴリズムと勾配降下法

はじめに

最適化問題は、科学、工学、経済学など、さまざまな分野で重要な役割を果たしています。多くの最適化問題は複雑であり、局所的最適解に陥りやすいという課題があります。本記事では、代表的な2つの最適化アルゴリズム、勾配降下法（Gradient Descent: GD）と遺伝的アルゴリズム（Genetic Algorithm: GA）を実装し、その性能を比較します。

この記事を書いたひと

デジタルリアクタ合同会社 代表
望月 優輝 / Yuki Mochizuki   

機械学習・統計、数値計算などの領域を軸としたソフトウェアエンジニアリングを専門としています。スタートアップからグローバル企業まで、さまざまなスケールの企業にて、事業価値に直結する計算システムを設計・構築してきました。主に機械学習の応用分野において、出版・発表、特許取得等の実績があり、また、IT戦略やデータベース技術等の情報処理に関する専門領域の国家資格を複数保有しています。九州大学理学部卒業。日本ITストラテジスト協会正会員。

対象読者:

• 最適化アルゴリズムの基礎を理解したい方
• 勾配降下法と遺伝的アルゴリズムの違いを知りたい方
• Pythonを用いた最適化問題の解法に興味がある方

記事のポイント:

• 勾配降下法と遺伝的アルゴリズムの原理と実装を解説
• 複数のベンチマーク関数を用いて、両アルゴリズムの性能を比較
• 実験結果に基づき、各アルゴリズムの長所と短所を考察

最適化アルゴリズムの概要

勾配降下法（Gradient Descent: GD）

勾配降下法は、目的関数の勾配（偏微分）を利用して、関数値が減少する方向に解を更新していく手法です。更新式は以下の通りです。

x_{t+1} = x_t - \eta \nabla f(x_t)

ここで、η は学習率、\nabla f(x_t) は点 x_t における目的関数 f の勾配です。

勾配降下法は、局所的には効率よく解を改善できますが、初期値によっては局所的最適解に陥りやすいです。

モーメンタム付き勾配降下法（Momentum GD）

モーメンタム付き勾配降下法は、通常の勾配降下法に「慣性」を導入し、過去の更新方向の情報を利用します。

v_{t+1} = \mu v_t - \eta \nabla f(x_t)
x_{t+1} = x_t + v_{t+1}

ここで、\mu はモーメンタム係数、v_t は時刻 t における速度ベクトルです。

遺伝的アルゴリズム（Genetic Algorithm: GA）

遺伝的アルゴリズムは、生物の進化を模倣した最適化手法です。複数の解候補（個体）からなる集団を維持し、以下の操作を繰り返します。

1. 選択（Selection）：適合度（目的関数値）に基づいて、次世代に残す個体を選択
2. 交叉（Crossover）：選択された個体同士を「交配」させ、新しい個体を生成
3. 突然変異（Mutation）：低確率で個体の一部をランダムに変化させる

GAは局所的最適解に陥りにくく、大域的な探索能力が高いですが、計算コストが高く、収束に時間がかかります。

実験設計

テスト関数

実験には、以下の4つの非凸最適化問題のベンチマーク関数を使用します。

1. Rastrigin関数：多数の局所的最小値を持つ関数

   f(x) = 10n + \sum_{i=1}^{n} [x_i^2 - 10\cos(2\pi x_i)]

2. Ackley関数：多くの局所的最小値と1つの大域的最小値を持つ関数

   f(x) = -20\exp\left(-0.2\sqrt{\frac{1}{n}\sum_{i=1}^{n}x_i^2}\right) - \exp\left(\frac{1}{n}\sum_{i=1}^{n}\cos(2\pi x_i)\right) + 20 + e

   

3. Schwefel関数：大域的最小値が探索空間の端に近い位置にある関数

   f(x) = 418.9829n - \sum_{i=1}^{n}x_i\sin(\sqrt{|x_i|})

   

4. Rosenbrock関数：バナナ形状の谷を持ち、最適解の発見が難しい関数

   f(x) = \sum_{i=1}^{n-1}[100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2]

   

実装したアルゴリズム

1. 勾配降下法（GD）
2. モーメンタム付き勾配降下法（Momentum）
3. 遺伝的アルゴリズム（GA）

実験パラメータ

• 次元数：2次元、10次元、30次元
• 実行回数：各設定で10回実行
• 評価指標：最終的な解の質（目的関数値）、計算時間、成功率（閾値以下の解を見つけられた割合）

実験結果と考察

2次元での可視化結果

各テスト関数の2次元での等高線と、各アルゴリズムが見つけた最良解を可視化しました。

Rastrigin関数

Ackley関数

Schwefel関数

Rosenbrock関数

収束曲線の比較

各アルゴリズムの収束曲線を比較しました。横軸は正規化された反復回数、縦軸は対数スケールでの目的関数値です。

勾配降下法は初期の収束が速いですが、局所解に陥りやすく、最終的な解の質はGAより劣る場合が多いです。モーメンタム付き勾配降下法は通常の勾配降下法よりも良いですが、今回のように複雑な関数ではGAの方が最終的に良い解を見つけられています。

次元数による影響

次元数が増加すると、問題の複雑さは指数関数的に増大します（次元の呪い）。各次元での成功率（閾値以下の解を見つけられた割合）を比較しました。

次元数が増加するにつれて、すべてのアルゴリズムの性能は低下しますが、GAは他のアルゴリズムよりも高い成功率を維持しています。ただし、GAでも、30次元のような複雑な問題では、最適解を求めるのが難しくなることもわかります。

計算時間の比較

各アルゴリズムの平均計算時間を比較しました。

GAは勾配法と比較して計算時間が長いですが、これは、GAが多数の解候補を同時に評価するためです。しかし、最終的な解の質を考慮すると、複雑な問題ではこの計算コストは正当化される場合が多いです。

結論

実験結果から、以下の結論が得られます。

• 複雑な非凸最適化問題において、遺伝的アルゴリズムは勾配降下法と比較して局所的最適解に陥りにくい。
• GAは広範囲の探索と良い解の活用のバランスが取れており、複雑な問題で良い性能を発揮する。
• 次元数が増加しても、GAは他のアルゴリズムより高い成功率を維持する。
• GAは計算コストが高いが、複雑な問題では最終的な解の質がこのコストを正当化する場合が多い。
• GAで大域的な探索を行った後、勾配法で局所的な改善を行うハイブリッドアプローチが有効である可能性がある。

コード

import numpy as np
import matplotlib.pyplot as plt
import time
import japanize_matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import random
from typing import Callable, List, Tuple, Dict, Any

# シード値を固定して再現性を確保
np.random.seed(42)
random.seed(42)

# テスト関数の定義
def rastrigin(x: np.ndarray) -> float:
    """
    Rastrigin関数: 多数の局所的最小値を持つ関数
    global minimum: f(0,0,...,0) = 0
    """
    n = len(x)
    return 10 * n + np.sum(x**2 - 10 * np.cos(2 * np.pi * x))

def ackley(x: np.ndarray) -> float:
    """
    Ackley関数: 多くの局所的最小値と1つの大域的最小値を持つ関数
    global minimum: f(0,0,...,0) = 0
    """
    a, b, c = 20, 0.2, 2 * np.pi
    n = len(x)
    sum1 = np.sum(x**2)
    sum2 = np.sum(np.cos(c * x))
    term1 = -a * np.exp(-b * np.sqrt(sum1 / n))
    term2 = -np.exp(sum2 / n)
    return term1 + term2 + a + np.exp(1)

def schwefel(x: np.ndarray) -> float:
    """
    Schwefel関数: 大域的最小値が探索空間の端に近い位置にある関数
    global minimum: f(420.9687,...,420.9687) = 0
    """
    n = len(x)
    return 418.9829 * n - np.sum(x * np.sin(np.sqrt(np.abs(x))))

def rosenbrock(x: np.ndarray) -> float:
    """
    Rosenbrock関数: バナナ形状の谷を持ち、最適解の発見が難しい関数
    global minimum: f(1,1,...,1) = 0
    """
    n = len(x)
    return np.sum(100.0 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)

# 勾配の計算（数値微分）
def numerical_gradient(f: Callable, x: np.ndarray, h: float = 1e-4) -> np.ndarray:
    """
    数値微分による勾配計算
    """
    grad = np.zeros_like(x)
    for i in range(len(x)):
        x_forward = x.copy()
        x_backward = x.copy()
        x_forward[i] += h
        x_backward[i] -= h
        grad[i] = (f(x_forward) - f(x_backward)) / (2 * h)
    return grad

# 勾配降下法の実装
class GradientDescent:
    def __init__(self, learning_rate: float = 0.01, max_iter: int = 1000, tol: float = 1e-6):
        self.learning_rate = learning_rate
        self.max_iter = max_iter
        self.tol = tol
        self.history = []

    def optimize(self, f: Callable, initial_point: np.ndarray, bounds: List[Tuple[float, float]]) -> Tuple[np.ndarray, float, List]:
        """
        勾配降下法による最適化
        """
        x = initial_point.copy()
        self.history = [{'x': x.copy(), 'f': f(x)}]

        for i in range(self.max_iter):
            grad = numerical_gradient(f, x)
            x_new = x - self.learning_rate * grad

            # 境界制約の適用
            for j in range(len(x)):
                x_new[j] = max(min(x_new[j], bounds[j][1]), bounds[j][0])

            # 収束判定
            if np.linalg.norm(x_new - x) < self.tol:
                x = x_new
                self.history.append({'x': x.copy(), 'f': f(x)})
                break

            x = x_new
            self.history.append({'x': x.copy(), 'f': f(x)})

        return x, f(x), self.history

# モーメンタム付き勾配降下法
class MomentumGD:
    def __init__(self, learning_rate: float = 0.01, momentum: float = 0.9, max_iter: int = 1000, tol: float = 1e-6):
        self.learning_rate = learning_rate
        self.momentum = momentum
        self.max_iter = max_iter
        self.tol = tol
        self.history = []

    def optimize(self, f: Callable, initial_point: np.ndarray, bounds: List[Tuple[float, float]]) -> Tuple[np.ndarray, float, List]:
        """
        モーメンタム付き勾配降下法による最適化
        """
        x = initial_point.copy()
        velocity = np.zeros_like(x)
        self.history = [{'x': x.copy(), 'f': f(x)}]

        for i in range(self.max_iter):
            grad = numerical_gradient(f, x)
            velocity = self.momentum * velocity - self.learning_rate * grad
            x_new = x + velocity

            # 境界制約の適用
            for j in range(len(x)):
                x_new[j] = max(min(x_new[j], bounds[j][1]), bounds[j][0])

            # 収束判定
            if np.linalg.norm(x_new - x) < self.tol:
                x = x_new
                self.history.append({'x': x.copy(), 'f': f(x)})
                break

            x = x_new
            self.history.append({'x': x.copy(), 'f': f(x)})

        return x, f(x), self.history

# 遺伝的アルゴリズムの実装
class GeneticAlgorithm:
    def __init__(self, 
                 pop_size: int = 100, 
                 crossover_rate: float = 0.8, 
                 mutation_rate: float = 0.1, 
                 max_iter: int = 100,
                 selection_method: str = 'tournament',
                 crossover_method: str = 'uniform',
                 elitism: bool = True):
        self.pop_size = pop_size
        self.crossover_rate = crossover_rate
        self.mutation_rate = mutation_rate
        self.max_iter = max_iter
        self.selection_method = selection_method
        self.crossover_method = crossover_method
        self.elitism = elitism
        self.history = []

    def initialize_population(self, dim: int, bounds: List[Tuple[float, float]]) -> np.ndarray:
        """
        初期集団の生成
        """
        population = np.zeros((self.pop_size, dim))
        for i in range(self.pop_size):
            for j in range(dim):
                population[i, j] = np.random.uniform(bounds[j][0], bounds[j][1])
        return population

    def evaluate_population(self, f: Callable, population: np.ndarray) -> np.ndarray:
        """
        集団の評価
        """
        fitness = np.zeros(self.pop_size)
        for i in range(self.pop_size):
            fitness[i] = f(population[i])
        return fitness

    def selection(self, population: np.ndarray, fitness: np.ndarray) -> np.ndarray:
        """
        選択操作
        """
        if self.selection_method == 'tournament':
            # トーナメント選択
            selected = np.zeros_like(population)
            for i in range(len(population)):
                idx1, idx2 = np.random.choice(len(population), 2, replace=False)
                if fitness[idx1] < fitness[idx2]:  # 最小化問題なので小さい方が良い
                    selected[i] = population[idx1].copy()
                else:
                    selected[i] = population[idx2].copy()
            return selected
        else:
            # ルーレット選択（最小化問題なので適合度を反転）
            fitness_inv = 1.0 / (fitness + 1e-10)  # ゼロ除算を避ける
            prob = fitness_inv / np.sum(fitness_inv)
            indices = np.random.choice(len(population), len(population), p=prob)
            return population[indices].copy()

    def crossover(self, parents: np.ndarray) -> np.ndarray:
        """
        交叉操作
        """
        offspring = parents.copy()

        for i in range(0, len(parents), 2):
            if i + 1 < len(parents) and np.random.random() < self.crossover_rate:
                if self.crossover_method == 'uniform':
                    # 一様交叉
                    mask = np.random.random(len(parents[i])) < 0.5
                    offspring[i, mask] = parents[i+1, mask]
                    offspring[i+1, mask] = parents[i, mask]
                elif self.crossover_method == 'one_point':
                    # 一点交叉
                    point = np.random.randint(1, len(parents[i]))
                    offspring[i, point:] = parents[i+1, point:]
                    offspring[i+1, point:] = parents[i, point:]
                else:
                    # 二点交叉
                    points = sorted(np.random.choice(len(parents[i]) - 1, 2, replace=False) + 1)
                    offspring[i, points[0]:points[1]] = parents[i+1, points[0]:points[1]]
                    offspring[i+1, points[0]:points[1]] = parents[i, points[0]:points[1]]

        return offspring

    def mutation(self, offspring: np.ndarray, bounds: List[Tuple[float, float]]) -> np.ndarray:
        """
        突然変異操作
        """
        for i in range(len(offspring)):
            for j in range(len(offspring[i])):
                if np.random.random() < self.mutation_rate:
                    # ランダムな値に変異
                    offspring[i, j] = np.random.uniform(bounds[j][0], bounds[j][1])
        return offspring

    def optimize(self, f: Callable, dim: int, bounds: List[Tuple[float, float]]) -> Tuple[np.ndarray, float, List]:
        """
        遺伝的アルゴリズムによる最適化
        """
        # 初期集団の生成と評価
        population = self.initialize_population(dim, bounds)
        fitness = self.evaluate_population(f, population)

        best_idx = np.argmin(fitness)
        best_solution = population[best_idx].copy()
        best_fitness = fitness[best_idx]

        self.history = [{'x': best_solution.copy(), 'f': best_fitness}]

        for generation in range(self.max_iter):
            # エリート保存
            if self.elitism:
                elite = population[best_idx].copy()
                elite_fitness = fitness[best_idx]

            # 選択
            selected = self.selection(population, fitness)

            # 交叉
            offspring = self.crossover(selected)

            # 突然変異
            offspring = self.mutation(offspring, bounds)

            # 評価
            population = offspring
            fitness = self.evaluate_population(f, population)

            # エリートの復元
            if self.elitism:
                worst_idx = np.argmax(fitness)
                if elite_fitness < fitness[worst_idx]:
                    population[worst_idx] = elite
                    fitness[worst_idx] = elite_fitness

            # 最良解の更新
            current_best_idx = np.argmin(fitness)
            if fitness[current_best_idx] < best_fitness:
                best_solution = population[current_best_idx].copy()
                best_fitness = fitness[current_best_idx]

            self.history.append({'x': best_solution.copy(), 'f': best_fitness})

        return best_solution, best_fitness, self.history

# 実験用関数
def run_experiment(function_name: str, 
                  function: Callable, 
                  dim: int, 
                  bounds: List[Tuple[float, float]], 
                  n_runs: int = 10,
                  visualize: bool = False) -> Dict[str, Any]:
    """
    各アルゴリズムで実験を実行し、結果を返す
    """
    results = {
        'function': function_name,
        'dimension': dim,
        'bounds': bounds,
        'gd': {'best_solutions': [], 'best_values': [], 'times': [], 'convergence': []},
        'momentum': {'best_solutions': [], 'best_values': [], 'times': [], 'convergence': []},
        'ga': {'best_solutions': [], 'best_values': [], 'times': [], 'convergence': []}
    }

    # 各アルゴリズムのパラメータ設定
    gd = GradientDescent(learning_rate=0.01, max_iter=1000)
    momentum = MomentumGD(learning_rate=0.01, momentum=0.9, max_iter=1000)
    ga = GeneticAlgorithm(pop_size=100, crossover_rate=0.8, mutation_rate=0.1, max_iter=100)

    for run in range(n_runs):
        print(f"Running experiment {run+1}/{n_runs} for {function_name} in {dim}D")

        # 初期点の生成（GDとMomentumで共通）
        initial_point = np.array([np.random.uniform(b[0], b[1]) for b in bounds])

        # 勾配降下法
        start_time = time.time()
        gd_solution, gd_value, gd_history = gd.optimize(function, initial_point, bounds)
        gd_time = time.time() - start_time

        results['gd']['best_solutions'].append(gd_solution)
        results['gd']['best_values'].append(gd_value)
        results['gd']['times'].append(gd_time)
        results['gd']['convergence'].append([h['f'] for h in gd_history])

        # モーメンタム付き勾配降下法
        start_time = time.time()
        momentum_solution, momentum_value, momentum_history = momentum.optimize(function, initial_point, bounds)
        momentum_time = time.time() - start_time

        results['momentum']['best_solutions'].append(momentum_solution)
        results['momentum']['best_values'].append(momentum_value)
        results['momentum']['times'].append(momentum_time)
        results['momentum']['convergence'].append([h['f'] for h in momentum_history])

        # 遺伝的アルゴリズム
        start_time = time.time()
        ga_solution, ga_value, ga_history = ga.optimize(function, dim, bounds)
        ga_time = time.time() - start_time

        results['ga']['best_solutions'].append(ga_solution)
        results['ga']['best_values'].append(ga_value)
        results['ga']['times'].append(ga_time)
        results['ga']['convergence'].append([h['f'] for h in ga_history])

    # 可視化（2次元の場合のみ）
    if visualize and dim == 2:
        visualize_results(function_name, function, bounds, results)

    return results

def visualize_results(function_name: str, 
                     function: Callable, 
                     bounds: List[Tuple[float, float]], 
                     results: Dict[str, Any]) -> None:
    """
    結果の可視化
    """
    # 関数の可視化
    x = np.linspace(bounds[0][0], bounds[0][1], 100)
    y = np.linspace(bounds[1][0], bounds[1][1], 100)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)

    for i in range(len(x)):
        for j in range(len(y)):
            Z[j, i] = function(np.array([X[j, i], Y[j, i]]))

    # 3Dプロット
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, alpha=0.8)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('f(x, y)')
    ax.set_title(f'{function_name} Function')
    plt.savefig(f'{function_name.lower()}_3d.png')
    plt.close()

    # 等高線プロット
    plt.figure(figsize=(12, 10))
    contour = plt.contour(X, Y, Z, 20, cmap='viridis')
    plt.colorbar(contour)

    # 各アルゴリズムの最良解をプロット
    best_gd_idx = np.argmin(results['gd']['best_values'])
    best_momentum_idx = np.argmin(results['momentum']['best_values'])
    best_ga_idx = np.argmin(results['ga']['best_values'])

    # 最良解の座標を取得
    gd_x, gd_y = results['gd']['best_solutions'][best_gd_idx][0], results['gd']['best_solutions'][best_gd_idx][1]
    momentum_x, momentum_y = results['momentum']['best_solutions'][best_momentum_idx][0], results['momentum']['best_solutions'][best_momentum_idx][1]
    ga_x, ga_y = results['ga']['best_solutions'][best_ga_idx][0], results['ga']['best_solutions'][best_ga_idx][1]

    # 異なるマーカー形状とサイズで表示
    plt.scatter(gd_x, gd_y, color='red', marker='o', s=150, label='GD', edgecolors='black', linewidths=1.5, zorder=10, alpha=0.3)
    plt.scatter(momentum_x, momentum_y, color='blue', marker='s', s=150, label='Momentum', edgecolors='black', linewidths=1.5, zorder=11, alpha=0.3)
    plt.scatter(ga_x, ga_y, color='green', marker='^', s=200, label='GA', edgecolors='black', linewidths=1.5, zorder=12, alpha=0.3)

    # 実際の最良解の位置を点線で表示
    for x_pos, y_pos, color, name in [
        (results['gd']['best_solutions'][best_gd_idx][0], results['gd']['best_solutions'][best_gd_idx][1], 'red', 'GD'),
        (results['momentum']['best_solutions'][best_momentum_idx][0], results['momentum']['best_solutions'][best_momentum_idx][1], 'blue', 'Momentum'),
        (results['ga']['best_solutions'][best_ga_idx][0], results['ga']['best_solutions'][best_ga_idx][1], 'green', 'GA')
    ]:
        if (x_pos != gd_x and y_pos != gd_y) or (x_pos != momentum_x and y_pos != momentum_y) or (x_pos != ga_x and y_pos != ga_y):
            plt.plot([x_pos, gd_x if name == 'GD' else (momentum_x if name == 'Momentum' else ga_x)], 
                     [y_pos, gd_y if name == 'GD' else (momentum_y if name == 'Momentum' else ga_y)], 
                     color=color, linestyle='--', linewidth=0.8, alpha=0.7)

    # 最良解の値をテキストで表示
    plt.annotate(f'GD: {results["gd"]["best_values"][best_gd_idx]:.2f}', 
                 xy=(gd_x, gd_y), xytext=(10, 10), 
                 textcoords='offset points', fontsize=9)
    plt.annotate(f'Momentum: {results["momentum"]["best_values"][best_momentum_idx]:.2f}', 
                 xy=(momentum_x, momentum_y), xytext=(10, 10), 
                 textcoords='offset points', fontsize=9)
    plt.annotate(f'GA: {results["ga"]["best_values"][best_ga_idx]:.2f}', 
                 xy=(ga_x, ga_y), xytext=(10, 10), 
                 textcoords='offset points', fontsize=9)

    plt.xlabel('x')
    plt.ylabel('y')
    plt.title(f'{function_name} Function Contour with Best Solutions')
    plt.legend(loc='upper right')
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.savefig(f'{function_name.lower()}_contour.png', dpi=300, bbox_inches='tight')
    plt.close()

    # 収束曲線の可視化
    plt.figure(figsize=(12, 8))

    # 各アルゴリズムの収束曲線（最も良かった実行のみ）
    best_gd_run = np.argmin(results['gd']['best_values'])
    best_momentum_run = np.argmin(results['momentum']['best_values'])
    best_ga_run = np.argmin(results['ga']['best_values'])

    gd_convergence = np.array(results['gd']['convergence'][best_gd_run])
    momentum_convergence = np.array(results['momentum']['convergence'][best_momentum_run])
    ga_convergence = np.array(results['ga']['convergence'][best_ga_run])

    # GAの反復回数が少ないため、x軸を正規化
    gd_x = np.linspace(0, 1, len(gd_convergence))
    momentum_x = np.linspace(0, 1, len(momentum_convergence))
    ga_x = np.linspace(0, 1, len(ga_convergence))

    plt.plot(gd_x, gd_convergence, label='GD', color='red')
    plt.plot(momentum_x, momentum_convergence, label='Momentum', color='blue')
    plt.plot(ga_x, ga_convergence, label='GA', color='green')

    plt.xlabel('Normalized Iteration')
    plt.ylabel('Function Value')
    plt.title(f'Convergence Curves for {function_name} Function')
    plt.legend()
    plt.yscale('log')
    plt.grid(True)
    plt.savefig(f'{function_name.lower()}_convergence.png')
    plt.close()

def analyze_results(all_results: Dict[str, Dict[str, Any]]) -> None:
    """
    実験結果の分析と可視化
    """
    # 各関数、各次元での最良値の比較
    for function_name, results_by_dim in all_results.items():
        for dim, results in results_by_dim.items():
            print(f"\n{function_name} Function ({dim}D):")

            gd_best = np.min(results['gd']['best_values'])
            momentum_best = np.min(results['momentum']['best_values'])
            ga_best = np.min(results['ga']['best_values'])

            print(f"  GD Best: {gd_best:.6f}")
            print(f"  Momentum Best: {momentum_best:.6f}")
            print(f"  GA Best: {ga_best:.6f}")

            gd_time = np.mean(results['gd']['times'])
            momentum_time = np.mean(results['momentum']['times'])
            ga_time = np.mean(results['ga']['times'])

            print(f"  GD Avg Time: {gd_time:.6f}s")
            print(f"  Momentum Avg Time: {momentum_time:.6f}s")
            print(f"  GA Avg Time: {ga_time:.6f}s")

    # 各関数での成功率の比較（最適解に十分近い解を見つけられた割合）
    plt.figure(figsize=(15, 10))

    functions = list(all_results.keys())
    dimensions = [2, 10, 30]

    for i, dim in enumerate(dimensions):
        success_rates = []

        for function_name in functions:
            results = all_results[function_name][dim]

            # 成功の閾値（関数によって異なる）
            if function_name == 'Rastrigin':
                threshold = 10.0
            elif function_name == 'Ackley':
                threshold = 5.0
            elif function_name == 'Schwefel':
                threshold = 100.0
            else:  # Rosenbrock
                threshold = 100.0

            gd_success = np.mean(np.array(results['gd']['best_values']) < threshold)
            momentum_success = np.mean(np.array(results['momentum']['best_values']) < threshold)
            ga_success = np.mean(np.array(results['ga']['best_values']) < threshold)

            success_rates.append([gd_success, momentum_success, ga_success])

        ax = plt.subplot(1, 3, i+1)
        success_rates = np.array(success_rates)

        x = np.arange(len(functions))
        width = 0.25

        # バーの作成
        gd_bars = ax.bar(x - width, success_rates[:, 0], width, label='GD', color='red')
        momentum_bars = ax.bar(x, success_rates[:, 1], width, label='Momentum', color='blue')
        ga_bars = ax.bar(x + width, success_rates[:, 2], width, label='GA', color='green')

        # バーの上に数値を表示
        def add_labels(bars):
            for bar in bars:
                height = bar.get_height()
                ax.annotate(f'{height:.2f}',
                            xy=(bar.get_x() + bar.get_width() / 2, height),
                            xytext=(0, 3),  # 3ポイント上
                            textcoords="offset points",
                            ha='center', va='bottom',
                            rotation=60,
                            fontsize=8)

        add_labels(gd_bars)
        add_labels(momentum_bars)
        add_labels(ga_bars)

        ax.set_ylabel('Success Rate')
        ax.set_title(f'{dim}D')
        ax.set_ylim(0, 1.1)
        ax.set_xticks(x)
        ax.set_xticklabels(functions)
        ax.legend()

    plt.suptitle('Success Rate by Algorithm, Function, and Dimension')
    plt.tight_layout()
    plt.savefig('success_rates.png')
    plt.close()

    # 各アルゴリズムの収束曲線の比較（2次元の場合）
    plt.figure(figsize=(15, 10))

    for i, function_name in enumerate(functions):
        results = all_results[function_name][2]  # 2次元の結果

        ax = plt.subplot(2, 2, i+1)

        # 最も良かった実行の収束曲線を使用
        best_gd_run = np.argmin(results['gd']['best_values'])
        best_momentum_run = np.argmin(results['momentum']['best_values'])
        best_ga_run = np.argmin(results['ga']['best_values'])

        gd_convergence = np.array(results['gd']['convergence'][best_gd_run])
        momentum_convergence = np.array(results['momentum']['convergence'][best_momentum_run])
        ga_convergence = np.array(results['ga']['convergence'][best_ga_run])

        # GAの反復回数が少ないため、x軸を正規化
        gd_x = np.linspace(0, 1, len(gd_convergence))
        momentum_x = np.linspace(0, 1, len(momentum_convergence))
        ga_x = np.linspace(0, 1, len(ga_convergence))

        ax.plot(gd_x, gd_convergence, label='GD', color='red')
        ax.plot(momentum_x, momentum_convergence, label='Momentum', color='blue')
        ax.plot(ga_x, ga_convergence, label='GA', color='green')

        ax.set_xlabel('Normalized Iteration')
        ax.set_ylabel('Function Value')
        ax.set_title(f'{function_name}')
        ax.legend()
        ax.set_yscale('log')
        ax.grid(True)

    plt.suptitle('Convergence Curves for 2D Functions')
    plt.tight_layout()
    plt.savefig('convergence_comparison.png')
    plt.close()

# メイン実行部分
if __name__ == "__main__":
    # 実験設定
    functions = {
        'Rastrigin': rastrigin,
        'Ackley': ackley,
        'Schwefel': schwefel,
        'Rosenbrock': rosenbrock
    }

    dimensions = [2, 10, 30]
    n_runs = 10  # 各設定での実行回数

    # 関数ごとの探索範囲
    bounds = {
        'Rastrigin': [(-5.12, 5.12)],
        'Ackley': [(-32.768, 32.768)],
        'Schwefel': [(-500, 500)],
        'Rosenbrock': [(-5, 10)]
    }

    # 結果を格納する辞書
    all_results = {}

    # 各関数、各次元で実験を実行
    for function_name, function in functions.items():
        all_results[function_name] = {}

        for dim in dimensions:
            # 境界条件を次元に合わせて拡張
            dim_bounds = bounds[function_name] * dim

            # 実験実行
            results = run_experiment(
                function_name=function_name,
                function=function,
                dim=dim,
                bounds=dim_bounds,
                n_runs=n_runs,
                visualize=(dim == 2)  # 2次元の場合のみ可視化
            )

            all_results[function_name][dim] = results

    # 結果の分析
    analyze_results(all_results)

    print("All experiments completed successfully!") 

関連投稿:

Graph Neural Networksを用いた異常検知入門 数値積分の実装と比較：台形法からガウス求積法まで 時系列データの予測：LSTMとTransformerの性能比較 シミュレーテッドアニーリングで解く巡回セールスマン問題