lagrangian 0.5.0.0 → 0.6.0.0

raw patch · 4 files changed

+161/−189 lines, 4 filesdep ~addep ~basedep ~hmatrixPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: ad, base, hmatrix

API changes (from Hackage documentation)

Files

• lagrangian.cabal +21/−67
• src/Numeric/AD/Lagrangian.hs +24/−16
• src/Numeric/AD/Lagrangian/Internal.hs +100/−94
• tests/Main.hs +16/−12

@@ -1,29 +1,14 @@--- Initial lagrangian.cabal generated by cabal init.  For further --- documentation, see http://haskell.org/cabal/users-guide/---- The name of the package. name:                lagrangian---- The package version.  See the Haskell package versioning policy (PVP) --- for standards guiding when and how versions should be incremented.--- http://www.haskell.org/haskellwiki/Package_versioning_policy--- PVP summary:      +-+------- breaking API changes---                   | | +----- non-breaking API additions---                   | | | +--- code changes with no API change-version:             0.5.0.0---- A short (one-line) description of the package.+version:             0.6.0.0 synopsis:            Solve Lagrange multiplier problems---- A longer description of the package.-description:      - Numerically solve convex Lagrange multiplier problems with conjugate gradient descent. +description:+ Numerically solve convex Lagrange multiplier problems with conjugate gradient descent.  .  For some background on the method of Lagrange multipliers checkout the wikipedia page  <http://en.wikipedia.org/wiki/Lagrange_multiplier>  .  Here is an example from the Wikipedia page on Lagrange multipliers- Maximize f(x, y) = x + y, subject to the constraint x^2 + y^2 = 1 + Maximize f(x, y) = x + y, subject to the constraint x^2 + y^2 = 1  .  @    \> maximize 0.00001 (\\[x, y] -> x + y) [(\\[x, y] -> x^2 + y^2) \<=\> 1] 2@@ -33,78 +18,47 @@  For more information look here: <http://en.wikipedia.org/wiki/Lagrange_multiplier#Example_1>  .  For example, to find the maximum entropy with the constraint that the probabilities sum- to one. + to one.  .  @    \> maximize 0.00001 (negate . sum . map (\\x -> x * log x)) [sum \<=\> 1] 3    Right ([0.33, 0.33, 0.33], [-0.09])  @  .- The first elements of the result pair are the arguments for the + The first elements of the result pair are the arguments for the  objective function at the maximum. The second elements are the Lagrange multipliers.  .--- URL for the project homepage or repository. homepage:            http://github.com/jfischoff/lagrangian---- The license under which the package is released. license:             BSD3---- The file containing the license text. license-file:        LICENSE---- The package author(s). author:              Jonathan Fischoff---- An email address to which users can send suggestions, bug reports, and --- patches. maintainer:          jonathangfischoff@gmail.com---- A copyright notice.--- copyright:           -+-- copyright: category:            Math- build-type:          Simple---- Constraint on the version of Cabal needed to build this package. cabal-version:       >=1.8 - library-  -- Modules exported by the library.   exposed-modules: Numeric.AD.Lagrangian-  -  -- Modules included in this library but not exported.-  other-modules: Numeric.AD.Lagrangian.Internal      -  -  -- Other library packages from which modules are imported.-  build-depends:    base ==4.6.*, -                    nonlinear-optimization ==0.3.*, -                    vector ==0.10.*, -                    ad ==3.4.*,-                    hmatrix == 0.14.*-  -  -- Directories containing source files.+  other-modules: Numeric.AD.Lagrangian.Internal+  ghc-options:          -Wall+  build-depends:    base >=4.5 && < 5,+                    nonlinear-optimization ==0.3.*,+                    vector ==0.10.*,+                    ad >= 4 && < 5,+                    hmatrix >= 0.14 && < 0.17   hs-source-dirs:      src  Test-Suite tests   Hs-Source-Dirs: src, tests   type:       exitcode-stdio-1.0   main-is:    Main.hs-  build-depends: base ==4.6.*,-                 nonlinear-optimization ==0.3.*, -                 vector ==0.10.*, -                 ad ==3.4.*,-                 hmatrix == 0.14.*, -                 test-framework ==0.8.*, -                 test-framework-hunit ==0.3.*, +  build-depends: base >=4.5 && < 5,+                 nonlinear-optimization ==0.3.*,+                 vector ==0.10.*,+                 ad >= 4 && <5,+                 hmatrix >= 0.14 && < 0.17,+                 test-framework ==0.8.*,+                 test-framework-hunit ==0.3.*,                  test-framework-quickcheck2 ==0.3.*,                  HUnit == 1.2.*-------- 

