自動微分と区間演算で非線形関数の根を全て探索するアルゴリズム #Haskell

自動微分について調べていたら『精度保証付き数値計算の基礎』という本に出会いました。その中に非常に面白い応用が載っていたので紹介したいと思います。

allsol :: (RealFloat a, Ord a)
       => (forall b. Floating b => b -> b) -- 根を求める非線形関数
       -> [Interval a]                     -- 探索する区間
       -> [Interval a]                     -- 根が含まれている区間

これが今回紹介する、与えられた非線形関数$f$の与えられた区間における根、すなわち$f(x)=0$を満たすような$x$をただ一つ含む区間を全て探索してくれる関数（の型）です。実装を見る前にその威力を実際に使って確かめてみましょう。

> f x = (x - 1) * (x - 2) * (x - 3)
> allsol f [-1e6...1e6]
[ 0.9837477320983485 ... 1.0146383074141294
, 1.9629261534755658 ... 2.0370738462121643
, 2.9031421559426076 ... 3.1018561522872043]

3根とも正解ですね区間の粗さが少し気になりますが、それは後で改善することにしましょう。

allsol の第一引数は Floating の関数であれば何でも良いので、もっと複雑な関数の根を求めることも出来ます。

みんな大好き連結だけど弧状連結じゃないアレです。これは原点も含めると根が無限個現れてしまうので探索範囲をいい感じに制限します。ちなみにグラフはdesmosを使って描きました。

> f x = sin (pi / x)
> allsol f [0.1...2]
[ 0.10936163341584652 ... 0.11286057988814575
, 0.12346379609469978 ... 0.12649565594828047
, 0.1414290559567462  ... 0.14436828896462503
, 0.1636713078757271  ... 0.16966180319537505
, 0.19999593237639138 ... 0.20000370455557456
, 0.24994360127228576 ... 0.2500533295003087
, 0.33124692229537783 ... 0.33531559720983745
, 0.4731280040565742  ... 0.5255002766863539
, 0.9495376830393973  ... 1.043355630992312]

期待通り$1/n$を含む区間を列挙してくれていますね¹！

さてそれではそろそろ allsol の実装を見てみましょう
こんなにすごい関数なんだからきっと実装はとても複雑なはず…

import Prelude hiding (null)
import Control.Monad (guard)

import Numeric.AD
import Numeric.Interval

allsol _ [] = []
allsol f (x:xs) =
  let result = do
        guard $ 0 `member` f x
        let c = midpoint x
            r = 1 / diff f c
            m = 1 - singleton r * diff f x
            k = singleton (c - r * f c) + m * (x - singleton c)
        guard $ not . null $ k `intersection` x
        if k `isSubsetOf` x
           then Just (Right k)
           else Just (Left $ bisect (k `intersection` x))
   in case result of
        Nothing            ->     allsol f xs
        Just (Right k)     -> k : allsol f xs
        Just (Left (y, z)) ->     allsol f (y : z : xs)

と思いきや、なななんとたったの16行で済んでしまいます！（本当にこれだけです）

その秘密はadとintervalsという自動微分と区間演算の素晴らしいライブラリがHaskellにはあるからなんですね²。

いやいや実装の簡単さもさることながら本当にすごいのはアルゴリズムの中で使われている数学なんです。だって考えても見て下さい。普通に考えたら適当に与えられた非線形関数の根がどこに何個あるのかを当てるのってめちゃくちゃ難しくないですか？

実はこれを可能にしているのはKrawczyk法³と呼ばれる強力な根の存在判定手法があるからなんです。この記事ではKrawczyk法について説明した後、allsolの実装について改めて解説したいと思います。

allsolの実装は柏木先生の"非線形方程式の全解探索"というサイトを大いに参考にさせていただいておりますこの記事では一変数のallsolしか扱いませんが本来は多変数の関数にも適用できる強力なアルゴリズムです。多変数のallsolは実際にWeb上で試すこともできるので是非試してみて下さい。
Verified Nonlinear Equation Solver (All Solution Finder)

区間演算

