diff --git a/sgd.cabal b/sgd.cabal
--- a/sgd.cabal
+++ b/sgd.cabal
@@ -1,5 +1,5 @@
 name:               sgd
-version:            0.3.7
+version:            0.4.0
 synopsis:           Stochastic gradient descent
 description:
     Implementation of a Stochastic Gradient Descent optimization method.
@@ -26,21 +26,22 @@
     build-depends:
         base            >= 4        && < 5
       , containers      >= 0.4      && < 0.6
-      , vector          >= 0.10     && < 0.11
-      , random          >= 1.0      && < 1.1
-      , primitive       >= 0.5      && < 0.6
-      , logfloat        >= 0.12     && < 0.13
+      , vector          >= 0.10     && < 0.13
+      , random          >= 1.0      && < 1.2
+      , primitive       >= 0.5      && < 0.7
+      , logfloat        >= 0.12     && < 0.14
       , monad-par       >= 0.3.4    && < 0.4
-      , deepseq         >= 1.3      && < 1.4
-      , binary          >= 0.5      && < 0.8
+      , deepseq         >= 1.3      && < 1.5
+      , binary          >= 0.5      && < 0.9
       , bytestring      >= 0.9      && < 0.11
       , mtl             >= 2.0      && < 2.3
-      , filepath        >= 1.3      && < 1.4
-      , temporary       >= 1.1      && < 1.2
+      , filepath        >= 1.3      && < 1.5
+      , temporary       >= 1.1      && < 1.3
       , lazy-io         >= 0.1      && < 0.2
 
     exposed-modules:
         Numeric.SGD
+      , Numeric.SGD.Momentum
       , Numeric.SGD.Dataset
       , Numeric.SGD.LogSigned
       , Numeric.SGD.Grad
diff --git a/src/Numeric/SGD.hs b/src/Numeric/SGD.hs
--- a/src/Numeric/SGD.hs
+++ b/src/Numeric/SGD.hs
@@ -21,7 +21,7 @@
 ) where
 
 
-import           Control.Monad (forM_)
+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
@@ -57,7 +57,7 @@
 
 
 -- | Vector of parameters.
-type Para       = U.Vector Double 
+type Para       = U.Vector Double
 
 
 -- | Type synonym for mutable vector with Double values.
@@ -75,17 +75,31 @@
     -> 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
+  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 :: 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 stdGen x
       | done k > iterNum = do
         frozen <- U.unsafeFreeze x
@@ -94,6 +108,20 @@
       | otherwise = do
         (batch, stdGen') <- sample stdGen 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
@@ -105,7 +133,7 @@
         scale (gain k) u
 
         x' <- U.unsafeThaw frozen
-        apply u x'
+        u `addTo` x'
         doIt u (k+1) stdGen' x'
 
 
@@ -128,8 +156,8 @@
 
 -- | Apply gradient to the parameters vector, that is add the first vector to
 -- the second one.
-apply :: MVect -> MVect -> IO ()
-apply w v = do 
+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
diff --git a/src/Numeric/SGD/Momentum.hs b/src/Numeric/SGD/Momentum.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/SGD/Momentum.hs
@@ -0,0 +1,201 @@
+{-# LANGUAGE RecordWildCards #-}
+
+
+-- | A version of `Numeric.SGD` extended with 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
+
+
+-- | 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 the gradient with
+      -- momentum becomes much 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 (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 :: 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)
+
+
+-- | 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)