share/lbfgs/maxima-lbfgs.lisp

   1 ; Copyright 2006 by Robert Dodier
   2 ; Released under the terms of the GNU General Public License.
   3
   4 ; Example 1:
   5 ; This is the same FOM as computed by FGCOMPUTE in the program sdrive.f
   6 #|
   7 load (lbfgs);
   8 t1[j] := 1 - u[j];
   9 t2[j] := 10*(u[j + 1] - u[j]^2);
  10 n : 8;
  11 FOM : sum (t1[2*j - 1]^2 + t2[2*j - 1]^2, j, 1, n/2);
  12 lbfgs_nfeval_max : 100;
  13 lbfgs_ncorrections : 25;
  14 lbfgs (FOM, '[u[1], u[2], u[3], u[4], u[5], u[6], u[7], u[8]],
  15        [-1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1], 1e-4, [1, 3]);
  16
  17   => [u[1] = 1.00000047321587,u[2] = 1.000000931806471,
  18       u[3] = 1.00000047321587,u[4] = 1.000000931806471,
  19       u[5] = 1.00000047321587,u[6] = 1.000000931806471,
  20       u[7] = 1.00000047321587,u[8] = 1.000000931806471]$
  21
  22  |#
  23
  24 ; Example 2:
  25 ; This computes the least-squares fit to a function A/(1 + exp(-B * (x - C)))
  26 #|
  27 load (lbfgs);
  28 FOM : '((1/length(X))*sum((F(X[i]) - Y[i])^2, i, 1, length(X)));
  29 X : [1, 2, 3, 4, 5];
  30 Y : [0, 0.5, 1, 1.25, 1.5];
  31 F(x) := A/(1 + exp(-B*(x - C)));
  32 ''FOM;
  33 estimates : lbfgs (FOM, '[A, B, C], [1, 1, 1], 1e-4, [1, 3]);
  34 plot2d ([F(x), [discrete, X, Y]], [x, -1, 6]), ''estimates;
  35  |#
  36
  37 ; estimates => [A = 1.461933911464101, B = 1.601593973254801, C = 2.528933072164855]
  38
  39 ; Example 3:
  40 ; Gradient of FOM specified directly (instead of constructing it via diff)
  41 #|
  42 load (lbfgs);
  43 F(a, b, c) := (a - 5)^2 + (b - 3)^4 + (c - 2)^6;
  44 F_grad : map (lambda ([x], diff (F(a, b, c), x)), [a, b, c]);
  45 estimates : lbfgs ([F(a, b, c), F_grad], [a, b, c], [0, 0, 0], 1e-4, [1, 0]);
  46  |#
  47
  48 ; estimates => [a = 5.000086823042934, b = 3.052395429705181, c = 1.927980629919583]
  49
  50 ; Example 4:
  51 ; same as Example 3, but FOM and gradient are computed by Lisp functions
  52 ; Construct Lisp functions suitable for this example via translate
  53 #|
  54 load (lbfgs);
  55 F(a, b, c) := (a - 5)^2 + (b - 3)^4 + (c - 2)^6;
  56 F1(a, b, c) := ''(diff(F(a, b, c), a));
  57 F2(a, b, c) := ''(diff(F(a, b, c), b));
  58 F3(a, b, c) := ''(diff(F(a, b, c), c));
  59 translate (F, F1, F2, F3);
  60 :lisp (mremprop '|$f| 'mexpr)
  61 :lisp (mremprop '|$f1| 'mexpr)
  62 :lisp (mremprop '|$f2| 'mexpr)
  63 :lisp (mremprop '|$f3| 'mexpr)
  64 estimates : lbfgs ('[F(a, b, c), [F1(a, b, c), F2(a, b, c), F3(a, b, c)]],
  65                    [a, b, c], [0, 0, 0], 1e-4, [1, 0]);
  66  |#
  67
  68 ; estimates => [a = 5.000086823042934, b = 3.052395429705181, c = 1.927980629919583]
  69
  70 (in-package :maxima)
  71
  72 (defmvar $lbfgs_nfeval_max 100)
  73 (defmvar $lbfgs_ncorrections 25)
  74
  75 (defmfun $lbfgs (FOM-expr x-list x-initial eps iprint-list)
  76
  77   (if
  78     (or (and (symbolp FOM-expr) (mfboundp FOM-expr))
  79         (and (consp FOM-expr) (eq (caar FOM-expr) 'lambda)))
  80     (merror "lbfgs: figure of merit, if specified alone, cannot be a function name or lambda."))
  81
  82   (let*
  83     ((n (length (cdr x-list)))
  84      (m $lbfgs_ncorrections)
  85      (nwork (+ (* n (+ (* 2 m) 1)) (* 2 m)))
  86
  87      (xtol +flonum-epsilon+)
  88      (iflag 0)
  89
  90      (scache (make-array n :element-type 'flonum))
  91      (w (make-array nwork :element-type 'flonum))
  92      (diag (make-array n :element-type 'flonum))
  93      (g (make-array n :element-type 'flonum))
  94      (x (make-array n :element-type 'flonum))
  95      (iprint (make-array 2 :element-type 'f2cl-lib:integer4))
  96      (diagco nil)
  97
  98      (return-value '((mlist)))
  99      (f 0.0)
 100
 101      FOM-function FOM-grad-lambda FOM-grad-expr FOM-grad-function)
 102
 103     (if ($listp FOM-expr)
 104       (progn
 105         (setq FOM-function (coerce-float-fun ($first FOM-expr) x-list))
 106         (setq FOM-grad-expr ($second FOM-expr))
 107         (setq FOM-grad-function (coerce-float-fun FOM-grad-expr x-list)))
 108       (progn
 109         (setq FOM-function (coerce-float-fun FOM-expr x-list))
 110         (setq FOM-grad-lambda `(lambda (x) (meval (list '($diff) ',FOM-expr x))))
 111         (setq FOM-grad-expr `((mlist) ,@(mapcar (coerce FOM-grad-lambda 'function) (cdr x-list))))
 112         (setq FOM-grad-function (coerce-float-fun FOM-grad-expr x-list))))
 113
 114     (setf (aref iprint 0) (nth 1 iprint-list))
 115     (setf (aref iprint 1) (nth 2 iprint-list))
 116     (setf diagco f2cl-lib:%false%)
 117
 118     (let (($numer t)) (setq x-initial ($float (meval x-initial))))
 119     ;; Fill X from the initial value, skipping over the Maxima mlist marker.
 120     (replace x x-initial :start2 1)
 121
 122     (common-lisp-user::/blockdata-lb2/)
 123
 124     (dotimes (nfeval $lbfgs_nfeval_max)
 125 ; (format t "hey nfeval ~S~%" nfeval)
 126 ; (format t "hey x ~S~%" x)
 127       (setf f (apply 'funcall `(,FOM-function ,@(coerce x 'list))))
 128       (let ((g-list (apply 'funcall `(,FOM-grad-function ,@(coerce x 'list)))))
 129 ; (format t "hey g-list ~S~%" g-list)
 130         (when (eq t g-list)
 131           (merror "Evaluation of gradient at ~M failed.  Bad initial point?~%"
 132                   (list* '(mlist) (coerce x 'list))))
 133         ;; Copy the answer from the grad function to the g array,
 134         ;; skipping over the mlist marker.
 135         (replace g g-list :start2 1))
 136 ; (format t "hey f ~S g ~S~%" f g)
 137
 138       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7
 139                                   var-8 var-9 var-10 var-11 var-12)
 140         (common-lisp-user::lbfgs n m x f g diagco diag iprint eps xtol w iflag scache)
 141         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7
 142                          var-8 var-9 var-10 var-12))
 143         (setf iflag var-11)
 144         ; Set return value to result of line search.
 145         ; That's what's returned if algorithm doesn't converge; better than nothing, I hope.
 146         (setq return-value (append '((mlist)) (mapcar #'(lambda (a b) `((mequal) ,a ,b)) (cdr x-list) (coerce scache 'list))))
 147         (cond
 148           ((eql iflag 0)
 149            (return)))))
 150
 151     return-value))