最終更新: 2023/3/10

非線形方程式の全解探索

柏木 雅英

1. はじめに

非線形方程式と探索領域を与え、その領域内の全ての解を計算する。 このプログラムが正常に終了すれば、見つかった解は間違いなく与えられた 領域内の全ての解を網羅しているはずである。

基本原理は単純だが、高速化するための様々なテクニックを使っていて プログラムはかなり複雑である。 まだまだ不完全だがとりあえず使用している手法の解説を 「非線形方程式の全解探索について」 (Version: 2023/3/10) にまとめた。

2. ファイル構成

allsol.hpp

下請けとして、

を使っている。

3. 使い方

3.1 簡単な使い方

とりあえず試すには、 KVライブラリのWebデモ の中の、 Verified Nonlinear Equation Solver (All Solution Finder) を使える。それで生成されるSource Codeも参考になるかも。

例えば、連立方程式

x2 + y2 - 1 = 0
x - y = 0
の、-10≤x≤10, -10≤y≤10における全ての解を求めることを考える。 次のようなファイルを作れば良い (allsol-sample1.cc)。 
#include <kv/allsol.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

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

    std::cout.precision(17);

    I.resize(2);
    I(0) = itv(-10., 10.);
    I(1) = itv(-10., 10.);

    kv::allsol(Func(), I);
}
解きたい問題を記述するには、このようにT型のベクトルを受け取り T型のベクトルを返すような関数オブジェクトを定義する。 更に、解きたい領域を表す区間ベクトルを用意し、allsolに解きたい問題とともに 渡す。これを実行すると、 
[2]([-0.9375,-0.625],[-0.85898437500000003,-0.625])(ex)
[2]([-0.70710678118654769,-0.70710678118654735],[-0.70710678118654769,-0.70710678118654735])(ex:improved)
[2]([0.625,0.9375],[0.625,0.85898437500000003])(ex)
[2]([0.70710678118654735,0.70710678118654769],[0.70710678118654735,0.70710678118654769])(ex:improved)
ne_test: 27, ex_test: 2, ne: 12, ex: 2, giveup: 0
のように表示されるはず。 (ex)は解が発見されたことを、直後の(ex:emproved)はその解を反復改良したものを 表す。最終行は、 を表す。

3.2 実行時に画面に出る情報

第3引数で画面にどのくらい情報を出すか制御できる。 
    kv::allsol(Func(), I, 2);
のようにすると (allsol-sample2.cc)、 最終行に、 
ne_test: 0, ex_test: 0, unknown: 1, ne: 0, ex: 0, giveup: 0    
ne_test: 1, ex_test: 0, unknown: 2, ne: 0, ex: 0, giveup: 0    
ne_test: 2, ex_test: 0, unknown: 3, ne: 0, ex: 0, giveup: 0    
(中略)
ne_test: 19, ex_test: 0, unknown: 8, ne: 6, ex: 0, giveup: 0    
ne_test: 20, ex_test: 0, unknown: 7, ne: 7, ex: 0, giveup: 0    
ne_test: 21, ex_test: 1, unknown: 6, ne: 7, ex: 1, giveup: 0    
ne_test: 22, ex_test: 1, unknown: 5, ne: 8, ex: 1, giveup: 0    
ne_test: 23, ex_test: 1, unknown: 4, ne: 9, ex: 1, giveup: 0    
ne_test: 24, ex_test: 1, unknown: 3, ne: 10, ex: 1, giveup: 0    
ne_test: 25, ex_test: 1, unknown: 2, ne: 11, ex: 1, giveup: 0    
ne_test: 26, ex_test: 2, unknown: 1, ne: 11, ex: 2, giveup: 0    
ne_test: 27, ex_test: 2, unknown: 0, ne: 12, ex: 2, giveup: 0    
ne_test: 27, ex_test: 2, ne: 12, ex: 2, giveup: 0
のように途中経過の区間数がリアルタイムに(CRのみでLFせずに)表示される (allsol-sample2.log)。 この例題では一瞬で終わってしまうだろうが、もう少し大きい問題で 進行状況を観察したいときに指定するとよい。

