最終更新: 2016/12/25

大域最適化

柏木 雅英

1. はじめに

n変数関数の与えられた超直方体領域(区間ベクトル)での最大化、最小化、 最大値、最小値を求める。

ひたすら区間を分割するという素朴なアルゴリズムなので、効率は良くない。 凸関数など条件がいい場合は、他の方法を考えた方が良さそう。 アルゴリズムは、 globalopt.pdf (区間演算による大域最適化、version: 2016/12/25) の通り。

分割数を少なくすれば計算時間の増大は軽減されるので、 少ない分割数で関数の最大値、最小値を計算させて 区間演算の過大評価を抑える目的で使用することも考えられる。

2. ファイル構成

optimize.hpp
(内部で区間演算と自動微分を使っているので、その関連ファイルも必要)

3. 使い方

test-optimize.cc に使用例がある。以下に詳細を説明する。

3.1 最小値を与える区間を求める

まず、
f(x,y) = 1/((x-2)2 + (y-3)2+1) - 1/((x-5)2+(y-6)2+1)
という2変数関数を([0,10] × [0,10])で最小化してみよう。グラフの概形は

のような感じ。 プログラムを示す(opt-sample1.cc)。

#include <kv/optimize.hpp>

namespace ub = boost::numeric::ublas;

struct Func {
    template <class T> T operator() (const ub::vector<T>& x){
        return 1/(pow(x(0)-2,2) + pow(x(1)-3,2)+1) - 1/(pow(x(0)-5,2)+pow(x(1)-6,2)+1);
    }
};

typedef kv::interval<double> itv;

int main()
{
    int i;
    ub::vector<itv> I;
    std::list< ub::vector<itv> > result;
    std::list< ub::vector<itv> >::iterator p;

    std::cout.precision(17);

    I.resize(2);
    for (i=0; i<I.size(); i++) I(i) = itv(0., 10.);

    result = minimize(I, Func(), 1e-5);

    p = result.begin();
    while (p != result.end()) {
        std::cout << *(p++) << "\n";
    }
}
与えられた区間Iを区間幅が1e-5以下になるまで分割し、最小値を与える 可能性のある区間のリストを返す。 実行すると、
[2]([5.0082302093505859,5.0082683563232422],[6.0082340240478515,6.0082626342773438])
と表示される。最小値を与える可能性のある区間は複数になることもある。 なお、複数の解のうち、隣り合った区間同士を自動的に結合して一つにする 機能がデフォルトでonになっている。
    result = minimize(I, Func(), 1e-5, false);
のように第4引数にfalseを与えるとこの機能をoffにすることができる。 (opt-sample2.cc)。 実行すると、
[2]([5.0082302093505859,5.00823974609375],[6.0082340240478515,6.0082435607910157])
[2]([5.0082302093505859,5.00823974609375],[6.0082435607910156,6.0082530975341797])
[2]([5.0082302093505859,5.00823974609375],[6.0082530975341796,6.0082626342773438])
[2]([5.00823974609375,5.0082492828369141],[6.0082340240478515,6.0082435607910157])
[2]([5.008249282836914,5.0082588195800782],[6.0082340240478515,6.0082435607910157])
[2]([5.00823974609375,5.0082492828369141],[6.0082435607910156,6.0082530975341797])
[2]([5.00823974609375,5.0082492828369141],[6.0082530975341796,6.0082626342773438])
[2]([5.008249282836914,5.0082588195800782],[6.0082435607910156,6.0082530975341797])
[2]([5.008249282836914,5.0082588195800782],[6.0082530975341796,6.0082626342773438])
[2]([5.0082588195800781,5.0082683563232422],[6.0082340240478515,6.0082435607910157])
[2]([5.0082588195800781,5.0082683563232422],[6.0082435607910156,6.0082530975341797])
[2]([5.0082588195800781,5.0082683563232422],[6.0082530975341796,6.0082626342773438])
のようになる。 また、分割の限界を区間幅でなく
    result = minimize(I, Func(), 10000);
