Optimization — 最適化

概要

sangi の最適化モジュールは、1次元の黄金分割法から、メタヒューリスティクス・線形計画・凸最適化 (FISTA/ADMM) まで、幅広い最適化アルゴリズムを提供する。ヘッダオンリーであり、*_impl.hpp を include するだけで利用できる。

• 1次元最小化 — 黄金分割法、Brent 法
• 多次元最小化 (勾配ベース) — 最急降下法、BFGS、Newton-CG、共役勾配法、L-BFGS
• 多次元最小化 (微分不要) — Nelder-Mead、Hooke-Jeeves
• メタヒューリスティクス — 遺伝的アルゴリズム、焼きなまし法、差分進化法、CMA-ES、jSO
• 非線形最小二乗 — Gauss-Newton (LM 正則化)、曲線フィッティング
• 制約付き最適化 — L-BFGS-B、ペナルティ法、拡張ラグランジュ法
• 線形計画 — 2フェーズ改訂シンプレックス法
• 凸最適化 — FISTA、ADMM

ビルド

// 1次元最小化
#include <math/optimization/optimization_1d.hpp>

// 多次元最小化 (勾配ベース + 微分不要 + メタヒューリスティクス + 最小二乗 + 制約付き)
#include <math/optimization/optimization_nd.hpp>

// 線形計画
#include <math/optimization/linear_programming.hpp>

// 凸最適化 (FISTA, ADMM)
#include <math/optimization/ConvexOptimization.hpp>

宣言/実装分離方式: 各ヘッダは関数宣言のみを含む。実装は *_impl.hpp に分離されており、 利用時に *_impl.hpp を include することでテンプレート実体化が行われる。

CMake

target_include_directories(my_app PRIVATE <sangi>/include)

直接コンパイル (MSVC)

テンプレート実体化のため、利用する側で *_impl.hpp を include する必要がある。

// 1次元最小化を使う場合
#include <math/optimization/optimization_1d_impl.hpp>
// 多次元最小化を使う場合
#include <math/optimization/optimization_nd_impl.hpp>

共通型

ヘッダ: optimization_base.hpp

OptimizationOptions<T>

OptimizationResult<P, V>

1次元問題では OptimizationResult1D<T> = OptimizationResult<T, T> が使われる。

LineSearchParams<T>

1次元最小化

ヘッダ: optimization_1d.hpp

golden_section_search

template<concepts::OrderedField T>
OptimizationResult1D<T> golden_section_search(
    const std::function<T(T)>& f,
    T a, T b,
    const OptimizationOptions<T>& options = {});

黄金分割法による区間 $[a, b]$ 上の最小化。区間を黄金比で反復的に縮小する。確実に収束するが線形収束である。

戻り値: OptimizationResult1D<T> — 最小点と最小値を返す。

brent_minimize

template<concepts::OrderedField T>
OptimizationResult1D<T> brent_minimize(
    const std::function<T(T)>& f,
    T a, T b,
    const OptimizationOptions<T>& options = {});

Brent 法による区間 $[a, b]$ 上の最小化。放物線補間と黄金分割法を組み合わせ、超線形収束を実現する。1次元最小化の実用上の最良選択である。

多次元最小化 (勾配ベース)

ヘッダ: optimization_nd.hpp

