diff --git a/lagrangian.cabal b/lagrangian.cabal --- a/lagrangian.cabal +++ b/lagrangian.cabal @@ -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 . 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: . 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.* - - - - - - - - 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,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 + ) 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,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 diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -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)