GSLLで最小値探索: one-dimensional minimizer

Posted on Jun 12, 2021

GNU Scientific Library for Lisp (GSLL)の軽い導入はGSLLの導入(Roswell)で行なったものとする。

基本的に使っている処理系はSBCL。

導入

GSLLには数多くの函数が含まれている。 GSLLのドキュメントとしてまとまったものは無いが、本家GSLのドキュメントはGSL - GNU Scientific LibraryとGSL reference manual Japanese translationとにあり、大変充実している。とくに日本語訳をしてくださっている富永大介さんに感謝。

上述のようにGSLのドキュメントは大変充実している一方でGSLLのドキュメントはほとんど無い。それゆえGSLLを使う際は基本的に

目的の操作・函数をGSLのドキュメントで確認
それに対応するCommon Lispでの函数をgsll:gsl-lookupを使うなどして調べる
GSLのドキュメントで使い方を確認し、Common Lispに読み替える

の手順を経ることになる、と思われる。

自分のようなプログラミングに不慣れな人からすると、ドキュメントはもちろん（サンプルコードも含めて）Cについて書かれているため、それを別の言語であるところのCommon Lispに置き換える手順3でわりと高いハードルを感じる。 gsll:gammaのような数値の引数を与えれば数値が返ってくるような普通の函数であれば容易に読み替えられるが、微分・積分・求根・最適化などになってくると把握に時間がかかる。

ドキュメントとにらめっこするにしても限界がある。ドキュメント以外にも直接的なヒントを得る方法がある。それはGSLLのソースコードにある例を見ることである。 Roswellを使って導入しているなら~/.roswell/lisp/quicklisp/dists/quicklisp/software/gsll-quicklisp-eeeda841-git/solve-minimize-fit/にあると思われる（人によってgsll-quicklisp-eeeda841-gitのバージョンが違うはず）。

GSLLを使う際は(ql:quickload :gsll)のあと忘れずに

(setq *default-read-float-format* 'double-float)

を評価。

最小値探索

GSLLに含まれる最小値探索に該当するものは、GSLのドキュメントの

28 Simulated Annealing （第25章シミュレーティド・アニーリング）
37 One Dimensional Minimization （第34章一次元関数の最適化）
39 Multidimensional Minimization （第36章多次元関数の最適化）

である。それぞれの詳しい説明はドキュメントで確認してもらうとして、以下ではGSLLでの使い方をソースコードの例に習って紹介する。

今回はone-dimensional minimizationのみ扱う。

One Dimensional Minimization

one-dimensional minimizationはminimization-one.lispに含まれている。

最小値を求める方法は、最小値が存在する区間を指定されたアルゴリズムに従って狭めていく方法が使われる。指定された回数か、指定された区間の幅に収まれば探索は終了する。

区間$[a, b] \ (a<b)$で函数$f$の最小値を探す際、初期値$c$は $$ f(a) > f(c) < f(b) $$ を満たす必要があることが注意点である（というよりもこの条件でちょっと手こずった）。初期値$c$はあくまで推定される最小値の位置であり、まずは境界値$f(a)$と$f(b)$とより小さい値$f(c)$（最小値の候補）を少なくとも1つ示さないと探索は始めない、ということらしい。このことからも感じるように、$f$として想定されるのはもちろん連続函数。

ソースコードにある例を若干修正して、

函数$f$について
区間$[a, b] \ (a<b)$で
$f(a) > f(c) < f(b)$を満たす最小値候補$c \ (a<c<b)$から出発し
Brent法を用いた

最小値探索を行なうコードを書けば、例えば

(defun minimize-one
    (a b f c &optional (minimizer-type gsll:+brent-fminimizer+))
  (let ((max-iter 100)
        (minimizer (gsll:make-one-dimensional-minimizer minimizer-type f c a b)))
    (loop :for iter :from 0
       :for min = (gsll:solution minimizer)
       :for lower = (gsll:fminimizer-x-lower minimizer)
       :for upper = (gsll:fminimizer-x-upper minimizer)
       :do
         (gsll:iterate minimizer)
       :while
         (and (< iter max-iter)
              (not (gsll:min-test-interval
                    lower upper
                    gsll:*default-absolute-error* gsll:*default-relative-error*)))
       :finally
         (return  (values min iter (- upper lower))))))

のように書ける。 gsll:make-one-dimensional-minimizerで最小値探索を行なうインスタンスminimizerを作り、それを用いてループ毎にminとlowerとupperとを更新する。ループ回数がmax-iterに達するか、区間の幅がgsll:min-test-intervalで指定される幅より小さくなったら終了し、最小値の位置（とループ回数と区間の幅と）を返す。

オプショナル引数のgsll:+brent-fminimizer+で最小値を求めるアルゴリズムを指定している。アルゴリズムは

gsll:+golden-section-fminimizer+: 黄金分割法のみを用いた方法、線形収束する。
gsll:+brent-fminimizer+: Brentの方法と呼ばれる方法、放物線による補間と黄金分割法とを組み合わせて収束性と安定性とを改良。
gsll:+quad-golden-fminimizer+: Brentの方法を改良した方法、ステップが区間から出ないように調整されている。

の3つがある。デフォルトではgsll:+brent-fminimizer+を用いるようにしている。

例として$f(x) = (x-\pi)^2$に適用すると

(flet ((test-f (x)
         (expt (- x pi) 2)))
  (minimize-one 0.0 4.0 #'test-f 3.0))
; 3.141592653589794
; 5
; 9.362675701396483e-8

(flet ((test-f (x)
         (expt (- x pi) 2)))
  (minimize-one 0.0 4.0 #'test-f 1.0))
; Evaluation aborted on #<GSLL:INVALID-ARGUMENT {10020D5683}>.

$f(0) \simeq 9.87, \ f(4) \simeq 0.737$より、最小値候補$f(3) \simeq 0.020$は条件を満たしている。一方後者は$f(1) \simeq 4.59$で境界値$f(4)$より大きいのでエラーが出る。

この初期値の設定がなかなかにやっかいなので、対処療法的に

(defun minimize-one& (a b f &optional (minimizer-type gsll:+brent-fminimizer+))
  (let* ((interval (- b a))
          (c (if (< (funcall f a) (funcall f b))
                 (+ a (* 1.0e-2 interval))
                 (- b (* 1.0e-2 interval)))))
    (minimize-one a b f c minimizer-type)))

とすれば手で与える必要はなくなる。

(flet ((test-f (x)
         (expt (- x pi) 2)))
  (minimize-one& 0.0 4.0 #'test-f))
; 3.141592653589793
; 5
; 9.362675701396483e-8

小さいほうの境界値からちょっとだけずらした点を使うことでエラーは出にくくなっている。こういう「ちょっと」で区間の長さの1.0e-2倍を使うのは行儀がよろしくないとは思うが、あくまで対処療法ということで……。