ad 4.2.4 → 4.3

raw patch · 12 files changed

+247/−34 lines, 12 filesdep +HerbiePluginPVP ok

version bump matches the API change (PVP)

Dependencies added: HerbiePlugin

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +6/−0
• README.markdown +7/−7
• ad.cabal +10/−1
• include/instances.h +1/−1
• src/Numeric/AD/Internal/Dense.hs +3/−0
• src/Numeric/AD/Internal/Forward.hs +3/−0
• src/Numeric/AD/Internal/Kahn.hs +4/−0
• src/Numeric/AD/Internal/Reverse.hs +3/−0
• src/Numeric/AD/Internal/Sparse.hs +84/−20
• src/Numeric/AD/Internal/Tower.hs +18/−1
• src/Numeric/AD/Newton.hs +101/−3
• src/Numeric/AD/Rank1/Halley.hs +7/−1

@@ -1,3 +1,9 @@+4.3+---+* Made drastic improvements in the performance of `Tower` and `Sparse` modes thanks to the help of Björn von Sydow.+* Added constrained convex optimization.+* Incorporated some suggestions from [herbie](http://herbie.uwplse.org/z) for improving floating point accuracy.+ 4.2.4 ----- * Added `Newton.Double` modules for performance.

@@ -31,10 +31,10 @@     Prelude Numeric.AD> diff' (exp . log) 2     (2.0,1.0) -You can compute the derivative of a function with a constant parameter using `auto` from Numeric.AD.Types:+You can compute the derivative of a function with a constant parameter using `auto`: -    Prelude Numeric.AD Numeric.AD.Types> let t = 2.0 :: Double-    Prelude Numeric.AD Numeric.AD.Types> diff (\ x -> auto t * sin x) 0+    Prelude Numeric.AD> let t = 2.0 :: Double+    Prelude Numeric.AD> diff (\ x -> auto t * sin x) 0     2.0  You can use a symbolic numeric type, like the one from `simple-reflect` to obtain symbolic derivatives:@@ -65,22 +65,22 @@  The answer: -    Prelude Numeric.AD Numeric.AD.Types> headJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD> headJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     7.38905609893065  The gradient: -    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD> headJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     [14.7781121978613,7.38905609893065]  The hessian (n * n matrix of 2nd derivatives) -    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     [[29.5562243957226,22.16716829679195],[22.16716829679195,7.38905609893065]]  Or even higher order tensors of derivatives as a jet. -    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     [[[59.1124487914452,44.3343365935839],[44.3343365935839,14.7781121978613]],[[44.3343365935839,14.7781121978613],[14.7781121978613,7.38905609893065]]]  Note the redundant values caused by the various symmetries in the tensors. The `ad` library is careful to compute

@@ -1,5 +1,5 @@ name:          ad-version:       4.2.4+version:       4.3 license:       BSD3 license-File:  LICENSE copyright:     (c) Edward Kmett 2010-2015,@@ -75,6 +75,10 @@   type: git   location: git://github.com/ekmett/ad.git +flag herbie+  default: False+  manual: True+ library   hs-source-dirs: src   include-dirs: include@@ -110,6 +114,11 @@    if impl(ghc < 7.8)     build-depends: tagged >= 0.7 && < 1++  if flag(herbie)+    build-depends: HerbiePlugin >= 0.1 && < 0.2+    cpp-options: -DHERBIE+    ghc-options: -fplugin=Herbie    exposed-modules:     Numeric.AD

@@ -21,7 +21,7 @@   fromInteger n = auto (fromInteger n)   (+)          = (<+>) -- binary (+) 1 1   (-)          = binary (-) (auto 1) (auto (-1)) -- TODO: <-> ? as it is, this might be pretty bad for Tower-  (*)          = lift2 (*) (\x y -> (y, x))+  (*)          = mul -- lift2 (*) (\x y -> (y, x))   negate       = lift1 negate (const (auto (-1)))   abs          = lift1 abs signum   signum a     = lift1 signum (const zero) a

@@ -185,6 +185,9 @@     (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)     productRule dbi dci = dadb * dbi + dci * dadc +mul :: (Traversable f, Num a) => Dense f a -> Dense f a -> Dense f a+mul = lift2 (*) (\x y -> (y, x))+ #define BODY1(x)   (Traversable f, x) #define BODY2(x,y) (Traversable f, x, y) #define HEAD Dense f a

@@ -214,3 +214,6 @@ transposeWith f as = snd . mapAccumL go xss0 where   go xss b = (tail <$> xss, f b (head <$> xss))   xss0 = toList <$> as++mul :: (Num a) => Forward a -> Forward a -> Forward a+mul = lift2 (*) (\x y -> (y, x))

@@ -158,6 +158,10 @@     a = f pb pc     (dadb, dadc) = df (Id a) (Id pb) (Id pc) ++mul :: Num a => Kahn a -> Kahn a -> Kahn a+mul = lift2 (*) (\x y -> (y, x))+ #define HEAD Kahn a #include <instances.h> 

@@ -199,6 +199,9 @@     a = f pb pc     (dadb, dadc) = df (Id a) (Id pb) (Id pc) +mul :: (Reifies s Tape, Num a) => Reverse s a -> Reverse s a -> Reverse s a+mul = lift2 (*) (\x y -> (y, x))+ #define BODY1(x) (Reifies s Tape,x) #define BODY2(x,y) (Reifies s Tape,x,y) #define HEAD Reverse s a

@@ -39,6 +39,9 @@   , vgrads   , Grad(..)   , Grads(..)+  , terms+  , deriv+  , primal   ) where  import Prelude hiding (lookup)@@ -48,7 +51,7 @@ import Control.Comonad.Cofree import Control.Monad (join) import Data.Data-import Data.IntMap (IntMap, mapWithKey, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)+import Data.IntMap (IntMap, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup) import qualified Data.IntMap as IntMap import Data.Number.Erf import Data.Traversable@@ -74,13 +77,17 @@ -- | We only store partials in sorted order, so the map contained in a partial -- will only contain partials with equal or greater keys to that of the map in -- which it was found. This should be key for efficiently computing sparse hessians.--- there are only (n + k - 1) choose k distinct nth partial derivatives of a+-- there are only (n + k - 1) choose (k - 1) distinct nth partial derivatives of a -- function with k inputs. data Sparse a   = Sparse !a (IntMap (Sparse a))   | Zero   deriving (Show, Data, Typeable) +{-++These functions are now unused.+ dropMap :: Int -> IntMap a -> IntMap a dropMap n = snd . IntMap.split (n - 1) {-# INLINE dropMap #-}@@ -93,6 +100,7 @@     (fmap (* b) (dropMap n da))     (fmap (a *) (dropMap n db)) {-# INLINE times #-}+-}  vars :: (Traversable f, Num a) => f a -> f (Sparse a) vars = snd . mapAccumL var 0 where@@ -169,45 +177,42 @@ a <+> Zero = a Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs +-- The instances for Jacobian for Sparse and Tower are almost identical;+-- could easily be made exactly equal by small changes. instance Num a => Jacobian (Sparse a) where   type D (Sparse a) = Sparse a   unary f _ Zero = auto (f 0)-  unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs+  unary f dadb (Sparse pb bs) = Sparse (f pb) $ IntMap.map (* dadb) bs    lift1 f _ Zero = auto (f 0)-  lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) bs+  lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ IntMap.map (* df b) bs    lift1_ f _  Zero = auto (f 0)   lift1_ f df b@(Sparse pb bs) = a where-    a = Sparse (f pb) $ mapWithKey (times (df a b)) bs+    a = Sparse (f pb) $ IntMap.map ((df a b) *) bs    binary f _    _    Zero           Zero           = auto (f 0 0)-  binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ mapWithKey (times dadc) dc-  binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ mapWithKey (times dadb) db+  binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ IntMap.map (dadc *) dc+  binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ IntMap.map (dadb *) db   binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $-    unionWith (<+>)-      (mapWithKey (times dadb) db)-      (mapWithKey (times dadc) dc)+    unionWith (<+>)  (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)    lift2 f _  Zero             Zero = auto (f 0 0)-  lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ mapWithKey (times dadc) dc where dadc = snd (df zero c)-  lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ mapWithKey (times dadb) db where dadb = fst (df b zero)+  lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ IntMap.map (dadc *) dc where dadc = snd (df zero c)+  lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ IntMap.map (* dadb) db where dadb = fst (df b zero)   lift2 f df b@(Sparse pb db) c@(Sparse pc dc) = Sparse (f pb pc) da where     (dadb, dadc) = df b c-    da = unionWith (<+>)-      (mapWithKey (times dadb) db)-      (mapWithKey (times dadc) dc)+    da = unionWith (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)    lift2_ f _  Zero             Zero = auto (f 0 0)-  lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (mapWithKey (times (fst (df a b zero))) db)-  lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (mapWithKey (times (snd (df a zero c))) dc)+  lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (IntMap.map (fst (df a b zero) *) db)+  lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (IntMap.map (* snd (df a zero c)) dc)   lift2_ f df b@(Sparse pb db) c@(Sparse pc dc) = a where     (dadb, dadc) = df a b c     a = Sparse (f pb pc) da-    da = unionWith (<+>)-      (mapWithKey (times dadb) db)-      (mapWithKey (times dadc) dc)+    da = unionWith (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc) + #define HEAD Sparse a #include "instances.h" @@ -254,3 +259,62 @@ vgrads i = unpacks (unsafeGrads (packs i)) where   unsafeGrads f as = ds as $ apply f as {-# INLINE vgrads #-}++isZero :: Sparse a -> Bool+isZero Zero = True+isZero _ = False++-- |+-- A monomial is used to indicate order of differentiation.+-- For a k-ary function, it represented as a list of k non-negative Ints.+-- MI [n_0,n_1...n_{k-1}] denotes differentiation n_0 times with respect+-- to variable 0, n_1 times to variable 1, etc.+-- Trailing zeros omitted for efficiency.+--+-- Add 1 to variable k (i.e.differentiate once more wrt variable k).+incMonomial :: Int -> [Int] -> [Int]+incMonomial k [] = replicate k 0 ++ [1]+incMonomial 0 (a:as) = a+1:as+incMonomial k (a:as) = a:incMonomial (k-1) as++-- deriv f mi is the derivative of f of order mi (including higher derivatives).+deriv :: Sparse a -> [Int] -> Sparse a+deriv f mi = indx 0 mi f where+  indx _ [] f = f+  indx _ _ Zero = Zero+  indx v (0:as) f = indx (v+1) as f+  indx v (a:as) (Sparse _ df) = maybe Zero (indx v (a-1 : as)) (lookup v df)++-- The value of the derivative of (f*g) of order mi is+--       sum [a*primal (deriv f b)*primal (deriv g c) | (a,b,c) <- terms mi ]+-- It is a bit more complicated in mul' below, since we build the whole tree of+-- derivatives and want to prune the tree with Zeros as much as possible.+-- The number of terms in the sum for order MI as of differentiation has+-- sum (map (+1) as) terms, so this is *much* more efficient+-- than the naive recursive differentiation with 2^(sum as) terms.+-- The coefficients a, which collect equivalent derivatives, are suitable products+-- of binomial coefficients.+terms :: [Int]-> [(Integer,[Int],[Int])]+terms [] = [(1,[],[])]+terms (a:as) = concatMap (f ps) (zip (bins!!a) [0..a]) where+  ps = terms as+  bins = iterate next [1]+  next xs@(_:ts) = 1 : zipWith (+) xs ts ++ [1]+  next [] = error "impossible"+  f ps (b,k) = map (\(w,ks,is) -> (w*b,(k:ks),(a-k:is))) ps++mul :: Num a => Sparse a -> Sparse a -> Sparse a+mul Zero _ = Zero+mul _ Zero = Zero+mul f@(Sparse _ am) g@(Sparse _ bm) = Sparse (primal f * primal g) (derivs 0 []) where+  derivs v mi = IntMap.unions (map fn [v..kMax]) where+    fn w+      | and zs = IntMap.empty+      | otherwise = IntMap.singleton w (Sparse (sum ds) (derivs w mi'))+      where+        mi' = incMonomial w mi+        (zs,ds) = unzip (map derVal (terms mi'))+        derVal (bin,mif,mig) = (isZero fder || isZero gder, fromIntegral bin * primal fder * primal gder) where+          fder = deriv f mif+          gder = deriv g mig+  kMax = max (maximum (-1:IntMap.keys am)) (maximum (-1:IntMap.keys bm))

@@ -34,7 +34,7 @@   , tower   ) where -import Prelude hiding (all)+import Prelude hiding (all, sum) #if __GLASGOW_HASKELL__ < 710 import Control.Applicative hiding ((<**>)) #endif@@ -163,6 +163,23 @@ _        <**> Tower []  = auto 1 x        <**> Tower [y] = lift1 (**y) (\z -> y *^ z <**> Tower [y-1]) x x        <**> y         = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y++-- mul xs ys = [ sum [xs!!j * ys!!(k-j)*bin k j | j <- [0..k]] | k <- [0..] ]+-- adapted for efficiency and to handle finite lists xs, ys+mul:: Num a => Tower a -> Tower a -> Tower a+mul (Tower []) _ = Tower []+mul (Tower (a:as)) (Tower bs) = Tower (convs' [1] [a] as bs)+  where convs' _ _ _ [] = []+        convs' ps ars as bs = sumProd3 ps ars bs :+              case as of+                 [] -> convs'' (next' ps) ars bs+                 a:as -> convs' (next ps) (a:ars) as bs+        convs'' _ _ [] = undefined -- convs'' never called with last argument empty+        convs'' _ _ [_] = []+        convs'' ps ars (_:bs) = sumProd3 ps ars bs : convs'' (next' ps) ars bs+        next xs = 1 : zipWith (+) xs (tail xs) ++ [1] -- next row in Pascal's triangle+        next' xs = zipWith (+) xs (tail xs) ++ [1] -- end part of next row in Pascal's triangle+        sumProd3 as bs cs = sum (zipWith3 (\x y z -> x*y*z) as bs cs)  #define HEAD Tower a #include <instances.h>

@@ -1,8 +1,14 @@-{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE GADTs #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE DeriveTraversable #-}+{-# LANGUAGE ParallelListComp #-} ----------------------------------------------------------------------------- -- | -- Copyright   :  (c) Edward Kmett 2010-2015@@ -21,14 +27,18 @@   , fixedPoint   , extremum   -- * Gradient Ascent/Descent (Reverse AD)-  , gradientDescent+  , gradientDescent, constrainedDescent, CC(..), eval   , gradientAscent   , conjugateGradientDescent   , conjugateGradientAscent   , stochasticGradientDescent   ) where +#if __GLASGOW_HASKELL__ < 710+import Data.Foldable (Foldable, all, sum)+#else import Data.Foldable (all, sum)+#endif import Data.Reflection (Reifies) import Data.Traversable import Numeric.AD.Internal.Combinators@@ -38,7 +48,7 @@ import Numeric.AD.Internal.Reverse (Reverse, Tape) import Numeric.AD.Internal.Type (AD(..)) import Numeric.AD.Mode-import Numeric.AD.Mode.Reverse as Reverse (gradWith, gradWith')+import Numeric.AD.Mode.Reverse as Reverse (gradWith, gradWith', grad') import Numeric.AD.Rank1.Kahn as Kahn (Kahn, grad) import qualified Numeric.AD.Rank1.Newton as Rank1 import Prelude hiding (all, mapM, sum)@@ -122,6 +132,94 @@         x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx         (fx1, xgx1) = Reverse.gradWith' (,) f x1 {-# INLINE gradientDescent #-}++data SEnv (f :: * -> *) a = SEnv { sValue :: a, origEnv :: f a }+  deriving (Functor, Foldable, Traversable)++-- | Convex constraint, CC, is a GADT wrapper that hides the existential+-- ('s') which is so prevalent in the rest of the API.  This is an+-- engineering convenience for managing the skolems.+data CC f a where+  CC :: forall f a. (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> CC f a++-- |@constrainedDescent obj fs env@ optimizes the convex function @obj@+-- subject to the convex constraints @f <= 0@ where @f `elem` fs@. This is+-- done using a log barrier to model constraints (i.e. Boyd, Chapter 11.3).+-- The returned optimal point for the objective function must satisfy @fs@,+-- but the initial environment, @env@, needn't be feasible.+constrainedDescent :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)+                   => (forall s. Reifies s Tape => f (Reverse s a)+                                                -> Reverse s a)+                   -> [CC f a]+                   -> f a+                   -> [(a,f a)]+constrainedDescent objF [] env =+  map (\x -> (eval objF x, x)) (gradientDescent objF env)+constrainedDescent objF cs env =+    let s0       = 1 + maximum [eval c env | CC c <- cs]+        -- ^ s0 = max ( f_i(0) )+        cs'      = [CC (\(SEnv sVal rest) -> c rest - sVal) | CC c <- cs]+        -- ^ f_i' = f_i - s0  and thus f_i' <= 0+        envS     = SEnv s0 env+        -- feasible point for f_i', use gd to find feasiblity for f_i+        cc       = constrainedConvex' (CC sValue) cs' envS ((<=0) . sValue)+    in case dropWhile ((0 <) . fst) (take (2^(20::Int)) cc) of+        []                  -> []+        (_,envFeasible) : _ ->+            constrainedConvex' (CC objF) cs (origEnv envFeasible) (const True)+{-# INLINE constrainedDescent #-}++eval :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> a+eval f e = fst (grad' f e)+{-# INLINE eval #-}++-- | Like 'constrainedDescent' except the initial point must be feasible.+constrainedConvex' :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)+                   => CC f a+                   -> [CC f a]+                   -> f a+                   -> (f a -> Bool)+                   -> [(a,f a)]+constrainedConvex' objF cs env term =+  -- 1. Transform cs using a log barrier with increasing t values.+  let os   = map (mkOpt objF cs) tValues+  -- 2. Iteratively run gradientDescent on each os.+      envs =  [(undefined,env)] :+              [gD (snd $ last e) o+                          | o  <- os+                          | e  <- limEnvs+                          ]+      -- Obtain a finite number of elements from the initial len tValues - 1 lists.+      limEnvs = zipWith id nrSteps envs+  in dropWhile (not . term . snd) (concat $ drop 1 limEnvs)+ where+  tValues = map realToFrac $ take 64 $ iterate (*2) (2 :: a)+  nrSteps = [take 20 | _ <- [1..length tValues]] ++ [id]+  -- | `gD f e` is gradient descent with the evaulated result+  gD e (CC f)  = (eval f e, e) :+                 map (\x -> (eval f x, x)) (gradientDescent f e)+{-# INLINE constrainedConvex' #-}++-- @mkOpt u fs t@ converts an inequality convex problem (@u,fs@) into an+-- unconstrained convex problem using log barrier @u + -(1/t)log(-f_i)@.+-- As @t@ increases the approximation is more accurate but the gradient+-- decreases, making the gradient descent more expensive.+mkOpt :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)+      => CC f a -> [CC f a]+      -> a -> CC f a+mkOpt (CC o) xs t = CC (\e -> o e + sum (map (\(CC c) -> iHat t c e) xs))+{-# INLINE mkOpt #-}++iHat :: forall a f. (Traversable f, RealFloat a, Floating a, Ord a)+     => a+     -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a)+     -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a)+iHat t c e =+   let r = c e+   in if r >= 0 || isNaN r+        then 1  / 0+        else (-1 / auto t) * log( - (c  e))+{-# INLINE iHat #-}  -- | The 'stochasticGradientDescent' function approximates -- the true gradient of the constFunction by a gradient at

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} ----------------------------------------------------------------------------- -- | -- Copyright   :  (c) Edward Kmett 2010-2015@@ -49,7 +50,12 @@ findZero f = go where   go x = x : if x == xn then [] else go xn where     (y:y':y'':_) = diffs0 f x-    xn = x - 2*y*y'/(2*y'*y'-y*y'')+    xn = x - 2*y*y'/(2*y'*y'-y*y'') -- 9.606671960457536 bits error+       -- = x - recip (y'/y - y''/ y') -- "improved error" = 6.640625e-2 bits+       -- = x - y' / (y'/y/y' - y''/2) -- "improved error" = 1.4+#ifdef HERBIE+{-# ANN findZero "NoHerbie" #-}+#endif {-# INLINE findZero #-}  -- | The 'inverse' function inverts a scalar function using