diff --git a/lagrangian.cabal b/lagrangian.cabal --- a/lagrangian.cabal +++ b/lagrangian.cabal @@ -10,37 +10,25 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.2.0.2 +version: 0.3.0.0 -- A short (one-line) description of the package. -synopsis: Solve lagrangian multiplier problems +synopsis: Solve lagrange multiplier problems -- A longer description of the package. description: Numerically solve convex lagrange multiplier problems with conjugate gradient descent. . - Convexity is key, otherwise the descent algorithm can return the wrong answer. - . - Convexity can be tested by assuring that the hessian of the lagrangian is positive - definite over region the function is defined in. - . - I have provided test that the hessian is positive definite at a point, which is something, - but not enough to ensure that the whole function is convex. - . - Be that as it may, if you know what the your lagrangian is convex you can use 'solve' to - find the minimum. - . - For example, find the maximum entropy with the constraint that the probabilities add - up to one. + For example, find the maximum entropy with the constraint that the probabilities sum + to one. . - @ - solve 0.00001 (negate . sum . map (\x -> x * log x), [(sum, 1)]) 3 @ - . - Gives the answer ([0.33, 0.33, 0.33], [-0.09]) + \> solve 0.00001 (negate . sum . map (\x -> x * log x)) [sum \<=\> 1] 3 + ([0.33, 0.33, 0.33], [-0.09]) + @ . - The first elements of the result pair are the arguments for the objective function at the minimum. - The second elements are the lagrange multipliers. + The first elements of the result pair are the arguments for the + objective function at the minimum. The second elements are the lagrange multipliers. . -- URL for the project homepage or repository. homepage: http://github.com/jfischoff/lagrangian diff --git a/src/Numeric/AD/Lagrangian.hs b/src/Numeric/AD/Lagrangian.hs --- a/src/Numeric/AD/Lagrangian.hs +++ b/src/Numeric/AD/Lagrangian.hs @@ -1,29 +1,21 @@ -- |Numerically solve convex lagrange multiplier problems with conjugate gradient descent. --- --- Convexity is key, otherwise the descent algorithm can return the wrong answer. -- --- Convexity can be tested by assuring that the hessian of the lagrangian is positive --- definite over region the function is defined in. --- --- I have provided test that the hessian is positive definite at a point, which is something, --- but not enough to ensure that the whole function is convex. --- --- Be that as it may, if you know what the your lagrangian is convex you can use 'solve' to --- find the minimum. --- -- For example, find the maximum entropy with the constraint that the probabilities add -- up to one. -- --- @ --- solve 0.00001 (negate . sum . map (\x -> x * log x), [(sum, 1)]) 3 --- @ -- --- Gives the answer ([0.33, 0.33, 0.33], [-0.09]) +-- >>> solve 0.00001 (negate . sum . map (\x -> x * log x)) [sum <=> 1] 3 +-- ([0.33, 0.33, 0.33], [-0.09]) -- -- The first elements of the result pair are the arguments for the objective function at the minimum. -- The second elements are the lagrange multipliers. module Numeric.AD.Lagrangian ( + -- *** Helper types + AD2, + (<=>), + Constraint, + -- ** Solver solve, - feasible, - Constraint) where -import Numeric.AD.Lagrangian.Internal (solve, feasible, Constraint) + -- *** Experimental features + feasible) where +import Numeric.AD.Lagrangian.Internal (AD2, (<=>), solve, feasible, Constraint) diff --git a/src/Numeric/AD/Lagrangian/Internal.hs b/src/Numeric/AD/Lagrangian/Internal.hs --- a/src/Numeric/AD/Lagrangian/Internal.hs +++ b/src/Numeric/AD/Lagrangian/Internal.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE Rank2Types, FlexibleContexts #-} module Numeric.AD.Lagrangian.Internal where import Numeric.Optimization.Algorithms.HagerZhang05 import qualified Data.Vector.Unboxed as U @@ -10,97 +10,94 @@ import Numeric.LinearAlgebra.Algorithms import qualified Data.Packed.Vector as V import qualified Data.Packed.Matrix as M - --- In general I am fighting against the lack of type inference rank two types. --- Hopefully some of the explicit type signatures can be removed. - +import Numeric.AD.Internal.Tower --- The type for the contraints. --- Given a constraint g(x, y, ...) = c, we would represent it as (g, c). +infixr 1 <=> +-- | This is just a little bit of sugar for (,) to make constraints look like +-- equals +(<=>) :: ([a] -> a) -> a -> Constraint a +g <=> c = (g,c) +-- | The type for the contraints. +-- Given a constraint @g(x, y, ...) = c@, we would represent it as @(g, c)@. +-- or with sugar @g@ '<=>' @c@ type Constraint a = ([a] -> a, a) --- | This is not a true feasibility test for the function. I am not sure exactly how to --- implement that. This just checks the feasiblility at point. If this ever returns --- false, 'solve' can fail. -feasible :: (forall a. Floating a => ([a] -> a, [Constraint a], [a])) - -> Bool -feasible params = result where - obj :: Floating a => [a] -> a - obj argsAndLams = squaredGrad lang argsAndLams - - lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a) - lang = lagrangian fAndGs (length point) - - fAndGs :: (forall a. Floating a => ([a] -> a, [Constraint a])) - fAndGs = (\(x, y, _) -> (x, y)) params - - point :: Floating a => [a] - point = (\(_, _, x) -> x) params - - h :: [[Double]] - h = hessian obj point - -- I want the hessian as a matrix - hessianMatrix = M.fromLists h - - -- make sure all of the eigenvalues are positive - result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix +type AD2 s r a = AD s (AD r Double) --- | This is the lagrangrain multiplier solver. It is assumed that the +-- | This is the lagrangian multiplier solver. It is assumed that the -- objective function and all of the constraints take in the --- same about of arguments. +-- same amount of arguments. solve :: Double - -> (forall a. Floating a => ([a] -> a, [Constraint a])) -- ^ A pair of the function to minimize and the constraints - -> Int -- ^ The arity of the objective function and the constraints. - -> Either (Result, Statistics) ([Double], [Double]) -- ^ Either an explaination of why the gradient descent failed or a pair of the arguments at the minimum and the lagrange multipliers -solve tolerance params argCount = result where - obj :: Floating a => [a] -> a - obj argsAndLams = squaredGrad lang argsAndLams - - lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a) - lang = lagrangian params argCount + -> (forall s r. (Mode s, Mode r) => [AD2 s r Double] -> AD2 s r Double) + -- ^ The function to minimize + -> (forall s r. (Mode s, Mode r) => [Constraint (AD2 s r Double)] ) + -- ^ The constraints as pairs @g \<=\> c@ which represent equations + -- of the form @g(x, y, ...) = c@ + -> Int + -- ^ The arity of the objective function which should equal the arity of + -- the constraints. + -> Either (Result, Statistics) (S.Vector Double, S.Vector Double) + -- ^ Either an explaination of why the gradient descent failed or a pair + -- containing the arguments at the minimum and the lagrange multipliers +solve tolerance toMin constraints argCount = result where + -- The function to minimize for the langrangian is the squared gradient + obj argsAndLams = + squaredGrad (lagrangian toMin constraints argCount) argsAndLams - constraintCount = length (snd params) + -- The mode does matter but I need to add annotation for the type checker + constraintCount = length (constraints :: [Constraint (AD Tower (AD Tower Double))]) - guess = U.fromList $ replicate (argCount + constraintCount) (1.0 :: Double) - - result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) - tolerance guess (toFunction obj) (toGradient obj) - Nothing) of - - (vs, ToleranceStatisfied, _) -> Right (take argCount . S.toList $ vs, - drop argCount . S.toList $ vs) + -- perhaps this should be exposed + guess = U.replicate (argCount + constraintCount) (1.0 :: Double) + + result = case unsafePerformIO $ + optimize + (defaultParameters { printFinal = False }) + tolerance + guess + (VFunction (lowerFU obj . U.toList)) + (VGradient (U.fromList . grad obj . U.toList)) + Nothing of + (vs, ToleranceStatisfied, _) -> Right (S.take argCount vs, + S.drop argCount vs) (_, x, y) -> Left (x, y) --- Convert a objective function and a list of constraints to a lagrangian -lagrangian :: Floating a - => ([a] -> a, [Constraint a]) - -> Int - -> [a] - -> a -lagrangian (f, constraints) argsLength argsAndLams = result where - -- L(x, y, ..., lam0, lam1, ...) = f(x, y, ...) + - result = f args + (sum $ zipWith (*) lams appliedConstraints) +lagrangian :: Num a + => ([a] -> a) + -> [Constraint a] + -> Int + -> [a] + -> a +lagrangian f constraints argCount argsAndLams = result where + args = take argCount argsAndLams + lams = drop argCount argsAndLams - -- Apply the arguments to the constraint function - -- and subtract to set equal to zero -- (g, c) <=> g(x, y, ...) = c <=> g(x, y, ...) - c = 0 - appliedConstraints = map (\(f, c) -> f args - c) constraints - - -- Split the input by args and lambdas. - -- It is assumed that the args for f and g's come before the - -- lambdas for the constraints - args = take argsLength argsAndLams - lams = drop argsLength argsAndLams - -sumMap f = sum . map f + appliedConstraints = fmap (\(f, c) -> f args - c) constraints + + -- L(x, y, ..., lam0, ...) = f(x, y, ...) + lam0 * (g0 - c0) ... + result = f args + (sum . zipWith (*) lams $ appliedConstraints) -squaredGrad :: Num a - => (forall s. Mode s => [AD s a] -> AD s a) -> [a] -> a -squaredGrad f vs = sumMap (\x -> x*x) (grad f vs) +squaredGrad :: Num a + => (forall s. Mode s => [AD s a] -> AD s a) + -> [a] -> a +squaredGrad f vs = sum . fmap (\x -> x*x) . grad f $ vs -toFunction :: (forall a. Floating a => [a] -> a) -> Function Simple -toFunction f = VFunction (f . U.toList) +-- | WARNING. Experimental. +-- This is not a true feasibility test for the function. I am not sure +-- exactly how to implement that. This just checks the feasiblility at point. +-- If this ever returns false, 'solve' can fail. +feasible :: (forall s r. (Mode s, Mode r) => [AD2 s r Double] -> AD2 s r Double) + -> (forall s r. (Mode s, Mode r) => [Constraint (AD2 s r Double)] ) + -> [Double] + -> Bool +feasible toMin constraints points = result where + obj argsAndLams = + squaredGrad (lagrangian toMin constraints $ length points) argsAndLams + + hessianMatrix = M.fromLists . hessian obj $ points + + -- make sure all of the eigenvalues are positive + result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix -toGradient :: (forall a. Floating a => [a] -> a) -> Gradient Simple -toGradient f = VGradient (U.fromList . grad f . U.toList) diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -5,6 +5,7 @@ import Test.Framework.Providers.QuickCheck2 (testProperty) import Numeric.AD.Lagrangian.Internal import Control.Applicative +import qualified Data.Vector.Storable as S main = defaultMain [ testGroup "trival test" [ @@ -15,17 +16,16 @@ noConstraints = (fst <$> actual) @?= Right expected where - actual = solve 0.00001 (f, []) 1 - expected = [1] + actual = solve 0.00001 f [] 1 + expected = S.fromList [1] f [x] = -(x - 1) ^2 --class Approximate a where -- x =~= y :: a -> a -> Bool - -entropyTest = (sum . map abs $ zipWith (-) actual expected) < 0.02 @?= True where - Right actual = fst <$> solve 0.00001 (f, [(\xs -> sum xs, 1.0)]) 3 - expected = [0.33, 0.33, 0.33] +entropyTest = (S.sum . S.map abs $ S.zipWith (-) actual expected) < 0.02 @?= True where + Right actual = fst <$> solve 0.00001 f [(\xs -> sum xs, 1)] 3 + expected = S.fromList [0.33, 0.33, 0.33] f :: Floating a => [a] -> a f = negate . sum . map (\x -> x * log x)