【C#】黄金分割法による最適化 - 旅行好きなソフトエンジニアの備忘録

今回は黄金分割法を使って一変数関数の最適化を行う方法をメモします。

public class GoldenSectionSearch
{
    private static readonly double _ratio = (3.0 - Math.Sqrt(5.0))/2.0;

    /// <summary>
    /// 黄金分割法による最適化を行うメソッド
    /// </summary>
    /// <param name="f">一変数関数</param>
    /// <param name="minX">探索区間最小値</param>
    /// <param name="maxX">探索区間最大値</param>
    /// <param name="iteration">計算回数</param>
    /// <returns>指定区間における関数の最小値を与えるx</returns>
    public static double Compute(Func<double, double> f, double minX, double maxX, int iteration)
    {
        ValidateArguments(minX, maxX, iteration);

        double xn1 = minX + _ratio*(maxX - minX);
        double xn2 = minX + (1.0 - _ratio)*(maxX - minX);
        double fxn1 = f(xn1);
        double fxn2 = f(xn2);
        double solution = xn1;
            
        for (int i = 0; i < iteration; ++i)
        {
            if (fxn1 <= fxn2)
            {
                maxX = xn2;
                xn2 = xn1;
                xn1 = minX + _ratio*(maxX - minX);
                solution = xn1;
                    
                fxn2 = fxn1;
                fxn1 = f(xn1);
            }
            else
            {
                minX = xn1;
                xn1 = xn2;
                xn2 = minX + (1.0 - _ratio)*(maxX - minX);
                solution = xn2;

                fxn1 = fxn2;
                fxn2 = f(xn2);
            }
        }

        return solution;
    }

    private static void ValidateArguments(double minX, double maxX, int iteration)
    {
        if (maxX <= minX)
        {
            throw new ArgumentException("maxX must be larger than minX", "maxX");
        }

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

    /// <summary>
    /// 許容誤差範囲に収めるための計算回数を返すメソッド
    /// </summary>
    /// <param name="minX">探索区間最小値</param>
    /// <param name="maxX">探索区間最大値</param>
    /// <param name="acceptableError">許容誤差範囲</param>
    /// <returns>計算回数</returns>
    public static int ComputeIterationCount(double minX, double maxX, double acceptableError)
    {
        ValidateArguments(minX, maxX, acceptableError);

        double range = maxX - minX;
        double count = (Math.Log10(range) - Math.Log10(acceptableError))/(Math.Log10(2.0) - Math.Log10(-1.0 + Math.Sqrt(5.0)));
        return (int)Math.Floor(count);
    }

    private static void ValidateArguments(double minX, double maxX, double acceptableError)
    {
        if (maxX <= minX)
        {
            throw new ArgumentException("maxX must be larger than minX", "maxX");
        }
            
        if (acceptableError <= 0.0)
        {
            throw new ArgumentOutOfRangeException("acceptableError", "acceptableError must be positive");
        }
    }
}

GoldenSectionSearchクラスのComputeIterationCountメソッドについてですが、黄金分割法では一回の探索で探索範囲が $\frac{(-1+√5)}{2}$ 倍になるため、この性質を利用して解の誤差を許容範囲内に収めるためには何回計算が必要かを調べるためのメソッドです。具体的には探索範囲をR、許容誤差をEとすると、 $R\times(\frac{-1+\sqrt{5}}{2})^n \leqq E$ を満たす $n$ を求めています。

次にGoldenSectionSearchクラスを利用して $f(x)=(x-10)^2+x$ の最小値を与える $x$ を探索範囲[0, 27]、許容誤差0.01として探索します。

// 探索範囲最小値
double minX = 0.0;
// 探索範囲最大値
double maxX = 27.0;
// 解の許容誤差
double acceptableError = 0.01;
// 解を許容誤差に収めるために必要な計算回数
int iteration = GoldenSectionSearch.ComputeIterationCount(minX, maxX, acceptableError);
// 黄金分割法で解を探索する
double answer = GoldenSectionSearch.Compute(x => Math.Pow(x - 10.0, 2) + x, minX, maxX, iteration);