sgd 0.4.0.1 → 0.5.0.0

raw patch · 15 files changed

+1569/−745 lines, 15 filesdep +data-defaultdep +hmatrixdep +pipesdep −lazy-iodep ~containersdep ~temporaryPVP ok

version bump matches the API change (PVP)

Dependencies added: data-default, hmatrix, pipes

Dependencies removed: lazy-io

Dependency ranges changed: containers, temporary

API changes (from Hackage documentation)

Files

• examples/example1.hs +0/−104
• sgd.cabal +18/−15
• src/Numeric/SGD.hs +205/−129
• src/Numeric/SGD/AdaDelta.hs +96/−0
• src/Numeric/SGD/DataSet.hs +133/−0
• src/Numeric/SGD/Dataset.hs +0/−103
• src/Numeric/SGD/Grad.hs +0/−133
• src/Numeric/SGD/LogSigned.hs +0/−85
• src/Numeric/SGD/Momentum.hs +56/−176
• src/Numeric/SGD/ParamSet.hs +465/−0
• src/Numeric/SGD/Sparse.hs +167/−0
• src/Numeric/SGD/Sparse/Grad.hs +131/−0
• src/Numeric/SGD/Sparse/LogSigned.hs +85/−0
• src/Numeric/SGD/Sparse/Momentum.hs +201/−0
• src/Numeric/SGD/Type.hs +12/−0

