added minimizeVectorBFGS2

4 files changed, 47 insertions, 12 deletions

diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index a1d009b..b95765f 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -2,14 +2,14 @@
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.GSL.Minimization
-Copyright   :  (c) Alberto Ruiz 2006
+Copyright   :  (c) Alberto Ruiz 2006-9
 License     :  GPL-style
 
 Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses ffi
 
-Minimization of a multidimensional function Minimization of a multidimensional function using some of the algorithms described in:
+Minimization of a multidimensional function using some of the algorithms described in:
 
 <http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
 
@@ -17,6 +17,7 @@ Minimization of a multidimensional function Minimization of a multidimensional f
 -----------------------------------------------------------------------------
 module Numeric.GSL.Minimization (
     minimizeConjugateGradient,
+    minimizeVectorBFGS2,
     minimizeNMSimplex
 ) where
 
@@ -132,7 +133,7 @@ The path to the solution can be graphically shown by means of:
 
 @'Graphics.Plot.mplot' $ drop 2 ('toColumns' p)@
 
--}     
+-}
 minimizeConjugateGradient ::
        Double        -- ^ initial step size
     -> Double        -- ^ minimization parameter   
@@ -142,7 +143,27 @@ minimizeConjugateGradient ::
     -> ([Double] -> [Double])      -- ^ gradient
     -> [Double]             -- ^ starting point
     -> ([Double], Matrix Double)        -- ^ solution vector, and the optimization trajectory followed by the algorithm
-minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ do
+minimizeConjugateGradient = minimizeWithDeriv 0
+
+{- | Taken from the GSL manual:
+
+The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. This is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.
+
+The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's Practical Methods of Optimization, Algorithms 2.6.2 and 2.6.4. It supercedes the original bfgs routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance tol corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches).
+-}
+minimizeVectorBFGS2 ::
+       Double        -- ^ initial step size
+    -> Double        -- ^ minimization parameter tol
+    -> Double        -- ^ desired precision of the solution (gradient test)
+    -> Int           -- ^ maximum number of iterations allowed
+    -> ([Double] -> Double) -- ^ function to minimize
+    -> ([Double] -> [Double])      -- ^ gradient
+    -> [Double]             -- ^ starting point
+    -> ([Double], Matrix Double)        -- ^ solution vector, and the optimization trajectory followed by the algorithm
+minimizeVectorBFGS2 = minimizeWithDeriv 1
+
+
+minimizeWithDeriv method istep minimpar tol maxit f df xi = unsafePerformIO $ do
     let xiv = fromList xi
         n = dim xiv
         f' = f . toList
@@ -151,7 +172,7 @@ minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ d
     dfp <- mkVecVecfun (aux_vTov df')
     rawpath <- withVector xiv $ \xiv' ->
                     createMIO maxit (n+2)
-                         (c_minimizeConjugateGradient fp dfp istep minimpar tol (fi maxit) // xiv')
+                         (c_minimizeWithDeriv method fp dfp istep minimpar tol (fi maxit) // xiv')
                          "minimizeDerivV"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
@@ -160,9 +181,8 @@ minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ d
     freeHaskellFunPtr dfp
     return (sol,path)
 
-
 foreign import ccall "gsl-aux.h minimizeWithDeriv"
-    c_minimizeConjugateGradient :: FunPtr (CInt -> Ptr Double -> Double)
+    c_minimizeWithDeriv :: CInt -> FunPtr (CInt -> Ptr Double -> Double)
                                 -> FunPtr (CInt -> Ptr Double -> Ptr Double -> IO ())
                                 -> Double -> Double -> Double -> CInt
                                 -> TVM
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 77e793e..3802574 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -441,7 +441,7 @@ void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g)
 }
 
 // conjugate gradient
-int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*), 
+int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*, double*), 
                       double initstep, double minimpar, double tolgrad, int maxit, 
                       KRVEC(xi), RMAT(sol)) {
     REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
@@ -463,7 +463,11 @@ int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*),
     // Starting point
     KDVVIEW(xi);
     // conjugate gradient fr
-    T = gsl_multimin_fdfminimizer_conjugate_fr;
+    switch(method) {
+        case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }
+        case 1 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }
+        default: ERROR(BAD_CODE);
+    }
     s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
     gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
     do {
diff --git a/lib/Numeric/GSL/gsl-aux.h b/lib/Numeric/GSL/gsl-aux.h
index cd17ef0..e88322c 100644
--- a/lib/Numeric/GSL/gsl-aux.h
+++ b/lib/Numeric/GSL/gsl-aux.h
@@ -39,7 +39,7 @@ int polySolve(KRVEC(a), CVEC(z));
 int minimize(double f(int, double*), double tolsize, int maxit, 
                  KRVEC(xi), KRVEC(sz), RMAT(sol));
 
-int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*),
+int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*, double*),
                       double initstep, double minimpar, double tolgrad, int maxit, 
                       KRVEC(xi), RMAT(sol));
 
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 7758268..278df78 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -1,9 +1,9 @@
 {-# LANGUAGE CPP #-}
-{-# OPTIONS_GHC -fno-warn-unused-imports #-}
+{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Tests
-Copyright   :  (c) Alberto Ruiz 2007
+Copyright   :  (c) Alberto Ruiz 2007-9
 License     :  GPL-style
 
 Maintainer  :  Alberto Ruiz (aruiz at um dot es)
@@ -105,6 +105,16 @@ expmTest2 = expm nd2 :~15~: (2><2)
  , 2.718281828459045
  , 2.718281828459045 ]
 
+---------------------------------------------------------------------
+
+minimizationTest = TestList [ utest "minimization conj grad" (minim1 f df [5,7] ~~ [1,2])
+                            , utest "minimization bg2"       (minim2 f df [5,7] ~~ [1,2])
+                            ]
+    where f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
+          df [x,y] = [20*(x-1), 40*(y-2)]
+          a ~~ b = fromList a |~| fromList b
+          minim1 g dg ini = fst $ minimizeConjugateGradient 1E-2 1E-4 1E-3 30 g dg ini
+          minim2 g dg ini = fst $ minimizeVectorBFGS2 1E-2 1E-2 1E-3 30 g dg ini
 
 ---------------------------------------------------------------------
 
@@ -206,6 +216,7 @@ runTests n = do
         , exponentialTest
         , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8)
         , utest "polySolve" (polySolveProp [1,2,3,4])
+        , minimizationTest
         ]
     return ()