@@ -1,23 +1,31 @@--- |Numerically solve convex lagrange multiplier problems with conjugate gradient descent. ---  ---  Here is an example from the Wikipedia page on Lagrange multipliers.---  Maximize f(x, y) = x + y, subject to the constraint x^2 + y^2 = 1 ---  ---  >>> maximize 0.00001 (\[x, y] -> x + y) [(\[x, y] -> x^2 + y^2) <=> 1] 2---  Right ([0.707,0.707], [-0.707])+-- | Numerically solve convex Lagrange-multiplier problems with conjugate+--   gradient descent.  --  ---  The first elements of the result pair are the arguments for the objective function at the minimum. ---  The second elements are the lagrange multipliers.+--   Consider an example from the Wikipedia page on Lagrange multipliers in+--   which we want to maximize the function f(x, y) = x + y, subject to the+--   constraint x^2 + y^2 = 1:+--   +--   >>> maximize (\[x, y] -> x + y) [(\[x, y] -> x^2 + y^2) <=> 1] 0.00001 2+--   Right ([0.707,0.707], [-0.707])+--   +--   The 'Right' indicates success; the first element of the pair is the+--   argument of the objective function at the maximum, and the second element+--   is a list of Lagrange multipliers.+ module Numeric.AD.Lagrangian (-    -- *** Helper types-    AD2,-    FU(..),-    (<=>),+    -- *** Constraint type     Constraint,-    -- ** Solver+    (<=>),+    -- ** Optimizers     maximize,     minimize,     -- *** Experimental features     feasible) where-import Numeric.AD.Lagrangian.Internal (AD2, FU(..), -    (<=>), maximize, minimize, feasible, Constraint)++import Numeric.AD.Lagrangian.Internal+    ( Constraint+    , (<=>)+    , maximize+    , minimize+    , feasible+    )

