Optimization with GNU Scientific Library for Lisp

Description

I have been itching to write a post on writing an optimization/minimization code using GNU Scientific Library for Lisp (GSLL). For this post, I had to go into the GSLL code itself which contains a good set of examples. I highly recommend reading it. So, here it goes:

How to use GSLL univariate minimization routines

Let me start with optimizing a real function with one variable. For a theoretical survey on the subject of one dimensional optimization problems I recommend Wotao Yin’s lecture notes.

First, let us load up GSLL:

(require :gsll)

(SB-CLTL2)

Next, let us define the function we are going to optimize:

(defun fn (x)
  (+ (expt x 4) (* -1.4d1 (expt x 3)) (* 6.0d1 (expt x 2)) (* -7.0d1 x)))

FN

In GSLL, I need to create an minimizer which takes a function, an interval and an initial point.

(defparameter minimizer
  (gsl:make-one-dimensional-minimizer
   gsl:+quad-golden-fminimizer+
   #'fn
   1d0 0.0d0 2.0d0
   (fn 1d0) (fn 0.0d0) (fn 2.0d0)))

MINIMIZER

I have used a specific minimizer but there are few others. Now, I will iterate the minimizer until it converges unless it takes more than 1000 steps:

(let ((absolute-error 1.0d-3)
      (size 1.0d0)
      (flag nil))
  (do ((i 0 (1+ i)))
      ((or flag (> i 1000))
       (cons (gsl:solution minimizer) i))
    (gsl:iterate minimizer)
    (setf size (- (gsl:fminimizer-x-upper minimizer)
                  (gsl:fminimizer-x-lower minimizer)))
    (setf flag (gsl:min-test-size size absolute-error))))

(0.7809263974476345d0 . 16)

How to use GSLL multivariate minimization routines

For the multivariate optimization function, I am going to re-implement ordinary least square regression as my running example. First, I am going to need data on which I am going to fit a regression line:

(defparameter data
  (let ((a 2.0d0)
        (b 1.0d0))
    (loop repeat 200 collect
         (let ((x (random 1.0d0)))
            (cons x (+ (* a x) b (* 0.1 (random 1.0d0))))))))

DATA

Next, the function we need to optimize. This is the ordinary euclidean distance between the data points and their interpolation:

(defun fn (u v)
  (sqrt (loop for point in data sum
             (let ((x (car point))
                   (y (cdr point)))
               (expt (- (+ (* u x) v) y) 2)))))

FN

Of course, we need a minimizer:

(defparameter minimizer
  (gsl:make-multi-dimensional-minimizer-f
   gsl:+simplex-nelder-mead+ ;; method
   2                         ;; dimension
   #'fn
   (grid:make-foreign-array 'double-float :initial-contents '(0.0d0 0.0d0))      ;; initial approximation
   (grid:make-foreign-array 'double-float :initial-contents '(1.0d-2 1.0d-2))))  ;; step-size

MINIMIZER

Let us roll the ball:

(let ((absolute-error 1.0d-3)
      (flag nil))
  (do ((i 0 (1+ i)))
      ((or flag (> i 1000))
       (cons (gsl:solution minimizer) i))
    (gsl:iterate minimizer)
    (setf flag (gsl:min-test-size (gsl:size minimizer) absolute-error))))

((1.998 1.052) . 46)

We get an approximate solution in 46 steps. Nice!

An analysis

In the past here and here, I used minimizers that myself wrote. But, as pleasurable as it was to write my own routines, it is safer and less error-prone to use a highly optimized well-written stable specialized library such as GSLL.