steepest_descent

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> steepest_descent(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const OptimizationOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

最急降下法。$-\nabla f$ 方向に Armijo バックトラッキングライン探索でステップする。実装が単純だが収束が遅いため、主にベースラインとして使用する。

bfgs_minimize

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
OptimizationResult<V, T> bfgs_minimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const OptimizationOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

BFGS 準ニュートン法。ヘッセ逆行列 $H$ を rank-2 更新で近似し、超線形収束を達成する。中規模問題の標準的選択である。$O(n^2)$ メモリが必要であることに注意すること。

newton_cg_minimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> newton_cg_minimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const std::function<V(const V&, const V&)>& d2f,
    const V& x0,
    const OptimizationOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

Newton 共役勾配法 (Truncated Newton)。内部 CG でヘッセ行列方程式 $H\mathbf{p} = -\mathbf{g}$ を近似的に解くため、ヘッセ行列を陽に保持しない。d2f(x, v) はヘッセ行列とベクトル $v$ の積 $H(x) v$ を返す関数である。

conjugateGradientMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> conjugateGradientMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const ConjugateGradientOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

非線形共役勾配法。Polak-Ribiere+ (デフォルト) または Fletcher-Reeves を選択可能。$n$ ステップごとに自動リスタートする。$O(n)$ メモリで大規模問題に適する。

lbfgsMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> lbfgsMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const LBFGSOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

L-BFGS 法 (Limited-memory BFGS)。過去 $m$ 個の勾配差分を保持し、二重ループ再帰で $H^{-1}\mathbf{g}$ を $O(mn)$ で計算する。大規模問題 (変数 1000+) の標準的選択である。

多次元最小化 (微分不要)

ヘッダ: optimization_nd.hpp

nelderMeadMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> nelderMeadMinimize(
    const std::function<T(const V&)>& f,
    const V& x0,
    const NelderMeadOptions<T>& options = {});

Nelder-Mead シンプレックス法。$n+1$ 頂点のシンプレックスを反射・拡大・縮小・全体縮小で変形し、目的関数値を改善する。微分不要で頑健だが、高次元では効率が低下する。収束判定は関数値スプレッドとシンプレックス直径の二重判定である。

hookeJeevesMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> hookeJeevesMinimize(
    const std::function<T(const V&)>& f,
    const V& x0,
    const HookeJeevesOptions<T>& options = {});

Hooke-Jeeves パターン探索法。各座標方向に $\pm\Delta x$ で試行し (探索移動)、改善方向にパターン移動する。ステップサイズを $\sqrt{2}$ 倍で拡大/縮小し、全変数の $\Delta x$ が min_step 以下になれば収束とみなす。

メタヒューリスティクス

ヘッダ: optimization_nd.hpp

大域的最適化のための確率的手法群。すべて微分不要であり、探索範囲 [lower, upper] を指定する。

geneticAlgorithmMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> geneticAlgorithmMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const GeneticAlgorithmOptions<T>& options = {});

実数コード型遺伝的アルゴリズム。トーナメント選択、BLX-$\alpha$ 交叉 ($\alpha = 0.5$)、ガウス変異、エリート保存を使用する。

simulatedAnnealingMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> simulatedAnnealingMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const SimulatedAnnealingOptions<T>& options = {});

焼きなまし法。Metropolis 判定で局所解からの脱出が可能。指数冷却スケジュール $T_{k+1} = r \cdot T_k$ を使用する。

differentialEvolutionMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> differentialEvolutionMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const DifferentialEvolutionOptions<T>& options = {});

差分進化法 (DE/rand/1/bin)。集団ベースの確率的最適化で、微分不要である。ランダム基底ベクトル、1差分ベクトル、二項交叉を使用する。

cmaesMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> cmaesMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const CMAESOptions<T>& options = {});

CMA-ES (共分散行列適応進化戦略)。共分散行列 $C$ を適応的に更新して探索分布の形状を学習する。ill-conditioned / non-separable 問題に強い。中規模問題 ($n <$ 数百) に適する。

jsoMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> jsoMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const JSOOptions<T>& options = {});

jSO (改良版適応的差分進化法)。L-SHADE の改良版に jSO 固有の $F/CR$ 初期化改良を組み合わせた手法である。current-to-pbest/1 変異、成功履歴ベースの適応、線形人口サイズ縮小を特徴とする。

非線形最小二乗

ヘッダ: optimization_nd.hpp

