【C#】共役勾配法による多変数関数の最適化 - 旅行好きなソフトエンジニアの備忘録

共役勾配法を利用して多変数関数の最適化を行うクラスを実装します。共役勾配法ではここ（【C#】勾配ベクトルを計算するクラス - 旅行好きなソフトエンジニアの備忘録）で作成したGradientというクラスを利用します。

public class ConjugateGradientMethod
{
    private readonly Func<double[], double> _f;
    private readonly double[] _xn;
    private readonly double _learningRate;
    private double[] _currentGradient;
    private double[] _direction;

    /// <summary>
    /// コンストラクタ
    /// </summary>
    /// <param name="f">多変数関数</param>
    /// <param name="initialX">初期値</param>
    /// <param name="learningRate">学習係数</param>
    public ConjugateGradientMethod(Func<double[], double> f, double[] initialX, double learningRate)
    {
        ValidateArguments(initialX, learningRate);

        _f = f;
        _xn = new double[initialX.Length];
        Array.Copy(initialX, _xn, _xn.Length);
        _learningRate = learningRate;
        // 一回目の探索方向は最急降下法と同じ
        _currentGradient = Gradient.Compute(_f, _xn);
        _direction = _currentGradient.Select(i => -i).ToArray();
    }

    private static void ValidateArguments(double[] initialX, double learningRate)
    {
        if (initialX == null)
        {
            throw new ArgumentNullException("initialX");
        }
        if (initialX.Length <= 1)
        {
            throw new ArgumentOutOfRangeException("initialX", "initialX length must be more than 1");
        }
        if (learningRate <= 0.0)
        {
            throw new ArgumentOutOfRangeException("learningRate", "learningRate must be positive");
        }
    }

    /// <summary>
    /// 現在の解を取得するためのプロパティ
    /// </summary>
    public double[] Xn
    {
        get { return _xn; }
    }

    /// <summary>
    /// 共役勾配法を利用して多変数関数の最小値を探索するメソッド
    /// </summary>
    /// <param name="f">多変数関数</param>
    /// <param name="initialX">初期値</param>
    /// <param name="iteration">計算回数</param>
    /// <param name="learningRate">学習係数</param>
    /// <returns>関数の最小値を与える変数</returns>
    public static double[] Compute(Func<double[], double> f, double[] initialX, int iteration, double learningRate)
    {
        ValidateArguments(initialX, iteration, learningRate);

        var xn = new double[initialX.Length];
        Array.Copy(initialX, xn, xn.Length);

        // 一回目の探索方向は最急降下法と同じ
        double[] currentGradient = Gradient.Compute(f, xn);
        double[] direction = currentGradient.Select(i => -i).ToArray();
        for (int i = 0; i < iteration; ++i)
        {
            for (int j = 0; j < xn.Length; ++j)
            {
                xn[j] += learningRate*direction[j];
            }

            // 探索方向を更新する
            double[] newGradient = Gradient.Compute(f, xn);
            double alpha = Length(newGradient)/Length(currentGradient);
            for (int j = 0; j < direction.Length; ++j)
            {
                direction[j] = -newGradient[j] + alpha*direction[j];
            }

            // 新しい探索方向の計算のため勾配を更新する
            currentGradient = newGradient;
        }

        return xn;
    }

    private static void ValidateArguments(double[] initialX, int iteration, double learningRate)
    {
        ValidateArguments(initialX, learningRate);

        if (iteration <= 0)
        {
            throw new ArgumentOutOfRangeException("iteration", "iteration must be positive");
        }
    }

    /// <summary>
    /// ベクトルの長さの二乗値を計算するメソッド
    /// </summary>
    /// <param name="x">ベクトル</param>
    /// <returns>ベクトルの長さの二乗値</returns>
    private static double Length(double[] x)
    {
        return x.Select(i => i*i).Sum();
    }

    /// <summary>
    /// 現在の解を更新するメソッド
    /// </summary>
    public void Update()
    {
        for (int i = 0; i < _xn.Length; ++i)
        {
            _xn[i] += _learningRate*_direction[i];
        }

        // 探索方向を更新する
        double[] newGradient = Gradient.Compute(_f, _xn);
        double alpha = Length(newGradient)/Length(_currentGradient);
        for (int i = 0; i < _direction.Length; ++i)
        {
            _direction[i] = -newGradient[i] + alpha*_direction[i];
        }

        // 新しい探索方向の計算のため勾配を保存しておく
        Array.Copy(newGradient, _currentGradient, _currentGradient.Length);
    }
}

次にConjugateGradientMethodクラスを利用して、 $f(x_0,x_1)=(-x_0+2x_1+4)^2+(-3x_0+x_1-2)^2$ の最小値を与える $(x_0,x_1)$ を求めます。

// 最適化する多変数関数
Func<double[], double> f = x => Math.Pow(-x[0] + 2.0*x[1] + 4.0, 2) + Math.Pow(-3.0*x[0] + x[1] - 2.0, 2);
// 初期値を(x0,x1)=(-2,4)とする
var initialX = new double[] { -2.0, 4.0 };
// 計算回数を33回としてみる
int iteration = 33;
// 学習係数
double learningRate = 0.01;
double[] answer = ConjugateGradientMethod.Compute(f, initialX, iteration, learningRate);

最適解は $(x_0,x_1)=(-1.6,-2.8)$ ですが、0.01程度の精度の値を得ることが出来ています。また、この例で示した関数は最急降下法（【C#】最急降下法による多変数関数の最適化 - 旅行好きなソフトエンジニアの備忘録）が苦手とする関数で、同程度の精度を得るためには共役勾配法の計算回数が33回に対して最急降下法の場合は159回の計算を必要とします。
なお、更新時の $x$ の遷移を様子を確認したい場合はstaticメソッドではなくConjugateGradientMethodクラスをインスタンス化して下記のように利用します。

// 最適化する多変数関数
Func<double[], double> f = x => Math.Pow(-x[0] + 2.0*x[1] + 4.0, 2) + Math.Pow(-3.0*x[0] + x[1] - 2.0, 2);
// 初期値を(x0,x1)=(-2,4)とする
var initialX = new double[] { -2.0, 4.0 };
// 学習係数
double learningRate = 0.01;
// ConjugateGradientMethodをインスタンス化して更新を繰り返す
var optimizer = new ConjugateGradientMethod(f, initialX, learningRate);
int iteration = 33;
var xHistory = new List<double[]>();
for (int i = 0; i < iteration; ++i)
{
    optimizer.Update();
    xHistory.Add(optimizer.Xn);
}