@@ -1,119 +1,125 @@ {-# LANGUAGE Rank2Types, FlexibleContexts #-}+ module Numeric.AD.Lagrangian.Internal where-import Numeric.Optimization.Algorithms.HagerZhang05+ import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Storable as S-import Numeric.AD-import GHC.IO                   (unsafePerformIO)-import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Numeric.LinearAlgebra.Algorithms import qualified Data.Packed.Vector as V import qualified Data.Packed.Matrix as M-import Numeric.AD.Internal.Tower -infixr 1 <=>--- | Build a 'Constraint' from a function and a constant-(<=>) :: (forall s r. (Mode s, Mode r) => [AD2 s r a] -> AD2 s r a) -> a -> Constraint a-g <=> c = (FU g,c)+import GHC.IO (unsafePerformIO) --- | A constraint of the form @g(x, y, ...) = c@. See '<=>' for constructing a 'Constraint'.-type Constraint a = (FU a, a)+import Numeric.AD+import Numeric.Optimization.Algorithms.HagerZhang05+import Numeric.LinearAlgebra.Algorithms -type AD2 s r a = AD s (AD r a)+-- | An equality constraint of the form @g(x, y, ...) = c@. Use '<=>' to+-- construct a 'Constraint'.+newtype Constraint = Constraint+  {unConstraint :: forall a. (Floating a) => ([a] -> a, a)} --- | A newtype wrapper for working with the rank 2 types constraint functions. -newtype FU a = FU {unFU :: forall s r. (Mode s, Mode r) => [AD2 s r a] -> AD2 s r a}+infixr 1 <=>+-- | Build a 'Constraint' from a function and a constant+(<=>) :: (forall a. (Floating a) => [a] -> a)+      -> (forall b. (Floating b) => b)+      -> Constraint+f <=> c = Constraint (f, c) --- | This is the lagrangian multiplier solver. It is assumed that the ---   objective function and all of the constraints take in the ---   same amount of arguments.-minimize :: Double-      -> (forall s r. (Mode s, Mode r) => [AD2 s r Double] -> AD2 s r Double) -        -- ^ The function to minimize-      -> [Constraint 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 explanation of why the gradient descent failed or a pair -      --   containing the arguments at the minimum and the lagrange multipliers-minimize 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 constraints -    -    -- perhaps this should be exposed-    guess = U.replicate (argCount + constraintCount) (1.0 :: Double) +-- | Numerically minimize the Langrangian. The objective function and each of+-- the constraints must take the same number of arguments.+minimize :: (forall a. (Floating a) => [a] -> a)+         -- ^ The objective function to minimize+         -> [Constraint]+         -- ^ A list of constraints @g \<=\> c@ corresponding to equations of+         -- the form @g(x, y, ...) = c@+         -> Double+         -- ^ Stop iterating when the largest component of the gradient is+         -- smaller than this value+         -> Int+         -- ^ The arity of the objective function, which must equal the arity of+         -- the constraints+         -> Either (Result, Statistics) (V.Vector Double, V.Vector Double)+         -- ^ Either a 'Right' containing the argmin and the Lagrange+         -- multipliers, or a 'Left' containing an explanation of why the+         -- gradient descent failed+minimize f constraints tolerance argCount = result where+    -- At a constrained minimum of `f`, the gradient of the Lagrangian must be+    -- zero. So we square the Lagrangian's gradient (making it non-negative) and+    -- minimize it.+    (sqGradLgn, gradSqGradLgn) = (fst . g, snd . g) where+        g = grad' $ squaredGrad $ lagrangian f constraints argCount     -    result = case unsafePerformIO $ -                    optimize -                        (defaultParameters { printFinal = False }) -                        tolerance -                        guess -                        (VFunction (lowerFU obj . U.toList)) -                        (VGradient (U.fromList . grad obj . U.toList))+    -- Perhaps this should be exposed. ...+    guess = U.replicate (argCount + length constraints) 1++    result = case unsafePerformIO $+                    optimize+                        (defaultParameters {printFinal = False})+                        tolerance+                        guess+                        (VFunction (sqGradLgn . U.toList))+                        (VGradient (U.fromList . gradSqGradLgn . U.toList))                         Nothing of-       (vs, ToleranceStatisfied, _) -> Right (S.take argCount vs, -                                              S.drop argCount vs) +       (vs, ToleranceStatisfied, _) -> Right (S.take argCount vs,+                                              S.drop argCount vs)        (_, x, y) -> Left (x, y)-       --- | Finding the maximum is the same as the minimum with the objective function inverted-maximize :: Double-      -> (forall s r. (Mode s, Mode r) => [AD2 s r Double] -> AD2 s r Double) -        -- ^ The function to maximize-      -> [Constraint 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 explanation of why the gradient descent failed or a pair -      --   containing the arguments at the minimum and the lagrange multipliers-maximize tolerance toMax constraints argCount = -    minimize tolerance (negate1 . toMax) constraints argCount -lagrangian :: (Num a, Mode s, Mode r)-           => (forall s r. (Mode s, Mode r) => [AD2 s r a] -> AD2 s r a) -           -> [Constraint a]+-- | Numerically maximize the Langrangian. The objective function and each of+-- the constraints must take the same number of arguments.+maximize :: (forall a. (Floating a) => [a] -> a)+         -- ^ The objective function to minimize+         -> [Constraint]+         -- ^ A list of constraints @g \<=\> c@ corresponding to equations of+         -- the form @g(x, y, ...) = c@+         -> Double+         -- ^ Stop iterating when the largest component of the gradient is+         -- smaller than this value+         -> Int+         -- ^ The arity of the objective function, which must equal the arity of+         -- the constraints+         -> Either (Result, Statistics) (V.Vector Double, V.Vector Double)+         -- ^ Either a 'Right' containing the argmax and the Lagrange+         -- multipliers, or a 'Left' containing an explanation of why the+         -- gradient ascent failed+maximize f = minimize $ negate . f++lagrangian :: (Floating a)+           => (forall b. (Floating b) => [b] -> b)+           -> [Constraint]            -> Int-           -> [AD2 s r a]  -           -> AD2 s r a+           -> [a]+           -> a lagrangian f constraints argCount argsAndLams = result where     args = take argCount argsAndLams     lams = drop argCount argsAndLams-    +     -- g(x, y, ...) = c <=> g(x, y, ...) - c = 0-    appliedConstraints = fmap (\(FU f, c) -> f args - (auto . auto) c) constraints-    -    -- L(x, y, ..., lam0, ...) = f(x, y, ...) + lam0 * (g0 - c0) ... -    result = f args + (sum . zipWith (*) lams $ appliedConstraints)+    appliedConstraints = fmap (\(Constraint (g, c)) -> g args - c) constraints -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+    -- L(x, y, ..., lam0, ...) = f(x, y, ...) + lam0 * (g0 - c0) ...+    result = (f args) + (sum . zipWith (*) lams $ appliedConstraints) --- | 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 a 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)-         -> [Constraint Double]-         -> [Double]+squaredGrad :: (Floating a)+            => (forall b. (Floating b) => [b] -> b)+            -> [a]+            -> a+squaredGrad f = sum . fmap square . grad f where+    square x = x * x++-- | 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 a+--   point. If this ever returns false, 'solve' can fail.+feasible :: (Floating a, Field a, M.Element a)+         => (forall b. (Floating b) => [b] -> b)+         ->[Constraint]+         -> [a]          -> 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+feasible f constraints points = result where+    sqGradLgn :: (Floating a) => [a] -> a+    sqGradLgn = squaredGrad $ lagrangian f constraints $ length points +    hessianMatrix = M.fromLists . hessian sqGradLgn $ points +    -- make sure all of the eigenvalues are positive+    result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix

@@ -1,12 +1,17 @@ module Main where++import Control.Applicative+ import Test.Framework (defaultMain, testGroup, defaultMainWithArgs) import Test.Framework.Providers.HUnit import Test.HUnit import Test.Framework.Providers.QuickCheck2 (testProperty)-import Numeric.AD.Lagrangian.Internal-import Control.Applicative+ import qualified Data.Vector.Storable as S +import Numeric.AD.Lagrangian.Internal++ main = defaultMain [         testGroup "trival test" [             testCase "noConstraints" noConstraints,@@ -16,22 +21,21 @@           noConstraints = (fst <$> actual) @?= Right expected where-    actual    = minimize 0.00001 f [] 1+    actual    = minimize f [] 0.00001 1     expected  = S.fromList [1]     f [x] = -(x - 1) ^2      --class Approximate a where --    x =~= y :: a -> a -> Bool ---entropyTest = (S.sum . S.map abs $ S.zipWith (-) actual expected) < 0.02 @?= True  where-    Right actual = fst <$> maximize 0.00001 f [sum <=> 1] 3+entropyTest = absDifference < 0.02 @?= True where+    absDifference = (S.sum . S.map abs $ S.zipWith (-) actual expected)+    Right actual = fst <$> maximize entropy [sum <=> 1] 0.00001 3     expected  = S.fromList [0.33, 0.33, 0.33]-    f :: Floating a => [a] -> a-    f = negate . sum . map (\x -> x * log x)     -    +--------------------------------------------------------------------------------+-- Objective functions to test+-------------------------------------------------------------------------------- -    -    +entropy :: (Floating a) => [a] -> a+entropy = negate . sum . fmap (\x -> x * log x)