diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,7 +6,17 @@ and this project adheres to the [Haskell Package Versioning Policy](https://pvp.haskell.org/). -## Unreleased +## 0.1.1.0 - 2023-06-21 + +* Support L-BFGS-B algorithm (when `with-lbfgsb` is enabled) +* Add some algorithm specific parameters +* Add instructions for installing dependent libraries +* Add `with-lbfgs` flag, which is `true` by default, but you can turn-off + the flag to build without L-BFGS. +* Add some instances of standard type classes: `Eq OptimizationException`, + `Show Result`, and `Show Statistics`. +* Return correct statistics for L-BFGS and L-BFGS-B. +* Fix many bugs ## 0.1.0.1 - 2023-06-03 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -6,6 +6,8 @@ Unified interface to various numerical optimization algorithms. +The aim of the package is to provide a convenient interface like Python's [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html). + Note that the package name is numeric-optimization and not numeri**cal**-optimization. The name `numeric-optimization` comes from the module name `Numeric.Optimization`. @@ -44,8 +46,69 @@ |---------|-------------------|---------------|-| |CG\_DESCENT|[CG_DESCENT-C](https://www.math.lsu.edu/~hozhang/SoftArchive/CG_DESCENT-C-3.0.tar.gz)|[nonlinear-optimization](https://hackage.haskell.org/package/nonlinear-optimization)|Requires `with-cg-descent` flag| |Limited memory BFGS (L-BFGS)|[liblbfgs](https://github.com/chokkan/liblbfgs)|[lbfgs](https://hackage.haskell.org/package/lbfgs)| +|Limited memory BFGS with bounds constraints (L-BFGS-B)|[L-BFGS-B](http://users.iems.northwestern.edu/~nocedal/lbfgsb.html)|[l-bfgs-b](https://hackage.haskell.org/package/l-bfgs-b)|Requires `with-lbfgsb` flag| |Newton's method|Pure Haskell implementation using [HMatrix](https://hackage.haskell.org/package/hmatrix)|-| +## Installation + +### Installing Prerequisites + +#### BLAS and LAPACK + +You may need to install BLAS and LAPACK for `hmatrix`. + +##### Windows (MSYS2): +``` +$ pacman -S mingw-w64-x86_64-lapack +``` + +or if you use MSYS2 installed by `stack` + +``` +$ stack exec -- pacman -S mingw-w64-x86_64-lapack +``` + +##### Debian and Ubuntu Linux: +``` +$ apt-get install libblas-dev liblapack-dev +``` + +`libblas-dev` and `liblapack-dev` are reference implementations. +You need to install optimized ones for better performance. +(See [DebianScience/LinearAlgebraLibraries](https://wiki.debian.org/DebianScience/LinearAlgebraLibraries)) + + +##### macOS + +By default `hmatrix` uses BLAS and LAPACK provided by Accelerate Framework provided by macOS. + +#### liblbfgsb + +If you want to use L-BFGS-B, you have to install the development package of `liblbfgsb`. + +##### Ubuntu Linux: +``` +$ apt-get install liblbfgsb-dev +``` + +##### Homebrew (macOS and Linux): +``` +$ brew install msakai/tap/liblbfgsb +``` + +##### Windows (MSYS2): +``` +$ wget https://github.com/msakai/mingw-w64-liblbfgsb/releases/download/v3.0-1/mingw-w64-x86_64-liblbfgsb-3.0-1-any.pkg.tar.zst +$ pacman -U mingw-w64-x86_64-liblbfgsb-3.0-1-any.pkg.tar.zst +``` + +or if you use MSYS2 installed by `stack` + +``` +$ wget https://github.com/msakai/mingw-w64-liblbfgsb/releases/download/v3.0-1/mingw-w64-x86_64-liblbfgsb-3.0-1-any.pkg.tar.zst +$ stack exec -- pacman -Sy +$ stack exec -- pacman -U mingw-w64-x86_64-liblbfgsb-3.0-1-any.pkg.tar.zst +``` ## Related Packages diff --git a/numeric-optimization.cabal b/numeric-optimization.cabal --- a/numeric-optimization.cabal +++ b/numeric-optimization.cabal @@ -5,10 +5,10 @@ -- see: https://github.com/sol/hpack name: numeric-optimization -version: 0.1.0.1 +version: 0.1.1.0 synopsis: Unified interface to various numerical optimization algorithms description: Please see the README on GitHub at -category: Math, Algorithms, Optimisation, Optimization +category: Math, Algorithms, Optimisation, Optimization, Numeric, Numerical homepage: https://github.com/msakai/nonlinear-optimization-ad#readme bug-reports: https://github.com/msakai/nonlinear-optimization-ad/issues author: Masahiro Sakai @@ -42,6 +42,16 @@ manual: True default: False +flag with-lbfgs + description: Enable L-BFGS (since 0.1.1.0) + manual: True + default: True + +flag with-lbfgsb + description: Enable L-BFGS-B (since 0.1.1.0) + manual: True + default: False + library exposed-modules: Numeric.Optimization @@ -55,7 +65,7 @@ , constraints , data-default-class >=0.1.2.0 && <0.2 , hmatrix >=0.20.0.0 - , lbfgs ==0.1.* + , numeric-limits ==0.1.* , primitive >=0.6.4.0 , vector >=0.12.0.2 && <0.14 default-language: Haskell2010 @@ -65,6 +75,18 @@ nonlinear-optimization >=0.3.7 && <0.4 else cpp-options: + if flag(with-lbfgs) + cpp-options: -DWITH_LBFGS + build-depends: + lbfgs ==0.1.* + else + cpp-options: + if flag(with-lbfgsb) + cpp-options: -DWITH_LBFGSB + build-depends: + l-bfgs-b >=0.1.0.1 && <0.2 + else + cpp-options: executable rosenbrock main-is: rosenbrock.hs @@ -98,6 +120,7 @@ , base >=4.12 && <5 , containers >=0.6.0.1 && <0.7 , data-default-class >=0.1.2.0 && <0.2 + , hmatrix , hspec >=2.7.1 && <3.0 , numeric-optimization , vector >=0.12.0.2 && <0.14 diff --git a/src/Numeric/Optimization.hs b/src/Numeric/Optimization.hs --- a/src/Numeric/Optimization.hs +++ b/src/Numeric/Optimization.hs @@ -16,7 +16,7 @@ -- Stability : provisional -- Portability : non-portable -- --- This module aims to provides unifined interface to various numerical +-- This module aims to provide unified interface to various numerical -- optimization, like [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html) in Python. -- -- In this module, you need to explicitly provide the function to calculate the @@ -63,7 +63,9 @@ , hasOptionalDict ) where +import Control.Applicative import Control.Exception +import Control.Monad import Control.Monad.Primitive import Control.Monad.ST import Data.Coerce @@ -78,12 +80,21 @@ import qualified Data.Vector.Generic.Mutable as VGM import qualified Data.Vector.Storable.Mutable as VSM import Foreign.C +#ifdef WITH_LBFGS import qualified Numeric.LBFGS.Vector as LBFGS +import qualified Numeric.LBFGS.Raw as LBFGS (unCLBFGSResult, lbfgserrCanceled) +#endif #ifdef WITH_CG_DESCENT import qualified Numeric.Optimization.Algorithms.HagerZhang05 as CG #endif +#ifdef WITH_LBFGSB +import qualified Numeric.LBFGSB as LBFGSB +import qualified Numeric.LBFGSB.Result as LBFGSB +#endif +import Numeric.Limits import Numeric.LinearAlgebra (Matrix) import qualified Numeric.LinearAlgebra as LA +import System.IO.Unsafe -- | Selection of numerical optimization algorithms @@ -117,8 +128,25 @@ -- * [2] -- -- * [3] + | LBFGSB + -- ^ Limited memory BFGS algorithm with bound constraints (L-BFGS-B) [1][2][3] + -- + -- The implementation is provided by l-bfgs-b package [4] + -- which is a bindign to L-BFGS-B fortran code [5]. + -- + -- * [1] R. H. Byrd, P. Lu and J. Nocedal. [A Limited Memory Algorithm for Bound Constrained Optimization](http://www.ece.northwestern.edu/~nocedal/PSfiles/limited.ps.gz), (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208. + -- + -- * [2] C. Zhu, R. H. Byrd and J. Nocedal. [L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization](http://www.ece.northwestern.edu/~nocedal/PSfiles/lbfgsb.ps.gz) (1997), ACM Transactions on Mathematical Software, Vol 23, Num. 4, pp. 550-560. + -- + -- * [3] J. L. Morales and J. Nocedal. [L-BFGS-B: Remark on Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization](http://www.ece.northwestern.edu/~morales/PSfiles/acm-remark.pdf) (2011), ACM Transactions on Mathematical Software, Vol 38, Num. 7, pp. 1–4 + -- + -- * [4] + -- + -- * [5] + -- + -- @since 0.1.1.0 | Newton - -- ^ Native implementation of Newton method + -- ^ Naïve implementation of Newton method in Haskell -- -- This method requires both gradient and hessian. deriving (Eq, Ord, Enum, Show, Bounded) @@ -126,12 +154,21 @@ -- | Whether a 'Method' is supported under the current environment. isSupportedMethod :: Method -> Bool +#ifdef WITH_LBFGS isSupportedMethod LBFGS = True +#else +isSupportedMethod LBFGS = False +#endif #ifdef WITH_CG_DESCENT isSupportedMethod CGDescent = True #else isSupportedMethod CGDescent = False #endif +#ifdef WITH_LBFGSB +isSupportedMethod LBFGSB = True +#else +isSupportedMethod LBFGSB = False +#endif isSupportedMethod Newton = True @@ -139,16 +176,65 @@ -- -- TODO: -- --- * How to pass algorithm specific parameters? +-- * Better way to pass algorithm specific parameters? -- --- * Separate 'callback' from other more concrete serializeable parameters? +-- * Separate 'paramsCallback' from other more concrete serializeable parameters? data Params a = Params { paramsCallback :: Maybe (a -> IO Bool) -- ^ If callback function returns @True@, the algorithm execution is terminated. , paramsTol :: Maybe Double - -- ^ Tolerance for termination. When 'tol' is specified, the selected algorithm sets - -- some relevant solver-specific tolerance(s) equal to 'tol'. + -- ^ Tolerance for termination. When @tol@ is specified, the selected algorithm sets + -- some relevant solver-specific tolerance(s) equal to @tol@. + -- + -- If specified, this value is used as defaults for 'paramsFTol' and 'paramsGTol'. + , paramsFTol :: Maybe Double + -- ^ 'LBFGS' stops iteration when delta-based convergence test + -- (see 'paramsPast') is enabled and the following condition is + -- met: + -- + -- \[ + -- \left|\frac{f' - f}{f}\right| < \mathrm{ftol}, + -- \] + -- + -- where @f'@ is the objective value of @past@ ('paramsPast') iterations ago, + -- and @f@ is the objective value of the current iteration. + -- The default value is @1e-5@. + -- + -- 'LBFGSB' stops iteration when the following condition is met: + -- + -- \[ + -- \frac{f^k - f^{k+1}}{\mathrm{max}\{|f^k|,|f^{k+1}|,1\}} \le \mathrm{ftol}. + -- \] + -- + -- The default value is @1e7 * ('epsilon' :: Double) = 2.220446049250313e-9@. + -- + -- @since 0.1.1.0 + , paramsGTol :: Maybe Double + -- ^ 'LBFGSB' stops iteration when \(\mathrm{max}\{|\mathrm{pg}_i| \mid i = 1, \ldots, n\} \le \mathrm{gtol}\) + -- where \(\mathrm{pg}_i\) is the i-th component of the projected gradient. + -- + -- @since 0.1.1.0 + , paramsMaxIters :: Maybe Int + -- ^ Maximum number of iterations. + -- + -- Currently only 'LBFGSB', 'CGDescent', and 'Newton' uses this. + -- + -- @since 0.1.1.0 + , paramsPast :: Maybe Int + -- ^ Distance for delta-based convergence test in 'LBFGS' + -- + -- This parameter determines the distance, in iterations, to compute + -- the rate of decrease of the objective function. If the value of this + -- parameter is @Nothing@, the library does not perform the delta-based + -- convergence test. The default value is @Nothing@. + -- + -- @since 0.1.1.0 + , paramsMaxCorrections :: Maybe Int + -- ^ The maximum number of variable metric corrections used in 'LBFGSB' + -- to define the limited memory matrix. + -- + -- @since 0.1.1.0 } instance Default (Params a) where @@ -156,6 +242,11 @@ Params { paramsCallback = Nothing , paramsTol = Nothing + , paramsFTol = Nothing + , paramsGTol = Nothing + , paramsMaxIters = Nothing + , paramsPast = Nothing + , paramsMaxCorrections = Nothing } instance Contravariant Params where @@ -185,6 +276,7 @@ , resultStatistics :: Statistics -- ^ Statistics of optimizaion process } + deriving (Show) instance Functor Result where fmap f result = @@ -203,18 +295,23 @@ -- ^ Total number of function evaluations. , gradEvals :: Int -- ^ Total number of gradient evaluations. + , hessianEvals :: Int + -- ^ Total number of hessian evaluations. , hessEvals :: Int -- ^ Total number of hessian evaluations. } + deriving (Show) +{-# DEPRECATED hessEvals "Use hessianEvals instead" #-} + -- | The bad things that can happen when you use the library. data OptimizationException = UnsupportedProblem String | UnsupportedMethod Method | GradUnavailable | HessianUnavailable - deriving (Show) + deriving (Show, Eq) instance Exception OptimizationException @@ -226,8 +323,8 @@ -- -- In the simplest case, @'VS.Vector' Double -> Double@ is a instance -- of 'IsProblem' class. It is enough if your problem does not have --- constraints and the selected algorithm does not further information --- (e.g. gradients and hessians), +-- constraints and the selected algorithm does not require further +-- information (e.g. gradients and hessians), -- -- You can equip a problem with other information using wrapper types: -- @@ -322,7 +419,7 @@ -- | Bounds for unconstrained problems, i.e. (-∞,+∞). boundsUnconstrained :: Int -> V.Vector (Double, Double) -boundsUnconstrained n = V.replicate n (-1/0, 1/0) +boundsUnconstrained n = V.replicate n (-infinity, infinity) -- | Whether all lower bounds are -∞ and all upper bounds are +∞. isUnconstainedBounds :: V.Vector (Double, Double) -> Bool @@ -374,10 +471,18 @@ Just Dict -> minimize_CGDescent Nothing -> \_ _ _ -> throwIO GradUnavailable #endif +#ifdef WITH_LBFGS minimize LBFGS = case optionalDict @(HasGrad prob) of Just Dict -> minimize_LBFGS Nothing -> \_ _ _ -> throwIO GradUnavailable +#endif +#ifdef WITH_LBFGSB +minimize LBFGSB = + case optionalDict @(HasGrad prob) of + Just Dict -> minimize_LBFGSB + Nothing -> \_ _ _ -> throwIO GradUnavailable +#endif minimize Newton = case optionalDict @(HasGrad prob) of Nothing -> \_ _ _ -> throwIO GradUnavailable @@ -399,6 +504,10 @@ cg_params = CG.defaultParameters { CG.printFinal = False + , CG.maxItersFac = + case paramsMaxIters params of + Nothing -> CG.maxItersFac CG.defaultParameters + Just m -> fromIntegral m / fromIntegral (VG.length x0) } mf :: forall m. PrimMonad m => CG.PointMVector m -> m Double @@ -456,12 +565,15 @@ , funcEvals = fromIntegral $ CG.funcEvals stat , gradEvals = fromIntegral $ CG.gradEvals stat , hessEvals = 0 + , hessianEvals = 0 } } #endif +#ifdef WITH_LBFGS + minimize_LBFGS :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double)) minimize_LBFGS _params prob _ | not (isNothing (bounds prob)) = throwIO (UnsupportedProblem "LBFGS does not support bounds") minimize_LBFGS _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "LBFGS does not support constraints") @@ -471,8 +583,8 @@ let lbfgsParams = LBFGS.LBFGSParameters - { LBFGS.lbfgsPast = Nothing - , LBFGS.lbfgsDelta = fromMaybe 0 $ paramsTol params + { LBFGS.lbfgsPast = paramsPast params + , LBFGS.lbfgsDelta = fromMaybe 1e-5 $ paramsFTol params <|> paramsTol params , LBFGS.lbfgsLineSearch = LBFGS.DefaultLineSearch , LBFGS.lbfgsL1NormCoefficient = Nothing } @@ -505,7 +617,7 @@ x <- VG.freeze (coerce xvec :: VSM.IOVector Double) #endif callback x - return $ if shouldStop then 1 else 0 + return $ if shouldStop then fromIntegral (LBFGS.unCLBFGSResult LBFGS.lbfgserrCanceled) else 0 (result, x_) <- LBFGS.lbfgs lbfgsParams evalFun progressFun instanceData (VG.toList x0) let x = VG.fromList x_ @@ -546,6 +658,7 @@ LBFGS.InvalidParameters -> (False, "A logic error (negative line-search step) occurred.") LBFGS.IncreaseGradient -> (False, "The current search direction increases the objective function value.") + iters <- readIORef iterRef nEvals <- readIORef evalCounter return $ @@ -559,54 +672,119 @@ , resultHessianInv = Nothing , resultStatistics = Statistics - { totalIters = undefined - , funcEvals = nEvals + 1 - , gradEvals = nEvals + 1 + { totalIters = iters + , funcEvals = nEvals + 1 -- +1 is for computing 'resultValue' + , gradEvals = nEvals , hessEvals = 0 + , hessianEvals = 0 } } +#endif + +#ifdef WITH_LBFGSB + +minimize_LBFGSB :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double)) +minimize_LBFGSB _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "LBFGSB does not support constraints") +minimize_LBFGSB params prob x0 = do + funcEvalRef <- newIORef (0::Int) + gradEvalRef <- newIORef (0::Int) + + let bounds' = + case bounds prob of + Nothing -> [] + Just vec -> map convertB (VG.toList vec) + convertB (lb, ub) = + ( if isInfinite lb && lb < 0 + then Nothing + else Just lb + , if isInfinite ub && ub > 0 + then Nothing + else Just ub + ) + func' x = unsafePerformIO $ do + modifyIORef' funcEvalRef (+1) + evaluate (func prob x) + grad' x = unsafePerformIO $ do + modifyIORef' gradEvalRef (+1) + evaluate (grad prob x) + + let -- | @m@: The maximum number of variable metric corrections used + -- to define the limited memory matrix. /Suggestion:/ @5@. + m = fromMaybe 5 (paramsMaxCorrections params) + + -- | @factr@: Iteration stops when the relative change in function value + -- is smaller than @factr*eps@, where @eps@ is a measure of machine precision + -- generated by the Fortran code. @1e12@ is low accuracy, @1e7@ is moderate, + -- and @1e1@ is extremely high. Must be @>=1@. /Suggestion:/ @1e7@. + factr = fromMaybe 1e7 $ (/ epsilon) <$> (paramsFTol params <|> paramsTol params) + + -- ^ @pgtol@: Iteration stops when the largest component of the projected + -- gradient is smaller than @pgtol@. Must be @>=0@. /Suggestion:/ @1e-5@. + pgtol = fromMaybe 1e-5 $ paramsGTol params <|> paramsTol params + + -- | @'Just' steps@ means the minimization is aborted if it has not converged after + -- @steps>0@ iterations. 'Nothing' signifies no limit. + steps = paramsMaxIters params + + result <- evaluate $ LBFGSB.minimize m factr pgtol steps bounds' x0 func' grad' + + let x = LBFGSB.solution result + (success, msg) = + case LBFGSB.stopReason result of + LBFGSB.Converged -> (True, "The solution converged.") + LBFGSB.StepCount -> (False, "The number of steps exceeded the user's request.") + LBFGSB.Other msg -> (False, msg) + + funcEvals <- readIORef funcEvalRef + gradEvals <- readIORef gradEvalRef + + return $ + Result + { resultSuccess = success + , resultMessage = msg + , resultSolution = x + , resultValue = func prob x + , resultGrad = Nothing + , resultHessian = Nothing + , resultHessianInv = Nothing + , resultStatistics = + Statistics + { totalIters = length (LBFGSB.backtrace result) + , funcEvals = funcEvals + , gradEvals = gradEvals + , hessEvals = 0 + , hessianEvals = 0 + } + } + +#endif + + minimize_Newton :: (HasGrad prob, HasHessian prob) => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double)) minimize_Newton _params prob _ | not (isNothing (bounds prob)) = throwIO (UnsupportedProblem "Newton does not support bounds") minimize_Newton _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "Newton does not support constraints") minimize_Newton params prob x0 = do let tol = fromMaybe 1e-6 (paramsTol params) - loop !x !y !g !h !n = do - shouldStop <- - case paramsCallback params of - Just callback -> callback x - Nothing -> return False - if shouldStop then do - return $ - Result - { resultSuccess = False - , resultMessage = "The minimization process has been canceled." - , resultSolution = x - , resultValue = y - , resultGrad = Just g - , resultHessian = Just h - , resultHessianInv = Nothing - , resultStatistics = - Statistics - { totalIters = n - , funcEvals = n - , gradEvals = n - , hessEvals = n - } - } - else do - let p = h LA.<\> g - x' = VG.zipWith (-) x p - if LA.norm_Inf (VG.zipWith (-) x' x) > tol then do - let (y', g') = grad' prob x' - h' = hessian prob x' - loop x' y' g' h' (n+1) - else do + + loop !x !y !g !h !iter = do + shouldStop <- msum <$> sequence + [ pure $ case paramsMaxIters params of + Just maxIter | maxIter <= iter -> Just "maximum number of iterations reached" + _ -> Nothing + , case paramsCallback params of + Nothing -> return Nothing + Just callback -> do + flag <- callback x + return $ if flag then Just "The minimization process has been canceled." else Nothing + ] + case shouldStop of + Just reason -> return $ Result - { resultSuccess = True - , resultMessage = "success" + { resultSuccess = False + , resultMessage = reason , resultSolution = x , resultValue = y , resultGrad = Just g @@ -614,15 +792,43 @@ , resultHessianInv = Nothing , resultStatistics = Statistics - { totalIters = n - , funcEvals = n - , gradEvals = n - , hessEvals = n + { totalIters = iter + , funcEvals = iter + 1 + , gradEvals = iter + 1 + , hessEvals = iter + 1 + , hessianEvals = iter + 1 } } + Nothing -> do + let p = h LA.<\> g + x' = VG.zipWith (-) x p + if LA.norm_Inf (VG.zipWith (-) x' x) > tol then do + let (y', g') = grad' prob x' + h' = hessian prob x' + loop x' y' g' h' (iter + 1) + else do + return $ + Result + { resultSuccess = True + , resultMessage = "success" + , resultSolution = x + , resultValue = y + , resultGrad = Just g + , resultHessian = Just h + , resultHessianInv = Nothing + , resultStatistics = + Statistics + { totalIters = iter + , funcEvals = iter + 1 + , gradEvals = iter + 1 + , hessEvals = iter + 1 + , hessianEvals = iter + 1 + } + } + let (y0, g0) = grad' prob x0 h0 = hessian prob x0 - loop x0 y0 g0 h0 1 + loop x0 y0 g0 h0 0 -- ------------------------------------------------------------------------ diff --git a/test/IsClose.hs b/test/IsClose.hs --- a/test/IsClose.hs +++ b/test/IsClose.hs @@ -1,4 +1,5 @@ {-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} module IsClose @@ -30,6 +31,7 @@ import qualified Data.Vector.Storable as VS import qualified Data.Vector.Unboxed as VU import GHC.Stack (HasCallStack) +import Numeric.LinearAlgebra as LA import Test.HUnit import Text.Printf @@ -121,6 +123,11 @@ instance (AllClose r v, VU.Unbox v) => AllClose r (VU.Vector v) where allCloseRaw tol xs ys | VG.length xs == VG.length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (VG.toList xs) (VG.toList ys)]) + | otherwise = Ap Nothing + +instance (AllClose r v, Num v, LA.Container Vector v) => AllClose r (LA.Matrix v) where + allCloseRaw tol xs ys + | LA.size xs == LA.size ys = allCloseRaw tol (flatten xs) (flatten ys) | otherwise = Ap Nothing -- ------------------------------------------------------------------------ diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,21 +1,224 @@ {-# LANGUAGE OverloadedLists #-} +{-# LANGUAGE LambdaCase #-} import Test.Hspec +import Control.Exception +import Control.Monad +import Data.IORef import Data.Vector.Storable (Vector) +import Numeric.LinearAlgebra (Matrix, (><)) import Numeric.Optimization import IsClose main :: IO () main = hspec $ do - describe "minimize" $ do - context "when given rosenbrock function" $ - it "returns the global optimum" $ do - result <- minimize LBFGS def (WithGrad rosenbrock rosenbrock') [-3,-4] - resultSuccess result `shouldBe` True - assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + describe "minimize CGDescent" $ do + when (isSupportedMethod CGDescent) $ do + context "when given rosenbrock function" $ + it "returns the global optimum" $ do + let prob = WithGrad rosenbrock rosenbrock' + result <- minimize CGDescent def prob [-3,-4] + resultSuccess result `shouldBe` True + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + assertAllClose (def :: Tol Double) (resultValue result) (func prob (resultSolution result)) + case resultGrad result of + Nothing -> return () + Just g -> assertAllClose (def :: Tol Double) g (grad prob (resultSolution result)) + resultHessian result `shouldBe` Nothing + resultHessianInv result `shouldBe` Nothing + let stat = resultStatistics result + totalIters stat `shouldSatisfy` (> 0) + funcEvals stat `shouldSatisfy` (> 0) + gradEvals stat `shouldSatisfy` (> 0) + hessianEvals stat `shouldBe` 0 + context "when given paramsMaxIters" $ + it "stops iterations early" $ do + let prob = WithGrad rosenbrock rosenbrock' + result <- minimize CGDescent def{ paramsMaxIters = Just 2 } prob [1000, 1000] + -- XXX: It seems that CG_DESCENT-C-3.0 report a number number of iterations that is 1 greater than the actual value + totalIters (resultStatistics result) `shouldSatisfy` (\i -> 0 < i && i <= 2+1) + resultSuccess result `shouldBe` False + context "when given a function without gradient" $ do + it "should throw GradUnavailable" $ do + minimize CGDescent def rosenbrock [-3,-4] `shouldThrow` (GradUnavailable ==) + + context "when given a problem with bounds" $ do + it "should throw UnsupportedProblem" $ do + minimize CGDescent def (rosenbrock `WithGrad` rosenbrock' `WithBounds` [(-4,2), (-5,2)]) [-3,-4] + `shouldThrow` (\case { UnsupportedProblem _ -> True; _ -> False }) + + describe "minimize LBFGS" $ do + when (isSupportedMethod LBFGS) $ do + context "when given rosenbrock function" $ + it "returns the global optimum" $ do + let prob = WithGrad rosenbrock rosenbrock' + result <- minimize LBFGS def prob [-3,-4] + resultSuccess result `shouldBe` True + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + assertAllClose (def :: Tol Double) (resultValue result) (func prob (resultSolution result)) + case resultGrad result of + Nothing -> return () + Just g -> assertAllClose (def :: Tol Double) g (grad prob (resultSolution result)) + resultHessian result `shouldBe` Nothing + resultHessianInv result `shouldBe` Nothing + let stat = resultStatistics result + totalIters stat `shouldSatisfy` (>0) + funcEvals stat `shouldSatisfy` (>0) + gradEvals stat `shouldSatisfy` (>0) + hessianEvals stat `shouldBe` 0 + + context "when given rosenbrock function with past" $ + it "returns the global optimum" $ do + let prob = WithGrad rosenbrock rosenbrock' + result <- minimize LBFGS def{ paramsPast = Just 1 } prob [-3,-4] + resultSuccess result `shouldBe` True + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + assertAllClose (def :: Tol Double) (resultValue result) (func prob (resultSolution result)) + + context "when given callback" $ + it "stops iterations early" $ do + let prob = rosenbrock `WithGrad` rosenbrock' + counter <- newIORef (0 :: Int) + let callback x = do + evaluate x + cnt <- readIORef counter + writeIORef counter (cnt + 1) + return (cnt >= 2) + result <- minimize LBFGS def{ paramsCallback = Just callback } prob [1000, 1000] + totalIters (resultStatistics result) `shouldBe` 3 -- ??? + resultSuccess result `shouldBe` False + + context "when given a function without gradient" $ do + it "should throw GradUnavailable" $ do + minimize LBFGS def rosenbrock [-3,-4] `shouldThrow` (GradUnavailable ==) + + context "when given a problem with bounds" $ do + it "should throw UnsupportedProblem" $ do + minimize LBFGS def (rosenbrock `WithGrad` rosenbrock' `WithBounds` [(-4,2), (-5,2)]) [-3,-4] + `shouldThrow` (\case { UnsupportedProblem _ -> True; _ -> False }) + + describe "minimize LBFGSB" $ do + when (isSupportedMethod LBFGSB) $ do + context "when given rosenbrock function" $ + it "returns the global optimum" $ do + let prob = rosenbrock `WithGrad` rosenbrock' + result <- minimize LBFGSB def prob [-3,-4] + resultSuccess result `shouldBe` True + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + assertAllClose (def :: Tol Double) (resultValue result) (func prob (resultSolution result)) + case resultGrad result of + Nothing -> return () + Just g -> assertAllClose (def :: Tol Double) g (grad prob (resultSolution result)) + evaluate $ resultHessian result + evaluate $ resultHessianInv result + let stat = resultStatistics result + totalIters stat `shouldSatisfy` (>0) + funcEvals stat `shouldSatisfy` (>0) + gradEvals stat `shouldSatisfy` (>0) + hessianEvals stat `shouldBe` 0 + + context "when given rosenbrock function with ftol (low accuracy)" $ + it "returns a solution not close enough to global optimum" $ do + let prob = rosenbrock `WithGrad` rosenbrock' + eps = 2.220446049250313e-16 + result <- minimize LBFGSB def{ paramsFTol = Just (1e12 * eps) } prob [-3,-4] + resultSuccess result `shouldBe` True + allClose (def :: Tol Double) (resultSolution result) [1,1] `shouldBe` False + + context "when given rosenbrock function with bounds" $ + it "returns the global optimum" $ do + let prob = rosenbrock `WithGrad` rosenbrock' `WithBounds` [(-4,2), (-5,2)] + result <- minimize LBFGSB def prob [-3,-4] + resultSuccess result `shouldBe` True + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + + context "when given rosenbrock function with bounds (-infinity, +infinity)" $ + it "returns the global optimum" $ do + let prob = rosenbrock `WithGrad` rosenbrock' `WithBounds` boundsUnconstrained 2 + result <- minimize LBFGSB def prob [-3,-4] + resultSuccess result `shouldBe` True + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + + context "when given paramsMaxIters" $ + it "stops iterations early" $ do + let prob = WithGrad rosenbrock rosenbrock' + result <- minimize LBFGSB def{ paramsMaxIters = Just 2 } prob [1000, 1000] + totalIters (resultStatistics result) `shouldSatisfy` (\i -> 0 < i && i <= 2) + resultSuccess result `shouldBe` False + + context "when given a function without gradient" $ do + it "should throw GradUnavailable" $ do + minimize LBFGSB def rosenbrock [-3,-4] `shouldThrow` (GradUnavailable ==) + + describe "minimize Newton" $ do + when (isSupportedMethod Newton) $ do + context "when given rosenbrock function" $ + it "returns the global optimum" $ do + let prob = rosenbrock `WithGrad` rosenbrock' `WithHessian` rosenbrock'' + result <- minimize Newton def prob [-3,-4] + resultSuccess result `shouldBe` True + totalIters (resultStatistics result) `shouldSatisfy` (> 0) + assertAllClose (def :: Tol Double) (resultSolution result) [1,1] + assertAllClose (def :: Tol Double) (resultValue result) (func prob (resultSolution result)) + case resultGrad result of + Nothing -> return () + Just g -> assertAllClose (def :: Tol Double) g (grad prob (resultSolution result)) + case resultHessian result of + Nothing -> return () + Just h -> assertAllClose (def :: Tol Double) h (hessian prob (resultSolution result)) + let stat = resultStatistics result + totalIters stat `shouldSatisfy` (>0) + funcEvals stat `shouldSatisfy` (>0) + gradEvals stat `shouldSatisfy` (>0) + hessianEvals stat `shouldSatisfy` (>0) + + context "when given paramsMaxIters" $ + it "stops iterations early" $ do + let prob = rosenbrock `WithGrad` rosenbrock' `WithHessian` rosenbrock'' + result <- minimize Newton def{ paramsMaxIters = Just 2 } prob [1000, 1000] + totalIters (resultStatistics result) `shouldSatisfy` (\i -> 0 < i && i <= 2) + resultSuccess result `shouldBe` False + assertAllClose (def :: Tol Double) (resultValue result) (func prob (resultSolution result)) + case resultGrad result of + Nothing -> return () + Just g -> assertAllClose (def :: Tol Double) g (grad prob (resultSolution result)) + case resultHessian result of + Nothing -> return () + Just h -> assertAllClose (def :: Tol Double) h (hessian prob (resultSolution result)) + + context "when given callback" $ + it "stops iterations early" $ do + let prob = rosenbrock `WithGrad` rosenbrock' `WithHessian` rosenbrock'' + counter <- newIORef (0 :: Int) + let callback x = do + evaluate x + cnt <- readIORef counter + writeIORef counter (cnt + 1) + return (cnt >= 2) + result <- minimize Newton def{ paramsCallback = Just callback } prob [1000, 1000] + resultSuccess result `shouldBe` False + let stat = resultStatistics result + totalIters stat `shouldBe` 2 + funcEvals stat `shouldSatisfy` (> 0) + gradEvals stat `shouldSatisfy` (> 0) + hessianEvals stat `shouldSatisfy` (> 0) + + context "when given a function without gradient" $ do + it "should throw GradUnavailable" $ do + minimize Newton def (rosenbrock `WithHessian` rosenbrock'') [-3,-4] `shouldThrow` (GradUnavailable ==) + + context "when given a function without Hessian" $ do + it "should throw HessianUnavailable" $ do + minimize Newton def (rosenbrock `WithGrad` rosenbrock') [-3,-4] `shouldThrow` (HessianUnavailable ==) + + context "when given a problem with bounds" $ do + it "should throw UnsupportedProblem" $ do + minimize Newton def (rosenbrock `WithGrad` rosenbrock' `WithHessian` rosenbrock'' `WithBounds` [(-4,2), (-5,2)]) [-3,-4] + `shouldThrow` (\case { UnsupportedProblem _ -> True; _ -> False }) + -- https://en.wikipedia.org/wiki/Rosenbrock_function rosenbrock :: Vector Double -> Double rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x) @@ -24,6 +227,13 @@ rosenbrock' [x,y] = [ 2 * (1 - x) * (-1) + 100 * 2 * (y - sq x) * (-2) * x , 100 * 2 * (y - sq x) + ] + +rosenbrock'' :: Vector Double -> Matrix Double +rosenbrock'' [x,y] = + (2><2) + [ 2 + 100 * 2 * (-2) * ((y - sq x) + (x * (-2) * x)), 100 * 2 * (-2) * x + , 100 * 2 * (-2) * x, 100 * 2 ] sq :: Floating a => a -> a