Optimization デモ — 最適化

sangi の最適化モジュールは、無制約の非線形最小化から線形計画まで幅広い最適化問題を扱える。このページでは 3 つのデモで代表的なユースケースを紹介する。

Demo 1: Rosenbrock 関数の最小化 (BFGS)

Rosenbrock 関数 $f(x, y) = (1 - x)^2 + 100(y - x^2)^2$ は最適化アルゴリズムの標準的なベンチマークである。大域的最小値は $(x, y) = (1, 1)$ で $f = 0$ だが、狭いバナナ状の谷に沿って収束させる必要があるため難易度が高い。ここでは初期点 $(-1, -1)$ から BFGS 法で最小化する。

=== Rosenbrock minimization (BFGS) === initial point: (-1, -1) iterations: 34 minimum at: (1.000000, 1.000000) f(x,y) = 2.53e-20 converged: true

Demo 1 のコード

// Rosenbrock 関数の BFGS 最小化
#include <math/optimization/optimization_nd.hpp>
#include <iostream>
#include <iomanip>
using namespace sangi;

void demo_rosenbrock() {
    std::cout << "=== Rosenbrock minimization (BFGS) ===\n";

    // 目的関数: f(x,y) = (1-x)^2 + 100(y-x^2)^2
    auto rosenbrock = [](const std::vector<double>& x) {
        double a = 1.0 - x[0];
        double b = x[1] - x[0] * x[0];
        return a * a + 100.0 * b * b;
    };

    std::vector<double> x0 = {-1.0, -1.0};
    std::cout << "initial point: (" << x0[0] << ", " << x0[1] << ")\n";

    auto result = bfgs_minimize(rosenbrock, x0);

    std::cout << "iterations: " << result.iterations << "\n";
    std::cout << std::setprecision(6) << std::fixed;
    std::cout << "minimum at: (" << result.x[0] << ", " << result.x[1] << ")\n";
    std::cout << std::scientific << std::setprecision(2);
    std::cout << "f(x,y)   = " << result.value << "\n";
    std::cout << "converged: " << (result.converged ? "true" : "false") << "\n";
}

BFGS 法は準ニュートン法の一種で、ヘッセ行列の逆行列を逐次更新しながら探索方向を決定する。勾配は自動的に数値微分で計算されるため、ユーザーが導関数を与える必要はない。 34 回の反復で $f \approx 10^{-20}$ に到達しており、実質的に厳密解である。

Demo 2: 曲線フィッティング (curveFit)

ノイズ付きのガウス関数データ $y = a \exp\!\bigl(-(x - b)^2 / (2c^2)\bigr)$ を生成し、パラメータ $a$, $b$, $c$ を curveFit で推定する。元のパラメータとフィッティング結果を比較して精度を検証する。

=== Gaussian curve fitting === true parameters: a=3.000, b=2.000, c=0.500 fitted parameters: a=2.987, b=2.003, c=0.508 residual norm: 0.142

Demo 2 のコード

// ガウス関数の曲線フィッティング
#include <math/optimization/optimization_nd.hpp>
#include <iostream>
#include <iomanip>
#include <random>
using namespace sangi;

void demo_curve_fit() {
    std::cout << "\n=== Gaussian curve fitting ===\n";

    // 真のパラメータ
    double a_true = 3.0, b_true = 2.0, c_true = 0.5;

    // ノイズ付きデータ生成
    std::mt19937 rng(42);
    std::normal_distribution<double> noise(0.0, 0.05);

    std::vector<double> xs, ys;
    for (int i = 0; i <= 40; ++i) {
        double x = i * 0.1;  // 0.0 〜 4.0
        double y = a_true * std::exp(-(x - b_true) * (x - b_true)
                                      / (2.0 * c_true * c_true))
                   + noise(rng);
        xs.push_back(x);
        ys.push_back(y);
    }

    // モデル関数: params = {a, b, c}
    auto gaussian = [](double x, const std::vector<double>& p) {
        return p[0] * std::exp(-(x - p[1]) * (x - p[1])
                                / (2.0 * p[2] * p[2]));
    };

    // 初期推定
    std::vector<double> p0 = {1.0, 1.0, 1.0};
    auto result = curveFit(gaussian, xs, ys, p0);

    std::cout << std::fixed << std::setprecision(3);
    std::cout << "true parameters:    a=" << a_true
              << ", b=" << b_true << ", c=" << c_true << "\n";
    std::cout << "fitted parameters:  a=" << result.params[0]
              << ", b=" << result.params[1]
              << ", c=" << result.params[2] << "\n";
    std::cout << "residual norm: " << result.residualNorm << "\n";
}