のように分割数で指定することもできる。 (opt-sample3.cc)。 第3引数の型がintだと、分割数の指定とみなされる。 実際には、領域の次元をn、分割数をmとして区間幅を
w = ((Π0≤i<nwidth(I(i))) / m)1/n
で計算して区間幅versionの関数を呼び出しているだけ。 この例題の場合は「10000分割=各成分を100分割=区間幅0.1」ということになる。 実行すると
[2]([4.921875,5.078125],[5.9375,6.09375])
となる。

なお、optimizeという関数もあるが、使い方はminimizeと全く同じ。 古いversionとの互換性のために残してある。

3.2 最小値を求める

最小値を与える区間は不要で、最小値のみ欲しいという場合は わずらわしいリストを使う必要がなくより簡単になる。 (opt-sample4.cc)。
#include <kv/optimize.hpp>

namespace ub = boost::numeric::ublas;

struct Func {
    template <class T> T operator() (const ub::vector<T>& x){
        return 1/(pow(x(0)-2,2) + pow(x(1)-3,2)+1) - 1/(pow(x(0)-5,2)+pow(x(1)-6,2)+1);
    }
};

typedef kv::interval<double> itv;

int main()
{
    int i;
    ub::vector<itv> I;

    std::cout.precision(17);

    I.resize(2);
    for (i=0; i<I.size(); i++) I(i) = itv(0., 10.);

    std::cout << minimize_value(I, Func(), 10000) << "\n";
}
関数名はminimize_value。 中では最小値を与える可能性のある区間のリストを受け取って、 その区間で関数値を再度計算している。 minimize関数と同様に、区間幅と区間数のどちらの指定も可能。 実行すると、
[-0.9498710098808304,-0.93640942001493898]
となる。

3.3 最大化問題

目的関数に"-"を付ければ最大化問題は最小化問題に直せるが、 minimize / minimize_valueの代わりにmaximize / maximize_valueを呼びだせば、 内部で自動的にその作業を行なう。

maximizeの例 (opt-sample5.cc) とその実行結果。

[2]([1.875,2.03125],[2.890625,3.046875])

maximize_valueの例 (opt-sample5.cc) とその実行結果。

[0.94167103624138703,0.94890324783426994]

3.4 1変数関数の最小化/最大化

1変数関数を最小化/最大化するのに、いちいち1変数のベクトルを受け取る 関数を作るのは面倒なので、1変数専用の minimize / minimize_value / maximize / maximize_valueを用意した。

以下は、関数

f(x) = sin(x) + log(x)
を区間[1,10]で、最小値を与える区間、最小値、 最大値を与える区間、最大値を計算させる例 (opt-sample7.cc)。
#include <kv/optimize.hpp>

struct Func {
    template <class T> T operator() (const T& x){
        return sin(x) + log(x);
    }
};

typedef kv::interval<double> itv;

int main()
{
    std::list<itv> result;
    std::list<itv>::iterator p;

    std::cout.precision(17);

    result = minimize(itv(1,10), Func(), 10000);
    p = result.begin();
    while (p != result.end()) {
        std::cout << *(p++) << "\n";
    }

    std::cout <<  minimize_value(itv(1,10), Func(), 10000) << "\n";

    result = maximize(itv(1,10), Func(), 10000);
    p = result.begin();
    while (p != result.end()) {
        std::cout << *(p++) << "\n";
    }

    std::cout <<  maximize_value(itv(1,10), Func(), 10000) << "\n";
}
関数定義、領域指定、受け取るリストなどの型がvectorで無いので、 すっきり書ける。中身は多変数版を呼び出しているだけ。 計算結果は次の通り。
[4.487060546875,4.48870849609375]
[0.52635445707781691,0.52659927176792277]
[7.97906494140625,7.98016357421875]
[3.0689398258148803,3.0690775057551321]

4. おわりに

まだ試作品で完成度は高くない。