非線型最適化テスト

【ミッション： 最急降下法を実装せよ】

新学期早々、ちょっと必要に迫られて、非線型最適化手法をてきとうに実装する事になった。

今回は簡単のために、以下のような多変数関数 f(x, y) = x² + y² を最小化する x 及び y の値を求めるプログラムを作ってみる。

---

ちなみに Java で表現すると、関数 f は以下のようにコーディングできる（static なのは main メソッドから直接テストできるように）。

view plainprint?

1. // f(x,y) = x^2 + y^2  
2. static double f(double x, double y) {  
3.     return x * x + y * y;  
4. }  

// f(x,y) = x^2 + y^2
static double f(double x, double y) {
    return x * x + y * y;
}

このくらいならば即座に x = 0, y = 0 という解が計算できるが、これをより一般の問題に敷衍できるよう、汎用性のある手法を探りたいと思う。

大学初年級程度の解析学の話になるが、関数 f(x, y) が最小値（より正確には極値）をとるとき、(x, y) における f の勾配 ∇f  = ( ∂f(x, y) / ∂x ,  ∂f(x, y) / ∂y ) は 0 となり、それ以外の場合、∇f は、点 (x, y) の近傍で関数値が最も急激に増大する方向のベクトルとなる。

そこで、適当な初期値 (x₀, y₀) を与え、各変数を以下のように更新していけば、理論上、関数値は極値へと漸近するはずである。

x_i+1 = x_i - λ ･ ∂f(x_i, y_i) / ∂x_i
y_i+1 = y_i - λ ･ ∂f(x_i, y_i) / ∂y_i

ただし λ は、任意の正数とする（ λ が大きい場合、更新量が大きくなるため、適正な解を飛び越えてしまう可能性がある。 λ が小さい場合、更新量も小さくなるため、なかなか解に収束しなくなる）。

上記の最適化法は、最急降下法と呼ばれている。

実はこの最急降下法、理論は簡単だがそれなりに強力な手法なので、計算効率を重視しない場合にはこれで充分な事も多い。

---

しかし、最急降下法には主に以下のような問題点がある。

1. ∇f を計算する必要があるため、そもそも微分できない関数には適さない。
2. 最大値・最小値以外の極値が存在する場合、初期値によっては局所解に陥る場合がある。
3. 関数によっては収束が遅い場合がある。

2番目に関しては、遺伝的アルゴリズムや焼きなまし法の採用による局所解の脱出が考えられる。

3番目は、更に効率の良い Gauss-Newton 法や Levenberg–Marquardt 法といったアルゴリズムが存在するので、そちらを使えばよい（とは言っても、Levenberg–Marquardt 法は最急降下法と Gauss-Newton 法を混ぜただけ）。

問題は1番目だが、一般的に多変数関数 f(x₁, …, x_i-1, x_i, x_i+1, …, x_n) の x_i における偏導関数が未知のとき、任意の点 (x₁, …, x_i-1, x_i, x_i+1, …, x_n) における x_i の偏微分係数を“数値的”に求める事ができる。

ただし、h は極めて微小な数値とする。

上式を用いて各 x_i について偏微分係数を求めれば、n 次元空間上の任意の点 (x₁, …, x_i-1, x_i, x_i+1, … x_n) における勾配ベクトルの各要素が得られる。

これを Java で書くと、以下のようになる。

view plainprint?

1. // ∇f(x0, y0)を求める      
2. static double[] getGradF(double x0, double y0) {  
3.     // 微小数  
4.     double h = 0.00000001;  
5.           
6.     // xで偏微分  
7.     double dx = (f(x0+h, y0) - f(x0, y0)) / h;  
8.           
9.     // yで偏微分  
10.     double dy = (f(x0, y0+h) - f(x0, y0)) / h;  
11.           
12.     // 戻り値: ∇f = (dx, dy)  
13.     double[] result = {dx, dy};  
14.           
15.     return result;  
16. }  

// ∇f(x0, y0)を求める    
static double[] getGradF(double x0, double y0) {
    // 微小数
    double h = 0.00000001;
        
    // xで偏微分
    double dx = (f(x0+h, y0) - f(x0, y0)) / h;
        
    // yで偏微分
    double dy = (f(x0, y0+h) - f(x0, y0)) / h;
        
    // 戻り値: ∇f = (dx, dy)
    double[] result = {dx, dy};
        
    return result;
}

---

上記のメソッドを用いる事で、適当な初期値から関数を最適化する事ができる。

以下のコードでは、無限ループの中で反復的に解ベクトルを更新するようコーディングしている。

更新後の解ベクトルと、更新前の解ベクトルのユークリッド距離が閾値（たとえば0.0000000000001）未満ならば、収束したと見なして反復処理を抜ける。

view plainprint?

1. public static void main(String[] args) {  
2.     // てきとうに初期値をセット  
3.     double x0 = 200;  
4.     double y0 = 200;  
5.           
6.     // 更新用  
7.     double x1;  
8.     double y1;  
9.           
10.     System.out.println("初期値は、" +  
11.             "( x = " + x0 + ",\t y = " + y0 +" ) です。" +  
12.             "\t f(x, y) = " + f(x0, y0) + "です。");  
13.                   
14.     while(true) {  
15.         double[] gradF = getGradF(x0, y0);  
16.         double lambda = 0.1;  
17.               
18.         // 更新  
19.         x1 = x0 - lambda*gradF[0];  
20.         y1 = y0 - lambda*gradF[1];  
21.               
22.         // 収束したら終了  
23.         double e = 0.0000000000001;   // 収束条件  
24.         double distance   
25.             = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));  
26.         if (distance < e) break;  
27.               
28.         x0 = x1;  
29.         y0 = y1;  
30.     }  
31.     System.out.println("最適化が終了しました。 " +  
32.             "( x = " + x1 + ",\t y = " + y1 +" ) です。" +  
33.             "\t f(x, y) = " + f(x1, y1) + "です。");  
34. }  

public static void main(String[] args) {
    // てきとうに初期値をセット
    double x0 = 200;
    double y0 = 200;
        
    // 更新用
    double x1;
    double y1;
        
    System.out.println("初期値は、" +
            "( x = " + x0 + ",\t y = " + y0 +" ) です。" +
            "\t f(x, y) = " + f(x0, y0) + "です。");
                
    while(true) {
        double[] gradF = getGradF(x0, y0);
        double lambda = 0.1;
            
        // 更新
        x1 = x0 - lambda*gradF[0];
        y1 = y0 - lambda*gradF[1];
            
        // 収束したら終了
        double e = 0.0000000000001;   // 収束条件
        double distance 
            = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
        if (distance < e) break;
            
        x0 = x1;
        y0 = y1;
    }
    System.out.println("最適化が終了しました。 " +
            "( x = " + x1 + ",\t y = " + y1 +" ) です。" +
            "\t f(x, y) = " + f(x1, y1) + "です。");
}

---

これだけではよく解らんので、実際に極値に収束していく様子を視覚化した図を示す。

うまくいったよ。やったね！