GSLLで最小値探索: multidimensional minimizer

Posted on Jun 19, 2021

GNU Scientific Library for Lisp (GSLL)の軽い導入はGSLLの導入(Roswell)で行なったものとする。また1次元方向における最小値探索は、前回の記事（GSLLで最小値探索: one-dimensional minimizer）を参考。

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

導入

前回に引き続き、~/.roswell/lisp/quicklisp/dists/quicklisp/software/gsll-quicklisp-eeeda841-git/solve-minimize-fit/（gsll-quicklisp-eeeda841-gitの名前はインストールしたバージョンによって異なるはず）にあるソースコードの例を微修正し、比較的使いやすい形にした多次元空間での最小値探索を行なう函数を紹介する。

GSLのドキュメント（本家、日本語訳）にあるように、GSLLに含まれる最小値探索に該当するものは

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

である。それぞれの詳しい説明はドキュメントで確認してください。

多次元空間における最小値探索という意味ではSimulated AnnealingとMultidimensional Minimizationの2つが該当する。今回は後者のみ扱う。

GSLLを使うときは(ql:quickload :gsll)のあとに

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

を評価することを忘れずに。

Multidimensional Minimization

多次元最小値探索とは$n$個の引数をもつ函数$f$が任意のパラメータの組$\{x_1, x_2, \dots, x_n\}$に対して $$ f(m_1, m_2, \dots, m_n) \leq f(x_1, x_2, \dots, x_n) $$ を満たすパラメータの組$\{m_1, m_2, \dots, m_n\}$を見つけることである。 $f$は下に有界であるとする。ソースコードの例で実際に評価する函数は $$ f(x_1, x_2) = 10(x_1-1)^2 + 20(x_2-2)^2 + 30 $$ である。最小値は$f(1,2)=30$である。

GSLLに含まれる多次元最小値探索は導函数を使うものと使わないものとに大別される。導函数使うものは共役勾配法（CG法）を用い、使わないものはシンプレックス法を用いる。順に紹介する。

minimize-multi-d
minimize-multi-de
minimize-multi-nd

導函数を用いる探索

導函数を用いる勾配法は

最小点の候補で勾配を計算する。
勾配とゼロとの差が許容範囲内に収まっているか確認する。
- （許容範囲に収まっていなければ）勾配に基づいて最小点の候補を更新して1に戻る。
- （許容範囲に収まっていれば）3へ進む。
候補を真の最小点とし、探索を終了する。

が基本的な流れであり、勾配法の一種として共役勾配法がある。実際には別途更新回数の上限が設けられ、その上限に達した時の候補を最小点として更新を打ち切ることもある。

GSLLに含まれる具体なアルゴリズムは

GSLL:+CONJUGATE-FLETCHER-REEVES+: フレッチャー・リーブズ法と呼ばれる共役勾配法。
GSLL:+CONJUGATE-POLAK-RIBIERE+: ポラック・リビエールによるの共役勾配法。
GSLL:+VECTOR-BFGS+: ブロイデン・フレッチャー・ゴールドファーブ・シャノによる共役勾配法。
GSLL:+VECTOR-BFGS2+: GSLL:+VECTOR-BFGS+を改良した方法。函数や導函数の評価回数が少なくて済む

の4つある。もともとのGSLには最急降下法アルゴリズムもあるが、実用性が弱いためGSLLには実装されていないようである。

GSLLでCG法を用いた最小値探索を行なうには、最小点を知りたい函数$f$以外に

導函数$\mathrm{d}f/\mathrm{d}x_i$ ($i=1, 2, \dots, n$)
函数$f$と導函数$\mathrm{d}f/\mathrm{d}x_i$とを両方計算する函数

を用意する必要がある。このときの書き方でも2パターンの書き方があり、函数$f$の引数としてパラメータの組$\{x_1, x_2, \dots, x_n\}$を

配列として渡す場合と
各要素を別々の引数として渡す場合と

がある。

まず1の配列として渡す場合、ソースコードにある例をほぼそのまま持ってきたものが

