lagrangian 0.2.0.2 → 0.3.0.0

raw patch · 4 files changed

+101/−124 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• lagrangian.cabal +9/−21
• src/Numeric/AD/Lagrangian.hs +10/−18
• src/Numeric/AD/Lagrangian/Internal.hs +76/−79
• tests/Main.hs +6/−6

@@ -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

@@ -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)

@@ -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) 

@@ -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)