なんと言っても主役は区間になるので、まずは区間演算について見ていきましょう。具体的には区間同士の足し算や掛け算、区間を引数に取る関数の計算について説明します。どうしてそういう事を考えるのかというと、例えば$\sqrt{2}$という数は無理数なので小数として正確な値を紙に書くことは出来ません。しかし、だいたい1.414から1.415の間にあるというように区間によって表すことは可能です。この表現はもちろん誤差を含んでいますが場合によっては、特に有限桁の値しか扱えない計算機でも扱えるという意味では非常に便利な表現です。それでは1.414から1.415の間にある$\sqrt{2}$という数と1.732から1.733の間にある$\sqrt{3}$という数を足し算した$\sqrt{2}+\sqrt{3}$という数はだいたいどれぐらいの数の間に存在するでしょうか？

区間の四則演算は以下のように定義することが出来ます⁴。

足し算: $[a, b] + [c, d] = [a + c, b + d]$
引き算: $[a, b] - [c, d] = [a - d, b - c]$
掛け算: $[a, b] \times [c, d] = [{\rm min}\{ac,ad,bc,bd\}, {\rm max}\{ac,ad,bc,bd\}]$
割り算: $[a, b] \div [c, d] = [a, b] \times \left[\frac{1}{d}, \frac{1}{c}\right]$

先程も触れたintervalsというライブラリを使って実際にこれらの区間演算を試してみましょう。invervalsでは ... という演算子を使って区間を構築することが出来ます。

> 1.414...1.415
1.414 ... 1.415

これを使って先程の$\sqrt{2}+\sqrt{3}$を計算してみましょう。

> (1.414...1.415) + (1.732...1.733)
3.146 ... 3.148

実際、$\sqrt{2}+\sqrt{3}$の値は$3.1462643..$という値なのでちゃんとこの値を含むような区間が計算されていることが分かります。 memberという関数を使えばその値が区間に含まれているかどうかを実際に計算することも出来ます。

> (sqrt 2 + sqrt 3) `member` ((1.414...1.415) + (1.732...1.733))
True

ちゃんと含まれてることが分かりましたね

他の演算も試してみましょう。

> (1.414...1.415) - (1.732...1.733)
-0.3190000000000002 ... -0.31699999999999995

> (1.414...1.415) * (1.732...1.733)
2.449048 ... 2.452195

> (1.414...1.415) / (1.732...1.733)
0.8159261396422388 ... 0.8169745958429562

うん、多分合ってるんでしょう。

次に区間上の関数について考えます。例えば

f(x) = x^2 + x + 1

という関数があった時に$f$に区間$[0,1]$を与えた結果、つまり$f$による区間$[0, 1]$の像はどうなるでしょうか。

グラフを見ると分かるように$f([0, 1])$は$[1, 3]$となることが分かります。しかしこの$f([0, 1])$をプログラムで直接計算することは困難です。そこで$f$を構成する各演算を区間演算だと思って計算するという方法が考えられます。このように$f$を区間上の関数として解釈することを区間拡張と呼び$f$の区間拡張を$f_{[\ ]}$という記号で表します。実際に$f_{[\ ]}([0, 1])$を評価してみましょう。

\begin{matrix}
f_{[\ ]}
&=& [0, 1] \times [0, 1] + [0, 1] + [1, 1] \\
&=& [0, 1] + [0, 1] + [1, 1] \\
&=& [1, 3]
\end{matrix}

となり、$f$による像と一致することが分かりました。

え、当たり前じゃないかって？それじゃあ今度は

f(x) = x^2 - x + 1

という関数を考えてみましょう。

グラフを見ると$f([0, 1])$は$[0.75, 1]$だということが分かります。ところが$f_{[\ ]}([0, 1])$を評価してみると、

\begin{matrix}
f_{[\ ]}
&=& [0, 1] \times [0, 1] - [0, 1] + [1, 1] \\
&=& [0, 1] - [0, 1] + [1, 1] \\
&=& [-1, 1] + [1, 1] \\
&=& [0, 2]
\end{matrix}

となり、$f$による像と違う結果が得られました。

以上の議論から分かる通り、関数$f$による区間の像と区間拡張された関数$f_{[\ ]}$による区間の計算結果が一致するとは限りません。ただし$f$が区間$I$上で連続であれば$f(I)\subset f_{[\ ]}(I)$が成り立つことが知られています。

Krawczyk法

それでは一変数の微分可能な非線形関数$f$が区間$I$の中に根を持つかどうかを判定する方法を考えましょう⁵。

天下り式に⁶

