等式制約あり最適化問題と拡張ラグランジュ乗数法

拡張ラグランジュ乗数法

等式制約あり最適化問題

\displaystyle{
\max_{\mathbf{x}} \: f(\mathbf{x}) \:\: \text{s.t.} \:\: g(\mathbf{x}) = 0
}

を解くためのアルゴリズムとして、拡張ラグランジュ乗数法 (Augmented Lagrangian algorithm) という手法があります。これは、同じく制約あり最適化問題を解くことができる手法であるラグランジュの未定乗数法 (the method of Lagrange multipliers) やペナルティ法 (Penalty method) とよく似ていますが、(おそらく)これらよりもより実用的な手法となっています。

ラグランジュの未定乗数法

ラグランジュの未定乗数法ではラグランジュ関数 (Lagrange function)

\displaystyle{
\mathcal{L}(\mathbf{x}, \lambda) = f(\mathbf{x}) - \lambda \cdot g(\mathbf{x})
}

を考え、これの \mathbf{x}\lambda に対する偏微分が全て 0 となる状態を探索する(例えばラグランジュ関数を目的関数とみなして制約なし最適化問題を解く)ことで等式制約あり最適化問題を解きます。

ペナルティ法

ペナルティ法では「等式制約がどれだけ満たされていないか」に対応する項(ペナルティ)を目的関数に足した制約なし最適化問題

\displaystyle{
\max_{\mathbf{x}} \: \left\{ \: f(\mathbf{x}) - c (g(\mathbf{x}))^2 \right\}
}

を考え、これを解くことで、元の等式制約あり最適化問題を解きます。一回ではなく、徐々に係数  c の値を大きくしていきながら、繰り返しこの問題を解くことで、等式制約が(ほぼ)満たされた解を得られることになります。

拡張ラグランジュ乗数法

拡張ラグランジュ乗数法では、ラグランジュ関数を上述のペナルティによって拡張したようなものを制約なし最適化問題の目的関数とみなし、これを繰り返し解くことで、元の等式制約あり最適化問題を解きます。

なお拡張ラグランジュ乗数法は "the original method of multipliers" と呼ばれることもあるようです*1*2

ライブラリを用いた利用例

NLopt というライブラリには拡張ラグランジュ乗数法が実装されており、これを使うと簡単に等式制約あり最適化問題を解くことができます。

以下の C++ による実装例では

\displaystyle{
\min_{x, y} \: x + y \:\: \text{s.t.} \:\: x^2 + y^2 = 1
}

という等式制約あり最適化問題を解いています。内部で繰り返し解く制約なし最適化問題に使用するアルゴリズムとして Nelder-Mead 法を指定しています。

#include <iostream>
#include <nlopt.hpp>

// f(x, y) = x + y
double objective_function(const std::vector<double> &x, std::vector<double>& /*grad*/, void* /*data*/)
{
    return x[0] + x[1];
}

// g(x, y) = x^2 + y^2 - 1
double equality_constraint(const std::vector<double> &x, std::vector<double>& /*grad*/, void* /*data*/)
{
    return x[0] * x[0] + x[1] * x[1] - 1.0;
}

int main()
{
    // Specify the augmented Lagrangian algorithm
    nlopt::opt optimizer(nlopt::LN_AUGLAG, 2);

    // Specify the objective function and the equality constraint
    optimizer.set_max_objective(objective_function, nullptr);
    optimizer.add_equality_constraint(equality_constraint, nullptr);

    // Specify the Nelder-Mead algorithm for solving each step
    nlopt::opt local_optimizer(nlopt::LN_NELDERMEAD, 2);
    optimizer.set_local_optimizer(local_optimizer);

    // Set initial solution
    std::vector<double> x = { 0.0, 0.0 };

    // Run optimization
    try
    {
        double value;
        optimizer.optimize(x, value);
    }
    catch (std::runtime_error e)
    {
        std::cerr << e.what() << std::endl;
    }

    // Print the result
    std::cout << "Solution: " << x[0] << ", " << x[1] << std::endl;

    return 0;
}

注意:本記事の言葉遣いは必ずしも厳密でないので、そのつもりで読んでください。