@@ -1,104 +0,0 @@-{-# LANGUAGE RecordWildCards #-}--import Control.Applicative ((<$>), (<*>))-import Control.Monad (replicateM)-import System.IO (hSetBuffering, stdout, BufferMode (NoBuffering))-import qualified System.Random as R-import qualified Data.Vector as V-import qualified Data.Vector.Unboxed as U-import qualified Numeric.SGD as S----------------------------------------------------------------------------------- Dataset generation----------------------------------------------------------------------------------- | Element of a dataset.-type Elem = [(Int, Double)]---- | Random dataset element.-elemR-    :: Int              -- ^ Maximum number of element items-    -> (Int, Int)       -- ^ Range for item's first component-    -> (Double, Double) -- ^ Range for item's second component-    -> IO Elem          -- ^ Result-elemR nMax xr yr = do-    n <- R.randomRIO (0, max 0 nMax)-    replicateM n ((,) <$> R.randomRIO xr <*> R.randomRIO yr)---- | Random dataset.-dataSetR-    :: Int              -- ^ Dataset size-    -> Int              -- ^ Number of model parameters-    -> Int              -- ^ Maximum number of items in data element-    -> (Double, Double) -- ^ Range for item's second component-    -> IO (V.Vector Elem)   -- ^ Result-dataSetR m n k yRan =-    V.fromList <$> replicateM m (elemR k (0, n-1) yRan)----------------------------------------------------------------------------------- Objective function and gradient----------------------------------------------------------------------------------- | An objective function. The SGD method can be used when--- the objective function is defined in a form of a sum.-goal :: S.Para -> [Elem] -> Double-goal para =-    sum . map perElem-  where-    perElem xs = sum-        [ (para U.! k - x) ^ (2 :: Int)-        | (k, x) <- xs ]---- | Since the goal function has a form of a sum, it is sufficient to define--- the gradient over one element only. The gradient with respect to the dataset--- is a sum of gradients over its individual elements.-grad :: S.Para -> Elem -> S.Grad-grad para xs = S.fromList-    -- [ (k, 2 * (x - para U.! k))-    [ (k, 2 * (para U.! k - x))-    | (k, x) <- xs ]---- | Negate gradient. We use it to find the minimum of the objective function.-negGrad :: (S.Para -> Elem -> S.Grad)-        -> (S.Para -> Elem -> S.Grad)-negGrad g para x = fmap negate (g para x)----------------------------------------------------------------------------------- SGD----------------------------------------------------------------------------------- | Notification run by the sgdM function every parameters update.-notify :: S.SgdArgs -> V.Vector Elem -> S.Para -> Int -> IO ()-notify S.SgdArgs{..} dataSet para k =-    if doneTotal k /= doneTotal (k - 1)-        then do-            let n = doneTotal k-                x = goal para (V.toList dataSet)-            putStrLn ("\n" ++ "[" ++ show n ++ "] f = " ++ show x)-        else-            putStr "."-  where-    doneTotal :: Int -> Int-    doneTotal = floor . done-    done :: Int -> Double-    done i-        = fromIntegral (i * batchSize)-        / fromIntegral (V.length dataSet)---- | Run the monadic version of SGD.-runSgdM-    :: Int              -- ^ Dataset size-    -> Int              -- ^ Number of model parameters-    -> Int              -- ^ Maximum number of items in data element-    -> S.SgdArgs        -- ^ SGD parameters-    -> IO S.Para-runSgdM m n k sgdArgs = do-    dataSet <- dataSetR m n k (-10, 10)-    let para = U.replicate n 0-    hSetBuffering stdout NoBuffering-    S.sgdM sgdArgs (notify sgdArgs dataSet) (negGrad grad) dataSet para---- | Run the monadic version of SGD with some default parameter values.-main = do-    let sgdArgs = S.sgdArgsDefault { S.iterNum = 50 }-    runSgdM 1000 1000000 10 sgdArgs

@@ -1,31 +1,28 @@ name:               sgd-version:            0.4.0.1+version:            0.5.0.0 synopsis:           Stochastic gradient descent description:-    Implementation of a Stochastic Gradient Descent optimization method.-    See examples directory in the source package for examples of usage.+    Stochastic gradient descent library.     .-    It is a preliminary implementation of the SGD method and API may change-    in future versions.+    Import "Numeric.SGD" to use the library. license:            BSD3 license-file:       LICENSE cabal-version:      >= 1.6-copyright:          Copyright (c) 2012 IPI PAN+copyright:          Copyright (c) 2012-2019 Jakub Waszczuk author:             Jakub Waszczuk maintainer:         waszczuk.kuba@gmail.com stability:          experimental-category:           Math, Algorithms+category:           Math homepage:           https://github.com/kawu/sgd build-type:         Simple -extra-source-files: examples/example1.hs- library     hs-source-dirs: src      build-depends:         base            >= 4        && < 5-      , containers      >= 0.4      && < 0.6+      , containers      >= 0.4      && < 0.7+      , pipes           >= 4.3      && < 4.4       , vector          >= 0.10     && < 0.13       , random          >= 1.0      && < 1.2       , primitive       >= 0.5      && < 0.7@@ -36,15 +33,21 @@       , bytestring      >= 0.9      && < 0.11       , mtl             >= 2.0      && < 2.3       , filepath        >= 1.3      && < 1.5-      , temporary       >= 1.1      && < 1.3-      , lazy-io         >= 0.1      && < 0.2+      , temporary       >= 1.1      && < 1.4+      , hmatrix         >= 0.19     && < 0.20+      , data-default    >= 0.7      && < 0.8      exposed-modules:         Numeric.SGD+      , Numeric.SGD.Type+      , Numeric.SGD.DataSet+      , Numeric.SGD.ParamSet+      , Numeric.SGD.AdaDelta       , Numeric.SGD.Momentum-      , Numeric.SGD.Dataset-      , Numeric.SGD.LogSigned-      , Numeric.SGD.Grad+      , Numeric.SGD.Sparse+      , Numeric.SGD.Sparse.Momentum+      , Numeric.SGD.Sparse.LogSigned+      , Numeric.SGD.Sparse.Grad      ghc-options: -Wall -O2 

@@ -1,164 +1,240 @@ {-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveGeneric #-}  --- | Stochastic gradient descent implementation using mutable--- vectors for efficient update of the parameters vector.--- A user is provided with the immutable vector of parameters--- so he is able to compute the gradient outside of the IO monad.--- Currently only the Gaussian priors are implemented.+-- | Main module of the stochastic gradient descent (SGD) library.  ----- This is a preliminary version of the SGD library and API may change--- in future versions.+-- SGD is a method for optimizing a global objective function defined as a sum+-- of smaller, differentiable functions.  The individual component functions+-- share the same set of parameters, represented by the `ParamSet` class.+--+-- To perform SGD, the gradients of the individual functions need to be+-- determined.  This can be done manually or automatically, using one of the+-- automatic differentiation libraries (ad, backprop) available in Haskell.+--+-- For instance, let's say we have a list of functions defined as:+--+-- > funs = [\x -> 0.3*x^2, \x -> -2*x, const 3, sin]+--+-- The global objective is then defined as:+--+-- > objective x = sum $ map ($x) funs+--+-- We can manually determine the individual derivatives:+--+-- > derivs = [\x -> 0.6*x, const (-2), const 0, cos]+--+-- or use an automatic differentiation library, for instance:+--+-- > import qualified Numeric.AD as AD+-- > derivs = map+-- >   (\k -> AD.diff (funs !! k))+-- >   [0..length funs-1]+--+-- Finally, `run` allows to approach a (potentially local) minimum of the+-- global objective function:+--+-- >>> run (momentum def id) (take 10000 $ cycle derivs) 0.0+-- 4.180177042912455+--+-- where:+-- +--     * @(take 10000 $ cycle derivs)@ is the stream of training examples+--     * @(momentum def id)@ is the selected SGD variant (`Mom.momentum`),+--     supplied with the default configuration (`def`) and the function (`id`)+--     for calculating the gradient from a training example+--     * @0.0@ is the initial parameter value   module Numeric.SGD-( SgdArgs (..)-, sgdArgsDefault-, Para-, sgd-, module Numeric.SGD.Grad-, module Numeric.SGD.Dataset-) where+  (+  -- * SGD variants+    Mom.momentum+  , Ada.adaDelta +  -- * Pure SGD+  , run -import           Control.Monad (forM_)-import qualified System.Random as R-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Unboxed.Mutable as UM-import qualified Control.Monad.Primitive as Prim+  -- * IO-based SGD+  , Config (..)+  , runIO -import           Numeric.SGD.Grad-import           Numeric.SGD.Dataset+  -- * Combinators+  , pipeSeq+  , pipeRan+  , result+  , every +  -- * Re-exports+  , def+  ) where --- | SGD parameters controlling the learning process.-data SgdArgs = SgdArgs-    { -- | Size of the batch-      batchSize :: Int-    -- | Regularization variance-    , regVar    :: Double-    -- | Number of iterations-    , iterNum   :: Double-    -- | Initial gain parameter-    , gain0     :: Double-    -- | After how many iterations over the entire dataset-    -- the gain parameter is halved-    , tau       :: Double } +import           GHC.Generics (Generic)+import           Numeric.Natural (Natural) --- | Default SGD parameter values.-sgdArgsDefault :: SgdArgs-sgdArgsDefault = SgdArgs-    { batchSize = 30-    , regVar    = 10-    , iterNum   = 10-    , gain0     = 1-    , tau       = 5 }+import qualified System.Random as R +import           Control.Monad (when, forM_) --- | Vector of parameters.-type Para       = U.Vector Double+import           Data.Functor.Identity (Identity(..))+import qualified Data.IORef as IO+import           Data.Default +import qualified Pipes as P+import qualified Pipes.Prelude as P+import           Pipes ((>->)) --- | Type synonym for mutable vector with Double values.-type MVect      = UM.MVector (Prim.PrimState IO) Double+import qualified Numeric.SGD.AdaDelta as Ada+import qualified Numeric.SGD.Momentum as Mom+import           Numeric.SGD.Type+import           Numeric.SGD.ParamSet+import           Numeric.SGD.DataSet  --- | A stochastic gradient descent method.--- A notification function can be used to provide user with--- information about the progress of the learning.-sgd-    :: SgdArgs                  -- ^ SGD parameter values-    -> (Para -> Int -> IO ())   -- ^ Notification run every update-    -> (Para -> x -> Grad)      -- ^ Gradient for dataset element-    -> Dataset x                -- ^ Dataset-    -> Para                     -- ^ Starting point-    -> IO Para                  -- ^ SGD result-sgd SgdArgs{..} notify mkGrad dataset x0 = do-  u <- UM.new (U.length x0)-  doIt u 0 (R.mkStdGen 0) =<< U.thaw x0-  where-    -- Gain in k-th iteration.-    gain k = (gain0 * tau) / (tau + done k)+-------------------------------+-- Pure SGD+------------------------------- -    -- Number of completed iterations over the full dataset.-    done :: Int -> Double-    done k-        = fromIntegral (k * batchSize)-        / fromIntegral (size dataset)-    -- doneTotal :: Int -> Int-    -- doneTotal = floor . done -    -- Regularization (Guassian prior)-    regularization k = regCoef-      where-        regCoef = (1.0 - gain k * iVar) ** coef-        iVar = 1.0 / regVar-        coef = fromIntegral batchSize-             / fromIntegral (size dataset)+-- | Traverse all the elements in the training data stream in one pass,+-- calculate the subsequent gradients, and apply them progressively starting+-- from the initial parameter values.+--+-- Consider using `runIO` if your training dataset is large.+run+  :: (ParamSet p)+  => SGD Identity e p+    -- ^ Selected SGD method+  -> [e]+    -- ^ Training data stream+  -> p+    -- ^ Initial parameters+  -> p+run sgd dataSet p0 = runIdentity $+  result p0 +    (P.each dataSet >-> sgd p0) ---     -- Regularization (Guassian prior) after a full dataset pass---     regularization k = 1.0 - (gain k / regVar) -    doIt u k stdGen x-      | done k > iterNum = do-        frozen <- U.unsafeFreeze x-        notify frozen k-        return frozen-      | otherwise = do-        (batch, stdGen') <- sample stdGen batchSize dataset+------------------------------- +-- Higher-level SGD+------------------------------- -        -- Regularization-        -- when (doneTotal (k - 1) /= doneTotal k) $ do-        --   <- we now apply regularization each step rather than each-        --      dataset pass-        let regParam = regularization k-        -- putStrLn $ "\nApplying regularization (params *= " ++ show regParam ++ ")"-        scale regParam x ---         -- Regularization---         when (doneTotal (k - 1) /= doneTotal k) $ do---           let regParam = regularization k---           putStrLn $ "\nApplying regularization (params *= " ++ show regParam ++ ")"---           scale regParam x+-- | High-level IO-based SGD configuration+data Config = Config+  { iterNum :: Natural+    -- ^ Number of iteration over the entire training dataset+  , batchRandom :: Bool+    -- ^ Should the mini-batch be selected at random?  If not, the subsequent+    -- training elements will be picked sequentially.  Random selection gives+    -- no guarantee of seeing each training sample in every epoch.+  , reportEvery :: Double+    -- ^ How often the value of the objective function should be reported (with+    -- `1` meaning once per pass over the training data)+  } deriving (Show, Eq, Ord, Generic) -        -- Freeze mutable vector of parameters. The frozen version is-        -- then supplied to external mkGrad function provided by user.-        frozen <- U.unsafeFreeze x-        notify frozen k+instance Default Config where+  def = Config+    { iterNum = 100+    , batchRandom = False+    , reportEvery = 1.0+    } -        -- let grad = M.unions (map (mkGrad frozen) batch)-        let grad = parUnions (map (mkGrad frozen) batch)-        addUp grad u-        scale (gain k) u -        x' <- U.unsafeThaw frozen-        u `addTo` x'-        doIt u (k+1) stdGen' x'+-- | Perform SGD in the IO monad, regularly reporting the value of the+-- objective function on the entire dataset.  A higher-level wrapper which+-- should be convenient to use when the training dataset is large.+--+-- An alternative is to use the simpler function `run`, or to build a custom+-- SGD pipeline based on lower-level combinators (`pipeSeq`, `Ada.adaDelta`,+-- `every`, `result`, etc.).+runIO+  :: (ParamSet p)+  => Config+    -- ^ SGD configuration+  -> SGD IO e p+    -- ^ Selected SGD method+  -> (e -> p -> Double)+    -- ^ Value of the objective function on a sample element (needed for model+    -- quality reporting)+  -> DataSet e+    -- ^ Training dataset+  -> p+    -- ^ Initial parameter values+  -> IO p+runIO Config{..} sgd quality0 dataSet net0 = do+  report net0+  result net0 $ pipeSeq dataSet+    >-> sgd net0+    >-> P.take realIterNum+    >-> every realReportPeriod report+  where+    -- Iteration scaling+    iterScale x = fromIntegral (size dataSet) * x+    -- Number of iterations and reporting period+    realIterNum = ceiling $ iterScale (fromIntegral iterNum :: Double)+    realReportPeriod = ceiling $ iterScale reportEvery+    -- Network quality over the entire training dataset+    report net = do+      putStr . show =<< quality net+      putStrLn $ " (norm_2 = " ++ show (norm_2 net) ++ ")"+    quality net = do+      res <- IO.newIORef 0.0+      forM_ [0 .. size dataSet - 1] $ \ix -> do+        x <- elemAt dataSet ix+        IO.modifyIORef' res (+ quality0 x net)+      IO.readIORef res  --- | Add up all gradients and store results in normal domain.-addUp :: Grad -> MVect -> IO ()-addUp grad v = do-    UM.set v 0-    forM_ (toList grad) $ \(i, x) -> do-        y <- UM.unsafeRead v i-        UM.unsafeWrite v i (x + y)+------------------------------- +-- Lower-level combinators+-------------------------------  --- | Scale the vector by the given value.-scale :: Double -> MVect -> IO ()-scale c v = do-    forM_ [0 .. UM.length v - 1] $ \i -> do-        y <- UM.unsafeRead v i-        UM.unsafeWrite v i (c * y)+-- | Pipe the dataset sequentially in a loop.+pipeSeq :: DataSet e -> P.Producer e IO ()+pipeSeq dataSet = do+  go (0 :: Int)+  where+    go k+      | k >= size dataSet = go 0+      | otherwise = do+          x <- P.lift $ elemAt dataSet k+          P.yield x+          go (k+1)  --- | Apply gradient to the parameters vector, that is add the first vector to--- the second one.-addTo :: MVect -> MVect -> IO ()-addTo w v = do-    forM_ [0 .. UM.length v - 1] $ \i -> do-        x <- UM.unsafeRead v i-        y <- UM.unsafeRead w i-        UM.unsafeWrite v i (x + y)+-- | Pipe the dataset randomly in a loop.+pipeRan :: DataSet e -> P.Producer e IO ()+pipeRan dataSet = do+  x <- P.lift $ do+    ix <- R.randomRIO (0, size dataSet - 1)+    elemAt dataSet ix+  P.yield x+  pipeRan dataSet+++-- | Extract the result of the SGD calculation (the last parameter+-- set flowing downstream).+result+  :: (Monad m)+  => p     +    -- ^ Default value (in case the stream is empty)+  -> P.Producer p m ()+    -- ^ Stream of parameter sets+  -> m p+result pDef = fmap (maybe pDef id) . P.last+++-- | Apply the given function every `k` param sets flowing downstream.+every :: (Monad m) => Int -> (p -> m ()) -> P.Pipe p p m x+every k f = do+  go (1 `mod` k)+  where+    go i = do+      paramSet <- P.await+      when (i == 0) $ do+        P.lift $ f paramSet+      P.yield paramSet+      go $ (i+1) `mod` k

@@ -0,0 +1,96 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveGeneric #-}+++-- | AdaDelta algorithm as described in the following paper:+--+--     * https://arxiv.org/pdf/1212.5701.pdf+++module Numeric.SGD.AdaDelta+  ( Config(..)+  , adaDelta+  ) where+++import           GHC.Generics (Generic)++import           Prelude hiding (div)+-- import           Control.Monad (when)++import           Data.Default++import qualified Pipes as P++import           Numeric.SGD.Type+import           Numeric.SGD.ParamSet+-- import           Numeric.SGD.Args+++-- | AdaDelta configuration+data Config = Config+  { decay :: Double+    -- ^ Exponential decay parameter+  , eps   :: Double+    -- ^ Epsilon value+  } deriving (Show, Eq, Ord, Generic)++instance Default Config where+  def = Config+    { decay = 0.9+    , eps = 1.0e-6+    }+++-- | Perform gradient descent using the AdaDelta algorithm.  +-- See "Numeric.SGD.AdaDelta" for more information.+adaDelta+  :: (Monad m, ParamSet p)+  => Config+    -- ^ AdaDelta configuration+  -> (e -> p -> p)+    -- ^ Gradient on a training element+  -> SGD m e p+adaDelta Config{..} gradient net0 =++  let zr = zero net0 +   in go (0 :: Integer) zr zr zr net0++  where++    go k expSqGradPrev expSqDeltaPrev deltaPrev net = do+      x <- P.await+      let grad = gradient x net+          expSqGrad = scale decay expSqGradPrev+                `add` scale (1-decay) (square grad)+          rmsGrad = squareRoot (pmap (+eps) expSqGrad)+          expSqDelta = scale decay expSqDeltaPrev+                 `add` scale (1-decay) (square deltaPrev)+          rmsDelta = squareRoot (pmap (+eps) expSqDelta)+          delta = (rmsDelta `mul` grad) `div` rmsGrad+          newNet = net `sub` delta+      P.yield newNet+      go (k+1) expSqGrad expSqDelta delta newNet+++-------------------------------+-- Utils+-------------------------------+++-- | Scaling+scale :: ParamSet p => Double -> p -> p+scale x = pmap (*x)+{-# INLINE scale #-}+++-- | Root square+squareRoot :: ParamSet p => p -> p+squareRoot = pmap sqrt+{-# INLINE squareRoot #-}+++-- | Square+square :: ParamSet p => p -> p+square x = x `mul` x+{-# INLINE square #-}

@@ -0,0 +1,133 @@+{-# LANGUAGE RecordWildCards #-}+++-- | Dataset abstraction.+++module Numeric.SGD.DataSet+( +-- * Dataset+  DataSet (..)+-- * Reading+, loadData+, randomSample+-- * Construction+, withVect+, withDisk+-- , withData+) where+++import           Control.Monad (forM_)+import qualified Control.Monad.State.Strict as S++import           System.IO.Temp (withTempDirectory)+import           System.IO.Unsafe (unsafeInterleaveIO)+import           System.FilePath ((</>))+import qualified System.Random as R++import           Data.Binary (Binary, encodeFile, decode)+import qualified Data.ByteString as B+import qualified Data.ByteString.Lazy as BL+import qualified Data.Vector as V+++------------------------------- +-- Type+-------------------------------+++-- | Dataset stored on a disk+data DataSet elem = DataSet+  { size :: Int +    -- ^ The size of the dataset; the individual indices are+    -- [0, 1, ..., size - 1]+  , elemAt :: Int -> IO elem+    -- ^ Get the dataset element with the given identifier+  }+++-------------------------------------------+-- Reading+-------------------------------------------+++-- | Lazily load the entire dataset from a disk.+loadData :: DataSet a -> IO [a]+loadData DataSet{..} = lazyMapM elemAt [0 .. size - 1]+++-- -- | A dataset sample of the given size.+-- sample :: R.RandomGen g => g -> Int -> DataSet a -> IO ([a], g)+-- sample g 0 _       = return ([], g)+-- sample g n dataset = do+--     (xs, g') <- sample g (n-1) dataset+--     let (i, g'') = R.next g'+--     x <- dataset `elemAt` (i `mod` size dataset)+--     return (x:xs, g'')+++-- | Random dataset sample with a specified number of elements (loaded eagerly)+randomSample :: Int -> DataSet a -> IO [a]+randomSample k dataSet+  | k <= 0 = return []+  | otherwise = do+      ix <- R.randomRIO (0, size dataSet - 1)+      x <- elemAt dataSet ix+      (x:) <$> randomSample (k-1) dataSet+++-------------------------------------------+-- Construction+-------------------------------------------+++-- | Construct dataset from a list of elements, store it as a vector, and run+-- the given handler.+withVect :: [a] -> (DataSet a -> IO b) -> IO b+withVect xs handler =+    handler dataset+  where+    v = V.fromList xs+    dataset = DataSet+        { size      = V.length v+        , elemAt    = \k -> return (v V.! k) }+++-- | Construct dataset from a list of elements, store it on a disk and run the+-- given handler.  Training elements must have the `Binary` instance for this+-- function to work.+withDisk :: Binary a => [a] -> (DataSet a -> IO b) -> IO b+withDisk xs handler = withTempDirectory "." ".sgd" $ \tmpDir -> do+    -- We use state monad to compute the number of dataset elements. +    n <- flip S.execStateT 0 $ forM_ (zip xs [0 :: Int ..]) $ \(x, ix) -> do+        S.lift $ encodeFile (tmpDir </> show ix) x+        S.modify (+1)++    -- Avoid decodeFile laziness when using some older versions of the binary+    -- library (as of year 2019, this could be probably simplified)+    let at ix = do+        cs <- B.readFile (tmpDir </> show ix)+        return . decode $ BL.fromChunks [cs]++    handler $ DataSet {size = n, elemAt = at}+++-------------------------------------------+-- Lazy IO Utils+-------------------------------------------+++-- | Lazily evaluate each action in the sequence from left to right,+-- and collect the results.+lazySequence :: [IO a] -> IO [a]+lazySequence (mx:mxs) = do+    x   <- mx+    xs  <- unsafeInterleaveIO (lazySequence mxs)+    return (x : xs)+lazySequence [] = return []+++-- | `lazyMapM` f is equivalent to `lazySequence` . `map` f.+lazyMapM :: (a -> IO b) -> [a] -> IO [b]+lazyMapM f = lazySequence . map f

@@ -1,103 +0,0 @@-{-# LANGUAGE RecordWildCards #-}----- | Dataset abstraction.---module Numeric.SGD.Dataset-( --- * Dataset-  Dataset (..)--- * Reading-, loadData-, sample--- * Construction-, withVect-, withDisk-, withData-) where---import           Control.Monad (forM_)-import           Data.Binary (Binary, encodeFile, decode)-import qualified Data.ByteString as B-import qualified Data.ByteString.Lazy as BL-import           System.IO.Temp (withTempDirectory)-import           System.FilePath ((</>))-import qualified System.Random as R-import qualified Data.Vector as V-import qualified Control.Monad.LazyIO as LazyIO-import qualified Control.Monad.State.Strict as S----- | A dataset with elements of type @a@.-data Dataset a = Dataset {-    -- | A size of the dataset.-      size      :: Int-    -- | Get dataset element with a given index.  The set of indices-    -- is of a {0, 1, .., size - 1} form.-    , elemAt    :: Int -> IO a }------------------------------------------------- Reading------------------------------------------------- | Lazily load dataset from a disk.-loadData :: Dataset a -> IO [a]-loadData Dataset{..} = LazyIO.mapM elemAt [0 .. size - 1]----- | A dataset sample of the given size.-sample :: R.RandomGen g => g -> Int -> Dataset a -> IO ([a], g)-sample g 0 _       = return ([], g)-sample g n dataset = do-    (xs, g') <- sample g (n-1) dataset-    let (i, g'') = R.next g'-    x <- dataset `elemAt` (i `mod` size dataset)-    return (x:xs, g'')------------------------------------------------- Construction------------------------------------------------- | Construct dataset from a vector of elements and run the--- given handler.-withVect :: [a] -> (Dataset a -> IO b) -> IO b-withVect xs handler =-    handler dataset-  where-    v = V.fromList xs-    dataset = Dataset-        { size      = V.length v-        , elemAt    = \k -> return (v V.! k) }----- | Construct dataset from a list of elements, store it on a disk--- and run the given handler.-withDisk :: Binary a => [a] -> (Dataset a -> IO b) -> IO b-withDisk xs handler = withTempDirectory "." ".sgd" $ \tmpDir -> do-    -- We use state monad to compute the number of dataset elements. -    n <- flip S.execStateT 0 $ forM_ (zip xs [0 :: Int ..]) $ \(x, ix) -> do-        S.lift $ encodeFile (tmpDir </> show ix) x-        S.modify (+1)--    -- We need to avoid decodeFile laziness when using some older-    -- versions of the binary library.-    let at ix = do-        cs <- B.readFile (tmpDir </> show ix)-        return . decode $ BL.fromChunks [cs]--    handler $ Dataset {size = n, elemAt = at}----- | Use disk or vector dataset representation depending on--- the first argument: when `True`, use `withDisk`, otherwise--- use `withVect`.-withData :: Binary a => Bool -> [a] -> (Dataset a -> IO b) -> IO b-withData x = case x of-    True    -> withDisk-    False   -> withVect

@@ -1,133 +0,0 @@-{-# LANGUAGE CPP #-}---- | A gradient is represented by an IntMap from gradient indices--- to values. Elements with no associated values in the gradient--- are assumed to have a 0 value assigned. Such elements are--- not interesting: when adding the gradient to the vector of--- parameters, only nonzero elements are taken into account.--- --- Each value associated with a gradient position is a pair of--- positive and negative components. They are stored separately--- to ensure high accuracy of computation results.--- Besides, both positive and negative components are stored--- in a logarithmic domain.--module Numeric.SGD.Grad-( Grad-, empty-, add-, addL-, fromList-, fromLogList-, toList-, parUnions-) where--import Data.List (foldl')-import Control.Applicative ((<$>), (<*>))-import Control.Monad.Par (Par, runPar, get)-#if MIN_VERSION_containers(0,4,2)-import Control.Monad.Par (spawn)-#else-import Control.DeepSeq (deepseq)-import Control.Monad.Par (spawn_)-#endif-#if MIN_VERSION_containers(0,5,0)-import qualified Data.IntMap.Strict as M-#else-import qualified Data.IntMap as M-#endif--import Numeric.SGD.LogSigned---- | Gradient with nonzero values stored in a logarithmic domain.--- Since values equal to zero have no impact on the update phase--- of the SGD method, it is more efficient to not to store those--- components in the gradient.-type Grad = M.IntMap LogSigned--{-# INLINE insertWith #-}-insertWith :: (a -> a -> a) -> M.Key -> a -> M.IntMap a -> M.IntMap a-#if MIN_VERSION_containers(0,5,0)-insertWith = M.insertWith-#elif MIN_VERSION_containers(0,4,1)-insertWith = M.insertWith'-#else-insertWith f k x m = -    M.alter g k m-  where-    g my = case my of-        Nothing -> Just x-        Just y  ->-            let z = f x y-            in  z `seq` Just z-#endif---- | Add normal-domain double to the gradient at the given position.-{-# INLINE add #-}-add :: Grad -> Int -> Double -> Grad-add grad i y = insertWith (+) i (logSigned y) grad ----- | Add log-domain, singed number to the gradient at the given position.-{-# INLINE addL #-}-addL :: Grad -> Int -> LogSigned -> Grad-addL grad i y = insertWith (+) i y grad ---- | Construct gradient from a list of (index, value) pairs.--- All values from the list are added at respective gradient--- positions.-{-# INLINE fromList #-}-fromList :: [(Int, Double)] -> Grad-fromList =-    let ins grad (i, y) = add grad i y-    in  foldl' ins empty---- | Construct gradient from a list of (index, signed, log-domain number)--- pairs.  All values from the list are added at respective gradient--- positions.-{-# INLINE fromLogList #-}-fromLogList :: [(Int, LogSigned)] -> Grad-fromLogList =-    let ins grad (i, y) = addL grad i y-    in  foldl' ins empty---- | Collect gradient components with values in normal domain.-{-# INLINE toList #-}-toList :: Grad -> [(Int, Double)]-toList =-    let unLog (i, x) = (i, toNorm x)-    in  map unLog . M.assocs---- | Empty gradient, i.e. with all elements set to 0.-{-# INLINE empty #-}-empty :: Grad-empty = M.empty---- | Perform parallel unions operation on gradient list. --- Experimental version.-parUnions :: [Grad] -> Grad-parUnions [] = error "parUnions: empty list"-parUnions xs = runPar (parUnionsP xs)---- | Parallel unions in the Par monad.-parUnionsP :: [Grad] -> Par Grad-parUnionsP [x] = return x-parUnionsP zs  = do-    let (xs, ys) = split zs-#if MIN_VERSION_containers(0,4,2)-    xsP <- spawn (parUnionsP xs)-    ysP <- spawn (parUnionsP ys)-    M.unionWith (+) <$> get xsP <*> get ysP-#else-    xsP <- spawn_ (parUnionsP xs)-    ysP <- spawn_ (parUnionsP ys)-    x <- M.unionWith (+) <$> get xsP <*> get ysP-    M.elems x `deepseq` return x-#endif-  where-    split []        = ([], [])-    split (x:[])    = ([x], [])-    split (x:y:rest)  =-        let (xs, ys) = split rest-        in  (x:xs, y:ys)

@@ -1,85 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}---- | Module provides data type for signed log-domain calculations.--module Numeric.SGD.LogSigned-( LogSigned (..)-, logSigned-, fromPos-, fromNeg-, toNorm-, toLogFloat-) where--import qualified Data.Number.LogFloat as L-import Data.Function (on)-import Control.DeepSeq (NFData(..))---- | Signed real value in the logarithmic domain.-data LogSigned = LogSigned-    { pos :: {-# UNPACK #-} !L.LogFloat     -- ^ Positive component-    , neg :: {-# UNPACK #-} !L.LogFloat     -- ^ Negative component-    } deriving Show--instance Eq LogSigned where-    (==) = (==) `on` toLogFloat--instance Ord LogSigned where-    compare = compare `on` toLogFloat---- All fields are strict and unpacked.-instance NFData LogSigned where-    rnf (LogSigned p q) = p `seq` q `seq` ()---- | Smart LogSigned constructor.-{-# INLINE logSigned #-}-logSigned :: Double -> LogSigned-logSigned x-    | x > 0     = LogSigned (L.logFloat x) zero-    | x < 0     = LogSigned zero (L.logFloat (-x))-    | otherwise = LogSigned zero zero---- | Make LogSigned from a positive, log-domain number.-{-# INLINE fromPos #-}-fromPos :: L.LogFloat -> LogSigned-fromPos x = LogSigned x zero---- | Make LogSigned from a negative, log-domain number.-{-# INLINE fromNeg #-}-fromNeg :: L.LogFloat -> LogSigned-fromNeg x = LogSigned zero x---- | Shift LogSigned to a normal domain.-{-# INLINE toNorm #-}-toNorm :: LogSigned -> Double-toNorm (LogSigned x y) = L.fromLogFloat x - L.fromLogFloat y---- | Change the 'LogSigned' to either negative 'Left' 'L.LogFloat'--- or positive 'Right' 'L.LogFloat'.-toLogFloat :: LogSigned -> Either L.LogFloat L.LogFloat-toLogFloat x = case signum x of-    -1  -> Left  $ neg x - pos x-    1   -> Right $ pos x - neg x-    _   -> Right $ L.logFloat (0 :: Double)--instance Num LogSigned where-    LogSigned x y + LogSigned x' y' =-        LogSigned (x + x') (y + y')-    LogSigned x y * LogSigned x' y' =-        LogSigned (x*x' + y*y') (x*y' + y*x')-    LogSigned x y - LogSigned x' y' =-        LogSigned (x + y') (y + x')-    negate  (LogSigned x y) = LogSigned y x-    abs     (LogSigned x y)-        | x >= y    = LogSigned x y-        | otherwise = LogSigned y x-    signum (LogSigned x y)-        | x > y     =  1-        | x < y     = -1-        | otherwise =  0-    fromInteger = logSigned . fromInteger--{-# INLINE zero #-}-zero :: L.LogFloat-zero = L.logFloat (0 :: Double)

@@ -1,201 +1,81 @@ {-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveGeneric #-}  --- | A version of `Numeric.SGD` extended with momentum.+-- | Stochastic gradient descent with momentum, following:+--+--     * http://ruder.io/optimizing-gradient-descent/index.html#momentum   module Numeric.SGD.Momentum-( SgdArgs (..)-, sgdArgsDefault-, Para-, sgd-, module Numeric.SGD.Grad-, module Numeric.SGD.Dataset-) where---import           Control.Monad (forM_, when)-import qualified System.Random as R-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Unboxed.Mutable as UM-import qualified Control.Monad.Primitive as Prim--import           Numeric.SGD.Grad-import           Numeric.SGD.Dataset+  ( Config(..)+  , momentum+  ) where  --- | SGD parameters controlling the learning process.-data SgdArgs = SgdArgs-    { -- | Size of the batch-      batchSize :: Int-    -- | Regularization variance-    , regVar    :: Double-    -- | Number of iterations-    , iterNum   :: Double-    -- | Initial gain parameter-    , gain0     :: Double-    -- | After how many iterations over the entire dataset-    -- the gain parameter is halved-    , tau       :: Double }-+import           GHC.Generics (Generic) --- | Default SGD parameter values.-sgdArgsDefault :: SgdArgs-sgdArgsDefault = SgdArgs-    { batchSize = 50-    , regVar    = 10-    , iterNum   = 10-    , gain0     = 0.25-      -- Without momentum I would rather go for '1', but with momentum the-      -- gradient gets significantly larger.-    , tau       = 5 }+import           Data.Default +import qualified Pipes as P --- | The gamma parameter which drives momentum.------ TODO: put in SgdArgs.----gamma :: Double-gamma = 0.9+import           Numeric.SGD.Type+import           Numeric.SGD.ParamSet  --- | Vector of parameters.-type Para       = U.Vector Double-+-- | Momentum configuration+data Config = Config+  { gain0 :: Double+  -- ^ Initial gain parameter, used to scale the gradient+  , tau :: Double+  -- ^ After how many gradient calculations the gain parameter is halved+  , gamma :: Double+  -- ^ Momentum term+  } deriving (Show, Eq, Ord, Generic) --- | Type synonym for mutable vector with Double values.-type MVect      = UM.MVector (Prim.PrimState IO) Double+instance Default Config where+  def = Config+    { gain0 = 0.01+    , gamma = 0.9+    , tau = 1000+    }  --- | A stochastic gradient descent method.--- A notification function can be used to provide user with--- information about the progress of the learning.-sgd-    :: SgdArgs                  -- ^ SGD parameter values-    -> (Para -> Int -> IO ())   -- ^ Notification run every update-    -> (Para -> x -> Grad)      -- ^ Gradient for dataset element-    -> Dataset x                -- ^ Dataset-    -> Para                     -- ^ Starting point-    -> IO Para                  -- ^ SGD result-sgd SgdArgs{..} notify mkGrad dataset x0 = do--  putStrLn $ "Running momentum!"--  -- A vector for the momentum gradient-  momentum <- UM.new (U.length x0)--  -- A worker vector for computing the actual gradients-  u <- UM.new (U.length x0)+-- | Stochastic gradient descent with momentum. See "Numeric.SGD.Momentum" for+-- more information.+momentum+  :: (Monad m, ParamSet p)+  => Config+    -- ^ Momentum configuration+  -> (e -> p -> p)+    -- ^ Gradient on a training element+  -> SGD m e p+momentum Config{..} gradient net0 = -  doIt momentum u 0 (R.mkStdGen 0) =<< U.thaw x0+  go (0 :: Integer) (zero net0) net0    where-    -- Gain in k-th iteration.-    gain k = (gain0 * tau) / (tau + done k) -    -- Number of completed iterations over the full dataset.-    done :: Int -> Double-    done k-        = fromIntegral (k * batchSize)-        / fromIntegral (size dataset)-    doneTotal :: Int -> Int-    doneTotal = floor . done--    -- Regularization (Guassian prior) parameter-    regularizationParam = regCoef-      where-        regCoef = iVar ** coef-        iVar = 1.0 / regVar-        coef = fromIntegral (size dataset)-             / fromIntegral batchSize--    doIt momentum u k stdGen x--      | done k > iterNum = do-        frozen <- U.unsafeFreeze x-        notify frozen k-        return frozen--      | otherwise = do--        -- Sample the dataset-        (batch, stdGen') <- sample stdGen batchSize dataset--        -- NEW: comment out-        -- -- Apply regularization to the parameters vector.-        -- scale (regularization k) x--        -- Freeze mutable vector of parameters. The frozen version is-        -- then supplied to external mkGrad function provided by user.-        frozen <- U.unsafeFreeze x-        notify frozen k--        -- Compute the gradient and put it in `u`-        let grad = parUnions (map (mkGrad frozen) batch)-        addUp grad u--        -- Apply regularization to `u`-        applyRegularization regularizationParam x u--        -- Scale the gradient-        scale (gain k) u--        -- Compute the new momentum-        updateMomentum gamma momentum u--        x' <- U.unsafeThaw frozen-        momentum `addTo` x'-        doIt momentum u (k+1) stdGen' x'----- | Compute the new momentum (gradient) vector.-applyRegularization-  :: Double -- ^ Regularization parameter-  -> MVect  -- ^ The parameters-  -> MVect  -- ^ The current gradient-  -> IO ()-applyRegularization regParam params grad = do-  forM_ [0 .. UM.length grad - 1] $ \i -> do-    x <- UM.unsafeRead grad i-    y <- UM.unsafeRead params i-    UM.unsafeWrite grad i $ x - regParam * y----- | Compute the new momentum (gradient) vector.-updateMomentum-  :: Double -- ^ The gamma parameter-  -> MVect  -- ^ The previous momentum-  -> MVect  -- ^ The scaled current gradient-  -> IO ()-updateMomentum gammaCoef momentum grad = do-  forM_ [0 .. UM.length momentum - 1] $ \i -> do-    x <- UM.unsafeRead momentum i-    y <- UM.unsafeRead grad i-    UM.unsafeWrite momentum i (gammaCoef * x + y)-+    -- Gain in the k-th iteration+    gain k+      = (gain0 * tau)+      / (tau + fromIntegral k) --- | Add up all gradients and store results in normal domain.-addUp :: Grad -> MVect -> IO ()-addUp grad v = do-    UM.set v 0-    forM_ (toList grad) $ \(i, x) -> do-        y <- UM.unsafeRead v i-        UM.unsafeWrite v i (x + y)+    go k moment net = do+      x <- P.await+      let grad = scale (gain k) (gradient x net)+          moment' = scale gamma moment `add` grad+          newNet = net `sub` moment'+      P.yield newNet+      go (k+1) moment' newNet  --- | Scale the vector by the given value.-scale :: Double -> MVect -> IO ()-scale c v = do-  forM_ [0 .. UM.length v - 1] $ \i -> do-    y <- UM.unsafeRead v i-    UM.unsafeWrite v i (c * y)+-------------------------------+-- Utils+-------------------------------  --- | Apply gradient to the parameters vector, that is add the first vector to--- the second one.-addTo :: MVect -> MVect -> IO ()-addTo w v = do-  forM_ [0 .. UM.length v - 1] $ \i -> do-    x <- UM.unsafeRead v i-    y <- UM.unsafeRead w i-    UM.unsafeWrite v i (x + y)+-- | Scaling+scale :: ParamSet p => Double -> p -> p+scale x = pmap (*x)+{-# INLINE scale #-}

@@ -0,0 +1,465 @@+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE EmptyCase #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveGeneric #-}+++module Numeric.SGD.ParamSet+  ( ParamSet(..)+  ) where+++import           GHC.Generics+import           GHC.TypeNats (KnownNat)++import           Prelude hiding (div)++import qualified Data.Map.Strict as M++import qualified Numeric.LinearAlgebra.Static as LA+++-- | Class of types that can be treated as parameter sets.  It provides basic+-- element-wise operations (addition, multiplication, mapping) which are+-- required to perform stochastic gradient descent.  Many of the operations+-- (`add`, `mul`, `sub`, `div`, etc.) have the same interpretation and follow+-- the same laws (e.g. associativity) as the corresponding operations in `Num`+-- and `Fractional`.  +-- +-- `zero` takes a parameter set as argument and "zero out"'s all its elements+-- (as in the backprop library).  This allows instances for `Maybe`, `M.Map`,+-- etc., where the structure of the parameter set is dynamic.  This leads to+-- the following property:+--+--     @add (zero x) x = x@+--+-- However, `zero` does not have to obey @(add (zero x) y = y)@.+--+-- A `ParamSet` can be also seen as a (structured) vector, hence `pmap` and+-- `norm_2`.  The latter is not strictly necessary to perform SGD, but it is+-- useful to control the training process.+--+-- `pmap` should obey the following law:+--+--     @pmap id x = x@+--+-- If you leave the body of an instance declaration blank, GHC Generics will be+-- used to derive instances if the type has a single constructor and each field+-- is an instance of `ParamSet`.+class ParamSet a where+  -- | Element-wise mapping+  pmap :: (Double -> Double) -> a -> a++  -- | Zero-out all elements+  zero :: a -> a+  zero = pmap (const 0.0)++--   -- | Element-wise negation+--   neg :: a -> a+--   neg = pmap (\x -> -x)++  -- | Element-wise addition+  add :: a -> a -> a+  -- | Elementi-wise substruction+  sub :: a -> a -> a++  -- | Element-wise multiplication+  mul :: a -> a -> a+  -- | Element-wise division+  div :: a -> a -> a++  -- | L2 norm+  norm_2 :: a -> Double++--   default zero :: (Generic a, GZero (Rep a)) => a -> a+--   zero = genericZero+--   {-# INLINE zero #-}++  default pmap+    :: (Generic a, GPMap (Rep a))+    => (Double -> Double) -> a -> a+  pmap = genericPMap+  {-# INLINE pmap #-}++  default add :: (Generic a, GAdd (Rep a)) => a -> a -> a+  add = genericAdd+  {-# INLINE add #-}++  default sub :: (Generic a, GSub (Rep a)) => a -> a -> a+  sub = genericSub+  {-# INLINE sub #-}++  default mul :: (Generic a, GMul (Rep a)) => a -> a -> a+  mul = genericMul+  {-# INLINE mul #-}++  default div :: (Generic a, GDiv (Rep a)) => a -> a -> a+  div = genericDiv+  {-# INLINE div #-}++  default norm_2 :: (Generic a, GNorm2 (Rep a)) => a -> Double+  norm_2 = genericNorm2+  {-# INLINE norm_2 #-}+++-- -- | 'add' using GHC Generics; works if all fields are instances of+-- -- 'ParamSet', but only for values with single constructors.+-- genericZero :: (Generic a, GZero (Rep a)) => a -> a+-- genericZero x = to $ gzero (from x)+-- {-# INLINE genericZero #-}+++-- | 'add' using GHC Generics; works if all fields are instances of+-- 'ParamSet', but only for values with single constructors.+genericAdd :: (Generic a, GAdd (Rep a)) => a -> a -> a+genericAdd x y = to $ gadd (from x) (from y)+{-# INLINE genericAdd #-}+++-- | 'sub' using GHC Generics; works if all fields are instances of+-- 'ParamSet', but only for values with single constructors.+genericSub :: (Generic a, GSub (Rep a)) => a -> a -> a+genericSub x y = to $ gsub (from x) (from y)+{-# INLINE genericSub #-}+++-- | 'div' using GHC Generics; works if all fields are instances of+-- 'ParamSet', but only for values with single constructors.+genericDiv :: (Generic a, GDiv (Rep a)) => a -> a -> a+genericDiv x y = to $ gdiv (from x) (from y)+{-# INLINE genericDiv #-}+++-- | 'mul' using GHC Generics; works if all fields are instances of+-- 'ParamSet', but only for values with single constructors.+genericMul :: (Generic a, GMul (Rep a)) => a -> a -> a+genericMul x y = to $ gmul (from x) (from y)+{-# INLINE genericMul #-}+++-- | 'norm_2' using GHC Generics; works if all fields are instances of+-- 'ParamSet', but only for values with single constructors.+genericNorm2 :: (Generic a, GNorm2 (Rep a)) => a -> Double+genericNorm2 x = gnorm_2 (from x)+{-# INLINE genericNorm2 #-}+++-- | 'pmap' using GHC Generics; works if all fields are instances of+-- 'ParamSet', but only for values with single constructors.+genericPMap :: (Generic a, GPMap (Rep a)) => (Double -> Double) -> a -> a+genericPMap f x = to $ gpmap f (from x)+{-# INLINE genericPMap #-}+++--------------------------------------------------+-- Generics+--+-- Partially borrowed from the backprop library+--------------------------------------------------+++-- -- | Helper class for automatically deriving 'add' using GHC Generics.+-- class GZero f where+--     gzero :: f t -> f t+-- +-- instance ParamSet p => GZero (K1 i p) where+--     gzero (K1 x) = K1 (zero x)+--     {-# INLINE gzero #-}+-- +-- instance (GZero f, GZero g) => GZero (f :*: g) where+--     gzero (x1 :*: y1) = x2 :*: y2+--       where+--         !x2 = gzero x1+--         !y2 = gzero y1+--     {-# INLINE gzero #-}+-- +-- instance GZero V1 where+--     gzero = \case {}+--     {-# INLINE gzero #-}+-- +-- instance GZero U1 where+--     gzero _ = U1+--     {-# INLINE gzero #-}+-- +-- instance GZero f => GZero (M1 i c f) where+--     gzero (M1 x) = M1 (gzero x)+--     {-# INLINE gzero #-}+-- +-- -- instance GZero f => GZero (f :.: g) where+-- --     gzero = Comp1 gzero+-- --     {-# INLINE gzero #-}+++-- | Helper class for automatically deriving 'add' using GHC Generics.+class GAdd f where+    gadd :: f t -> f t -> f t++instance ParamSet a => GAdd (K1 i a) where+    gadd (K1 x) (K1 y) = K1 (add x y)+    {-# INLINE gadd #-}++instance (GAdd f, GAdd g) => GAdd (f :*: g) where+    gadd (x1 :*: y1) (x2 :*: y2) = x3 :*: y3+      where+        !x3 = gadd x1 x2+        !y3 = gadd y1 y2+    {-# INLINE gadd #-}++instance GAdd V1 where+    gadd = \case {}+    {-# INLINE gadd #-}++instance GAdd U1 where+    gadd _ _ = U1+    {-# INLINE gadd #-}++instance GAdd f => GAdd (M1 i c f) where+    gadd (M1 x) (M1 y) = M1 (gadd x y)+    {-# INLINE gadd #-}++-- instance GAdd f => GAdd (f :.: g) where+--     gadd (Comp1 x) (Comp1 y) = Comp1 (gadd x y)+--     {-# INLINE gadd #-}+++-- | Helper class for automatically deriving 'sub' using GHC Generics.+class GSub f where+    gsub :: f t -> f t -> f t++instance ParamSet a => GSub (K1 i a) where+    gsub (K1 x) (K1 y) = K1 (sub x y)+    {-# INLINE gsub #-}++instance (GSub f, GSub g) => GSub (f :*: g) where+    gsub (x1 :*: y1) (x2 :*: y2) = x3 :*: y3+      where+        !x3 = gsub x1 x2+        !y3 = gsub y1 y2+    {-# INLINE gsub #-}++instance GSub V1 where+    gsub = \case {}+    {-# INLINE gsub #-}++instance GSub U1 where+    gsub _ _ = U1+    {-# INLINE gsub #-}++instance GSub f => GSub (M1 i c f) where+    gsub (M1 x) (M1 y) = M1 (gsub x y)+    {-# INLINE gsub #-}++-- instance GSub f => GSub (f :.: g) where+--     gsub (Comp1 x) (Comp1 y) = Comp1 (gsub x y)+--     {-# INLINE gsub #-}+++-- | Helper class for automatically deriving 'mul' using GHC Generics.+class GMul f where+    gmul :: f t -> f t -> f t++instance ParamSet a => GMul (K1 i a) where+    gmul (K1 x) (K1 y) = K1 (mul x y)+    {-# INLINE gmul #-}++instance (GMul f, GMul g) => GMul (f :*: g) where+    gmul (x1 :*: y1) (x2 :*: y2) = x3 :*: y3+      where+        !x3 = gmul x1 x2+        !y3 = gmul y1 y2+    {-# INLINE gmul #-}++instance GMul V1 where+    gmul = \case {}+    {-# INLINE gmul #-}++instance GMul U1 where+    gmul _ _ = U1+    {-# INLINE gmul #-}++instance GMul f => GMul (M1 i c f) where+    gmul (M1 x) (M1 y) = M1 (gmul x y)+    {-# INLINE gmul #-}++-- instance GMul f => GMul (f :.: g) where+--     gmul (Comp1 x) (Comp1 y) = Comp1 (gmul x y)+--     {-# INLINE gmul #-}+++-- | Helper class for automatically deriving 'div' using GHC Generics.+class GDiv f where+    gdiv :: f t -> f t -> f t++instance ParamSet a => GDiv (K1 i a) where+    gdiv (K1 x) (K1 y) = K1 (div x y)+    {-# INLINE gdiv #-}++instance (GDiv f, GDiv g) => GDiv (f :*: g) where+    gdiv (x1 :*: y1) (x2 :*: y2) = x3 :*: y3+      where+        !x3 = gdiv x1 x2+        !y3 = gdiv y1 y2+    {-# INLINE gdiv #-}++instance GDiv V1 where+    gdiv = \case {}+    {-# INLINE gdiv #-}++instance GDiv U1 where+    gdiv _ _ = U1+    {-# INLINE gdiv #-}++instance GDiv f => GDiv (M1 i c f) where+    gdiv (M1 x) (M1 y) = M1 (gdiv x y)+    {-# INLINE gdiv #-}++-- instance GDiv f => GDiv (f :.: g) where+--     gdiv (Comp1 x) (Comp1 y) = Comp1 (gdiv x y)+--     {-# INLINE gdiv #-}+++-- | Helper class for automatically deriving 'norm_2' using GHC Generics.+class GNorm2 f where+    gnorm_2 :: f t -> Double++instance ParamSet a => GNorm2 (K1 i a) where+    gnorm_2 (K1 x) = norm_2 x+    {-# INLINE gnorm_2 #-}++instance (GNorm2 f, GNorm2 g) => GNorm2 (f :*: g) where+    gnorm_2 (x1 :*: y1) =+      sqrt ((x2 ^ (2 :: Int)) + (y2 ^ (2 :: Int)))+      where+        !x2 = gnorm_2 x1+        !y2 = gnorm_2 y1+    {-# INLINE gnorm_2 #-}++instance GNorm2 V1 where+    gnorm_2 = \case {}+    {-# INLINE gnorm_2 #-}++instance GNorm2 U1 where+    gnorm_2 _ = 0+    {-# INLINE gnorm_2 #-}++instance GNorm2 f => GNorm2 (M1 i c f) where+    gnorm_2 (M1 x) = gnorm_2 x+    {-# INLINE gnorm_2 #-}++-- -- TODO: Make sure this makes sense+-- instance GNorm2 f => GNorm2 (f :.: g) where+--     gnorm_2 (Comp1 x) = gnorm_2 x+--     {-# INLINE gnorm_2 #-}+++-- | Helper class for automatically deriving 'pmap' using GHC Generics.+class GPMap f where+    gpmap :: (Double -> Double) -> f t -> f t++instance ParamSet a => GPMap (K1 i a) where+    gpmap f (K1 x) = K1 (pmap f x)+    {-# INLINE gpmap #-}++instance (GPMap f, GPMap g) => GPMap (f :*: g) where+    gpmap f (x1 :*: y1) = x2 :*: y2+      where+        !x2 = gpmap f x1+        !y2 = gpmap f y1+    {-# INLINE gpmap #-}++instance GPMap V1 where+    gpmap _ = \case {}+    {-# INLINE gpmap #-}++instance GPMap U1 where+    gpmap _ _ = U1+    {-# INLINE gpmap #-}++instance GPMap f => GPMap (M1 i c f) where+    gpmap f (M1 x) = M1 (gpmap f x)+    {-# INLINE gpmap #-}++-- instance GPMap f => GPMap (f :.: g) where+--     gpmap f (Comp1 x) = Comp1 (gpmap f x)+--     {-# INLINE gpmap #-}+++--------------------------------------------------+-- Basic instances+--------------------------------------------------+++instance ParamSet Double where+  zero = const 0+  pmap = id+  add = (+)+  sub = (-)+  mul = (*)+  div = (/)+  norm_2 = abs+++instance (KnownNat n) => ParamSet (LA.R n) where+  zero = const 0+  pmap = LA.dvmap+  add = (+)+  sub = (-)+  mul = (*)+  div = (/)+  norm_2 = LA.norm_2+++instance (KnownNat n, KnownNat m) => ParamSet (LA.L n m) where+  zero = const 0+  pmap = LA.dmmap+  add = (+)+  sub = (-)+  mul = (*)+  div = (/)+  norm_2 = LA.norm_2+++-- | `Nothing` represents a deactivated parameter set component. If `Nothing`+-- is given as an argument to one of the `ParamSet` operations, the result is+-- `Nothing` as well.+--+-- This differs from the corresponding instance in the backprop library, where+-- `Nothing` is equivalent to `Just 0`.  However, the implementation below+-- seems to correspond adequately enough to the notion that a particular+-- component is either active or not in both the parameter set and the+-- gradient, hence it doesn't make sense to combine `Just` with `Nothing`.+instance (ParamSet a) => ParamSet (Maybe a) where+  zero = fmap zero+  pmap = fmap . pmap++  add (Just x) (Just y) = Just (add x y)+  add _ _ = Nothing++  sub (Just x) (Just y) = Just (sub x y)+  sub _ _ = Nothing++  mul (Just x) (Just y) = Just (mul x y)+  mul _ _ = Nothing++  div (Just x) (Just y) = Just (div x y)+  div _ _ = Nothing++  norm_2 = maybe 0 norm_2+++-- | A map with different parameter sets (of the same type) assigned to the+-- individual keys.+--+-- When combining two maps with different sets of keys, only their intersection+-- is preserved.+instance (Ord k, ParamSet a) => ParamSet (M.Map k a) where+  zero = fmap zero+  pmap f = fmap (pmap f)+  add = M.intersectionWith add+  sub = M.intersectionWith sub+  mul= M.intersectionWith mul+  div= M.intersectionWith div+  norm_2 = sqrt . sum . map ((^(2::Int)) . norm_2)  . M.elems

@@ -0,0 +1,167 @@+{-# LANGUAGE RecordWildCards #-}+++-- | Stochastic gradient descent using mutable vectors for efficient parameter+-- update.  This module is intended for use with sparse features.  If you use+-- dense feature vectors (as arise e.g. in deep learning), have a look at+-- "Numeric.SGD".+-- +-- Currently only the Gaussian regularization is implemented.+--+-- SGD with momentum is known to converge faster than vanilla SGD.  It's+-- implementation can be found in "Numeric.SGD.Sparse.Momentum".+++module Numeric.SGD.Sparse+( SgdArgs (..)+, sgdArgsDefault+, Para+, sgd+, module Numeric.SGD.Sparse.Grad+, module Numeric.SGD.DataSet+) where+++import           Control.Monad (forM_)+-- import qualified System.Random as R+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM+import qualified Control.Monad.Primitive as Prim++import           Numeric.SGD.Sparse.Grad+import           Numeric.SGD.DataSet+++-- | SGD parameters controlling the learning process.+data SgdArgs = SgdArgs+    { -- | Size of the batch+      batchSize :: Int+    -- | Regularization variance+    , regVar    :: Double+    -- | Number of iterations+    , iterNum   :: Double+    -- | Initial gain parameter+    , gain0     :: Double+    -- | After how many iterations over the entire dataset+    -- the gain parameter is halved+    , tau       :: Double }+++-- | Default SGD parameter values.+sgdArgsDefault :: SgdArgs+sgdArgsDefault = SgdArgs+    { batchSize = 30+    , regVar    = 10+    , iterNum   = 10+    , gain0     = 1+    , tau       = 5 }+++-- | Vector of parameters.+type Para       = U.Vector Double+++-- | Type synonym for mutable vector with Double values.+type MVect      = UM.MVector (Prim.PrimState IO) Double+++-- | A stochastic gradient descent method.+-- A notification function can be used to provide user with+-- information about the progress of the learning.+sgd+    :: SgdArgs                  -- ^ SGD parameter values+    -> (Para -> Int -> IO ())   -- ^ Notification run every update+    -> (Para -> x -> Grad)      -- ^ Gradient for dataset element+    -> DataSet x                -- ^ Dataset+    -> Para                     -- ^ Starting point+    -> IO Para                  -- ^ SGD result+sgd SgdArgs{..} notify mkGrad dataset x0 = do+  u <- UM.new (U.length x0)+  -- doIt u 0 (R.mkStdGen 0) =<< U.thaw x0+  doIt u 0 =<< U.thaw x0+  where+    -- Gain in k-th iteration.+    gain k = (gain0 * tau) / (tau + done k)++    -- Number of completed iterations over the full dataset.+    done :: Int -> Double+    done k+        = fromIntegral (k * batchSize)+        / fromIntegral (size dataset)+    -- doneTotal :: Int -> Int+    -- doneTotal = floor . done++    -- Regularization (Guassian prior)+    regularization k = regCoef+      where+        regCoef = (1.0 - gain k * iVar) ** coef+        iVar = 1.0 / regVar+        coef = fromIntegral batchSize+             / fromIntegral (size dataset)++--     -- Regularization (Guassian prior) after a full dataset pass+--     regularization k = 1.0 - (gain k / regVar)++    doIt u k x+      | done k > iterNum = do+        frozen <- U.unsafeFreeze x+        notify frozen k+        return frozen+      | otherwise = do+        -- (batch, stdGen') <- sample stdGen batchSize dataset+        batch <- randomSample batchSize dataset++        -- Regularization+        -- when (doneTotal (k - 1) /= doneTotal k) $ do+        --   <- we now apply regularization each step rather than each+        --      dataset pass+        let regParam = regularization k+        -- putStrLn $ "\nApplying regularization (params *= " ++ show regParam ++ ")"+        scale regParam x++--         -- Regularization+--         when (doneTotal (k - 1) /= doneTotal k) $ do+--           let regParam = regularization k+--           putStrLn $ "\nApplying regularization (params *= " ++ show regParam ++ ")"+--           scale regParam x++        -- Freeze mutable vector of parameters. The frozen version is+        -- then supplied to external mkGrad function provided by user.+        frozen <- U.unsafeFreeze x+        notify frozen k++        -- let grad = M.unions (map (mkGrad frozen) batch)+        let grad = parUnions (map (mkGrad frozen) batch)+        addUp grad u+        scale (gain k) u++        x' <- U.unsafeThaw frozen+        u `addTo` x'+        doIt u (k+1) x'+++-- | Add up all gradients and store results in normal domain.+addUp :: Grad -> MVect -> IO ()+addUp grad v = do+    UM.set v 0+    forM_ (toList grad) $ \(i, x) -> do+        y <- UM.unsafeRead v i+        UM.unsafeWrite v i (x + y)+++-- | Scale the vector by the given value.+scale :: Double -> MVect -> IO ()+scale c v = do+    forM_ [0 .. UM.length v - 1] $ \i -> do+        y <- UM.unsafeRead v i+        UM.unsafeWrite v i (c * y)+++-- | Apply gradient to the parameters vector, that is add the first vector to+-- the second one.+addTo :: MVect -> MVect -> IO ()+addTo w v = do+    forM_ [0 .. UM.length v - 1] $ \i -> do+        x <- UM.unsafeRead v i+        y <- UM.unsafeRead w i+        UM.unsafeWrite v i (x + y)

@@ -0,0 +1,131 @@+{-# LANGUAGE CPP #-}++-- | A gradient is represented by an IntMap from gradient indices to values.+-- Elements with no associated values in the gradient are assumed to have a 0+-- value assigned.  Such elements are of no interest: when adding the gradient+-- to the vector of parameters, only non-zero elements are taken into account.+-- +-- Each value associated with a gradient position is a pair of positive and+-- negative components.  They are stored separately to ensure high accuracy of+-- computation results.  Besides, both positive and negative components are+-- stored in a logarithmic domain.++module Numeric.SGD.Sparse.Grad+( Grad+, empty+, add+, addL+, fromList+, fromLogList+, toList+, parUnions+) where++import Data.List (foldl')+import Control.Applicative ((<$>), (<*>))+import Control.Monad.Par (Par, runPar, get)+#if MIN_VERSION_containers(0,4,2)+import Control.Monad.Par (spawn)+#else+import Control.DeepSeq (deepseq)+import Control.Monad.Par (spawn_)+#endif+#if MIN_VERSION_containers(0,5,0)+import qualified Data.IntMap.Strict as M+#else+import qualified Data.IntMap as M+#endif++import Numeric.SGD.Sparse.LogSigned++-- | Gradient with nonzero values stored in a logarithmic domain.+-- Since values equal to zero have no impact on the update phase+-- of the SGD method, it is more efficient to not to store those+-- components in the gradient.+type Grad = M.IntMap LogSigned++{-# INLINE insertWith #-}+insertWith :: (a -> a -> a) -> M.Key -> a -> M.IntMap a -> M.IntMap a+#if MIN_VERSION_containers(0,5,0)+insertWith = M.insertWith+#elif MIN_VERSION_containers(0,4,1)+insertWith = M.insertWith'+#else+insertWith f k x m = +    M.alter g k m+  where+    g my = case my of+        Nothing -> Just x+        Just y  ->+            let z = f x y+            in  z `seq` Just z+#endif++-- | Add normal-domain double to the gradient at the given position.+{-# INLINE add #-}+add :: Grad -> Int -> Double -> Grad+add grad i y = insertWith (+) i (logSigned y) grad +++-- | Add log-domain, singed number to the gradient at the given position.+{-# INLINE addL #-}+addL :: Grad -> Int -> LogSigned -> Grad+addL grad i y = insertWith (+) i y grad ++-- | Construct gradient from a list of (index, value) pairs.+-- All values from the list are added at respective gradient+-- positions.+{-# INLINE fromList #-}+fromList :: [(Int, Double)] -> Grad+fromList =+    let ins grad (i, y) = add grad i y+    in  foldl' ins empty++-- | Construct gradient from a list of (index, signed, log-domain number)+-- pairs.  All values from the list are added at respective gradient+-- positions.+{-# INLINE fromLogList #-}+fromLogList :: [(Int, LogSigned)] -> Grad+fromLogList =+    let ins grad (i, y) = addL grad i y+    in  foldl' ins empty++-- | Collect gradient components with values in normal domain.+{-# INLINE toList #-}+toList :: Grad -> [(Int, Double)]+toList =+    let unLog (i, x) = (i, toNorm x)+    in  map unLog . M.assocs++-- | Empty gradient, i.e. with all elements set to 0.+{-# INLINE empty #-}+empty :: Grad+empty = M.empty++-- | Perform parallel unions operation on gradient list. +-- Experimental version.+parUnions :: [Grad] -> Grad+parUnions [] = error "parUnions: empty list"+parUnions xs = runPar (parUnionsP xs)++-- | Parallel unions in the Par monad.+parUnionsP :: [Grad] -> Par Grad+parUnionsP [x] = return x+parUnionsP zs  = do+    let (xs, ys) = split zs+#if MIN_VERSION_containers(0,4,2)+    xsP <- spawn (parUnionsP xs)+    ysP <- spawn (parUnionsP ys)+    M.unionWith (+) <$> get xsP <*> get ysP+#else+    xsP <- spawn_ (parUnionsP xs)+    ysP <- spawn_ (parUnionsP ys)+    x <- M.unionWith (+) <$> get xsP <*> get ysP+    M.elems x `deepseq` return x+#endif+  where+    split []        = ([], [])+    split (x:[])    = ([x], [])+    split (x:y:rest)  =+        let (xs, ys) = split rest+        in  (x:xs, y:ys)

@@ -0,0 +1,85 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++-- | Module provides data type for signed log-domain calculations.++module Numeric.SGD.Sparse.LogSigned+( LogSigned (..)+, logSigned+, fromPos+, fromNeg+, toNorm+, toLogFloat+) where++import qualified Data.Number.LogFloat as L+import Data.Function (on)+import Control.DeepSeq (NFData(..))++-- | Signed real value in the logarithmic domain.+data LogSigned = LogSigned+    { pos :: {-# UNPACK #-} !L.LogFloat     -- ^ Positive component+    , neg :: {-# UNPACK #-} !L.LogFloat     -- ^ Negative component+    } deriving Show++instance Eq LogSigned where+    (==) = (==) `on` toLogFloat++instance Ord LogSigned where+    compare = compare `on` toLogFloat++-- All fields are strict and unpacked.+instance NFData LogSigned where+    rnf (LogSigned p q) = p `seq` q `seq` ()++-- | Smart LogSigned constructor.+{-# INLINE logSigned #-}+logSigned :: Double -> LogSigned+logSigned x+    | x > 0     = LogSigned (L.logFloat x) zero+    | x < 0     = LogSigned zero (L.logFloat (-x))+    | otherwise = LogSigned zero zero++-- | Make LogSigned from a positive, log-domain number.+{-# INLINE fromPos #-}+fromPos :: L.LogFloat -> LogSigned+fromPos x = LogSigned x zero++-- | Make LogSigned from a negative, log-domain number.+{-# INLINE fromNeg #-}+fromNeg :: L.LogFloat -> LogSigned+fromNeg x = LogSigned zero x++-- | Shift LogSigned to a normal domain.+{-# INLINE toNorm #-}+toNorm :: LogSigned -> Double+toNorm (LogSigned x y) = L.fromLogFloat x - L.fromLogFloat y++-- | Change the 'LogSigned' to either negative 'Left' 'L.LogFloat'+-- or positive 'Right' 'L.LogFloat'.+toLogFloat :: LogSigned -> Either L.LogFloat L.LogFloat+toLogFloat x = case signum x of+    -1  -> Left  $ neg x - pos x+    1   -> Right $ pos x - neg x+    _   -> Right $ L.logFloat (0 :: Double)++instance Num LogSigned where+    LogSigned x y + LogSigned x' y' =+        LogSigned (x + x') (y + y')+    LogSigned x y * LogSigned x' y' =+        LogSigned (x*x' + y*y') (x*y' + y*x')+    LogSigned x y - LogSigned x' y' =+        LogSigned (x + y') (y + x')+    negate  (LogSigned x y) = LogSigned y x+    abs     (LogSigned x y)+        | x >= y    = LogSigned x y+        | otherwise = LogSigned y x+    signum (LogSigned x y)+        | x > y     =  1+        | x < y     = -1+        | otherwise =  0+    fromInteger = logSigned . fromInteger++{-# INLINE zero #-}+zero :: L.LogFloat+zero = L.logFloat (0 :: Double)

@@ -0,0 +1,201 @@+{-# LANGUAGE RecordWildCards #-}+++-- | A version of `Numeric.SGD.Sparse` extended with momentum.+++module Numeric.SGD.Sparse.Momentum+( SgdArgs (..)+, sgdArgsDefault+, Para+, sgd+, module Numeric.SGD.Sparse.Grad+, module Numeric.SGD.DataSet+) where+++import           Control.Monad (forM_)+-- import qualified System.Random as R+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM+import qualified Control.Monad.Primitive as Prim++import           Numeric.SGD.Sparse.Grad+import           Numeric.SGD.DataSet+++-- | SGD parameters controlling the learning process.+data SgdArgs = SgdArgs+    { -- | Size of the batch+      batchSize :: Int+    -- | Regularization variance+    , regVar    :: Double+    -- | Number of iterations+    , iterNum   :: Double+    -- | Initial gain parameter+    , gain0     :: Double+    -- | After how many iterations over the entire dataset+    -- the gain parameter is halved+    , tau       :: Double +    }+++-- | Default SGD parameter values.+sgdArgsDefault :: SgdArgs+sgdArgsDefault = SgdArgs+    { batchSize = 50+    , regVar    = 10+    , iterNum   = 10+    , gain0     = 0.25+      -- Without momentum I would rather go for '1', but with momentum the+      -- gradient gets significantly larger.+    , tau       = 5 +    }+++-- | The gamma parameter which drives momentum.+--+-- TODO: put in SgdArgs.+--+gamma :: Double+gamma = 0.9+++-- | Vector of parameters.+type Para       = U.Vector Double+++-- | Type synonym for mutable vector with Double values.+type MVect      = UM.MVector (Prim.PrimState IO) Double+++-- | A stochastic gradient descent method.+-- A notification function can be used to provide user with+-- information about the progress of the learning.+sgd+    :: SgdArgs                  -- ^ SGD parameter values+    -> (Para -> Int -> IO ())   -- ^ Notification run every update+    -> (Para -> x -> Grad)      -- ^ Gradient for dataset element+    -> DataSet x                -- ^ DataSet+    -> Para                     -- ^ Starting point+    -> IO Para                  -- ^ SGD result+sgd SgdArgs{..} notify mkGrad dataset x0 = do++  putStrLn $ "Running momentum!"++  -- A vector for the momentum gradient+  momentum <- UM.new (U.length x0)++  -- A worker vector for computing the actual gradients+  u <- UM.new (U.length x0)++  doIt momentum u 0 =<< U.thaw x0++  where+    -- Gain in k-th iteration.+    gain k = (gain0 * tau) / (tau + done k)++    -- Number of completed iterations over the full dataset.+    done :: Int -> Double+    done k+        = fromIntegral (k * batchSize)+        / fromIntegral (size dataset)++    -- Regularization (Guassian prior) parameter+    regularizationParam = regCoef+      where+        regCoef = iVar ** coef+        iVar = 1.0 / regVar+        coef = fromIntegral (size dataset)+             / fromIntegral batchSize++    doIt momentum u k x++      | done k > iterNum = do+        frozen <- U.unsafeFreeze x+        notify frozen k+        return frozen++      | otherwise = do++        -- Sample the dataset+        batch <- randomSample batchSize dataset++        -- NEW: comment out+        -- -- Apply regularization to the parameters vector.+        -- scale (regularization k) x++        -- Freeze mutable vector of parameters. The frozen version is+        -- then supplied to external mkGrad function provided by user.+        frozen <- U.unsafeFreeze x+        notify frozen k++        -- Compute the gradient and put it in `u`+        let grad = parUnions (map (mkGrad frozen) batch)+        addUp grad u++        -- Apply regularization to `u`+        applyRegularization regularizationParam x u++        -- Scale the gradient+        scale (gain k) u++        -- Compute the new momentum+        updateMomentum gamma momentum u++        x' <- U.unsafeThaw frozen+        momentum `addTo` x'+        doIt momentum u (k+1) x'+++-- | Compute the new momentum (gradient) vector.+applyRegularization+  :: Double -- ^ Regularization parameter+  -> MVect  -- ^ The parameters+  -> MVect  -- ^ The current gradient+  -> IO ()+applyRegularization regParam params grad = do+  forM_ [0 .. UM.length grad - 1] $ \i -> do+    x <- UM.unsafeRead grad i+    y <- UM.unsafeRead params i+    UM.unsafeWrite grad i $ x - regParam * y+++-- | Compute the new momentum (gradient) vector.+updateMomentum+  :: Double -- ^ The gamma parameter+  -> MVect  -- ^ The previous momentum+  -> MVect  -- ^ The scaled current gradient+  -> IO ()+updateMomentum gammaCoef momentum grad = do+  forM_ [0 .. UM.length momentum - 1] $ \i -> do+    x <- UM.unsafeRead momentum i+    y <- UM.unsafeRead grad i+    UM.unsafeWrite momentum i (gammaCoef * x + y)+++-- | Add up all gradients and store results in normal domain.+addUp :: Grad -> MVect -> IO ()+addUp grad v = do+    UM.set v 0+    forM_ (toList grad) $ \(i, x) -> do+        y <- UM.unsafeRead v i+        UM.unsafeWrite v i (x + y)+++-- | Scale the vector by the given value.+scale :: Double -> MVect -> IO ()+scale c v = do+  forM_ [0 .. UM.length v - 1] $ \i -> do+    y <- UM.unsafeRead v i+    UM.unsafeWrite v i (c * y)+++-- | Apply gradient to the parameters vector, that is add the first vector to+-- the second one.+addTo :: MVect -> MVect -> IO ()+addTo w v = do+  forM_ [0 .. UM.length v - 1] $ \i -> do+    x <- UM.unsafeRead v i+    y <- UM.unsafeRead w i+    UM.unsafeWrite v i (x + y)

@@ -0,0 +1,12 @@+module Numeric.SGD.Type+  ( SGD+  ) where+++import Pipes as P+++-- | SGD is a pipe which, given the initial parameter values, consumes training+-- elements of type `e` and outputs the subsequently calculated parameter sets+-- of type `p`.+type SGD m e p = p -> P.Pipe e p m ()