第3引数のdefault値(省略時の値)は1である。 0にすると、全くメッセージを出さなくなる。

3.3 見つかった解を受け取る

前の例では見つかった解の情報は画面に表示されるだけだったが、 当然見つかった解を受け取って何か処理を行いたいこともあるだろう。 また、 
    kv::allsol(Func(), I, 0);
のようにメッセージを抑制すると見つかった解も見られなくなってしまう。

関数allsolは、実は見つかった解のリストを返している。 これを受け取るには、次の ように区間ベクトルのリストを入れるための変数を作ってやり、それに 受け取ればよい (allsol-sample3.cc)。 

#include <kv/allsol.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

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

    std::cout.precision(17);

    I.resize(2);
    I(0) = itv(-10., 10.);
    I(1) = itv(-10., 10.);

    result = kv::allsol(Func(), I, 0);

    p = result.begin();
    while (p != result.end()) {
        std::cout << *p << "\n";
        p++;
    }
}

4. 重解を持つ場合

全解探索アルゴリズムは、残念ながら重解を持つ場合はその近傍で 解の存在も非存在も示すことが出来ず、無限ループに陥ってしまう。 例えば、連立方程式
x2 + y2 - 1 = 0
2 x2 - y -1 = 0
は2つの解と1つの重解を持つが、 次のプログラムを動かすと2つの解を見つけた後無限ループに陥ってしまう (allsol-sample4.cc)。

このような問題を扱うときには、分割区間がある程度小さくなったら 存在判定または非存在判定をあきらめる必要がある。 第4引数 (defaultは0) にその閾値を指定することが出来る。 

    kv::allsol(Func(), I, 1, 1e-8);
のように指定すると、区間幅が10-8以下であきらめるようになる。 次のプログラムを動かすと、2つの解と4つの不明領域を見つけてちゃんと停止する。 (allsol-sample5.cc)。 
#include <kv/allsol.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = 2. * x(0) * x(0) - x(1) - 1.;

        return y;
    }
};

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

    std::cout.precision(17);

    I.resize(2);
    I(0) = itv(-10., 10.);
    I(1) = itv(-10., 10.);

    kv::allsol(Func(), I, 1, 1e-8);
}
実行結果は次の通り。 
[2]([-0.89003295328093813,-0.81066597662619577],[0.44480715259119962,0.56696428571428593])(ex)
[2]([-0.86602540378443882,-0.86602540378443848],[0.49999999999999983,0.50000000000000023])(ex:improved)
[2]([0.81066597662619577,0.89003295328093813],[0.44480715259119962,0.56696428571428593])(ex)
[2]([0.86602540378443848,0.86602540378443882],[0.49999999999999983,0.50000000000000023])(ex:improved)
ne_test: 185, ex_test: 114, ne: 84, ex: 2, giveup: 4

なお、この不明領域は捨てられてしまうが、これを受け取って何らかの 処理をするには、第5引数に「区間ベクトルのリストへのポインタ」を 与えることが出来る (defaultはNULL) 。第5引数が非NULLであれば、 そこに不明領域を書き込む。 例えば、次のようにすればよい。 (allsol-sample6.cc)。 

