sgd (empty) → 0.1.0
raw patch · 7 files changed
+505/−0 lines, 7 filesdep +basedep +containersdep +logfloatsetup-changed
Dependencies added: base, containers, logfloat, monad-par, primitive, random, vector
Files
- LICENSE +26/−0
- Numeric/SGD.hs +163/−0
- Numeric/SGD/Grad.hs +96/−0
- Numeric/SGD/LogSigned.hs +70/−0
- Setup.lhs +4/−0
- examples/example1.hs +104/−0
- sgd.cabal +42/−0
+ LICENSE view
@@ -0,0 +1,26 @@+Copyright (c) 2012, IPI PAN+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Numeric/SGD.hs view
@@ -0,0 +1,163 @@+{-# LANGUAGE RecordWildCards #-}++-- | Stochastic gradient descent implementation using mutable+-- vectors for efficient update of the parameters vector.+-- A user is provided with the immutable version of parameters vector+-- so he is able to compute the gradient outside the IO/ST monad.+-- Currently only the Gaussian priors are implemented.+--+-- This is a preliminary version of the SGD library and API may change+-- in future versions.++module Numeric.SGD+( SgdArgs (..)+, sgdArgsDefault+, Dataset+, Para+, sgd+, sgdM+, module Numeric.SGD.Grad+) where++import Control.Applicative (Applicative)+import Control.Monad (forM_)+import Control.Monad.ST (ST, runST)+import qualified System.Random as R+import qualified Data.Vector as V+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++-- | 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 }++-- | Dataset with elements of x type.+type Dataset x = V.Vector x++-- | Vector of parameters.+type Para = U.Vector Double ++-- | Type synonym for mutable vector with Double values.+type MVect m = UM.MVector (Prim.PrimState m) Double++-- | Pure version of the stochastic gradient descent method.+sgd :: SgdArgs -- ^ SGD parameter values+ -> (Para -> x -> Grad) -- ^ Gradient for dataset element+ -> Dataset x -- ^ Dataset+ -> Para -- ^ Starting point+ -> Para -- ^ SGD result+sgd sgdArgs mkGrad dataset x0 =+ let dummy _ _ = return ()+ in runST $ sgdM sgdArgs dummy mkGrad dataset x0++-- | Monadic version of the stochastic gradient descent method.+-- A notification function can be used to provide user with+-- information about the progress of the learning.+{-# SPECIALIZE sgdM :: SgdArgs+ -> (Para -> Int -> IO ())+ -> (Para -> x -> Grad)+ -> Dataset x -> Para -> IO Para #-}+{-# SPECIALIZE sgdM :: SgdArgs+ -> (Para -> Int -> ST s ())+ -> (Para -> x -> Grad)+ -> Dataset x -> Para -> ST s Para #-}+sgdM+ :: (Applicative m, Prim.PrimMonad m)+ => SgdArgs -- ^ SGD parameter values+ -> (Para -> Int -> m ()) -- ^ Notification run every update+ -> (Para -> x -> Grad) -- ^ Gradient for dataset element+ -> Dataset x -- ^ Dataset+ -> Para -- ^ Starting point+ -> m Para -- ^ SGD result+sgdM 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)+ -- | Number of completed iterations over the full dataset.+ done k+ = fromIntegral (k * batchSize)+ / fromIntegral (V.length dataset) ++ doIt u k stdGen x+ | done k > iterNum = do+ frozen <- U.unsafeFreeze x+ notify frozen k+ return frozen+ | otherwise = do+ let (batch, stdGen') = sample stdGen batchSize dataset++ -- 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.unionsWith (<+>) (map (mkGrad frozen) batch)+ let grad = parUnions (map (mkGrad frozen) batch)+ addUp grad u+ scale (gain k) u++ x' <- U.unsafeThaw frozen+ apply u x'+ doIt u (k+1) stdGen' x'++-- | Add up all gradients and store results in normal domain.+{-# SPECIALIZE addUp :: Grad -> MVect IO -> IO () #-}+{-# SPECIALIZE addUp :: Grad -> MVect (ST s) -> ST s () #-}+addUp :: Prim.PrimMonad m => Grad -> MVect m -> m ()+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.+{-# SPECIALIZE scale :: Double -> MVect IO -> IO () #-}+{-# SPECIALIZE scale :: Double -> MVect (ST s) -> ST s () #-}+scale :: Prim.PrimMonad m => Double -> MVect m -> m ()+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.+{-# SPECIALIZE apply :: MVect IO -> MVect IO -> IO () #-}+{-# SPECIALIZE apply :: MVect (ST s) -> MVect (ST s) -> ST s () #-}+apply :: Prim.PrimMonad m => MVect m -> MVect m -> m ()+apply 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)++sample :: R.RandomGen g => g -> Int -> Dataset x -> ([x], g)+sample g 0 _ = ([], g)+sample g n dataset =+ let (xs, g') = sample g (n-1) dataset+ (i, g'') = R.next g'+ x = dataset V.! (i `mod` V.length dataset)+ in (x:xs, g'')
+ Numeric/SGD/Grad.hs view
@@ -0,0 +1,96 @@+-- | 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 Control.Applicative ((<$>), (<*>))+import Data.List (foldl')+import qualified Data.IntMap as M+import Control.Monad.Par.Scheds.Direct (Par, runPar, spawn, get)++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++-- | Add normal-domain double to the gradient at the given position.+{-# INLINE add #-}+add :: Grad -> Int -> Double -> Grad+add grad i y = M.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 = M.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 unoins in the Par monad.+parUnionsP :: [Grad] -> Par Grad+parUnionsP [x] = return x+parUnionsP zs = do+ let (xs, ys) = split zs+ xsP <- spawn (parUnionsP xs)+ ysP <- spawn (parUnionsP ys)+ M.unionWith (+) <$> get xsP <*> get ysP+ where+ split [] = ([], [])+ split (x:[]) = ([x], [])+ split (x:y:rest) =+ let (xs, ys) = split rest+ in (x:xs, y:ys)
+ Numeric/SGD/LogSigned.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++-- | Module provides data type for signed log-domain calculations.++module Numeric.SGD.LogSigned+( LogSigned (..)+, logSigned+, fromPos+, fromNeg+, toNorm+) where++import qualified Data.Number.LogFloat as L+import Control.Monad.Par (NFData)++instance NFData L.LogFloat++-- | Signed real value in the logarithmic domain.+data LogSigned = LogSigned+ { pos :: {-# UNPACK #-} !L.LogFloat -- ^ Positive component+ , neg :: {-# UNPACK #-} !L.LogFloat -- ^ Negative component+ }++-- All fields are strict and unpacked.+instance NFData LogSigned++-- | 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++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)
+ Setup.lhs view
@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain
+ examples/example1.hs view
@@ -0,0 +1,104 @@+{-# 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
+ sgd.cabal view
@@ -0,0 +1,42 @@+name: sgd+version: 0.1.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.+ .+ It is a preliminary implementation of the SGD method and API may change+ in future versions.+license: BSD3+license-file: LICENSE+cabal-version: >= 1.6+copyright: Copyright (c) 2012 IPI PAN+author: Jakub Waszczuk+maintainer: waszczuk.kuba@gmail.com+stability: experimental+category: Math, Algorithms+homepage: https://github.com/kawu/sgd+build-type: Simple++extra-source-files: examples/example1.hs++library+ build-depends:+ base >= 4 && < 5+ , containers+ , vector+ , random+ , primitive+ , logfloat+ , monad-par++ exposed-modules:+ Numeric.SGD+ , Numeric.SGD.LogSigned+ , Numeric.SGD.Grad++ ghc-options: -Wall -O2++source-repository head+ type: git+ location: git://github.com/kawu/sgd.git