{-# OPTIONS_GHC -Wall #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.Optimization -- Copyright : (c) Masahiro Sakai 2023 -- License : BSD-style -- -- Maintainer : masahiro.sakai@gmail.com -- Stability : provisional -- Portability : non-portable -- -- 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 -- gradient, but you can use -- [numeric-optimization-ad](https://hackage.haskell.org/package/numeric-optimization-ad) or -- [numeric-optimization-backprop](https://hackage.haskell.org/package/numeric-optimization-backprop) -- to define it using automatic differentiation. -- ----------------------------------------------------------------------------- module Numeric.Optimization ( -- * Main function minimize -- * Problem specification -- -- $problemDefinition , IsProblem (..) , HasGrad (..) , HasHessian (..) , Constraint (..) , boundsUnconstrained , isUnconstainedBounds -- ** Wrapper types , WithGrad (..) , WithHessian (..) , WithBounds (..) , WithConstraints (..) -- * Algorithm selection , Method (..) , isSupportedMethod , Params (..) -- * Result , Result (..) , Statistics (..) , OptimizationException (..) -- * Utilities and Re-export , Default (..) , Optionally (..) , hasOptionalDict ) where import Control.Applicative import Control.Exception import Control.Monad import Control.Monad.Primitive import Control.Monad.ST import Data.Coerce import Data.Constraint (Dict (..)) import Data.Default.Class import Data.Functor.Contravariant import Data.IORef import Data.Maybe import qualified Data.Vector as V import Data.Vector.Storable (Vector) import qualified Data.Vector.Generic as VG 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 data Method = CGDescent -- ^ Conjugate gradient method based on Hager and Zhang [1]. -- -- The implementation is provided by nonlinear-optimization package [3] -- which is a binding library of [2]. -- -- This method requires gradient but does not require hessian. -- -- * [1] Hager, W. W. and Zhang, H. /A new conjugate gradient/ -- /method with guaranteed descent and an efficient line/ -- /search./ Society of Industrial and Applied Mathematics -- Journal on Optimization, 16 (2005), 170-192. -- -- * [2] -- -- * [3] | LBFGS -- ^ Limited memory BFGS (L-BFGS) algorithm [1] -- -- The implementtion is provided by lbfgs package [2] -- which is a binding of liblbfgs [3]. -- -- This method requires gradient but does not require hessian. -- -- * [1] -- -- * [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 -- ^ Naïve implementation of Newton method in Haskell -- -- This method requires both gradient and hessian. deriving (Eq, Ord, Enum, Show, Bounded) -- | 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 -- | Parameters for optimization algorithms -- -- TODO: -- -- * Better way to pass algorithm specific 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@. -- -- 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 def = Params { paramsCallback = Nothing , paramsTol = Nothing , paramsFTol = Nothing , paramsGTol = Nothing , paramsMaxIters = Nothing , paramsPast = Nothing , paramsMaxCorrections = Nothing } instance Contravariant Params where contramap f params = params { paramsCallback = fmap ((. f)) (paramsCallback params) } -- | Optimization result data Result a = Result { resultSuccess :: Bool -- ^ Whether or not the optimizer exited successfully. , resultMessage :: String -- ^ Description of the cause of the termination. , resultSolution :: a -- ^ Solution , resultValue :: Double -- ^ Value of the function at the solution. , resultGrad :: Maybe a -- ^ Gradient at the solution , resultHessian :: Maybe (Matrix Double) -- ^ Hessian at the solution; may be an approximation. , resultHessianInv :: Maybe (Matrix Double) -- ^ Inverse of Hessian at the solution; may be an approximation. , resultStatistics :: Statistics -- ^ Statistics of optimizaion process } deriving (Show) instance Functor Result where fmap f result = result { resultSolution = f (resultSolution result) , resultGrad = fmap f (resultGrad result) } -- | Statistics of optimizaion process data Statistics = Statistics { totalIters :: Int -- ^ Total number of iterations. , funcEvals :: Int -- ^ 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, Eq) instance Exception OptimizationException -- $problemDefinition -- -- Problems are specified by types of 'IsProblem' type class. -- -- 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 require further -- information (e.g. gradients and hessians), -- -- You can equip a problem with other information using wrapper types: -- -- * 'WithBounds' -- -- * 'WithConstraints' -- -- * 'WithGrad' -- -- * 'WithHessian' -- -- If you need further flexibility or efficient implementation, you can -- define instance of 'IsProblem' by yourself. -- | Optimization problems class IsProblem prob where -- | Objective function -- -- It is called @fun@ in @scipy.optimize.minimize@. func :: prob -> Vector Double -> Double -- | Bounds -- bounds :: prob -> Maybe (V.Vector (Double, Double)) bounds _ = Nothing -- | Constraints constraints :: prob -> [Constraint] constraints _ = [] {-# MINIMAL func #-} -- | Optimization problem equipped with gradient information class IsProblem prob => HasGrad prob where -- | Gradient of a function computed by 'func' -- -- It is called @jac@ in @scipy.optimize.minimize@. grad :: prob -> Vector Double -> Vector Double grad prob = snd . grad' prob -- | Pair of 'func' and 'grad' grad' :: prob -> Vector Double -> (Double, Vector Double) grad' prob x = runST $ do gret <- VGM.new (VG.length x) y <- grad'M prob x gret g <- VG.unsafeFreeze gret return (y, g) -- | Similar to 'grad'' but destination passing style is used for gradient vector grad'M :: PrimMonad m => prob -> Vector Double -> VSM.MVector (PrimState m) Double -> m Double grad'M prob x gvec = do let y = func prob x VG.imapM_ (VGM.write gvec) (grad prob x) return y {-# MINIMAL grad | grad' | grad'M #-} -- | Optimization problem equipped with hessian information class IsProblem prob => HasHessian prob where -- | Hessian of a function computed by 'func' -- -- It is called @hess@ in @scipy.optimize.minimize@. hessian :: prob -> Vector Double -> Matrix Double -- | The product of the hessian @H@ of a function @f@ at @x@ with a vector @x@. -- -- It is called @hessp@ in @scipy.optimize.minimize@. -- -- See also . hessianProduct :: prob -> Vector Double -> Vector Double -> Vector Double hessianProduct prob x v = hessian prob x LA.#> v {-# MINIMAL hessian #-} -- | Optional constraint class Optionally c where optionalDict :: Maybe (Dict c) -- | Utility function to define 'Optionally' instances hasOptionalDict :: c => Maybe (Dict c) hasOptionalDict = Just Dict -- | Type of constraint -- -- Currently, no constraints are supported. data Constraint -- | Bounds for unconstrained problems, i.e. (-∞,+∞). boundsUnconstrained :: Int -> V.Vector (Double, Double) boundsUnconstrained n = V.replicate n (-infinity, infinity) -- | Whether all lower bounds are -∞ and all upper bounds are +∞. isUnconstainedBounds :: V.Vector (Double, Double) -> Bool isUnconstainedBounds = V.all p where p (lb, ub) = isInfinite lb && lb < 0 && isInfinite ub && ub > 0 -- | Minimization of scalar function of one or more variables. -- -- This function is intended to provide functionality similar to Python's @scipy.optimize.minimize@. -- -- Example: -- -- > {-# LANGUAGE OverloadedLists #-} -- > -- > import Data.Vector.Storable (Vector) -- > import Numeric.Optimization -- > -- > main :: IO () -- > main = do -- > (x, result, stat) <- minimize LBFGS def (WithGrad rosenbrock rosenbrock') [-3,-4] -- > print (resultSuccess result) -- True -- > print (resultSolution result) -- [0.999999999009131,0.9999999981094296] -- > print (resultValue result) -- 1.8129771632403013e-18 -- > -- > -- https://en.wikipedia.org/wiki/Rosenbrock_function -- > rosenbrock :: Vector Double -> Double -- > rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x) -- > -- > rosenbrock' :: Vector Double -> Vector Double -- > rosenbrock' [x,y] = -- > [ 2 * (1 - x) * (-1) + 100 * 2 * (y - sq x) * (-2) * x -- > , 100 * 2 * (y - sq x) -- > ] -- > -- > sq :: Floating a => a -> a -- > sq x = x ** 2 minimize :: forall prob. (IsProblem prob, Optionally (HasGrad prob), Optionally (HasHessian prob)) => Method -- ^ Numerical optimization algorithm to use -> Params (Vector Double) -- ^ Parameters for optimization algorithms. Use 'def' as a default. -> prob -- ^ Optimization problem to solve -> Vector Double -- ^ Initial value -> IO (Result (Vector Double)) #ifdef WITH_CG_DESCENT minimize CGDescent = case optionalDict @(HasGrad prob) of 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 Just Dict -> case optionalDict @(HasHessian prob) of Nothing -> \_ _ _ -> throwIO HessianUnavailable Just Dict -> minimize_Newton minimize method = \_ _ _ -> throwIO (UnsupportedMethod method) #ifdef WITH_CG_DESCENT minimize_CGDescent :: HasGrad prob => Params (Vector Double) -> prob -> Vector Double -> IO (Result (Vector Double)) minimize_CGDescent _params prob _ | not (isNothing (bounds prob)) = throwIO (UnsupportedProblem "CGDescent does not support bounds") minimize_CGDescent _params prob _ | not (null (constraints prob)) = throwIO (UnsupportedProblem "CGDescent does not support constraints") minimize_CGDescent params prob x0 = do let grad_tol = fromMaybe 1e-6 $ paramsTol params 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 mf mx = do x <- VG.unsafeFreeze mx return $ func prob x mg :: forall m. PrimMonad m => CG.PointMVector m -> CG.GradientMVector m -> m () mg mx mret = do x <- VG.unsafeFreeze mx _ <- grad'M prob x mret return () mc :: forall m. PrimMonad m => CG.PointMVector m -> CG.GradientMVector m -> m Double mc mx mret = do x <- VG.unsafeFreeze mx grad'M prob x mret (x, result, stat) <- CG.optimize cg_params grad_tol x0 (CG.MFunction mf) (CG.MGradient mg) (Just (CG.MCombined mc)) let (success, msg) = case result of CG.ToleranceStatisfied -> (True, "convergence tolerance satisfied") CG.FunctionChange -> (True, "change in func <= feps*|f|") CG.MaxTotalIter -> (False, "total iterations exceeded maxit") CG.NegativeSlope -> (False, "slope always negative in line search") CG.MaxSecantIter -> (False, "number secant iterations exceed nsecant") CG.NotDescent -> (False, "search direction not a descent direction") CG.LineSearchFailsInitial -> (False, "line search fails in initial interval") CG.LineSearchFailsBisection -> (False, "line search fails during bisection") CG.LineSearchFailsUpdate -> (False, "line search fails during interval update") CG.DebugTol -> (False, "debugger is on and the function value increases") CG.FunctionValueNaN -> (False, "function value became nan") CG.StartFunctionValueNaN -> (False, "starting function value is nan") return $ Result { resultSuccess = success , resultMessage = msg , resultSolution = x , resultValue = CG.finalValue stat , resultGrad = Nothing , resultHessian = Nothing , resultHessianInv = Nothing , resultStatistics = Statistics { totalIters = fromIntegral $ CG.totalIters stat , 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") minimize_LBFGS params prob x0 = do evalCounter <- newIORef (0::Int) iterRef <- newIORef (0::Int) let lbfgsParams = LBFGS.LBFGSParameters { LBFGS.lbfgsPast = paramsPast params , LBFGS.lbfgsDelta = fromMaybe 1e-5 $ paramsFTol params <|> paramsTol params , LBFGS.lbfgsLineSearch = LBFGS.DefaultLineSearch , LBFGS.lbfgsL1NormCoefficient = Nothing } instanceData :: () instanceData = () evalFun :: () -> VSM.IOVector CDouble -> VSM.IOVector CDouble -> CInt -> CDouble -> IO CDouble evalFun _inst xvec gvec _n _step = do modifyIORef' evalCounter (+1) #if MIN_VERSION_vector(0,13,0) x <- VG.unsafeFreeze (VSM.unsafeCoerceMVector xvec :: VSM.IOVector Double) y <- grad'M prob x (VSM.unsafeCoerceMVector gvec :: VSM.IOVector Double) #else x <- VG.unsafeFreeze (coerce xvec :: VSM.IOVector Double) y <- grad'M prob x (coerce gvec :: VSM.IOVector Double) #endif return (coerce y) progressFun :: () -> VSM.IOVector CDouble -> VSM.IOVector CDouble -> CDouble -> CDouble -> CDouble -> CDouble -> CInt -> CInt -> CInt -> IO CInt progressFun _inst xvec _gvec _fx _xnorm _gnorm _step _n iter _nev = do writeIORef iterRef $! fromIntegral iter shouldStop <- case paramsCallback params of Nothing -> return False Just callback -> do #if MIN_VERSION_vector(0,13,0) x <- VG.freeze (VSM.unsafeCoerceMVector xvec :: VSM.IOVector Double) #else x <- VG.freeze (coerce xvec :: VSM.IOVector Double) #endif callback x 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_ (success, msg) = case result of LBFGS.Success -> (True, "Success") LBFGS.Stop -> (True, "Stop") LBFGS.AlreadyMinimized -> (True, "The initial variables already minimize the objective function.") LBFGS.UnknownError -> (False, "Unknown error.") LBFGS.LogicError -> (False, "Logic error.") LBFGS.OutOfMemory -> (False, "Insufficient memory.") LBFGS.Canceled -> (False, "The minimization process has been canceled.") LBFGS.InvalidN -> (False, "Invalid number of variables specified.") LBFGS.InvalidNSSE -> (False, "Invalid number of variables (for SSE) specified.") LBFGS.InvalidXSSE -> (False, "The array x must be aligned to 16 (for SSE).") LBFGS.InvalidEpsilon -> (False, "Invalid parameter lbfgs_parameter_t::epsilon specified.") LBFGS.InvalidTestPeriod -> (False, "Invalid parameter lbfgs_parameter_t::past specified.") LBFGS.InvalidDelta -> (False, "Invalid parameter lbfgs_parameter_t::delta specified.") LBFGS.InvalidLineSearch -> (False, "Invalid parameter lbfgs_parameter_t::linesearch specified.") LBFGS.InvalidMinStep -> (False, "Invalid parameter lbfgs_parameter_t::max_step specified.") LBFGS.InvalidMaxStep -> (False, "Invalid parameter lbfgs_parameter_t::max_step specified.") LBFGS.InvalidFtol -> (False, "Invalid parameter lbfgs_parameter_t::ftol specified.") LBFGS.InvalidWolfe -> (False, "Invalid parameter lbfgs_parameter_t::wolfe specified.") LBFGS.InvalidGtol -> (False, "Invalid parameter lbfgs_parameter_t::gtol specified.") LBFGS.InvalidXtol -> (False, "Invalid parameter lbfgs_parameter_t::xtol specified.") LBFGS.InvalidMaxLineSearch -> (False, "Invalid parameter lbfgs_parameter_t::max_linesearch specified.") LBFGS.InvalidOrthantwise -> (False, "Invalid parameter lbfgs_parameter_t::orthantwise_c specified.") LBFGS.InvalidOrthantwiseStart-> (False, "Invalid parameter lbfgs_parameter_t::orthantwise_start specified.") LBFGS.InvalidOrthantwiseEnd -> (False, "Invalid parameter lbfgs_parameter_t::orthantwise_end specified.") LBFGS.OutOfInterval -> (False, "The line-search step went out of the interval of uncertainty.") LBFGS.IncorrectTMinMax -> (False, "A logic error occurred; alternatively, the interval of uncertainty became too small.") LBFGS.RoundingError -> (False, "A rounding error occurred; alternatively, no line-search step satisfies the sufficient decrease and curvature conditions.") LBFGS.MinimumStep -> (False, "The line-search step became smaller than lbfgs_parameter_t::min_step.") LBFGS.MaximumStep -> (False, "The line-search step became larger than lbfgs_parameter_t::max_step.") LBFGS.MaximumLineSearch -> (False, "The line-search routine reaches the maximum number of evaluations.") LBFGS.MaximumIteration -> (False, "The algorithm routine reaches the maximum number of iterations.") LBFGS.WidthTooSmall -> (False, "Relative width of the interval of uncertainty is at most lbfgs_parameter_t::xtol.") 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 $ Result { resultSuccess = success , resultMessage = msg , resultSolution = x , resultValue = func prob x , resultGrad = Nothing , resultHessian = Nothing , resultHessianInv = Nothing , resultStatistics = Statistics { 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 !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 = False , resultMessage = reason , 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 } } 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 0 -- ------------------------------------------------------------------------ instance IsProblem (Vector Double -> Double) where func f = f instance Optionally (HasGrad (Vector Double -> Double)) where optionalDict = Nothing instance Optionally (HasHessian (Vector Double -> Double)) where optionalDict = Nothing -- ------------------------------------------------------------------------ -- | Wrapper type for adding gradient function to a problem data WithGrad prob = WithGrad prob (Vector Double -> Vector Double) instance IsProblem prob => IsProblem (WithGrad prob) where func (WithGrad prob _g) = func prob bounds (WithGrad prob _g) = bounds prob constraints (WithGrad prob _g) = constraints prob instance IsProblem prob => HasGrad (WithGrad prob) where grad (WithGrad _prob g) = g instance HasHessian prob => HasHessian (WithGrad prob) where hessian (WithGrad prob _g) = hessian prob hessianProduct (WithGrad prob _g) = hessianProduct prob instance IsProblem prob => Optionally (HasGrad (WithGrad prob)) where optionalDict = hasOptionalDict instance Optionally (HasHessian prob) => Optionally (HasHessian (WithGrad prob)) where optionalDict = case optionalDict @(HasHessian prob) of Just Dict -> hasOptionalDict Nothing -> Nothing -- ------------------------------------------------------------------------ -- | Wrapper type for adding hessian to a problem data WithHessian prob = WithHessian prob (Vector Double -> Matrix Double) instance IsProblem prob => IsProblem (WithHessian prob) where func (WithHessian prob _hess) = func prob bounds (WithHessian prob _hess) = bounds prob constraints (WithHessian prob _hess) = constraints prob instance HasGrad prob => HasGrad (WithHessian prob) where grad (WithHessian prob _) = grad prob instance IsProblem prob => HasHessian (WithHessian prob) where hessian (WithHessian _prob hess) = hess instance Optionally (HasGrad prob) => Optionally (HasGrad (WithHessian prob)) where optionalDict = case optionalDict @(HasGrad prob) of Just Dict -> hasOptionalDict Nothing -> Nothing instance IsProblem prob => Optionally (HasHessian (WithHessian prob)) where optionalDict = hasOptionalDict -- ------------------------------------------------------------------------ -- | Wrapper type for adding bounds to a problem data WithBounds prob = WithBounds prob (V.Vector (Double, Double)) instance IsProblem prob => IsProblem (WithBounds prob) where func (WithBounds prob _bounds) = func prob bounds (WithBounds _prob bounds) = Just bounds constraints (WithBounds prob _bounds) = constraints prob instance HasGrad prob => HasGrad (WithBounds prob) where grad (WithBounds prob _bounds) = grad prob grad' (WithBounds prob _bounds) = grad' prob grad'M (WithBounds prob _bounds) = grad'M prob instance HasHessian prob => HasHessian (WithBounds prob) where hessian (WithBounds prob _bounds) = hessian prob hessianProduct (WithBounds prob _bounds) = hessianProduct prob instance Optionally (HasGrad prob) => Optionally (HasGrad (WithBounds prob)) where optionalDict = case optionalDict @(HasGrad prob) of Just Dict -> hasOptionalDict Nothing -> Nothing instance Optionally (HasHessian prob) => Optionally (HasHessian (WithBounds prob)) where optionalDict = case optionalDict @(HasHessian prob) of Just Dict -> hasOptionalDict Nothing -> Nothing -- ------------------------------------------------------------------------ -- | Wrapper type for adding constraints to a problem data WithConstraints prob = WithConstraints prob [Constraint] instance IsProblem prob => IsProblem (WithConstraints prob) where func (WithConstraints prob _constraints) = func prob bounds (WithConstraints prob _constraints) = bounds prob constraints (WithConstraints _prob constraints) = constraints instance HasGrad prob => HasGrad (WithConstraints prob) where grad (WithConstraints prob _constraints) = grad prob grad' (WithConstraints prob _constraints) = grad' prob grad'M (WithConstraints prob _constraints) = grad'M prob instance HasHessian prob => HasHessian (WithConstraints prob) where hessian (WithConstraints prob _constraints) = hessian prob hessianProduct (WithConstraints prob _constraints) = hessianProduct prob instance Optionally (HasGrad prob) => Optionally (HasGrad (WithConstraints prob)) where optionalDict = case optionalDict @(HasGrad prob) of Just Dict -> hasOptionalDict Nothing -> Nothing instance Optionally (HasHessian prob) => Optionally (HasHessian (WithConstraints prob)) where optionalDict = case optionalDict @(HasHessian prob) of Just Dict -> hasOptionalDict Nothing -> Nothing -- ------------------------------------------------------------------------