diff --git a/sgd.cabal b/sgd.cabal --- a/sgd.cabal +++ b/sgd.cabal @@ -1,5 +1,5 @@ name: sgd -version: 0.6.0.0 +version: 0.7.0.0 synopsis: Stochastic gradient descent description: Stochastic gradient descent library. @@ -25,6 +25,7 @@ , pipes >= 4.3 && < 4.4 , vector >= 0.10 && < 0.13 , random >= 1.0 && < 1.2 + , random-shuffle >= 0.0.4 && < 0.1 , primitive >= 0.5 && < 0.7 , logfloat >= 0.12 && < 0.14 , monad-par >= 0.3.4 && < 0.4 @@ -36,6 +37,7 @@ , temporary >= 1.1 && < 1.4 , hmatrix >= 0.19 && < 0.20 , data-default >= 0.7 && < 0.8 + , parallel >= 3.2 && < 3.3 exposed-modules: Numeric.SGD diff --git a/src/Numeric/SGD.hs b/src/Numeric/SGD.hs --- a/src/Numeric/SGD.hs +++ b/src/Numeric/SGD.hs @@ -61,10 +61,19 @@ , runIO -- * Combinators + -- ** Input , pipeSeq , pipeRan + -- ** Batch + , batch + , batchGradSeq + , batchGradPar + , batchGradPar' + -- ** Output , result - , every + -- ** Misc + , keepEvery + , decreasingBy -- * Re-exports , def @@ -72,13 +81,19 @@ import GHC.Generics (Generic) +-- import GHC.Conc (numCapabilities) import Numeric.Natural (Natural) -import qualified System.Random as R +-- import qualified System.Random as R -import Control.Monad (when, forM_) +import Control.Monad (when, forM_, forever) +import Control.Parallel.Strategies (parMap, rseq, rdeepseq, Strategy) +import Control.DeepSeq (NFData) +import qualified Control.Monad.State.Strict as State import Data.Functor.Identity (Identity(..)) +import Data.List (foldl1') -- , transpose) + import qualified Data.IORef as IO import Data.Default @@ -127,6 +142,10 @@ data Config = Config { iterNum :: Natural -- ^ Number of iteration over the entire training dataset + , batchSize :: Natural + -- ^ Mini-batch size + , batchOverlap :: Natural + -- ^ The number of overlapping elements in subsequent mini-batches , batchRandom :: Bool -- ^ Should the mini-batch be selected at random? If not, the subsequent -- training elements will be picked sequentially. Random selection gives @@ -139,6 +158,8 @@ instance Default Config where def = Config { iterNum = 100 + , batchSize = 1 + , batchOverlap = 0 , batchRandom = False , reportEvery = 1.0 } @@ -149,16 +170,16 @@ -- 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.). +-- SGD pipeline based on lower-level combinators (`pipeSeq`, `batch`, +-- `Adam.adam`, `keepEvery`, `result`, etc.). runIO :: (ParamSet p) => Config -- ^ SGD configuration - -> SGD IO e p - -- ^ Selected SGD method + -> SGD IO [e] p + -- ^ SGD pipe consuming mini-batches of dataset elements -> (e -> p -> Double) - -- ^ Value of the objective function on a sample element (needed for model + -- ^ Value of the objective function on a dataset element (used for model -- quality reporting) -> DataSet e -- ^ Training dataset @@ -166,21 +187,40 @@ -- ^ Initial parameter values -> IO p runIO Config{..} sgd quality0 dataSet net0 = do - report net0 - result net0 $ pipeSeq dataSet + _ <- report net0 + result net0 $ pipeData dataSet + >-> batch (fromIntegral batchSize) + >-> batchFilter >-> sgd net0 - >-> P.take realIterNum - >-> every realReportPeriod report + >-> keepEvery realReportPeriod + >-> P.take (fromIntegral iterNum) + >-> decreasingBy report where - -- Iteration scaling - iterScale x = fromIntegral (size dataSet) * x + -- Data streaming function + pipeData = forever . + if batchRandom + then pipeRan + else pipeSeq + -- Number of new elements in each subsequent batch + batchNew = max 1 + ( fromIntegral batchSize + - fromIntegral batchOverlap ) + -- Batch stream filter + batchFilter = do + P.await >>= P.yield + keepEvery batchNew + -- Iteration (epoch) scaling + iterScale x = fromIntegral (size dataSet) * x / fromIntegral batchNew -- Number of iterations and reporting period - realIterNum = ceiling $ iterScale (fromIntegral iterNum :: Double) + -- realIterNum = ceiling $ iterScale (fromIntegral iterNum :: Double) + -- realReportPeriod = floor $ iterScale reportEvery realReportPeriod = ceiling $ iterScale reportEvery -- Network quality over the entire training dataset report net = do - putStr . show =<< quality net + q <- quality net + putStr $ show q putStrLn $ " (norm_2 = " ++ show (norm_2 net) ++ ")" + return q quality net = do res <- IO.newIORef 0.0 forM_ [0 .. size dataSet - 1] $ \ix -> do @@ -194,29 +234,82 @@ ------------------------------- --- | Pipe the dataset sequentially in a loop. +-- | Pipe all the elements in the dataset sequentially. pipeSeq :: DataSet e -> P.Producer e IO () pipeSeq dataSet = do go (0 :: Int) where go k - | k >= size dataSet = go 0 + | k >= size dataSet = return () | otherwise = do x <- P.lift $ elemAt dataSet k P.yield x go (k+1) --- | Pipe the dataset randomly in a loop. +-- | Pipe all the elements in the dataset in a random order. 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 +pipeRan dataSet0 = do + dataSet <- P.lift $ shuffle dataSet0 + pipeSeq dataSet +-- | Group dataset elements into (mini-)batches of the given size. +batch :: (Monad m) => Int -> P.Pipe e [e] m () +batch k = flip State.evalStateT [] . forever $ do + x <- P.lift P.await + xs <- State.get + let xs' = take k (x:xs) + when (length xs' == k) $ do + P.lift (P.yield xs') + State.put xs' + + +-- | Adapt the gradient function to handle (mini-)batches. Relies on the @p@'s +-- `NFData` instance to efficiently calculate gradients in parallel. +batchGradPar + :: (ParamSet p, NFData p) + => (e -> p -> p) + -> ([e] -> p -> p) +batchGradPar = batchGradWith rdeepseq + + +-- | A version of `batchGradPar` with no `NFData` constraint. Evaluates the +-- sub-gradients calculated in parallel to weak head normal form. +batchGradPar' + :: (ParamSet p) + => (e -> p -> p) + -> ([e] -> p -> p) +batchGradPar' = batchGradWith rseq + + +-- | Adapt the gradient function to handle (mini-)batches. The sub-gradients +-- of the individual batch elements are evaluated in parallel based on the +-- given `Strategy`. +batchGradWith + :: (ParamSet p) + => Strategy p + -> (e -> p -> p) + -> ([e] -> p -> p) +batchGradWith strategy grad xs param = + case parMap strategy (\e -> grad e param) xs of + [] -> param + -- TODO: the fold is sequential, we could try to parallize it as well. + ps -> foldl1' add ps + + +-- | Adapt the gradient function to handle (mini-)batches. The function +-- calculates the individual sub-gradients sequentially. +batchGradSeq + :: (ParamSet p) + => (e -> p -> p) + -> ([e] -> p -> p) +batchGradSeq grad xs param = + case map (flip grad param) xs of + [] -> param + ps -> foldl1' add ps + + -- | Extract the result of the SGD calculation (the last parameter -- set flowing downstream). result @@ -229,14 +322,77 @@ 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) +-- -- | Apply the given monadic function to every @k@-th value 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 + + +-- | Keep every @k@-th element flowing downstream and discard all the others. +keepEvery :: (Monad m) => Int -> P.Pipe a a m x +keepEvery k = forever $ do + sequence_ $ replicate (k-1) P.await + P.await >>= P.yield +-- keepEvery k = do +-- go (1 `mod` k) +-- where +-- go i = do +-- x <- P.await +-- when (i == 0) $ do +-- P.yield x +-- go $ (i+1) `mod` k + + +-- -- | Keep the elements with the corresponding `True` in the argument list. +-- -- +-- -- TODO: (=) or (==) in the following example? And is this example correct? +-- -- @ +-- -- keep (forever True) = P.id +-- -- @ +-- keep :: (Monad m) => [Bool] -> P.Pipe a a m () +-- keep [] = return () +-- keep (b:bs) = do +-- x <- P.await +-- when b (P.yield x) +-- keep bs +-- +-- +-- -- | Create the mask to `keep` each @k@-th element flowing downstream. +-- every :: Int -> [Bool] +-- every k = cycle $ replicate (k-1) False ++ [True] + + +-- | Make the stream decreasing in the given (monadic) function by discarding +-- elements with values higher than those already seen. +decreasingBy :: (Monad m, Ord a) => (p -> m a) -> P.Pipe p p m x +decreasingBy f = do + x <- P.await + v <- P.lift (f x) + P.yield x + go v where - go i = do - paramSet <- P.await - when (i == 0) $ do - P.lift $ f paramSet - P.yield paramSet - go $ (i+1) `mod` k + go w = do + x <- P.await + v <- P.lift (f x) + when (v < w) (P.yield x) + go (min v w) + + +------------------------------- +-- Utils +------------------------------- + + +-- partition :: Int -> [a] -> [[a]] +-- partition n = +-- transpose . group n +-- where +-- group _ [] = [] +-- group k xs = take k xs : (group k $ drop k xs) diff --git a/src/Numeric/SGD/DataSet.hs b/src/Numeric/SGD/DataSet.hs --- a/src/Numeric/SGD/DataSet.hs +++ b/src/Numeric/SGD/DataSet.hs @@ -9,6 +9,7 @@ ( -- * Dataset DataSet (..) +, shuffle -- * Reading , loadData , randomSample @@ -26,11 +27,13 @@ import System.IO.Unsafe (unsafeInterleaveIO) import System.FilePath (()) import qualified System.Random as R +import System.Random.Shuffle (shuffleM) 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 +import qualified Data.Map.Strict as M ------------------------------- @@ -66,6 +69,18 @@ -- let (i, g'') = R.next g' -- x <- dataset `elemAt` (i `mod` size dataset) -- return (x:xs, g'') + + +-- | Shuffle the dataset. +shuffle :: DataSet a -> IO (DataSet a) +shuffle DataSet{..} = do + let ixs = [0 .. size - 1] + ixs' <- shuffleM ixs + let m = M.fromList (zip ixs ixs') + return $ DataSet + { size = size + , elemAt = elemAt . (m M.!) + } -- | Random dataset sample with a specified number of elements (loaded eagerly)