numeric-optimization 0.1.0.1 → 0.1.1.0
raw patch · 6 files changed
+583/−64 lines, 6 filesdep +l-bfgs-bdep +numeric-limitsdep ~hmatrix
Dependencies added: l-bfgs-b, numeric-limits
Dependency ranges changed: hmatrix
Files
- CHANGELOG.md +11/−1
- README.md +63/−0
- numeric-optimization.cabal +26/−3
- src/Numeric/Optimization.hs +260/−54
- test/IsClose.hs +7/−0
- test/Spec.hs +216/−6
CHANGELOG.md view
@@ -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
README.md view
@@ -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
numeric-optimization.cabal view
@@ -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 <https://github.com/msakai/nonlinear-optimization-ad/tree/master/numeric-optimization#readme>-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
src/Numeric/Optimization.hs view
@@ -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] <https://hackage.haskell.org/package/lbfgs> -- -- * [3] <https://github.com/chokkan/liblbfgs>+ | 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] <https://hackage.haskell.org/package/l-bfgs-b>+ --+ -- * [5] <http://users.iems.northwestern.edu/~nocedal/lbfgsb.html>+ --+ -- @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 -- ------------------------------------------------------------------------
test/IsClose.hs view
@@ -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 -- ------------------------------------------------------------------------
test/Spec.hs view
@@ -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