LBFGSを使ってみる

はじめに

f(x)=x^4 - 4x^3の最小値を、LBFGSという手法を用いたCommon Lispのプログラムを書いて求めます。
使うライブラリはcl-liblbfgs。liblbfgsというLBFGS法のCによる実装をCommon Lispから呼べるようにしたもので、私が書いたものです^^;

LBFGSは、従属変数に制約のない関数の最小値を求める問題を解いてくれます：minimize f(x), x = (x1, x2, ..., xN)
今回の場合、f(x)=x^4 - 4x^3で、N=1です。

cl-liblbfgsのインストール

1. liblbfgsのビルド
2. https://github.com/mhkoji/cl-liblbfgs をロード
3. liblbfgsをビルドしてliblbfgs.soというファイルが生成されたことを確認したら、(push <path/to/liblbfgs.so> cffi:*foreign-library-directories*)

最適化のコード

(use-package :cl-liblbfgs)

(cffi:defcallback progress :int
    ((instance  :pointer)
     (xs        :pointer)
     (g         :pointer)
     (fx        lbfgsfloatval-t)
     (xnorm     lbfgsfloatval-t)
     (gnorm     lbfgsfloatval-t)
     (step      lbfgsfloatval-t)
     (n         :int)
     (k         :int)
     (ls        :int))
  (declare (ignore instance xs g n ls))
  (format t "Iteration ~A:~%" k)
  (format t "  fx = ~A~%" fx)
  (format t "  xnorm = ~A, gnorm = ~A, step = ~A%" xnorm gnorm step)
  (format t "~%")
  0)


;; f(x) = x^4 - 4x^3
;; g(x) = 4x^3 - 12x^2
(cffi:defcallback eval-1 lbfgsfloatval-t
    ((instance :pointer)
     (xs       :pointer)
     (g        :pointer)
     (n        :int)
     (step     lbfgsfloatval-t))
  (declare (ignore instance step n))
  (let ((x (cffi:mem-ref xs 'lbfgsfloatval-t 0)))
    (let ((fx (- (* 1.0d0 x x x x) (* 4 x x x)))
          (gx (- (* 1.0d0 4 x x x) (* 12 x x))))
      (setf (cffi:mem-aref g 'lbfgsfloatval-t 0) gx)
      fx)))

(defun main-1 ()
  (cffi:with-foreign-objects ((xs 'lbfgsfloatval-t 1)
                              (param '(:struct lbfgs-parameter-t)))
    ;; 0.0d0だと死ぬ
    (setf (cffi:mem-aref xs 'lbfgsfloatval-t 0) 100.0d0)
    (lbfgs-parameter-init param)
    (let ((ret (lbfgs 1
                      xs
                      (cffi:null-pointer)
                      (cffi:callback eval-1)
                      (cffi:callback progress)
                      (cffi:null-pointer)
                      param)))
      (format t "L-BFGS optimization terminated with status code = ~A~%"
              ret)
      (format t "  x = ~A~%" (cffi:mem-aref xs 'lbfgsfloatval-t 0)))))

eval-1がライブラリから繰り返し呼ばれ、各xでのf(x)の値と、f(x)の微分であるg(x)の値を求めます。
LBFGSはこれらの値をもとに次のxの値を決めるようです。f(x)は関数の返り値となり、g(x)はgに直接代入します。

これを実行すると、最終的にxが3付近の値になります。
f(x)=x^4 - 4x^3はx=3のときに最小となるので、だいたいあっているわけです。