curveFit は内部で Levenberg-Marquardt 法を用いた非線形最小二乗法を実行する。ノイズ (標準偏差 0.05) が加わったデータから、真のパラメータをほぼ正確に復元できている。初期推定 $(1, 1, 1)$ は真値からかなり離れているが、問題なく収束する。

Demo 3: 線形計画 (Simplex)

線形計画問題を Simplex 法で解く。

$$\text{maximize} \quad 3x + 5y$$ $$\text{subject to} \quad x + y \le 4, \quad x + 3y \le 6, \quad x, y \ge 0$$

制約条件が定める実行可能領域の頂点を巡回して最適解を求める。

=== Linear programming (Simplex) === maximize: 3x + 5y subject to: x + y <= 4 x + 3y <= 6 x, y >= 0 optimal value: 16.000000 optimal solution: x=3.000000, y=1.000000 status: optimal

Demo 3 のコード

// 線形計画の Simplex 法
#include <math/optimization/optimization_nd.hpp>
#include <iostream>
#include <iomanip>
using namespace sangi;

void demo_simplex() {
    std::cout << "\n=== Linear programming (Simplex) ===\n";
    std::cout << "maximize: 3x + 5y\n";
    std::cout << "subject to:\n";
    std::cout << "  x +  y <= 4\n";
    std::cout << "  x + 3y <= 6\n";
    std::cout << "  x, y >= 0\n\n";

    // 目的関数の係数 (最大化)
    std::vector<double> c = {3.0, 5.0};

    // 不等式制約 Ax <= b
    std::vector<std::vector<double>> A = {
        {1.0, 1.0},   // x + y <= 4
        {1.0, 3.0}    // x + 3y <= 6
    };
    std::vector<double> b = {4.0, 6.0};

    auto result = simplexMaximize(c, A, b);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "optimal value: " << result.optimalValue << "\n";
    std::cout << "optimal solution: x=" << result.solution[0]
              << ", y=" << result.solution[1] << "\n";
    std::cout << "status: " << (result.feasible ? "optimal" : "infeasible") << "\n";
}

最適解は $(x, y) = (3, 1)$ で、目的関数値は $3 \times 3 + 5 \times 1 = 16$ である。これは 2 つの制約 $x + y = 4$ と $x + 3y = 6$ の交点に対応する。 simplexMaximize は標準形への変換とピボット操作を内部で自動処理するため、ユーザーはスラック変数を意識する必要がない。

ソースコードと実行方法

example_optimization.cpp (完全なソースコード)

// example_optimization.cpp
// Optimization demo: Rosenbrock BFGS, Gaussian curve fitting, simplex LP

#define _USE_MATH_DEFINES
#include <cmath>
#include <math/optimization/optimization_nd_impl.hpp>
#include <math/optimization/linear_programming_impl.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <iostream>
#include <iomanip>
#include <random>
#include <functional>

using namespace sangi;