gaussNewtonMinimize

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
OptimizationResult<V, T> gaussNewtonMinimize(
    const std::function<V(const V&)>& residuals,
    const std::function<M(const V&)>& jacobian,
    const V& x0,
    const GaussNewtonOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

Gauss-Newton 法 (Levenberg-Marquardt 正則化付き)。残差関数 $\mathbf{r}(\mathbf{x})$ の二乗和 $\sum r_i^2$ を最小化する。正規方程式 $(J^T J + \lambda I)\delta = -J^T\mathbf{r}$ を解き、Armijo ライン探索でステップする。value は残差の二乗和を返す。

leastSquaresFit

template<concepts::OrderedField T>
LeastSquaresFitResult<T> leastSquaresFit(
    const std::function<Vector<T>(const Vector<T>&)>& residuals,
    const std::function<Matrix<T>(const Vector<T>&)>& jacobian,
    const Vector<T>& p0,
    const LeastSquaresOptions<T>& options = {});

Levenberg-Marquardt 法による非線形最小二乗フィッティング。gain ratio ベースの適応的 $\lambda$ 制御を行い、パラメータ共分散行列・標準誤差等の統計量を返す。

LeastSquaresFitResult<T>

curveFit

template<concepts::OrderedField T>
LeastSquaresFitResult<T> curveFit(
    const std::function<T(T, const Vector<T>&)>& model,
    const std::vector<T>& xdata,
    const std::vector<T>& ydata,
    const Vector<T>& p0,
    const LeastSquaresOptions<T>& options = {});

高レベル曲線フィッティング (scipy.optimize.curve_fit 相当)。モデル関数 $y = f(x, \mathbf{p})$ を観測データにフィットする。残差とヤコビアンは内部で中心差分により自動構築される。

制約付き最適化

ヘッダ: optimization_nd.hpp

lbfgsbMinimize (境界制約)

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> lbfgsbMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const V& lower, const V& upper,
    const LBFGSBOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

L-BFGS-B 法。各変数に上下限制約 $l_i \le x_i \le u_i$ を持つ最適化。勾配射影法で実行可能領域に投影し、自由変数のみで L-BFGS ステップを計算する。無制限の場合は lower[i] = -max(), upper[i] = max() とする。

penaltyMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> penaltyMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const std::vector<std::function<T(const V&)>>& constraints,
    const V& x0,
    const ConstrainedOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

ペナルティ法。不等式制約 $c_i(\mathbf{x}) \le 0$ をペナルティ項に変換し、$\Phi(\mathbf{x}) = f(\mathbf{x}) + \mu \sum \max(0, c_i)^2$ を無制約最適化する。$\mu$ を段階的に増加させて制約違反を抑制する。

augmentedLagrangianMinimize

template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> augmentedLagrangianMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const std::vector<std::function<T(const V&)>>& eq_constraints,
    const std::vector<std::function<T(const V&)>>& ineq_constraints,
    const V& x0,
    const ConstrainedOptions<T>& options = {},
    const LineSearchParams<T>& line_search_params = {});

拡張ラグランジュ法。等式制約 $h_i(\mathbf{x}) = 0$ と不等式制約 $g_i(\mathbf{x}) \le 0$ の両方を処理する。ラグランジュ乗数を推定・更新することで、純粋ペナルティ法の ill-conditioning を回避する。

ConstrainedOptions<T>

線形計画

ヘッダ: linear_programming.hpp

simplex

template<typename T>
LPResult<T> simplex(
    const LinearProgram<T>& lp,
    T eps = T(1e-10),
    size_t maxIterations = 10000);

2フェーズ改訂シンプレックス法。$\min \mathbf{c}^T \mathbf{x}$ subject to 線形制約を解く。等式 (=)、不等式 (<=, >=) すべての制約型に対応する。

simplexMaximize

template<typename T>
LPResult<T> simplexMaximize(
    const LinearProgram<T>& lp,
    T eps = T(1e-10),
    size_t maxIterations = 10000);

最大化ラッパー。目的関数の符号を反転して simplex を呼び出す。

LinearProgram<T>

LinearConstraint<T>

LPResult<T>

凸最適化

ヘッダ: ConvexOptimization.hpp

fista

template<typename T>
FISTAResult<T> fista(
    std::function<Vector<T>(const Vector<T>&)> gradF,
    std::function<Vector<T>(const Vector<T>&, T)> proxG,
    const Vector<T>& x0,
    T L,
    size_t maxIter = 1000, T tol = T{1e-8});