#include <kv/allsol.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = 2. * x(0) * x(0) - x(1) - 1.;

        return y;
    }
};

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

    std::cout.precision(17);

    I.resize(2);
    I(0) = itv(-10., 10.);
    I(1) = itv(-10., 10.);

    kv::allsol(Func(), I, 1, 1e-8, &rest);

    p = rest.begin();
    while(p != rest.end()) {
        std::cout << *p << "\n";
        p++;
    }
}
実行すると、 
[2]([-0.89003295328093813,-0.81066597662619577],[0.44480715259119962,0.56696428571428593])(ex)
[2]([-0.86602540378443882,-0.86602540378443848],[0.49999999999999983,0.50000000000000023])(ex:improved)
[2]([0.81066597662619577,0.89003295328093813],[0.44480715259119962,0.56696428571428593])(ex)
[2]([0.86602540378443848,0.86602540378443882],[0.49999999999999983,0.50000000000000023])(ex:improved)
ne_test: 185, ex_test: 114, ne: 84, ex: 2, giveup: 4    
[2]([-1.2603898843129483e-08,1.2603898843129483e-08],[-1.0000000000000003,-0.99999999999999977])
のように、(0, -1)の重解の存在領域が除外されていることが分かる。 なお、あきらめた領域のうち互いに隣り合う領域は自動的に結合されるため、 ここで見つかった4つの領域は結合されて1つになって返されている。

5. カスタマイズ

allsol.hppは、どういうアルゴリズムが最善か筆者の中でも固まっていないため、 アルゴリズムの挙動を変更するための多数のコンパイルオプションがある。 この内容は今後変更されたり削除されたりする可能性があるが、 現状の一覧と簡単な説明を書いておく。

マクロ名 デフォルト値 機能
EDGE_RATIO 0.9 K(I)とIの共通部分が存在する状態で、K(I)の大きさがIのEDGE_RATIO倍以下ならIを分割せずにK(I)に置き換える。分割の境界に解が乗った場合対策。
USE_TRIM 3 1以上なら、f'(I)の情報を使ってIのうち解の存在しない部分を削るアルゴリズム(トリミング)を使う。0なら使わない。1-3は内部で使っている部分和計算のアルゴリズムの違い。
USE_ZERODIVIDE 2 トリミングで、f'(I)の成分が0を含む場合0除算が起きてしまうが、そのとき動作を規定する。0なら何もしない。1ならdivision_part1, division_part2を駆使してトリミングを実行する。2なら、更にトリミングによってIが2つに分かれてしまう場合でもトリミングを実行する。
USE_MULTITRIM 0 1なら、トリミングを行った後ある条件を満たしたら更に連続してトリミングを行う。
USE_AFFINEABS 0 1なら、f(c)とf'(I)を用いてfをaffine化して解の非存在を判定するアルゴリズムを使う。
USE_MIGSKIP 1 1なら、mig(f'(I))を利用して解の存在判定が決して通らない場合を察知し、存在判定をスキップする。
RECOVER_RATIO 0.1 トリミング等で解の特定の成分の区間幅が極端に小さくなると存在判定の性能に悪影響を及ぼすので、一度の分割またはトリミングでRECOVER_RATIO以下の大きさにはならないように制御する。
USE_SLOWDIVIDE 1 1なら、トリミング等で区間幅が半分以下に縮んでいたら、区間の分割を行わない。
USE_FI 1 1なら、f'(I)を計算する前にf(I)による解の非存在判定を行う。通常はf'(I)の副産物としてf(I)が得られるのでf'(I)を計算してからそれを行うが、f(I)とf'(I)の計算コストが大きく違う問題の場合は、先にf(I)で除外できた方が有利なこともありそう。
USE_SUPERLINEAR 1 1なら、解が見つかったときの区間反復において、f'(I_n)を再計算させて2次収束させる。0なら、f'(I_0)を使い続ける1次収束。
UNIFY_REST 1 1なら、重解等であきらめた領域をlistで返却する場合、その区間のうち隣接するものを結合する。
ENABLE_INFINITY 1 (2015/8/6追加) 1なら、入力区間として[0,∞]や[-∞,∞]などの無限区間を入力できる。
WEIGHTED_MAX 1 (2015/8/6追加) 1なら、区間の分割時に重み付き最大値を使う。

これらはお互いに関連しており、もはや筆者にもベストなパラメータを 調整する元気がない。誰か頑張ってくれないだろうか。

6. double以外を両端に持つ区間を使う

allsolはテンプレート関数になっており、内部に使う区間の両端に使う 数の型はカスタマイズ出来るようになっている。 第2引数に渡された区間の両端の型を使って、全体の計算を行うことになる。 詳細はinterval型の説明を見て欲しい。

例えば、 (allsol-sample1.cc) に対して、

変更しただけのプログラム (allsol-sample7.cc) で、4倍精度の全解探索を実行できる。 
#include <kv/allsol.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<kv::dd> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

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

    std::cout.precision(34);

    I.resize(2);
    I(0) = itv(-10., 10.);
    I(1) = itv(-10., 10.);

    kv::allsol(Func(), I);
}
以下は実行結果。解が高精度に計算されている。 
[2]([-0.9375,-0.625],[-0.8589843750000000277555756156289154,-0.625])(ex)
[2]([-0.7071067811865475799119955933626987,-0.7071067811865475799119955933626555],[-0.7071067811865475799119955933626987,-0.7071067811865475799119955933626555])(ex:improved)
[2]([0.625,0.9375],[0.625,0.8589843750000000277555756156289154])(ex)
[2]([0.7071067811865475799119955933626555,0.7071067811865475799119955933626987],[0.7071067811865475799119955933626555,0.7071067811865475799119955933626987])(ex:improved)
ne_test: 27, ex_test: 2, ne: 12, ex: 2, giveup: 0