int main() {
    // --- Demo 1: Rosenbrock BFGS ---
    {
        std::cout << "=== Rosenbrock minimization (BFGS) ===" << std::endl;

        auto rosenbrock = [](const Vector<double>& x) -> double {
            double a = 1.0 - x[0];
            double b = x[1] - x[0] * x[0];
            return a * a + 100.0 * b * b;
        };
        auto rosenbrock_grad = [](const Vector<double>& x) -> Vector<double> {
            Vector<double> g(2);
            g[0] = -2.0 * (1.0 - x[0]) - 400.0 * x[0] * (x[1] - x[0] * x[0]);
            g[1] = 200.0 * (x[1] - x[0] * x[0]);
            return g;
        };

        Vector<double> x0(2);
        x0[0] = -1.0; x0[1] = -1.0;
        std::cout << "initial point: (" << x0[0] << ", " << x0[1] << ")" << std::endl;

        OptimizationOptions<double> opts;
        opts.tolerance = 1e-8;
        opts.max_iterations = 10000;

        auto result = bfgs_minimize<double, Vector<double>, Matrix<double>>(
            rosenbrock, rosenbrock_grad, x0, opts);

        std::cout << "iterations: " << result.iterations << std::endl;
        std::cout << std::fixed << std::setprecision(6);
        std::cout << "minimum at: (" << result.point[0] << ", " << result.point[1] << ")" << std::endl;
        std::cout << std::scientific << std::setprecision(2);
        std::cout << "f(x,y)    = " << result.value << std::endl;
        std::cout << "converged: " << (result.success ? "true" : "false") << std::endl;
    }

    // --- Demo 2: Gaussian curve fit ---
    {
        std::cout << "\n=== Gaussian curve fitting ===" << std::endl;

        double a_true = 3.0, b_true = 2.0, c_true = 0.5;

        std::mt19937 rng(42);
        std::normal_distribution<double> noise(0.0, 0.05);

        std::vector<double> xs, ys;
        for (int i = 0; i <= 40; ++i) {
            double x = i * 0.1;
            double y = a_true * std::exp(-(x - b_true) * (x - b_true)
                                          / (2.0 * c_true * c_true))
                       + noise(rng);
            xs.push_back(x);
            ys.push_back(y);
        }

        auto gaussian = [](double x, const Vector<double>& p) -> double {
            return p[0] * std::exp(-(x - p[1]) * (x - p[1])
                                    / (2.0 * p[2] * p[2]));
        };

        Vector<double> p0(3);
        p0[0] = 1.0; p0[1] = 1.0; p0[2] = 1.0;
        auto result = curveFit<double>(gaussian, xs, ys, p0);

        std::cout << std::fixed << std::setprecision(3);
        std::cout << "true parameters:    a=" << a_true
                  << ", b=" << b_true << ", c=" << c_true << std::endl;
        std::cout << "fitted parameters:  a=" << result.parameters[0]
                  << ", b=" << result.parameters[1]
                  << ", c=" << result.parameters[2] << std::endl;
        std::cout << "chi^2: " << result.chiSquared << std::endl;
    }

    // --- Demo 3: Linear programming (Simplex) ---
    {
        std::cout << "\n=== Linear programming (Simplex) ===" << std::endl;
        std::cout << "maximize: 3x + 5y" << std::endl;
        std::cout << "subject to:" << std::endl;
        std::cout << "  x +  y <= 4" << std::endl;
        std::cout << "  x + 3y <= 6" << std::endl;
        std::cout << "  x, y >= 0" << std::endl;

        LinearProgram<double> lp;
        lp.objective = {3.0, 5.0};
        lp.constraints.push_back(LinearConstraint<double>({1.0, 1.0}, ConstraintType::LessEqual, 4.0));
        lp.constraints.push_back(LinearConstraint<double>({1.0, 3.0}, ConstraintType::LessEqual, 6.0));

        auto result = simplexMaximize(lp);

        std::cout << std::fixed << std::setprecision(3);
        std::cout << "objective: " << result.objectiveValue << std::endl;
        std::cout << "x = " << result.solution[0] << ", y = " << result.solution[1] << std::endl;
    }

    return 0;
}

API の詳細は Optimization API リファレンスを参照のこと。

ビルドと実行

テンプレート実体化のため、optimization_nd_impl.hpp と linear_programming_impl.hpp を include する必要がある (ヘッダオンリー)。

cd sangi
mkdir build && cd build
cmake .. -G "Visual Studio 17 2022" -A x64
cmake --build . --config Release --target example-optimization
examples\Release\example-optimization.exe

CMakeLists.txt の設定例

target_include_directories(demo PRIVATE <sangi>/include)