lagrangian (empty) → 0.1.0.0
raw patch · 6 files changed
+284/−0 lines, 6 filesdep +HUnitdep +addep +basesetup-changed
Dependencies added: HUnit, ad, base, hmatrix, nonlinear-optimization, test-framework, test-framework-hunit, test-framework-quickcheck2, vector
Files
- LICENSE +30/−0
- Setup.hs +2/−0
- lagrangian.cabal +85/−0
- src/Numeric/AD/Lagrangian.hs +27/−0
- src/Numeric/AD/Lagrangian/Internal.hs +105/−0
- tests/Main.hs +35/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2013, Jonathan Fischoff++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Jonathan Fischoff nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ lagrangian.cabal view
@@ -0,0 +1,85 @@+-- 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.1.0.0++-- A short (one-line) description of the package.+synopsis: Solve lagrangian multiplier problems++-- A longer description of the package.+-- description: ++-- 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: ++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.3.*,+ hmatrix == 0.14.*+ + -- Directories containing source files.+ 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.3.*,+ hmatrix == 0.14.*, + test-framework ==0.6.*, + test-framework-hunit ==0.2.*, + test-framework-quickcheck2 ==0.2.*,+ HUnit == 1.2.*++++++++
+ src/Numeric/AD/Lagrangian.hs view
@@ -0,0 +1,27 @@+-- |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 (negate . sum . map (\x -> x * log x), [(sum, 1)]) 3+-- @+-- +-- Gives the answer ([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 (+ solve) where+import Numeric.AD.Lagrangian.Internal (solve, feasible)
+ src/Numeric/AD/Lagrangian/Internal.hs view
@@ -0,0 +1,105 @@+{-# LANGUAGE Rank2Types #-}+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++-- In general I am fighting against the lack of type inference rank two types.+-- Hopefully some of the explicit type signatures can be removed.+++-- The type for the contraints.+-- Given a constraint g(x, y, ...) = c, we would represent it as (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 ++-- | This is the lagrangrain multiplier solver. It is assumed that the +-- objective function and all of the constraints take in the +-- same about of arguments.+solve :: (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 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+ + constraintCount = length (snd params)+ + guess = U.fromList $ replicate (argCount + constraintCount) (1.0 :: Double) ++ result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) + 0.00001 guess (toFunction obj) (toGradient obj)+ Nothing) of+ + (vs, ToleranceStatisfied, _) -> Right (take argCount . S.toList $ vs, + drop argCount . S.toList $ 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)+ + -- 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 ++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)++toFunction :: (forall a. Floating a => [a] -> a) -> Function Simple+toFunction f = VFunction (f . U.toList)++toGradient :: (forall a. Floating a => [a] -> a) -> Gradient Simple+toGradient f = VGradient (U.fromList . grad f . U.toList)+
+ tests/Main.hs view
@@ -0,0 +1,35 @@+module Main where+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++main = defaultMain [+ testGroup "trival test" [+ testCase "noConstraints" noConstraints,+ testCase "entropyTest" entropyTest+ ]+ ] + + +noConstraints = (fst <$> actual) @?= Right expected where+ actual = solve (f, []) 1+ expected = [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 (f, [(\xs -> sum xs, 1.0)]) 3+ expected = [0.33, 0.33, 0.33]+ f :: Floating a => [a] -> a+ f = negate . sum . map (\x -> x * log x)+ + ++ +