7. 無限区間の全解探索

(2015/8/6に機能追加) 試験的に、無限区間の全解探索を出来るようにしてみた。

例えば、連立方程式

x2 + y2 - 1 = 0
x - y = 0
の、-∞≤x≤∞, -∞≤y≤∞における全ての解を求めるには 次のようにする (allsol-sample8.cc)。 
#include <kv/allsol.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

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

    std::cout.precision(17);

    I.resize(2);
    I(0) = itv::whole();
    I(1) = itv::whole();

    kv::allsol(Func(), I);
}
ただ区間に無限大を含む区間を与えるだけでよい。 実行すると、 
[2]([-1.0000000000000005,-0.50000000000000022],[-1.0000000000000005,-0.50000000000000022])(ex)
[2]([-0.70710678118654769,-0.70710678118654735],[-0.70710678118654769,-0.70710678118654735])(ex:improved)
[2]([0.50000000000000022,1.0000000000000005],[0.50000000000000022,1.0000000000000005])(ex)
[2]([0.70710678118654735,0.70710678118654769],[0.70710678118654735,0.70710678118654769])(ex:improved)
ne_test: 4739, ex_test: 2, ne: 2062, ex: 2, giveup: 0
のようになる。

ENABLE_INFINITY=1である必要がある。また、WEIGHTED_MAX=1がお勧め。いずれもデフォルト値がそうなっている。

まだ試験的な機能で、実用的な時間では止まらない例題も多いと思われる。今後改良していくつもり。

8. 1変数バージョン

(2016/12/30に機能追加) 1変数関数の全解探索で、いちいち入出力の 変数を「1変数のvector」にしなくてよい別のallsol関数を用意した。 (allsol-sample9.cc)。 
#include <kv/allsol.hpp>

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

typedef kv::interval<double> itv;

int main()
{
    std::cout.precision(17);
    kv::allsol(Func(), itv(-10, 10));
}
引数の型 (区間ベクトルかただの区間か) でどちらが呼び出されるかは区別される。 中身は多変数版のallsolを呼び出しているだけ。

9. その他

allsol.hppにはallsol_listという関数もあり、こちらは第2引数の探索領域を 区間ベクトルでなく区間ベクトルのリストで与え、そのリスト中の全ての 区間に対して全解探索を行う。それ以外はallsolと同じ。

allsol.hppは大変複雑なので、見通し良く簡単にした(性能は低い) allsol-simple.hppもある。 関数名はallsol_simpleとallsol_list_simple。 ソースコードを読んで大体何をしているか理解するためのもので、 実用性は無い。

allsol-affine.hppというのもあるが、まだ試作品であり十分にテストされていない。 affine arithmeticを用いてfを評価して全解探索の性能向上を狙ったもの。 関数名はallsol_affineとallsol_list_affine。