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

拡張ラグランジュ乗数法

$\displaystyle{ \min_{\mathbf{x}} \: f(\mathbf{x}) \:\: \text{s.t.} \:\: g(\mathbf{x}) = 0 }$ *1

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

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

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

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

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

Lagrange multiplier - Wikipedia

ペナルティ法

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

$\displaystyle{ \min_{\mathbf{x}} \: \left\{ \: f(\mathbf{x}) + \frac{1}{2} \mu \{ g(\mathbf{x}) \} ^2 \right\} }$

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

Penalty method - Wikipedia

この手法は、 $\mu$ の値が大きくなったときに数値的に不安定な問題となってしまうという問題があります。

拡張ラグランジュ乗数法

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

$\displaystyle{ \min_{\mathbf{x}} \: \left\{ \: f(\mathbf{x}) + \frac{1}{2} \mu \{ g(\mathbf{x}) \} ^2 - \lambda g(\mathbf{x}) \right\} }$

という問題を繰り返し解いていきます。その際、

$\displaystyle{ \lambda \leftarrow \lambda - \mu g(\mathbf{x}) \\ \mu \leftarrow \alpha \mu \:\: (\alpha > 1) }$

という更新をしていきます*4。ペナルティ法同様に $\mu$ の値は増加していきますが、ラグランジュ乗数 $\lambda$ が更新されていく分だけ制約 $g(\mathbf{x}) = 0$ が満たされた状態に近づくという性質があり*5、ペナルティ法に比べずっと小さい $\mu$ でも有用な解に辿り着くことができます。

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

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

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

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

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

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

Nelder–Mead method - Wikipedia

#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;
}

C++による実装例 (2018/2/13追記)

等式制約あり最適化問題を、ペナルティ法および拡張ラグランジュ乗数法で解くサンプルプログラムを公開しました。

github.com

C++11 で記述されており、制約なし最適化問題を解くために NLopt が使用されています。ソースコードは CMake で管理されています。

*1:ここでは簡単のため制約条件は1つとしていますが、任意の個数の制約条件がある場合でも基本的な考え方は変わりません。

*2:2018/02/12追記: 英語版 Wikipedia によると、 $\lambda$ に対する偏微分が 0 になる点はラグランジュ関数の極小値や極大値ではなく停留点 (saddle point) となるらしく、通常の数値最適化アルゴリズムでラグランジュ関数を目的関数として計算しても解を発見できないようです。そこで、ラグランジュ関数そのものではなく、例えばラグランジュ関数の勾配の大きさを最小化するように数値最適化問題を解くことで、元々の等式制約あり最適化問題の解を得られるそうです。

*3:ペナルティ法では等式制約だけでなく不等式制約も扱うことができます。

*4:ただしパラメタ更新方法には様々な亜種があるようです。

*5:証明は教科書等を参照。

*6:http://www.mit.edu/~dimitrib/lagr_mult.html

*7:Interactive High-Quality Green-Screen Keying via Color Unmixing – Yağız Aksoy