FISTA (Fast Iterative Shrinkage-Thresholding Algorithm)。$\min f(\mathbf{x}) + g(\mathbf{x})$ を解く。$f$ は微分可能な凸関数 (Lipschitz 定数 $L$)、$g$ は近接作用素が計算可能な凸関数。Nesterov 加速により $O(1/k^2)$ 収束を達成する。

lassoFISTA

template<typename T>
FISTAResult<T> lassoFISTA(
    const Matrix<T>& A, const Vector<T>& b,
    T lambda, size_t maxIter = 1000, T tol = T{1e-8});

Lasso 回帰の FISTA 実装。$\min \frac{1}{2}\|A\mathbf{x} - \mathbf{b}\|^2 + \lambda\|\mathbf{x}\|_1$ を解く。Lipschitz 定数はパワーメソッドで自動推定される。

admmLasso

template<typename T>
ADMMResult<T> admmLasso(
    const Matrix<T>& A, const Vector<T>& b,
    T lambda, T rho = T{1},
    size_t maxIter = 1000, T tol = T{1e-6});

ADMM (Alternating Direction Method of Multipliers) による Lasso 回帰。変数分離 $\mathbf{x} = \mathbf{z}$ とペナルティ $\rho$ で分解する。

コード例

Rosenbrock 関数の BFGS 最小化

#include <math/optimization/optimization_nd.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // Rosenbrock: f(x,y) = (1-x)^2 + 100*(y-x^2)^2
    auto f = [](const Vector<double>& x) {
        double a = 1.0 - x[0];
        double b = x[1] - x[0] * x[0];
        return a * a + 100.0 * b * b;
    };

    auto df = [](const Vector<double>& x) -> Vector<double> {
        double dx = -2.0 * (1.0 - x[0]) - 400.0 * x[0] * (x[1] - x[0] * x[0]);
        double dy = 200.0 * (x[1] - x[0] * x[0]);
        return {dx, dy};
    };

    Vector<double> x0 = {-1.0, 1.0};
    auto result = bfgs_minimize<double, Vector<double>, Matrix<double>>(f, df, x0);

    std::cout << "最小点: (" << result.point[0] << ", " << result.point[1] << ")\n";
    std::cout << "最小値: " << result.value << "\n";
    std::cout << "反復回数: " << result.iterations << "\n";
    // 最小点: (1, 1), 最小値: 0
}

曲線フィッティング

#include <math/optimization/optimization_nd.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // モデル: y = a * exp(-b * x)
    auto model = [](double x, const Vector<double>& p) {
        return p[0] * std::exp(-p[1] * x);
    };

    // 観測データ
    std::vector<double> xdata = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5};
    std::vector<double> ydata = {5.0, 3.1, 1.8, 1.2, 0.7, 0.4};

    Vector<double> p0 = {4.0, 0.5};  // 初期推定
    auto result = curveFit<double>(model, xdata, ydata, p0);

    std::cout << "a = " << result.parameters[0]
              << " +/- " << result.standardErrors[0] << "\n";
    std::cout << "b = " << result.parameters[1]
              << " +/- " << result.standardErrors[1] << "\n";
    std::cout << "chi^2/dof = " << result.reducedChiSquared << "\n";
}

線形計画

#include <math/optimization/linear_programming.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // maximize x + 2y  subject to  x + y <= 4,  x <= 3,  y <= 3
    LinearProgram<double> lp;
    lp.objective = {1.0, 2.0};  // simplexMaximize なので最大化
    lp.constraints = {
        {{1.0, 1.0}, ConstraintType::LessEqual, 4.0},
        {{1.0, 0.0}, ConstraintType::LessEqual, 3.0},
        {{0.0, 1.0}, ConstraintType::LessEqual, 3.0}
    };
    auto result = simplexMaximize(lp);

    if (result.status == LPStatus::Optimal) {
        std::cout << "最適値: " << result.objectiveValue << "\n";
        std::cout << "x = " << result.solution[0]
                  << ", y = " << result.solution[1] << "\n";
        // 最適値: 7, x = 1, y = 3
    }
}