(defparameter *paraboloid-center* #(1.0d0 2.0d0))

(defun paraboloid-vector (xy)
  (let ((x (grid:aref xy 0))
        (y (grid:aref xy 1))
        (dp0 (aref *paraboloid-center* 0))
        (dp1 (aref *paraboloid-center* 1)))
    (+ (* 10 (expt (- x dp0) 2))
       (* 20 (expt (- y dp1) 2))
       30)))

(defun paraboloid-derivative (xy output)
  (let ((x (grid:aref xy 0))
        (y (grid:aref xy 1))
        (dp0 (aref *paraboloid-center* 0))
        (dp1 (aref *paraboloid-center* 1)))
    (setf (grid:aref output 0)
          (* 20 (- x dp0))
          (grid:aref output 1)
          (* 40 (- y dp1)))))

(defun paraboloid-and-derivative
    (arguments-gv-pointer value-pointer derivative-gv-pointer)
  (prog1
      (setf (grid:aref value-pointer 0)
            (paraboloid-vector arguments-gv-pointer))
    (paraboloid-derivative
     arguments-gv-pointer derivative-gv-pointer)))

である。これに対し、ソースコードの例から比較的使いやすいように書き換えた最小値を求める函数を

(defun minimize-multi-d
    (init-farray funcs &key (method gsll:+conjugate-fletcher-reeves+)
                         (step-size 0.01e0) (tolerance 1.0e-4))
  (let* ((minimizer
          (gsll:make-multi-dimensional-minimizer-fdf
           method (gsll:dim0 init-farray)
           funcs init-farray step-size tolerance nil)))
    (loop :with status = t
       :for iter :from 0 :below 100
       :while status
       :do
         (gsll:iterate minimizer)
         (setf status (not (gsll:min-test-gradient
                            (gsll:mfdfminimizer-gradient minimizer)
                            gsll:*default-absolute-error*)))
       :finally
         (return
           (values (gsll:solution minimizer)
                   (gsll:function-value minimizer))))))

としてみた。実際に計算させると

(minimize-multi-d (grid:make-foreign-array 'double-float
                                           :initial-contents '(5.0d0 7.0d0))
                  '(paraboloid-vector
                    paraboloid-derivative
                    paraboloid-and-derivative))
;; => #m(1.000000000000000d0 2.000000000000000d0)
;; => 30.0

となる。初期値を配列として渡しているので最小点も配列で返すようにした。 timeを使ってみると

Evaluation took:
  0.002 seconds of real time
  0.001767 seconds of total run time (0.001237 user, 0.000530 system)
  100.00% CPU
  3 lambdas converted
  4,061,614 processor cycles
  491,024 bytes consed

という結果が返ってきた。

次に各要素を別々の引数として与える場合は、こちらもソースコードの例をほぼそのまま持ってきたものが

(defun paraboloid-scalar (x y)
  (let ((dp0 (aref *paraboloid-center* 0))
        (dp1 (aref *paraboloid-center* 1)))
    (+ (* 10 (expt (- x dp0) 2))
       (* 20 (expt (- y dp1) 2))
       30)))

(defun paraboloid-derivative-scalar (x y)
  (let ((dp0 (aref *paraboloid-center* 0))
        (dp1 (aref *paraboloid-center* 1)))
    (values
     (* 20 (- x dp0))
     (* 40 (- y dp1)))))

(defun paraboloid-and-derivative-scalar (x y)
  (values-list
   (cons (paraboloid-scalar x y)
         (multiple-value-list (paraboloid-derivative-scalar x y)))))

である。こちらも比較的使いやすいであろう形に書き換えた最小値を求める函数を

(defun minimize-multi-de
    (init funcs &key (method gsll:+conjugate-fletcher-reeves+)
                  (step-size 0.01e0) (tolerance 1.0e-4))
  (labels ((cl-list1d (grid-1d)
             "one-dimensional grid:foreign-array -> list"
             (loop :for i :from 0 :below (grid:dim0 grid-1d)
                :collect (grid:aref grid-1d i))))
    (let* ((dim (length init))
           (initial
            (grid:make-foreign-array
             'double-float :initial-contents init))
           (minimizer
            (gsll:make-multi-dimensional-minimizer-fdf
             method dim
             funcs initial step-size tolerance t)))
      (declare (inline cl-list1d))
      (loop :with status = t
         :for iter :from 0 :below 100
         :while status
         :do
           (gsll:iterate minimizer)
           (setf status
                 (not (gsll:min-test-gradient
                       (gsll:mfdfminimizer-gradient minimizer)
                       gsll:*default-absolute-error*)))
         :finally
           (return
             (values (cl-list1d (gsll:solution minimizer))
                     (gsll:function-value minimizer)))))))

としてみた。計算してみると

(minimize-multi-de '(5.0d0 7.0d0)
                   '(paraboloid-scalar
                     paraboloid-derivative-scalar paraboloid-and-derivative-scalar))
;; => (0.9999999999999997 2.0)
;; => 30.0

となる。こちらは初期値として普通のリストを渡しているので最小点もリストで返すようにした。 timeを使ってみると

Evaluation took:
  0.016 seconds of real time
  0.015679 seconds of total run time (0.015679 user, 0.000000 system)
  100.00% CPU
  4 forms interpreted
  173 lambdas converted
  36,101,477 processor cycles
  2,584,384 bytes consed

となった。

両者を比較すると前者の方がパフォーマンスは良さそうである。後者の方は素のCommon Lispを扱うように利用できる利点がある。

導函数を用いない探索

導関数を用いないシンプレックス法は（$n$次元空間上で最小値探索を行なうとして）

（適当な長さ、適当な向きの）線形独立な$n$個の線分を生成する。
生成した$n$個の線分を最小点の候補を始点とした有効線分と見直し、$n$個の終点を得る。
始点と$n$個の終点とを結んで$n$次元単体（$n+1$個の頂点を結んで得られる多面体: 例えば$2$次元なら$3$つの頂点を結んで得られる三角形であり、$3$次元なら$4$つの頂点を結んで得られる三角錐）を考える。
単体の各頂点において函数値を計算し、値が最も大きい頂点を適当な変換により更新する（その過程で単体は拡大・縮小を繰り返す）。
単体の大きさが許容範囲に収まれば探索を終了する。
探索が終了した時、函数値が最も小さい頂点を最小点として返す。

の流れで行われる。 GSLLに含まれる具体的なアルゴリズムは

GSLL:+SIMPLEX-NELDER-MEAD-ON2+: 線形独立な線分を空間の各軸に平行な線分にとり、頂点を更新する適当な変換として鏡映変換（更新する頂点以外の頂点の位置とそこでの函数値とで重心を与え、その重心に対する対称変換）が用いられる。
GSLL:+SIMPLEX-NELDER-MEAD+: GSLL:+SIMPLEX-NELDER-MEAD-ON2+の計算量が$O(N^2)$であるのに対して$O(N)$に改良した方法。メモリの使用量は変わらずどちらも$O(N^2)$。
GSLL:+SIMPLEX-NELDER-MEAD-RANDOM+: GSLL:+SIMPLEX-NELDER-MEAD-ON2+で線形独立な線分を生成する際、乱数で生成した線分を用いる。単純な乱数を用いるため、規則的に線分が生成される恐れがあるそうである。

の3つである。

ソースコードの例にある函数は

(defparameter *paraboloid-center* #(1.0d0 2.0d0))

(defun paraboloid-scalar (x y)
  (let ((dp0 (aref *paraboloid-center* 0))
        (dp1 (aref *paraboloid-center* 1)))
    (+ (* 10 (expt (- x dp0) 2))
       (* 20 (expt (- y dp1) 2))
       30)))

である。これは導函数を用いる方法で配列を渡さない時の函数と全く同じである。

ソースコードの例を若干変更して比較的扱いやすいようにした最小値を求める函数を

(defun minimize-multi-nd
    (init func &key (method GSLL:+SIMPLEX-NELDER-MEAD+) (step-size 1.0e0))
  (labels ((cl-list1d (grid-1d)
             "one-dimensional grid:foreign-array -> list"
             (loop :for i :from 0 :below (grid:dim0 grid-1d)
                :collect (grid:aref grid-1d i))))
    (let* ((dim (length init))
           (step-size* (grid:make-foreign-array 'double-float
                                                :dimensions dim
                                                :initial-element step-size))
           (minimizer
            (gsll:make-multi-dimensional-minimizer-f
             method dim func
             (grid:make-foreign-array 'double-float :initial-contents init)
             step-size*)))
      (declare (inline cl-list1d))
      (loop :with status = t :and size
         :for iter :from 0 :below 100
         :while status
         :do
           (gsll:iterate minimizer)
           (setf size (gsll:size minimizer)
                 status (not (gsll:min-test-size size gsll:*default-absolute-error*)))
         :finally
           (return (values (cl-list1d (gsll:solution minimizer))
                           (gsll:function-value minimizer)))))))

としてみた。計算してみると

(minimize-multi-nd '(5.0d0 7.0d0) 'paraboloid-scalar)
;; => (0.9999963612837839 1.9999944070239908)
;; => 30.000000000758032

となる。 timeを使ってみると

Evaluation took:
  0.005 seconds of real time
  0.005048 seconds of total run time (0.005048 user, 0.000000 system)
  100.00% CPU
  33 lambdas converted
  11,603,590 processor cycles
  457,936 bytes consed

となる。