【C#】ラグランジュの未定乗数法による制約条件付き最適化 - 旅行好きなソフトエンジニアの備忘録

$g(\bf x)=0$ という制約条件付きで、関数 $f(\bf x)$ を最小化する、という問題をラグランジュの未定乗数法を使って解きます。ここでは例として、以下の制約条件下で目的関数を最小化する $(x_0,x_1,x_2)$ を探す、という問題を解きます。
目的関数　： $f(x_0,x_1,x_2)=x_0^2+x_1^2+x_2^2$
制約条件１： $g_1(x_0,x_1,x_2)=(x_0-1)^2+(x_1-1)^2+(x_2-1)^2-1=0$
制約条件２： $g_2(x_0,x_1,x_2)=x_2-1=0$
解き方は以下の手順に従います。

1. 関数と制約条件から評価関数を作成する

$\Pi =f-\sum_{j=1}^{m}{\lambda_j}{g_j}$
これを今回の例題に当てはめると評価関数は以下のように書くことができます。
$\Pi ={x_0^2+x_1^2+x_2^2}-{\lambda_1}{((x_0-1)^2+(x_1-1)^2+(x_2-1)^2-1)}-{\lambda_2}{(x_2-1)}$
この評価関数を作る部分をプログラムで書くと、以下のように書くことができます。

// fWithConstraints:作成する評価関数
// f:目的関数
// variableCount:目的関数に含まれる変数の数
// constraints:制約条件群
Func<double[], double> fWithConstraints = x =>
    {
        double v = f(x);
        for (int i = 0; i < constraints.Length; ++i)
        {
            v -= x[variableCount + i]*constraints[i](x);
        }
        return v;
    };

2. 作成した評価関数を偏微分する

ステップ1で作成した評価関数を $x_i, \lambda_j$ で偏微分し、 $\frac{\partial \Pi}{\partial x_i}=0, \frac{\partial \Pi}{\partial \lambda_j}=0$ を連立方程式で解きます。例題に当てはめると以下の連立方程式を解くことになります。
$\frac{\partial \Pi}{\partial x_0}=0$
$\frac{\partial \Pi}{\partial x_1}=0$
$\frac{\partial \Pi}{\partial x_2}=0$
$\frac{\partial \Pi}{\partial \lambda_0}=0$
$\frac{\partial \Pi}{\partial \lambda_1}=0$
ただ、悲しいことにこの連立方程式をプログラムでどのように解けば良いのか現時点で分からず、問題を言い換えようと思います。上記連立方程式が成り立つのであれば、上記の式を全て二乗して足し算しても下記が成り立つはずです。
$(\frac{\partial \Pi}{\partial x_0})^2+(\frac{\partial \Pi}{\partial x_1})^2+(\frac{\partial \Pi}{\partial x_2})^2+(\frac{\partial \Pi}{\partial \lambda_0})^2+(\frac{\partial \Pi}{\partial \lambda_1})^2=0$
このことから、以下の目的関数 $F(\bf x)$ の最適化問題と見なすことができます。
$F(x_0,x_1,x_2,\lambda_0,\lambda_1)=(\frac{\partial \Pi}{\partial x_0})^2+(\frac{\partial \Pi}{\partial x_1})^2+(\frac{\partial \Pi}{\partial x_2})^2+(\frac{\partial \Pi}{\partial \lambda_0})^2+(\frac{\partial \Pi}{\partial \lambda_1})^2=0$
この目的関数を作る部分をプログラムで書くと、以下のように書くことができます。ちなみにプログラム中のPartialDerivativeというクラスは偏微分を行うためのクラスで、ここ（【C#】偏微分を計算するクラス - 旅行好きなソフトエンジニアの備忘録）で作成したものを使っています。

// objectiveFunction:作成する目的関数
// fWithConstraints:評価関数
// variableCount:目的関数に含まれる変数の数
// constraints:制約条件群
Func<double[], double> objectiveFunction = x =>
    {
        double v = 0.0;
        for (int i = 0; i < variableCount + constraints.Length; ++i)
        {
            v += Math.Pow(PartialDerivative.Compute(fWithConstraints, x, i), 2);
        }
        return v;
    };

ここでステップ1とステップ2をまとめ、目的関数と制約条件から、最適化問題に適用するための新しい目的関数を作成するクラスを作成しておきます。

public class LagrangeMultiplierMethod
{
    /// <summary>
    /// ラグランジュの未定乗数法で最適化問題を解くための目的関数を作成するメソッド
    /// </summary>
    /// <param name="f">目的関数</param>
    /// <param name="variableCount">目的関数に含まれる変数の数</param>
    /// <param name="constraints">制約条件</param>
    /// <returns>目的関数</returns>
    public static Func<double[], double> CreateObjectiveFunction(Func<double[], double> f, int variableCount, Func<double[], double>[] constraints)
    {
        Func<double[], double> fWithConstraints = x =>
            {
                double v = f(x);
                for (int i = 0; i < constraints.Length; ++i)
                {
                    v -= x[variableCount + i]*constraints[i](x);
                }
                return v;
            };

        Func<double[], double> objectiveFunction = x =>
            {
                double v = 0.0;
                for (int i = 0; i < variableCount + constraints.Length; ++i)
                {
                    v += Math.Pow(PartialDerivative.Compute(fWithConstraints, x, i), 2);
                }
                return v;
            };

        return objectiveFunction;
    }
}

3. 連立方程式を解く

ステップ2で説明したように、今回は連立方程式を解くのではなく、ステップ2で作成した目的関数の最小値を与える $\bf x$ を求める最適化問題と見なします。目的関数は作成済みなので、後は例えばここ（【C#】共役勾配法による多変数関数の最適化 - 旅行好きなソフトエンジニアの備忘録）で作成したクラスを利用して問題を解くだけです。以上をまとめ、冒頭の例題を解いてみます。

// 目的関数
Func<double[], double> f = x => x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
// 制約条件１
Func<double[], double> g1 = x => Math.Pow(x[0] - 1.0, 2) + Math.Pow(x[1] - 1.0, 2) + Math.Pow(x[2] - 1.0, 2) - 1.0;
// 制約条件２
Func<double[], double> g2 = x => x[2] - 1.0;
// 制約条件を配列でまとめる
var constraints = new Func<double[], double>[] { g1, g2 };
// 目的関数に含まれる未知数はx0,x1,x2の3つ
int variableCount = 3;
// 最適化問題に適用する新たな目的関数を作成する
Func<double[], double> objectiveFunction = LagrangeMultiplierMethod.CreateObjectiveFunction(f, variableCount, constraints);

// x0,x1,lamda0,lamda1の初期値を0とする。x2=1は制約条件より明らか
var initialX = new double[] { 0.0, 0.0, 1.0, 0.0, 0.0 };
// 計算回数
int iteration = 1000;
// 学習係数
double learningRate = 0.001;
// 最適化問題を解く
double[] answer = ConjugateGradientMethod.Compute(objectiveFunction, initialX, iteration, learningRate);