diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Numeric/SGD.hs b/Numeric/SGD.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/SGD.hs
@@ -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'')
diff --git a/Numeric/SGD/Grad.hs b/Numeric/SGD/Grad.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/SGD/Grad.hs
@@ -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)
diff --git a/Numeric/SGD/LogSigned.hs b/Numeric/SGD/LogSigned.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/SGD/LogSigned.hs
@@ -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)
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/examples/example1.hs b/examples/example1.hs
new file mode 100644
--- /dev/null
+++ b/examples/example1.hs
@@ -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
diff --git a/sgd.cabal b/sgd.cabal
new file mode 100644
--- /dev/null
+++ b/sgd.cabal
@@ -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