g(x) = x - \frac{f(x)}{f'(c)}, c \in I

という形の関数を考えると（ちなみに$c$は$I$の中から適当に選んで構いません）

$f(x)=0$ であること

と

$g(x) = x$ となること

が同値になります。つまり$f$の根を見つけたければ$g$の不動点を見つければ良いわけです。不動点があることを証明するために縮小写像の原理を利用しましょう。これは完備距離空間のある写像$f: T \rightarrow T$が、ある1未満の正の定数$q$が存在して任意の$x, y \in T$について

\|f(x) - f(y)\| \leq q \|x - y\|

を満たすならば、$T$の中にただ一つ不動点を持つという定理です。

実は先程定義した$g$はこの不等式は満たすのですが⁷、問題は$I$から$I$への写像とみなせるかどうか、つまり$g(I) \subset I$を満たすとは限らないというところです⁸。これを満たさなければ$g$を$I$から$I$への写像とみなせないので不動点の存在を示すことも出来ません。

$g(I) \subset I$を確かめる一つの方法として$g$を区間拡張した$g_{[\ ]}$を使って

g_{[\ ]}(I) = I - \frac{f_{[\ ]}(I)}{f'(c)}

という区間を計算することを考えてみましょう。$g(I)\subset g_{[\ ]}(I)$なのでもし$g_{[\ ]}(I) \subset I$となることが分かれば$g(I)\subset g_{[\ ]}(I) \subset I$が分かるという算段です。しかし残念ながら区間演算の引き算の性質から⁹これは必ず$I$より大きな区間となってしまい$I$に含まれることはありません。なのでより精密な$g(I)$の評価方法が必要になります。そこで登場するのが平均値形式と呼ばれる方法です。

平均値形式

平均値の定理を思い出してください。そう、高校生の時に習ったであろうあれです。

区間$[a, b]$で連続であり$(a, b)$で微分可能な関数$f(x)$について、ある$c \in [a, b]$が存在して

\frac{f(b) - f(a)}{b - a} = f'(c)

を満たす。これが平均値の定理でした。

この式を少し変形すると

f(b) = f(a) + f'(c)(b - a)

となります。さらに$c$は区間$[a, b]$のある数だったので思い切って$c$のところに$[a, b]$を放り込むと

f(b) \in f(a) + f'([a, b])(b - a)

が分かります。さらに$f'$の区間拡張を考えると$f'([a, b]) \subset f'_{[\ ]}([a, b])$ が成り立つので

f(b) \in f(a) + f'_{[\ ]}([a, b])(b - a)

が成り立つことが分かりました。

これを踏まえて$f(x)$の$c \in I$における平均値形式を

f(c) + f'_{[\ ]}(I)(x - c)

と定めます。すると先程の平均値の定理に関する考察から、この平均値形式は常に$f(x)$自身を要素として持っていることが分かり、

f(x) \in f(c) + f'_{[\ ]}(I)(x - c)

$x$を区間$I$全体で動かして考えると、

f(I) \subset f(c) + f'_{[\ ]}(I)(I - c)

を満たすことが分かります。

これを使って先程の$g(x)$を評価してみましょう。$g(x)$の平均値形式を使えば上の関係式から

g(I) \subset c - \frac{f(c)}{f'(c)} + \left(1 - \frac{f'_{[\ ]}(I)}{f'(c)}\right)(I - c)=K(I)

となることが分かり、単純な区間拡張の時と違って、右辺の$K(I)$は必ずしも$I$より大きくなるとは言えません（やったー！）。
そこで３つの場合に分けて考えてみましょう。

$K(I) \subset I$ となる場合
この場合は$g(I) \subset K(I) \subset I$となり縮小写像の原理が満たされるので$f$は$I$の中にただ一つ根を持つことが分かります。

$K(I)\cap I = \emptyset$ となる場合
この場合は$g(x)$は$I$の中に不動点を持ち得ず、これは$f$が$I$の中に根を持つことと同値な条件だったので、$f$は$I$の中に根を持たないということが分かります。

$K(I)\cap I \neq \emptyset$ となる場合
この場合はなんとも言えません。なので区間を分割するなどして更に細かく調べていく必要が出てきます。

このように$K(I)$と$I$の交わり方を調べることで区間$I$の中にただ一つ根が存在することを判定する方法はKrawczyk法と呼ばれています。そしてこのKrawczyk法がallsolの中で中心的に使われているアルゴリズムなのです。

ところでKrawczyk法に出てくる$K(I)$はKrawczyk写像と呼ばれていて、これを計算するためには$f$の導関数の区間拡張$f'_{[\ ]}$を計算する必要があります。数値微分では導関数の値までは求められてもそれを区間拡張した関数がどのように振る舞うのかまで求めることは出来ません。しかし自動微分であれば導関数のプログラム自体を導出することができるので、そこに現れる演算子を区間演算と見なすことで導関数の区間拡張まで自動的に求めることが出来ます。精度保証付き数値計算の中で自動微分が重宝される理由も分かりますね。

allsolの解説

もう一度 allsol の実装を見てみましょう。

 1: allsol :: (RealFloat a, Ord a)
 2:        => (forall b. Floating b => b -> b) -- 根を求める非線形関数
 3:        -> [Interval a]                     -- 探索する区間
 4:        -> [Interval a]                     -- 根が含まれている区間
 5: allsol _ [] = []
 6: allsol f (x:xs) =
 7:   let result = do
 8:         guard $ 0 `member` f x
 9:         let c = midpoint x
10:             r = 1 / diff f c
11:             m = 1 - singleton r * diff f x
12:             k = singleton (c - r * f c) + m * (x - singleton c)
13:         guard $ not . null $ k `intersection` x
14:         if k `isSubsetOf` x
15:            then Just (Right k)
16:            else Just (Left $ bisect (k `intersection` x))
17:    in case result of
18:         Nothing            ->     allsol f xs
19:         Just (Right k)     -> k : allsol f xs
20:         Just (Left (y, z)) ->     allsol f (y : z : xs)

説明しやすいように行番号もつけました。

まず5行目は再帰関数の基底部だから良いとして6行目の再帰部から順番に見ていきましょう。7行目のresultはMaybe (Either a b) という型をしています。なので7行目からのdo構文の文脈はMaybeです。なぜdo構文を使っているのかというと与えられた区間に根がないことが分かったらその時点で次の区間の探索に進みたいのですがifで分岐するとネストが深くなって見にくくなってしまうので guard を使って早期リターンのような動作を実現しています。こうしたことでソースコードがフラットになり分岐の仕方もデータ構造に陽に現れるようになったので分かりやすくなりました。result を使った分岐は17行目から始まります。ここではcase式を使って、Nothingの場合は区間に根が無かったので次の区間の探索に進む、Just (Right k)の場合は正解が見つかったので結果に追加する、Just (Left (y, z))の場合は根が無いとは言えなかったので更に区間を分割して探索候補に追加するというようにに処理が分岐しています。

それではアルゴリズムのコアの部分である7行目からの処理を見ていきましょう。まず区間xのfによる像が0を含んでいるか確認します。もし含んでいなければそもそも根は存在しないのでここで区間xに関する処理をやめて次の区間の探索に進みます。もし含んでいたとしても複数の根が含まれている可能性があるのでまだなんとも言えない状況です。ちなみにfはforall b. Floating b => b -> bという型の関数でありaがRealFloatとOrdのインスタンスであれば Interval a も Floating のインスタンスになるので f は区間Interval a上の関数として振る舞うことが出来ます。9から12行目まではKrawczyk写像$K(I)$を計算しています。cは区間の要素なのでただの値ですがfは多相関数なので10行目ではdiffに適用された結果fの値上の導関数として振る舞っています。そして11行目ではfはdiffに適用された結果、区間上の導関数として振る舞います。

なんとfは8~11行目までの間に

f :: Interval a                  -> Interval a
f :: AD s (Forward a)            -> AD s (Forward a)
f :: AD s (Forward (Interval a)) -> AD s (Forward (Interval a))

という3種類の使われ方をしているのです¹⁰。そしてfはaではなくInterval aの関数とみなされるだけで区間拡張が行われているのです。多相関数恐るべしですね。

続きを見ていきましょう。13行目ではKrawczyk写像と元の区間の積集合が空かどうかを判定しています。もしこれらが交わりを持っていなければこの区間に根は含まれていないので次の区間に探索を進めます。そうでなければ最後にKrawczyk写像が元の区間にすっぽり包まれているかを確認します。もし部分集合になっていれば縮小写像の原理を満たすのでこの区間の中にただ一つ根が存在します。もしそうでなければなんとも言えない状況なので区間を二分割して探索候補の中に入れて探索候補が無くなるまで再帰的に処理を繰り返していきます。

以上が allsol の動作となります。数学を駆使した見事なアルゴリズムですね

おまけ: 区間ニュートン法

Krawczyk法の$g(x)$の形を見てニュートン法を思い出した人もいるかと思います。実はニュートン法の区間拡張である区間ニュートン法は別にあって以下のようなアルゴリズムになっています。

適当な初期区間$x_0$を決める
$N(x_i) = {\rm mid}(x_i) - \frac{f({\rm mid}x_i)}{f'(x_i)}$を計算する。ただし${\rm mid}$は区間の中点をとる操作であり、$0\notin f'(x_i)$とする
$x_{i+1} = N(x_i) \cap x_i$ とする
終了条件を満たすまで2,3を繰り返す

参考: 区間ニュートン法

これをHaskellで実装するとこんな感じです。

inm :: (RealFloat a, Ord a)
    => (forall b. Floating b => b -> b) -- 根を求める非線形関数
    -> Interval a                       -- 区間の初期値
    -> Interval a                       -- 計算結果の区間
inm f x
  | width x < 1e-6 = x
  | otherwise      =
    let c  = midpoint x
        nx = singleton c - (singleton (f c) / diff f x)
     in inm f (x `intersection` nx)

定義そのままですね。

実際に使ってみましょう。

> f x = x^2 - 2
> inm f (1...2)
1.4142135592945242 ... 1.4142135659471788

見事に$\sqrt{2}$が含まれてる区間を高い精度で求められています。これをallsolと組み合わせて用いると

> f x = (x - 1) * (x - 2) * (x - 3)
> map (inm f) $ allsol f [0...5]
[ 0.9999999855896213 ... 1.00000001438773
, 1.99999999832797   ... 2.0000000016716304
, 2.9999998953191924 ... 3.0000001327935335]

となりより精度の高い全ての根が含まれる区間を求めることが出来ました。

この記事のコードは以下のgistに公開しています。ぜひ自分でも試して遊んでみて下さい。
https://gist.github.com/lotz84/afa0ed91735bde23cc55220104c8fefa

あとがき

アドベントカレンダーの時期はみんなが大量に記事を書くから一つ一つの記事は短いほうが読まれやすいですよねという話を昨日のHaskellもくもく会で @igrep さんとしていた次の日にこの分量書いてしまった… ここまで来たみなさん、頑張って読んでいただいてありがとうございました精度保証付き数値計算は工学的な側面が強いのかと思いきや「ローレンツアトラクターは本当にストレンジアトラクターなのか？」という14番目のスメイルの問題と呼ばれる純粋な数学の問題を解く¹¹のにも使われていて非常に面白い分野です。実は今回紹介したallsolは重根が苦手だったりもするんですが、最初見たときは非常に魔法のようなアルゴリズムに見えました。こんなところにも自動微分が活用されているなんてますます自動微分の方も面白くなってきちゃいますね。それでは今回はこれで

おっと$1/10$が端点だから含まれていないことに気づいたそこのアナタ！鋭いですね〜 ↩
どちらもekmett先生作のライブラリです。すごい。 ↩
クラフチックと読むようです ↩
これは区間を集合としてみた時に全ての要素同士の演算結果を集めた集合として定義するという考え方と整合性があります ↩
本当は多変数でも考えられますがここでは簡単のため一変数に限ります。 ↩
ニュートン法を思い出すとなんとなくやりたいことが分かると思います ↩
証明は『精度保証付き数値計算の基礎』の定理6.3を参照 ↩
正確には$g(I)$は$I$の内部の部分集合になっている必要があります ↩
$A = [a, b], B = [c, d]$とし$C = [a, b] - [c, d] = [a - d, b - c]$とする。もし$c$が負の値であれば$b - c > b$より$C\nsubseteq A$、もし$d$が負の値であれば$a - d < a$より$C\nsubseteq A$、$c$が正の値である時$d$が負であることはなく、$d$が負の値である時$c$が正になることはないので一般に$A - B = C \subset A$は成り立たない。ただし$B = [0, 0]$であれば$C = A$。 ↩
ADやForwardは自動微分のライブラリであるadの型です ↩
精度保証付き数値計算の力学系への応用について ↩