packages feed

optimization (empty) → 0.1

raw patch · 24 files changed

+914/−0 lines, 24 filesdep +addep +basedep +directorybuild-type:Customsetup-changed

Dependencies added: ad, base, directory, distributive, doctest, filepath, linear, semigroupoids, vector

Files

+ .ghci view
@@ -0,0 +1,1 @@+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
+ .gitignore view
@@ -0,0 +1,13 @@+dist+docs+wiki+TAGS+tags+wip+.DS_Store+.*.swp+.*.swo+*.o+*.hi+*~+*#
+ .travis.yml view
@@ -0,0 +1,25 @@+language: haskell+before_install:+  # Uncomment whenever hackage is down.+  # - mkdir -p ~/.cabal && cp travis/config ~/.cabal/config && cabal update++  # Try installing some of the build-deps with apt-get for speed.+  - travis/cabal-apt-install $mode++install:+  - cabal configure $mode+  - cabal build++script:+  - $script && hlint src --cpp-define HLINT++notifications:+  irc:+    channels:+      - "irc.freenode.org#haskell-lens"+    skip_join: true+    template:+      - "\x0313foo\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}"++env:+  - mode="--enable-tests" script="cabal test"
+ .vim.custom view
@@ -0,0 +1,31 @@+" Add the following to your .vimrc to automatically load this on startup++" if filereadable(".vim.custom")+"     so .vim.custom+" endif++function StripTrailingWhitespace()+  let myline=line(".")+  let mycolumn = col(".")+  silent %s/  *$//+  call cursor(myline, mycolumn)+endfunction++" enable syntax highlighting+syntax on++" search for the tags file anywhere between here and /+set tags=TAGS;/++" highlight tabs and trailing spaces+set listchars=tab:‗‗,trail:‗+set list++" f2 runs hasktags+map <F2> :exec ":!hasktags -x -c --ignore src"<CR><CR>++" strip trailing whitespace before saving+" au BufWritePre *.hs,*.markdown silent! cal StripTrailingWhitespace()++" rebuild hasktags after saving+au BufWritePost *.hs silent! :exec ":!hasktags -x -c --ignore src"
+ CHANGELOG.markdown view
@@ -0,0 +1,3 @@+0.1+---+* Repository initialized
+ HLint.hs view
@@ -0,0 +1,1 @@+import "hint" HLint.Default
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright 2011 Edward Kmett++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. 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.++3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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.
+ README.markdown view
@@ -0,0 +1,27 @@+optimization+===++These are a set of implementations of various numerical optimization+methods in Haskell. Note that these implementations were originally+written as part of a class project; while at one point they worked+no attention has been given to numerical stability or the many other+potential difficulties of implementing robust numerical+methods. That being said, they should serve to succinctly illustrate+a number of optimization techniques from the modern optimization+literature.++Those seeking a high-level overview of some of these methods are+referred to Stephen Wright's excellent+[tutorial](http://videolectures.net/nips2010_wright_oaml/) from NIPS+2010. A deeper introduction can be found in Boyd and Vandenberghe's+*Complex Optimization* which available freely online.+++Contact Information+-------------------++Contributions and bug reports are welcome!++Please feel free to contact me through github or on the #haskell IRC channel on irc.freenode.net.++- Ben Gamari
+ Setup.lhs view
@@ -0,0 +1,44 @@+#!/usr/bin/runhaskell+\begin{code}+{-# OPTIONS_GHC -Wall #-}+module Main (main) where++import Data.List ( nub )+import Data.Version ( showVersion )+import Distribution.Package ( PackageName(PackageName), PackageId, InstalledPackageId, packageVersion, packageName )+import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) )+import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks )+import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose )+import Distribution.Simple.BuildPaths ( autogenModulesDir )+import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), fromFlag )+import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) )+import Distribution.Verbosity ( Verbosity )+import System.FilePath ( (</>) )++main :: IO ()+main = defaultMainWithHooks simpleUserHooks+  { buildHook = \pkg lbi hooks flags -> do+     generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi+     buildHook simpleUserHooks pkg lbi hooks flags+  }++generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO ()+generateBuildModule verbosity pkg lbi = do+  let dir = autogenModulesDir lbi+  createDirectoryIfMissingVerbose verbosity True dir+  withLibLBI pkg lbi $ \_ libcfg -> do+    withTestLBI pkg lbi $ \suite suitecfg -> do+      rewriteFile (dir </> "Build_" ++ testName suite ++ ".hs") $ unlines+        [ "module Build_" ++ testName suite ++ " where"+        , "deps :: [String]"+        , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg))+        ]+  where+    formatdeps = map (formatone . snd)+    formatone p = case packageName p of+      PackageName n -> n ++ "-" ++ showVersion (packageVersion p)++testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)]+testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys++\end{code}
+ optimization.cabal view
@@ -0,0 +1,85 @@+name:          optimization+category:      Math+version:       0.1+license:       BSD3+cabal-version: >= 1.10+license-file:  LICENSE+author:        Ben Gamari+maintainer:    Ben Gamari <bgamari@gmail.com>+stability:     experimental+homepage:      http://github.com/bgamari/optimization+bug-reports:   http://github.com/bgamari/optimization/issues+copyright:     Copyright (C) 2013 Ben Gamari+synopsis:      Numerical optimization+description:+  These are a set of implementations of various numerical optimization+  methods in Haskell. Note that these implementations were originally+  written as part of a class project; while at one point they worked+  no attention has been given to numerical stability or the many other+  potential difficulties of implementing robust numerical+  methods. That being said, they should serve to succinctly illustrate+  a number of optimization techniques from the modern optimization+  literature.+  .+  Those seeking a high-level overview of some of these methods are+  referred to Stephen Wright's excellent tutorial from NIPS 2010+  <http://videolectures.net/nips2010_wright_oaml/>. A deeper+  introduction can be found in Boyd and Vandenberghe's *Complex+  Optimization* which available freely online.++build-type:    Custom++extra-source-files:+  .ghci+  .gitignore+  .travis.yml+  .vim.custom+  CHANGELOG.markdown+  HLint.hs+  README.markdown+  travis/cabal-apt-install+  travis/config++source-repository head+  type: git+  location: git://github.com/bgamari/optimization.git++library+  hs-source-dirs: src+  default-language: Haskell2010+  ghc-options: -Wall -fno-warn-type-defaults+  build-depends:+    base                >= 4.4          && < 5,+    vector              >= 0.10         && < 1.0,+    ad                  >= 3.4          && < 4.0,+    linear              >= 1.0          && < 2.0,+    semigroupoids       >= 3.0          && < 4.0,+    distributive        >= 0.3          && < 0.4++  exposed-modules:+    Optimization.LineSearch+    Optimization.LineSearch.ConjugateGradient+    Optimization.LineSearch.BarzilaiBorwein+    Optimization.LineSearch.SteepestDescent+    Optimization.LineSearch.MirrorDescent+    Optimization.LineSearch.BFGS+    Optimization.TrustRegion.Nesterov1983+    Optimization.TrustRegion.Fista+    Optimization.TrustRegion.Newton+    Optimization.Constrained.Penalty+    Optimization.Constrained.ProjectedSubgradient+++test-suite doctests+  type:    exitcode-stdio-1.0+  main-is: doctests.hs+  default-language: Haskell2010+  build-depends:+    base,+    directory >= 1.0,+    doctest >= 0.9.1,+    filepath+  ghc-options: -Wall -threaded+  if impl(ghc<7.6.1)+    ghc-options: -Werror+  hs-source-dirs: tests
+ src/Optimization/Constrained/Penalty.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveTraversable, DeriveGeneric,+    FlexibleInstances, FlexibleContexts, TypeFamilies,+    KindSignatures, DataKinds, TypeOperators, RankNTypes, ExistentialQuantification #-}++module Optimization.Constrained.Penalty+  ( -- * Building the problem+    Opt+  , FU(..)+  , optimize+  , constrainEQ+  , constrainLT+  , constrainGT+    -- * Optimizing the problem+  , minimize+  , maximize+    -- * Finding the Lagrangian+  , lagrangian+  ) where++import           Numeric.AD.Types++import qualified Data.Vector as V++newtype FU f a = FU { runFU :: forall s. Mode s => f (AD s a) -> AD s a }++type V = V.Vector++-- | @Opt d f gs hs@ is a Lagrangian optimization problem with objective @f@+-- equality (@g(x) == 0@) constraints @gs@ and less-than (@h(x) < 0@)+-- constraints @hs@+data Opt f a = Opt (FU f a) (V (FU f a)) (V (FU f a))++optimize :: (forall s. Mode s => f (AD s a) -> AD s a) -> Opt f a+optimize f = Opt (FU f) V.empty V.empty++augment :: a -> V a -> V a+augment a xs = V.cons a xs++constrainEQ :: (forall s. Mode s => f (AD s a) -> AD s a)+            -> Opt f a -> Opt f a+constrainEQ g (Opt f gs hs) = Opt f (augment (FU g) gs) hs++constrainLT :: (forall s. Mode s => f (AD s a) -> AD s a)+            -> Opt f a -> Opt f a+constrainLT h (Opt f gs hs) = Opt f gs (augment (FU h) hs)++constrainGT :: (Num a) => (forall s. Mode s => f (AD s a) -> AD s a)+            -> Opt f a -> Opt f a+constrainGT h (Opt f gs hs) = Opt f gs (augment (FU $ negate . h) hs)++-- | Minimize the given constrained optimization problem+-- This is a basic penalty method approach+minimize :: (Functor f, Num a, Ord a, g ~ V)+         => (FU f a -> f a -> [f a])   -- ^ Primal minimizer+         -> Opt f a                    -- ^ The optimization problem of interest+         -> a                          -- ^ The penalty increase factor+         -> f a                        -- ^ The primal starting value+         -> g a                        -- ^ The dual starting value+         -> [f a]                      -- ^ Optimizing iterates+minimize minX opt alpha = go+  where go x0 l0 = let l1 = fmap (*alpha) l0+                       x1 = head $ drop 100 $ minX (FU $ \x -> augLagrangian opt x (fmap auto l1)) x0+                   in x1 : go x1 l1++-- | Maximize the given constrained optimization problem+maximize :: (Functor f, Num a, Ord a, g ~ V)+         => (FU f a -> f a -> [f a])   -- ^ Primal minimizer+         -> Opt f a                    -- ^ The optimization problem of interest+         -> a                          -- ^ The penalty increase factor+         -> f a                        -- ^ The primal starting value+         -> g a                        -- ^ The dual starting value+         -> [f a]                      -- ^ Optimizing iterates+maximize minX (Opt (FU f) gs hs) alpha =+    minimize minX (Opt (FU $ negate . f) gs hs) alpha++-- | The Lagrangian for the given constrained optimization+lagrangian :: (Num a) => Opt f a+           -> (forall s. Mode s => f (AD s a) -> V (AD s a) -> AD s a)+lagrangian (Opt (FU f) gs hs) x l =+    f x - V.sum (V.zipWith (\lamb (FU g)->lamb * g x) l gs)++-- | The augmented Lagrangian for the given constrained optimization+augLagrangian :: (Num a, Ord a) => Opt f a+           -> (forall s. Mode s => f (AD s a) -> V (AD s a) -> AD s a)+augLagrangian (Opt (FU f) gs hs) x l =+    f x + V.sum (V.zipWith (*) l $ V.concat [gs', hs'])+  where gs' = V.map (\(FU g) -> (g x)^2) gs+        hs' = V.map (\(FU h) -> (max 0 $ h x)^2) hs
+ src/Optimization/Constrained/ProjectedSubgradient.hs view
@@ -0,0 +1,114 @@+module Optimization.Constrained.ProjectedSubgradient+    ( -- * Projected subgradient method+      projSubgrad+    , linearProjSubgrad+      -- * Step schedules+    , StepSched+    , optimalStepSched+    , constStepSched+    , sqrtKStepSched+    , invKStepSched+      -- * Linear constraints+    , Constraint(..)+    , linearProjection+    ) where++import Linear+import Data.Traversable+import Data.Function (on)+import Data.List (maximumBy)++-- | A step size schedule+-- A list of functions (one per step) which, given a function's+-- gradient and value, will provide a size for the next step+type StepSched f a = [f a -> a -> a]++-- | @projSubgrad stepSizes proj a b x0@ minimizes the objective @A+-- x - b@ potentially projecting iterates into a feasible space with+-- @proj@ with the given step-size schedule+projSubgrad :: (Additive f, Traversable f, Metric f, Ord a, Fractional a)+            => StepSched f a  -- ^ A step size schedule+            -> (f a -> f a)   -- ^ Function projecting into the feasible space+            -> (f a -> f a)   -- ^ Gradient of objective function+            -> (f a -> a)     -- ^ The objective function+            -> f a            -- ^ Initial solution+            -> [f a]+projSubgrad stepSizes proj df f = go stepSizes+  where go (alpha:rest) x0 =+            let p = negated $ df x0+                step = alpha (df x0) (f x0)+                x1 = proj $ x0 ^+^ step *^ p+            in x1 : go rest x1+        go [] _ = []++-- | @linearProjSubgrad stepSizes proj a b x0@ minimizes the objective @A+-- x - b@ potentially projecting iterates into a feasible space with+-- @proj@ with the given step-size schedule+linearProjSubgrad :: (Additive f, Traversable f, Metric f, Ord a, Fractional a)+                  => StepSched f a  -- ^ A step size schedule+                  -> (f a -> f a)   -- ^ Function projecting into the feasible space+                  -> f a            -- ^ Coefficient vector @A@ of objective+                  -> a              -- ^ Constant @b@ of objective+                  -> f a            -- ^ Initial solution+                  -> [f a]+linearProjSubgrad stepSizes proj a b = go stepSizes+  where go (alpha:rest) x0 =+            let p = negated $ df x0+                step = alpha a (f x0)+                x1 = proj $ x0 ^+^ step *^ p+            in x1 : go rest x1+        go [] _ = []+        df _ = a+        f x = a `dot` x - b++-- | The optimal step size schedule when the optimal value of the+-- objective is known+optimalStepSched :: (Fractional a, Metric f)+                 => a    -- ^ The optimal value of the objective+                 -> StepSched f a+optimalStepSched fStar =+    repeat $ \gk fk->(fk - fStar) / quadrance gk++-- | Constant step size schedule+constStepSched :: a    -- ^ The step size+               -> StepSched f a+constStepSched gamma =+    repeat $ \_ _ -> gamma++-- | A non-summable step size schedule+sqrtKStepSched :: Floating a+               => a       -- ^ The size of the first step+               -> StepSched f a+sqrtKStepSched gamma =+    map (\k _ _ -> gamma / sqrt (fromIntegral k)) [0..]++-- | A square-summable step size schedule+invKStepSched :: Fractional a+              => a        -- ^ The size of the first step+              -> StepSched f a+invKStepSched gamma =+    map (\k _ _ -> gamma / fromIntegral k) [0..]++-- | A linear constraint. For instance, @Constr LT 2 (V2 1 3)@ defines+-- the constraint @x_1 + 3 x_2 <= 2@+data Constraint f a = Constr Ordering a (f a)+                    deriving (Show)++-- | Project onto a the space of feasible solutions defined by a set+-- of linear constraints+linearProjection :: (Fractional a, Ord a, RealFloat a, Metric f)+                 => [Constraint f a] -- ^ A set of linear constraints+                 -> f a -> f a+linearProjection constraints x =+    case unmet of+      []   -> x+      _    -> linearProjection constraints $ fixConstraint x+              $ maximumBy (flip compare `on` (`ap` x)) unmet+  where unmet = filter (not . met x) constraints+        ap (Constr _ b a) c = a `dot` c - b+        met c (Constr t a constr) = let y = constr `dot` c - a+                                    in case t of+                                       EQ -> abs y < 1e-4+                                       GT -> y >= 0 || abs y < 1e-4+                                       LT -> y <= 0 || abs y < 1e-4+        fixConstraint c (Constr _ b a) = c ^-^ (a `dot` c - b) *^ a ^/ quadrance a
+ src/Optimization/LineSearch.hs view
@@ -0,0 +1,100 @@+-- |+-- Module      : Optimization.LineSearch+-- Copyright   : (c) 2012-2013 Ben Gamari+-- License     : BSD-style (see the file LICENSE)+-- Maintainer  : Ben Gamari <bgamari@gmail.com>+-- Stability   : provisional+-- Portability : portable+--+-- Line search algorithms are a class of iterative optimization+-- methods. These methods are distinguished by the characteristic of,+-- starting from a point @x0@, choosing a direction @d@ (by some method)+-- to advance and then finding an optimal distance @a@ (known as the+-- step-size) to advance in this direction.+--+-- Here we provide several methods for determining this optimal+-- distance. These can be used with any of line-search optimization+-- algorithms found in this namespace.++module Optimization.LineSearch+    ( -- * Line search methods+      LineSearch+    , backtrackingSearch+    , armijoSearch+    , wolfeSearch+    , newtonSearch+    , secantSearch+    , constantSearch+    ) where++import Prelude hiding (pred)+import Linear++-- | A 'LineSearch' method 'search df p x' determines a step size+-- in direction 'p' from point 'x' for function 'f' with gradient 'df'+type LineSearch f a = (f a -> f a) -> f a -> f a -> a++-- | Armijo condition+--+-- The Armijo condition captures the intuition that step should+-- move far enough from its starting point to change the function enough,+-- as predicted by its gradient. This often finds its place as a criterion+-- for line-search+armijo :: (Num a, Additive f, Ord a, Metric f)+       => a -> (f a -> a) -> (f a -> f a) -> f a -> f a -> a -> Bool+armijo c1 f df x p a =+    f (x ^+^ a *^ p) <= f x + c1 * a * (df x `dot` p)++-- | Curvature condition+curvature :: (Num a, Ord a, Additive f, Metric f)+          => a -> (f a -> f a) -> f a -> f a -> a -> Bool+curvature c2 df x p a =+    df (x ^+^ a *^ p) `dot` p >= c2 * (df x `dot` p)++-- | Backtracking line search algorithm+--+-- @backtrackingSearch gamma alpha pred@ starts with the given step+-- size @alpha@ and reduces it by a factor of @gamma@ until the given+-- condition is satisfied.+backtrackingSearch :: (Num a, Ord a, Metric f)+                   => a -> a -> (a -> Bool) -> LineSearch f a+backtrackingSearch gamma alpha pred _ _ _ =+    head $ dropWhile (not . pred) $ nonzero $ iterate (*gamma) alpha+  where nonzero (x:xs) | not $ x > 0 = error "Backtracking search failed: alpha=0" -- FIXME+                       | otherwise   = x : nonzero xs+        nonzero [] = error "Backtracking search failed: no more iterates"++-- | Armijo backtracking line search algorithm+--+-- @armijoSearch gamma alpha c1@ starts with the given step size @alpha@+-- and reduces it by a factor of @gamma@ until the Armijo condition+-- is satisfied.+armijoSearch :: (Num a, Ord a, Metric f)+             => a -> a -> a -> (f a -> a) -> LineSearch f a+armijoSearch gamma alpha c1 f df p x =+    backtrackingSearch gamma alpha (armijo c1 f df x p) df p x++-- | Wolfe backtracking line search algorithm+--+-- @wolfeSearch gamma alpha c1@ starts with the given step size @alpha@+-- and reduces it by a factor of @gamma@ until both the Armijo and+-- curvature conditions is satisfied.+wolfeSearch :: (Num a, Ord a, Metric f)+             => a -> a -> a -> a -> (f a -> a) -> LineSearch f a+wolfeSearch gamma alpha c1 c2 f df p x =+    backtrackingSearch gamma alpha wolfe df p x+  where wolfe a = armijo c1 f df p x a && curvature c2 df x p a++-- | Line search by Newton's method+newtonSearch :: (Num a) => LineSearch f a+newtonSearch = undefined++-- | Line search by secant method with given tolerance+secantSearch :: (Num a, Fractional a) => a -> LineSearch f a+secantSearch = undefined++-- | Constant line search+--+-- @constantSearch c@ always chooses a step-size @c@.+constantSearch :: a -> LineSearch f a+constantSearch c _ _ _ = c
+ src/Optimization/LineSearch/BFGS.hs view
@@ -0,0 +1,31 @@+{-# LANGUAGE ScopedTypeVariables #-}++module Optimization.LineSearch.BFGS (bfgs) where++import Linear+import Optimization.LineSearch+import Control.Applicative+import Data.Traversable+import Data.Distributive+import Data.Foldable++-- | Broyden–Fletcher–Goldfarb–Shanno (BFGS) method+-- @bfgs search df b0 x0@ where @b0@ is the inverse Hessian (the+-- identity is often a good initial value).+bfgs :: ( Metric f, Additive f, Distributive f, Foldable f, Traversable f, Applicative f+        , Fractional a, Epsilon a)+     => LineSearch f a -> (f a -> f a) -> f (f a) -> f a -> [f a]+bfgs search df = go+    where go b0 x0 = let p1 = negated $ b0 !* df x0+                         alpha = search df p1 x0+                         s = alpha *^ p1+                         x1 = x0 ^+^ s+                         y = df x1 ^-^ df x0+                         -- Sherman-Morrison update of inverse Hessian+                         sy = s `dot` y+                         rho = if nearZero sy then 1000 else 1 / sy+                         i = kronecker (pure 1)+                         u = i !-! rho *!! outer y s+                         v = i !-! rho *!! outer s y+                         b1 = u !*! b0 !*! v !+! rho *!! outer s s+                     in x1 : go b1 x1
+ src/Optimization/LineSearch/BarzilaiBorwein.hs view
@@ -0,0 +1,17 @@+module Optimization.LineSearch.BarzilaiBorwein+    ( barzilaiBorwein+    ) where++import Linear++-- | Barzilai-Borwein 1988 is a non-monotonic optimization method+barzilaiBorwein :: (Additive f, Metric f, Functor f, Fractional a, Epsilon a)+                => (f a -> f a) -> f a -> f a -> [f a]+barzilaiBorwein df = go+  where go x0 x1 = let s = x1 ^-^ x0+                       z = df x1 ^-^ df x0+                       alpha = (s `dot` z) / (z `dot` z)+                       x2 = x1 ^-^ alpha *^ df x1+                   in if nearZero (z `dot` z)+                        then [x2]+                        else x2 : go x1 x2
+ src/Optimization/LineSearch/ConjugateGradient.hs view
@@ -0,0 +1,54 @@+module Optimization.LineSearch.ConjugateGradient+    ( -- * Conjugate gradient methods+      conjGrad+      -- * General line search+    , module Optimization.LineSearch+      -- * Beta expressions+    , Beta+    , fletcherReeves+    , polakRibiere+    , hestenesStiefel+    ) where++import Optimization.LineSearch+import Linear++-- | A beta expression 'beta df0 df1 p' is an expression for the+-- conjugate direction contribution given the derivative 'df0' and+-- direction 'p' for iteration 'k', 'df1' for iteration 'k+1'+type Beta f a = f a -> f a -> f a -> a++-- | Conjugate gradient method with given beta and line search method+--+-- The conjugate gradient method avoids the trouble encountered by the+-- steepest descent method on poorly conditioned problems (e.g. those with+-- a wide range of eigenvalues). It does this by choosing directions which+-- satisfy a condition of @A@ orthogonality, ensuring that steps in the+-- "unstretched" search space are orthogonal.+-- TODO: clarify explanation+{-# INLINEABLE conjGrad #-}+conjGrad :: (Num a, RealFloat a, Additive f, Metric f)+         => LineSearch f a -> Beta f a+         -> (f a -> f a) -> f a -> [f a]+conjGrad search beta df x0 = go (negated $ df x0) x0+  where go p x = let a = search df p x+                     x' = x ^+^ a *^ p+                     b = beta (df x) (df x') p+                     p' = negated (df x') ^+^ b *^ p+                 in x' : go p' x'++-- | Fletcher-Reeves expression for beta+{-# INLINEABLE fletcherReeves #-}+fletcherReeves :: (Num a, RealFloat a, Metric f) => Beta f a+fletcherReeves df0 df1 _ = norm df1 / norm df0++-- | Polak-Ribiere expression for beta+{-# INLINEABLE polakRibiere #-}+polakRibiere :: (Num a, RealFloat a, Metric f) => Beta f a+polakRibiere df0 df1 _ = df1 `dot` (df1 ^-^ df0) / norm df0++-- | Hestenes-Stiefel expression for beta+{-# INLINEABLE hestenesStiefel #-}+hestenesStiefel :: (Num a, RealFloat a, Metric f) => Beta f a+hestenesStiefel df0 df1 p0 =+    - (df1 `dot` (df1 ^-^ df0)) / (p0 `dot` (df1 ^-^ df0))
+ src/Optimization/LineSearch/MirrorDescent.hs view
@@ -0,0 +1,23 @@+module Optimization.LineSearch.MirrorDescent+    ( mirrorDescent ) where++import Optimization.LineSearch+import Linear++-- | Mirror descent method.+--+-- Originally described by Nemirovsky and Yudin and later elucidated+-- by Beck and Teboulle, the mirror descent method is a generalization of+-- the projected subgradient method for convex optimization.+-- Mirror descent requires the gradient of a strongly+-- convex function @psi@ (and its dual) as well as a way to get a+-- subgradient for each point of the objective function @f@.+mirrorDescent :: (Num a, Additive f)+              => LineSearch f a -> (f a -> f a) -> (f a -> f a)+              -> (f a -> f a) -> f a -> [f a]+mirrorDescent search dPsi dPsiStar df = go+  where go y0 = let x0 = dPsiStar y0+                    t0 = search df (df x0) x0+                    y1 = dPsi x0 ^-^ t0 *^ df x0+                    x1 = dPsiStar y1+                in x1 : go y1
+ src/Optimization/LineSearch/SteepestDescent.hs view
@@ -0,0 +1,22 @@+module Optimization.LineSearch.SteepestDescent+    ( -- * Steepest descent method+      steepestDescent+    ) where++import Optimization.LineSearch+import Linear++-- | Steepest descent method+--+-- @steepestDescent search f df x0@ optimizes a function @f@ with gradient @df@+-- with step size schedule @search@ starting from initial point @x0@+--+-- The steepest descent method chooses the negative gradient of the function+-- as its step direction.+{-# INLINEABLE steepestDescent #-}+steepestDescent :: (Num a, Ord a, Additive f, Metric f)+                => LineSearch f a -> (f a -> f a) -> f a -> [f a]+steepestDescent search df x0 = iterate go x0+  where go x = let p = negated (df x)+                   a = search df p x+               in x ^+^ a *^ p
+ src/Optimization/TrustRegion/Fista.hs view
@@ -0,0 +1,17 @@+module Optimization.TrustRegion.Fista+    ( -- * Fast Iterative Shrinkage-Thresholding Algorithm+      fista+    ) where++import Linear++-- | Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) with+-- constant stepsize+{-# INLINEABLE fista #-}+fista :: (Additive f, Fractional a, Floating a)+      => a -> (f a -> f a) -> f a -> [f a]+fista l df x0' = go x0' x0' 1+  where go x0 y1 t1 = let x1 = y1 ^-^ df y1 ^/ l+                          t2 = (1 + sqrt (1 + 4 * t1^2)) / 2+                          y2 = x1 ^+^ (t1-1) / t2 *^ (x1 ^-^ x0)+                      in x1 : go x1 y2 t2
+ src/Optimization/TrustRegion/Nesterov1983.hs view
@@ -0,0 +1,35 @@+module Optimization.TrustRegion.Nesterov1983+    ( -- * Nesterov's Optimal Gradient method+      optimalGradient+    ) where++import Linear++-- | Nesterov 1983+-- @optimalGradient kappa l df alpha0 x0@ is Nesterov's optimal+-- gradient method, first described in 1983. This method requires+-- knowledge of the Lipschitz constant @l@ of the gradient, the condition+-- number @kappa@, as well as an initial step size @alpha0@ in @(0,1)@.+{-# INLINEABLE optimalGradient #-}+optimalGradient :: (Additive f, Functor f, Ord a, Floating a, Epsilon a)+                => a -> a -> (f a -> f a) -> a -> f a -> [f a]+optimalGradient kappa l df a0' x0' = go a0' x0' x0'+  where go a0 x0 y0 = let x1 = y0 ^-^ df y0 ^/ l+                          alphas = quadratic 1 (a0^2 - 1/kappa) (-a0^2)+                          a1 = case filter (\x->x >= 0 && x <= 1) alphas of+                                 a:_  -> a+                                 []   -> error "No solution for alpha_{k+1}"+                          b1 = a0 * (1 - a0) / (a0^2 + a1)+                          y1 = x1 ^+^ b1 *^ (x1 ^-^ x0)+                      in x1 : go a0 x1 y1++-- | 'quadratic a b c' is the real solutions to a quadratic equation+-- 'a x^2 + b x + c == 0'+quadratic :: (Ord a, Floating a, Epsilon a)+          => a -> a -> a -> [a]+quadratic a b c+    | discr < 0      = []+    | nearZero discr = [-b / 2 / a]+    | otherwise      = [ (-b + sqrt discr) / 2 / a+                       , (-b - sqrt discr) / 2 / a ]+  where discr = b^2 - 4*a*c
+ src/Optimization/TrustRegion/Newton.hs view
@@ -0,0 +1,37 @@+module Optimization.TrustRegion.Newton+    ( -- * Newton's method+      newton+      -- * Matrix inversion methods+    , bicInv+    , bicInv'+    ) where++import Control.Applicative+import Data.Distributive (Distributive)+import Data.Functor.Bind (Apply)+import Data.Foldable (Foldable)+import Linear++-- | Newton's method+{-# INLINEABLE newton #-}+newton :: (Num a, Ord a, Additive f, Metric f, Foldable f)+       => (f a -> f a) -> (f a -> f (f a)) -> f a -> [f a]+newton df ddfInv x0 = iterate go x0+  where go x = x ^-^ ddfInv x !* df x++-- | Inverse by iterative method of Ben-Israel and Cohen+-- with given starting condition+bicInv' :: (Functor m, Distributive m, Additive m,+            Applicative m, Apply m, Foldable m, Conjugate a)+        => m (m a) -> m (m a) -> [m (m a)]+bicInv' a0 a = iterate go a0+  where go ak = 2 *!! ak !-! ak !*! a !*! ak++-- | Inverse by iterative method of Ben-Israel and Cohen+-- starting from 'alpha A^T'. Alpha should be set such that+-- 0 < alpha < 2/sigma^2 where sigma denotes the largest singular+-- value of A+bicInv :: (Functor m, Distributive m, Additive m,+           Applicative m, Apply m, Foldable m, Conjugate a)+       => a -> m (m a) -> [m (m a)]+bicInv alpha a = bicInv' (alpha *!! adjoint a) a
+ tests/doctests.hsc view
@@ -0,0 +1,73 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE ForeignFunctionInterface #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Main (doctests)+-- Copyright   :  (C) 2012-13 Edward Kmett+-- License     :  BSD-style (see the file LICENSE)+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>+-- Stability   :  provisional+-- Portability :  portable+--+-- This module provides doctests for a project based on the actual versions+-- of the packages it was built with. It requires a corresponding Setup.lhs+-- to be added to the project+-----------------------------------------------------------------------------+module Main where++import Build_doctests (deps)+import Control.Applicative+import Control.Monad+import Data.List+import System.Directory+import System.FilePath+import Test.DocTest++##ifdef mingw32_HOST_ARCH+##ifdef i386_HOST_ARCH+##define USE_CP+import Control.Applicative+import Control.Exception+import Foreign.C.Types+foreign import stdcall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool+foreign import stdcall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt+##elif defined(x86_64_HOST_ARCH)+##define USE_CP+import Control.Applicative+import Control.Exception+import Foreign.C.Types+foreign import ccall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool+foreign import ccall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt+##endif+##endif++-- | Run in a modified codepage where we can print UTF-8 values on Windows.+withUnicode :: IO a -> IO a+##ifdef USE_CP+withUnicode m = do+  cp <- c_GetConsoleCP+  (c_SetConsoleCP 65001 >> m) `finally` c_SetConsoleCP cp+##else+withUnicode m = m+##endif++main :: IO ()+main = withUnicode $ getSources >>= \sources -> doctest $+    "-isrc"+  : "-idist/build/autogen"+  : "-optP-include"+  : "-optPdist/build/autogen/cabal_macros.h"+  : "-hide-all-packages"+  : map ("-package="++) deps ++ sources++getSources :: IO [FilePath]+getSources = filter (isSuffixOf ".hs") <$> go "src"+  where+    go dir = do+      (dirs, files) <- getFilesAndDirectories dir+      (files ++) . concat <$> mapM go dirs++getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])+getFilesAndDirectories dir = do+  c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir+  (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c
+ travis/cabal-apt-install view
@@ -0,0 +1,27 @@+#! /bin/bash+set -eu++APT="sudo apt-get -q -y"+CABAL_INSTALL_DEPS="cabal install --only-dependencies --force-reinstall"++$APT update+$APT install dctrl-tools++# Find potential system packages to satisfy cabal dependencies+deps()+{+	local M='^\([^ ]\+\)-[0-9.]\+ (.*$'+	local G=' -o ( -FPackage -X libghc-\L\1\E-dev )'+	local E="$($CABAL_INSTALL_DEPS "$@" --dry-run -v 2> /dev/null \+		| sed -ne "s/$M/$G/p" | sort -u)"+	grep-aptavail -n -sPackage \( -FNone -X None \) $E | sort -u+}++$APT install $(deps "$@") libghc-quickcheck2-dev # QuickCheck is special+$CABAL_INSTALL_DEPS "$@" # Install the rest via Hackage++if ! $APT install hlint ; then+	$APT install $(deps hlint)+	cabal install hlint+fi+
+ travis/config view
@@ -0,0 +1,16 @@+-- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix+--+-- This is particularly useful for travis-ci to get it to stop complaining+-- about a broken build when everything is still correct on our end.+--+-- This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead+--+-- To enable this, uncomment the before_script in .travis.yml++remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive+remote-repo-cache: ~/.cabal/packages+world-file: ~/.cabal/world+build-summary: ~/.cabal/logs/build.log+remote-build-reporting: anonymous+install-dirs user+install-dirs global