ひたすら区間を分割するという素朴なアルゴリズムなので、効率は良くない。 凸関数など条件がいい場合は、他の方法を考えた方が良さそう。 アルゴリズムは、 globalopt.pdf (区間演算による大域最適化、version: 2021/3/30) の通り。
分割数を少なくすれば計算時間の増大は軽減されるので、 少ない分割数で関数の最大値、最小値を計算させて 区間演算の過大評価を抑える目的で使用することも考えられる。
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との互換性のために残してある。
#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]となる。
maximizeの例 (opt-sample5.cc) とその実行結果。
[2]([1.875,2.03125],[2.890625,3.046875])
maximize_valueの例 (opt-sample5.cc) とその実行結果。
[0.94167103624138703,0.94890324783426994]
以下は、関数
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]