diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,11 @@
+# Changelog for `numeric-optimization-backprop`
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to the
+[Haskell Package Versioning Policy](https://pvp.haskell.org/).
+
+## Unreleased
+
+## 0.1.0.0 - YYYY-MM-DD
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Masahiro Sakai (c) 2023
+
+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.
+
+    * Neither the name of Masahiro Sakai nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,28 @@
+# numeric-optimization-backprop
+
+Wrapper of [numeric-optimization](https://hackage.haskell.org/package/numeric-optimization) package for using with [backprop](https://hackage.haskell.org/package/backprop) package.
+
+## Example Usage
+
+```haskell
+{-# LANGUAGE FlexibleContexts #-}
+import Numeric.Optimization.Backprop
+import Lens.Micro
+
+main :: IO ()
+main = do
+  result <- minimize LBFGS def rosenbrock Nothing [] (-3,-4)
+  print (resultSuccess result)  -- True
+  print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]
+  print (resultValue result)  -- 1.8129771632403013e-18
+
+-- https://en.wikipedia.org/wiki/Rosenbrock_function
+rosenbrock :: Reifies s W => BVar s (Double, Double) -> BVar s Double
+rosenbrock t = sq (1 - x) + 100 * sq (y - sq x)
+  where
+    x = t ^^. _1
+    y = t ^^. _2
+
+sq :: Floating a => a -> a
+sq x = x ** 2
+```
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/examples/rosenbrock.hs b/examples/rosenbrock.hs
new file mode 100644
--- /dev/null
+++ b/examples/rosenbrock.hs
@@ -0,0 +1,20 @@
+{-# LANGUAGE FlexibleContexts #-}
+import Numeric.Optimization.Backprop
+import Lens.Micro
+
+main :: IO ()
+main = do
+  result <- minimize LBFGS def rosenbrock Nothing [] (-3,-4)
+  print (resultSuccess result)  -- True
+  print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]
+  print (resultValue result)  -- 1.8129771632403013e-18
+
+-- https://en.wikipedia.org/wiki/Rosenbrock_function
+rosenbrock :: Reifies s W => BVar s (Double, Double) -> BVar s Double
+rosenbrock t = sq (1 - x) + 100 * sq (y - sq x)
+  where
+    x = t ^^. _1
+    y = t ^^. _2
+
+sq :: Floating a => a -> a
+sq x = x ** 2
diff --git a/numeric-optimization-backprop.cabal b/numeric-optimization-backprop.cabal
new file mode 100644
--- /dev/null
+++ b/numeric-optimization-backprop.cabal
@@ -0,0 +1,91 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.35.1.
+--
+-- see: https://github.com/sol/hpack
+
+name:           numeric-optimization-backprop
+version:        0.1.0.0
+synopsis:       Wrapper of numeric-optimization package for using with backprop package
+description:    Please see the README on GitHub at <https://github.com/msakai/nonlinear-optimization-ad/tree/master/numeric-optimization-backprop#readme>
+category:       Math, Algorithms, Optimisation, Optimization
+homepage:       https://github.com/msakai/numeric-optimization-backprop#readme
+bug-reports:    https://github.com/msakai/numeric-optimization-backprop/issues
+author:         Masahiro Sakai
+maintainer:     masahiro.sakai@gmail.com
+copyright:      Masahiro Sakai &lt;masahiro.sakai@gmail.com&gt;
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+tested-with:
+    GHC == 9.4.5
+  , GHC == 9.2.7
+  , GHC == 9.0.2
+  , GHC == 8.10.7
+  , GHC == 8.8.4
+  , GHC == 8.6.5
+extra-source-files:
+    README.md
+    CHANGELOG.md
+
+source-repository head
+  type: git
+  location: https://github.com/msakai/numeric-optimization-backprop
+
+library
+  exposed-modules:
+      Numeric.Optimization.Backprop
+      Numeric.Optimization.Backprop.ToVector
+  other-modules:
+      Paths_numeric_optimization_backprop
+  hs-source-dirs:
+      src
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints
+  build-depends:
+      backprop >=0.2.6.3 && <0.3
+    , base >=4.12 && <5
+    , containers >=0.6.0.1 && <0.7
+    , data-default-class >=0.1.2.0 && <0.2
+    , mono-traversable >=1.0.15.1 && <1.1
+    , mtl >=2.2.2 && <2.4
+    , numeric-optimization >=0.1.0.0 && <0.2.0.0
+    , primitive >=0.6.4.0
+    , reflection >=2.1.5
+    , vector >=0.12.0.2 && <0.14
+  default-language: Haskell2010
+
+executable rosenbrock-backprop
+  main-is: rosenbrock.hs
+  other-modules:
+      Paths_numeric_optimization_backprop
+  hs-source-dirs:
+      examples
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      base >=4.12 && <5
+    , containers >=0.6.0.1 && <0.7
+    , data-default-class >=0.1.2.0 && <0.2
+    , microlens >=0.4.10 && <0.5
+    , numeric-optimization >=0.1.0.0 && <0.2.0.0
+    , numeric-optimization-backprop
+  default-language: Haskell2010
+
+test-suite numeric-optimization-backprop-test
+  type: exitcode-stdio-1.0
+  main-is: Spec.hs
+  other-modules:
+      IsClose
+      Paths_numeric_optimization_backprop
+  hs-source-dirs:
+      test
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      HUnit >=1.6.0.0 && <1.7
+    , base >=4.12 && <5
+    , containers >=0.6.0.1 && <0.7
+    , data-default-class >=0.1.2.0 && <0.2
+    , hspec >=2.7.1 && <3.0
+    , microlens >=0.4.10 && <0.5
+    , numeric-optimization >=0.1.0.0 && <0.2.0.0
+    , numeric-optimization-backprop
+  default-language: Haskell2010
diff --git a/src/Numeric/Optimization/Backprop.hs b/src/Numeric/Optimization/Backprop.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/Backprop.hs
@@ -0,0 +1,138 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.Backprop
+-- Copyright   :  (c) Masahiro Sakai 2023
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-- This module is a wrapper of "Numeric.Optimization" that uses
+-- [backprop](https://hackage.haskell.org/package/backprop)'s automatic differentiation.
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.Backprop
+  (
+  -- * Main function
+    minimize
+
+  -- * Problem specification
+  , Constraint (..)
+
+  -- * Algorithm selection
+  , Method (..)
+  , isSupportedMethod
+  , Params (..)
+
+  -- * Result
+  , Result (..)
+  , Statistics (..)
+  , OptimizationException (..)
+
+  -- * Utilities and Re-exports
+  , Default (..)
+  , ToVector
+  , module Numeric.Backprop
+  ) where
+
+
+import Data.Default.Class
+import Data.Functor.Contravariant
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as VG
+import qualified Data.Vector.Storable as VS
+import Numeric.Backprop
+import qualified Numeric.Optimization as Opt
+import Numeric.Optimization hiding (minimize, IsProblem (..))
+import Numeric.Optimization.Backprop.ToVector
+
+
+data Problem a
+  = Problem
+      (forall s. Reifies s W => BVar s a -> BVar s Double)
+      (Maybe (V.Vector (Double, Double)))
+      [Constraint]
+      a
+
+
+instance (ToVector a) => Opt.IsProblem (Problem a) where
+  func (Problem f _bounds _constraints x0) x = evalBP f (updateFromVector x0 x)
+
+  bounds (Problem _f bounds _constraints _template) = bounds
+
+  constraints (Problem _f _bounds constraints _template) = constraints
+
+
+instance (Backprop a, ToVector a) => Opt.HasGrad (Problem a) where
+  grad (Problem f _bounds _constraints x0) x = toVector $ gradBP f (updateFromVector x0 x)
+
+  grad'M (Problem f _bounds _constraints x0) x gvec = do
+    case backprop f (updateFromVector x0 x) of
+      (y, g) -> do
+        writeToMVector g gvec
+        return y
+
+
+instance (Backprop a, ToVector a) => Opt.Optionally (Opt.HasGrad (Problem a)) where
+  optionalDict = hasOptionalDict
+
+
+instance Opt.Optionally (Opt.HasHessian (Problem a)) where
+  optionalDict = Nothing
+
+
+-- | Minimization of scalar function of one or more variables.
+--
+-- This is a wrapper of 'Opt.minimize' and use "Numeric.Backprop" to compute gradient.
+--
+-- Example:
+--
+-- > {-# LANGUAGE FlexibleContexts #-}
+-- > import Numeric.Optimization.Backprop
+-- > import Lens.Micro
+-- > 
+-- > main :: IO ()
+-- > main = do
+-- >   (x, result, stat) <- minimize LBFGS def rosenbrock Nothing [] (-3,-4)
+-- >   print (resultSuccess result)  -- True
+-- >   print (resultSolution result)  -- [0.999999999009131,0.9999999981094296]
+-- >   print (resultValue result)  -- 1.8129771632403013e-18
+-- > 
+-- > -- https://en.wikipedia.org/wiki/Rosenbrock_function
+-- > rosenbrock :: Reifies s W => BVar s (Double, Double) -> BVar s Double
+-- > rosenbrock t = sq (1 - x) + 100 * sq (y - sq x)
+-- >   where
+-- >     x = t ^^. _1
+-- >     y = t ^^. _2
+-- > 
+-- > sq :: Floating a => a -> a
+-- > sq x = x ** 2
+minimize
+  :: forall a. (Backprop a, ToVector a)
+  => Method  -- ^ Numerical optimization algorithm to use
+  -> Params a  -- ^ Parameters for optimization algorithms. Use 'def' as a default.
+  -> (forall s. Reifies s W => BVar s a -> BVar s Double)  -- ^ Function to be minimized.
+  -> Maybe [(Double, Double)]  -- ^ Bounds
+  -> [Constraint]  -- ^ Constraints
+  -> a -- ^ Initial value
+  -> IO (Result a)
+minimize method params f bounds constraints x0 = do
+  let bounds' :: Maybe (V.Vector (Double, Double))
+      bounds' = fmap VG.fromList bounds
+
+      prob :: Problem a
+      prob = Problem f bounds' constraints x0
+
+      params' :: Params (VS.Vector Double)
+      params' = contramap (updateFromVector x0) params
+
+  result <- Opt.minimize method params' prob (toVector x0)
+  return $ fmap (updateFromVector x0) result
diff --git a/src/Numeric/Optimization/Backprop/ToVector.hs b/src/Numeric/Optimization/Backprop/ToVector.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/Backprop/ToVector.hs
@@ -0,0 +1,301 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE DefaultSignatures #-}
+{-# LANGUAGE EmptyCase #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.Backprop.ToVector
+-- Copyright   :  (c) Masahiro Sakai 2023
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-- Conversion between a type and 'VS.Vector' 'Double'.
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.Backprop.ToVector
+  (
+  -- * ToVector class
+    ToVector (..)
+  , toVector
+
+  -- * Utilities for defining ToVector class
+
+  -- ** Generics
+  , GToVector (..)
+
+  -- ** @Foldable@/@Traversable@-based definition
+  , dimFoldable
+  , writeToMVectorFoldable
+  , updateFromVectorTraversable
+
+  -- ** @MonoFoldable@/@MonoTraversable@-based definition
+  , dimMonoFoldable
+  , writeToMVectorMonoFoldable
+  , updateFromVectorMonoTraversable
+  ) where
+
+import Control.Monad.Primitive
+import Control.Monad.State
+import qualified Data.MonoTraversable as MT
+import Data.Traversable (mapAccumL)
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as VG
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Storable.Mutable as VSM
+import qualified Data.Vector.Unboxed as VU
+import GHC.Generics
+
+import qualified Data.Functor.Identity as Functor
+import qualified Data.Functor.Compose as Functor
+import qualified Data.Functor.Const as Functor
+import qualified Data.Functor.Product as Functor
+import qualified Data.Functor.Sum as Functor
+import Data.IntMap (IntMap)
+import Data.List.NonEmpty (NonEmpty)
+import Data.Map (Map)
+import Data.Monoid
+import qualified Data.Semigroup as SG
+import Data.Sequence (Seq)
+import Data.Void
+
+-- ------------------------------------------------------------------------
+
+-- | Type that can be converted to @'VS.Vector' 'Double'@ and back.
+--
+-- Laws that should be satisfied:
+--
+-- * @'VS.length' . 'toVector' = dim@
+--
+-- * @updateFromVector a ('toVector' a) = a@
+--
+-- * @updateFromVector (updateFromVector a v1) v2 = updateFromVector a v2@
+class ToVector a where
+  -- | Dimention of the resulting vector.
+  dim :: a -> Int
+
+  -- | Destination passing style version of 'toVector'.
+  writeToMVector :: PrimMonad m => a -> VSM.MVector (PrimState m) Double -> m ()
+
+  -- | Converting @'VS.Vector' 'Double'@ back to a value
+  updateFromVector :: a -> VS.Vector Double -> a
+
+  default dim :: (Generic a, GToVector (Rep a)) => a -> Int
+  dim x = gDim (from x)
+
+  default writeToMVector :: (Generic a, GToVector (Rep a), PrimMonad m) => a -> VSM.MVector (PrimState m) Double -> m ()
+  writeToMVector x vec = gWriteToMVector (from x) vec
+
+  default updateFromVector :: (Generic a, GToVector (Rep a)) => a -> VS.Vector Double -> a
+  updateFromVector x v = to (gUpdateFromVector (from x) v)
+
+-- | Converting a value to @'VS.Vector' 'Double'@.
+toVector :: ToVector a => a -> VS.Vector Double
+toVector x = VS.create $ do
+  vec <- VSM.new (dim x)
+  writeToMVector x vec
+  return vec
+
+-- ------------------------------------------------------------------------
+
+-- | Implementation of 'dim' for the type of the form @f a@ for @'Foldable' f@.
+dimFoldable :: (Foldable f, ToVector a) => f a -> Int
+dimFoldable = getSum . foldMap (Sum . dim)
+
+-- | Implementation of 'writeToMVector' for the type of the form @f a@ for @'Foldable' f@.
+writeToMVectorFoldable :: (Foldable f, ToVector a, PrimMonad m) => f a -> VSM.MVector (PrimState m) Double -> m ()
+writeToMVectorFoldable xs vec = foldM_ f vec xs
+  where
+    f vec' x =
+      case VSM.splitAt (dim x) vec' of
+        (vec1, vec2) -> do
+          writeToMVector x vec1
+          return vec2
+
+-- | Implementation of 'updateFromVectorTraversable' for the type of the form @f a@ for @'Traversable' f@.
+updateFromVectorTraversable :: (Traversable f, ToVector a) => f a -> VS.Vector Double -> f a
+updateFromVectorTraversable xs v0 = flip evalState v0 $ do
+  forM xs $ \x -> do
+    v <- get
+    case VS.splitAt (dim x) v of
+      (v1, v2) -> do
+        put v2
+        return (updateFromVector x v1)
+
+-- ------------------------------------------------------------------------
+
+-- | Implementation of 'dim' for a 'MT.MonoFoldable' type
+dimMonoFoldable :: (MT.MonoFoldable a, ToVector (MT.Element a)) => a -> Int
+dimMonoFoldable = getSum . MT.ofoldMap (Sum . dim)
+
+-- | Implementation of 'writeToMVector' for a 'MT.MonoFoldable' type
+writeToMVectorMonoFoldable :: (MT.MonoFoldable a, ToVector (MT.Element a), PrimMonad m) => a -> VSM.MVector (PrimState m) Double -> m ()
+writeToMVectorMonoFoldable xs vec = MT.ofoldM f vec xs >> return ()
+  where
+    f vec' x =
+      case VSM.splitAt (dim x) vec' of
+        (vec1, vec2) -> do
+          writeToMVector x vec1
+          return vec2
+
+-- | Implementation of 'updateFromVector' for a 'MT.MonoTraversable' type
+updateFromVectorMonoTraversable :: (MT.MonoTraversable a, ToVector (MT.Element a)) => a -> VS.Vector Double -> a
+updateFromVectorMonoTraversable xs v0 = flip evalState v0 $ do
+  MT.oforM xs $ \x -> do
+    v <- get
+    case VS.splitAt (dim x) v of
+      (v1, v2) -> do
+        put v2
+        return (updateFromVector x v1)
+
+-- ------------------------------------------------------------------------
+
+-- | Class of generic representation types that can be converted to/from 'VS.Vector' 'Double'.
+class GToVector f where
+  gDim :: f p -> Int
+  gWriteToMVector :: PrimMonad m => f p -> VSM.MVector (PrimState m) Double -> m ()
+  gUpdateFromVector :: f p -> VS.Vector Double -> f p
+
+instance GToVector V1 where
+  gDim x = case x of { }
+  gWriteToMVector _x _vec = return ()
+  gUpdateFromVector x _v = case x of { }
+
+instance GToVector U1 where
+  gDim _ = 0
+  gWriteToMVector _x _vec = return ()
+  gUpdateFromVector x _v = x
+
+instance (GToVector f, GToVector g) => GToVector (f :+: g) where
+  gDim (L1 x) = gDim x
+  gDim (R1 x) = gDim x
+  gWriteToMVector (L1 x) vec = gWriteToMVector x vec
+  gWriteToMVector (R1 x) vec = gWriteToMVector x vec
+  gUpdateFromVector (L1 x) v = L1 (gUpdateFromVector x v)
+  gUpdateFromVector (R1 x) v = R1 (gUpdateFromVector x v)
+
+instance (GToVector f, GToVector g) => GToVector (f :*: g) where
+  gDim (a :*: b) = gDim a + gDim b
+  gWriteToMVector (a :*: b) vec =
+    case VSM.splitAt (gDim a) vec of
+      (vec1, vec2) -> do
+        gWriteToMVector a vec1
+        gWriteToMVector b vec2
+  gUpdateFromVector (a :*: b) v =
+    case VS.splitAt (gDim a) v of
+      (vec1, vec2) -> (gUpdateFromVector a vec1 :*: gUpdateFromVector b vec2)
+
+instance (ToVector c) => GToVector (K1 i c) where
+  gDim (K1 x) = dim x
+  gWriteToMVector (K1 x) vec = writeToMVector x vec
+  gUpdateFromVector (K1 x) v = K1 (updateFromVector x v)
+
+instance (GToVector f) => GToVector (M1 i t f) where
+  gDim (M1 x) = gDim x
+  gWriteToMVector (M1 x) vec = gWriteToMVector x vec
+  gUpdateFromVector (M1 x) v = M1 (gUpdateFromVector x v)
+
+-- ------------------------------------------------------------------------
+
+instance ToVector Double where
+  dim _ = 1
+  writeToMVector x vec = VSM.write vec 0 x
+  updateFromVector _x v = v VS.! 0
+
+instance (a ~ Double) => ToVector (VS.Vector a) where
+  dim x = VS.length x
+#if MIN_VERSION_vector(0,12,2)
+  writeToMVector x vec = VS.imapM_ (VSM.write vec) x
+#else
+  writeToMVector x vec = flip evalStateT 0 $ VS.mapM_ (\e -> do{ i <- get; VSM.write vec i e; put (i+1) }) x
+#endif
+  updateFromVector _x v = v
+
+instance (a ~ Double) => ToVector (VU.Vector a) where
+  dim x = VU.length x
+#if MIN_VERSION_vector(0,12,2)
+  writeToMVector x vec = VU.imapM_ (VSM.write vec) x
+#else
+  writeToMVector x vec = flip evalStateT 0 $ VU.mapM_ (\e -> do{ i <- get; VSM.write vec i e; put (i+1) }) x
+#endif
+  updateFromVector _x v = VG.convert v
+
+instance (ToVector a) => ToVector (V.Vector a) where
+  dim xs = V.sum (V.map dim xs)
+  writeToMVector xs vec = V.foldM_ f vec xs
+    where
+      f vec' x =
+        case VSM.splitAt (dim x) vec' of
+          (vec1, vec2) -> do
+            writeToMVector x vec1
+            return vec2
+  updateFromVector xs v = snd $ mapAccumL f v xs
+    where
+      f v' x =
+        case VS.splitAt (dim x) v' of
+          (v1, v2) -> (v2, updateFromVector x v1)
+
+instance ToVector Void
+
+instance ToVector ()
+instance (ToVector a, ToVector b) => ToVector (a, b)
+instance (ToVector a, ToVector b, ToVector c) => ToVector (a, b, c)
+instance (ToVector a, ToVector b, ToVector c, ToVector d) => ToVector (a, b, c, d)
+instance (ToVector a, ToVector b, ToVector c, ToVector d, ToVector e) => ToVector (a, b, c, d, e)
+
+instance (ToVector a) => ToVector (Maybe a)
+
+instance ToVector a => ToVector (SG.Min a)
+instance ToVector a => ToVector (SG.Max a)
+instance ToVector a => ToVector (SG.First a)
+instance ToVector a => ToVector (SG.Last a)
+instance ToVector a => ToVector (SG.WrappedMonoid a)
+#if !MIN_VERSION_base(4,16,0)
+instance ToVector a => ToVector (SG.Option a)
+#endif
+instance (ToVector a, ToVector b) => ToVector (SG.Arg a b)
+
+instance ToVector a => ToVector (Dual a)
+instance ToVector a => ToVector (Sum a)
+instance ToVector a => ToVector (Product a)
+instance ToVector a => ToVector (First a)
+instance ToVector a => ToVector (Last a)
+instance ToVector (f a) => ToVector (Alt f a)
+instance ToVector (f a) => ToVector (Ap f a)
+
+instance ToVector a => ToVector (Functor.Identity a)
+instance ToVector (f (g a)) => ToVector (Functor.Compose f g a)
+instance ToVector w => ToVector (Functor.Const w a)
+instance (ToVector (f a), ToVector (g a)) => ToVector (Functor.Product f g a)
+instance (ToVector (f a), ToVector (g a)) => ToVector (Functor.Sum f g a)
+
+instance ToVector a => ToVector [a] where
+  dim = dimFoldable
+  writeToMVector = writeToMVectorFoldable
+  updateFromVector = updateFromVectorTraversable
+
+instance ToVector a => ToVector (NonEmpty a) where
+  dim = dimFoldable
+  writeToMVector = writeToMVectorFoldable
+  updateFromVector = updateFromVectorTraversable
+
+instance ToVector a => ToVector (Map k a) where
+  dim = dimFoldable
+  writeToMVector = writeToMVectorFoldable
+  updateFromVector = updateFromVectorTraversable
+
+instance ToVector a => ToVector (IntMap a) where
+  dim = dimFoldable
+  writeToMVector = writeToMVectorFoldable
+  updateFromVector = updateFromVectorTraversable
+
+instance ToVector a => ToVector (Seq a) where
+  dim = dimFoldable
+  writeToMVector = writeToMVectorFoldable
+  updateFromVector = updateFromVectorTraversable
+
+-- ------------------------------------------------------------------------
diff --git a/test/IsClose.hs b/test/IsClose.hs
new file mode 100644
--- /dev/null
+++ b/test/IsClose.hs
@@ -0,0 +1,137 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module IsClose
+  (
+  -- Tolerance type
+    Tol (..)
+
+  -- AllClose class
+  , AllClose (..)
+  , allCloseRawUnit
+  , allCloseRawRealFrac
+  , allCloseRawRealFloat
+
+  -- * Re-exports
+  , Default (..)
+
+  -- * HUnit
+  , assertAllClose
+  ) where
+
+import Data.Default.Class
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.Map (Map)
+import qualified Data.Map as Map
+import Data.Monoid
+import Data.Semigroup
+import GHC.Stack (HasCallStack)
+import Test.HUnit
+import Text.Printf
+
+-- ------------------------------------------------------------------------
+
+-- | Tolerance
+--
+-- Values @a@ and @b@ are considered /close/ if @abs (a - b) <= atol + rtol * abs b@.
+data Tol a
+  = Tol
+  { rtol :: a -- ^ The relative tolerance parameter (default: @1e-05@)
+  , atol :: a -- ^ The absolute tolerance parameter (default: @1e-08@)
+  , equalNan :: Bool -- ^ Whether to compare NaN’s as equal (default: @False@)
+  } deriving (Show)
+
+instance RealFrac a => Default (Tol a) where
+  def = Tol
+    { rtol = 1e-05
+    , atol = 1e-08
+    , equalNan = False
+    }
+
+-- ------------------------------------------------------------------------
+
+class Real r => AllClose r a where
+  -- | Returns number of mismatches, number of elements, maximal absolute difference, and maximal relative difference.
+  -- Returns @'Ap' 'Nothing'@ if given values are incomparable.
+  allCloseRaw :: Tol r -> a -> a -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+
+  -- | Returns 'True' if the two arrays are equal within the given tolerance; 'False' otherwise.
+  allClose :: Tol r -> a -> a -> Bool
+  allClose tol x y =
+    case getAp (allCloseRaw tol x y) of
+      Nothing -> False
+      Just (Sum numMismatched, _, _, _) -> numMismatched == 0
+
+allCloseRawRealFrac :: RealFrac r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+allCloseRawRealFrac t a b = Ap $ Just $
+  ( Sum $ if abs (a - b) <= atol t + rtol t * abs b then 0 else 1
+  , Sum 1
+  , Max (abs (a - b))
+  , Max (abs (a - b) / abs b)
+  )
+
+allCloseRawRealFloat :: RealFloat r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+allCloseRawRealFloat t a b
+  | isNaN a /= isNaN b = Ap Nothing
+  | otherwise = Ap $ Just $
+      ( Sum $ if (equalNan t && isNaN a && isNaN b) || a == b || abs (a - b) <= atol t + rtol t * abs b then 0 else 1
+      , Sum 1
+      , Max (abs (a - b))
+      , Max (abs (a - b) / abs b)
+      )
+
+allCloseRawUnit :: Num r => Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+allCloseRawUnit = Ap (Just (Sum 0, Sum 0, Max 0, Max 0))
+
+instance AllClose Rational Rational where
+  allCloseRaw = allCloseRawRealFrac
+
+instance AllClose Double Double where
+  allCloseRaw = allCloseRawRealFloat
+
+instance (AllClose r a) => AllClose r (Maybe a) where
+  allCloseRaw tol (Just a) (Just b) = allCloseRaw tol a b
+  allCloseRaw _ Nothing Nothing = allCloseRawUnit
+  allCloseRaw _ _ _ = Ap Nothing
+
+instance (AllClose r v) => AllClose r [v] where
+  allCloseRaw tol xs ys
+    | length xs == length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip xs ys])
+    | otherwise = Ap Nothing
+
+instance (Ord k, AllClose r v) => AllClose r (Map k v) where
+  allCloseRaw tol m1 m2
+    | Map.keys m1 == Map.keys m2 = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (Map.elems m1) (Map.elems m2)])
+    | otherwise = Ap Nothing
+
+instance (AllClose r v1, AllClose r v2) => AllClose r (v1, v2)  where
+  allCloseRaw tol (x1,y1) (x2,y2) = allCloseRaw tol x1 x2 <> allCloseRaw tol y1 y2
+
+-- ------------------------------------------------------------------------
+
+-- | Assert that two objects are equal up to desired tolerance.
+assertAllClose
+  :: (HasCallStack, AllClose r a, Show r, Show a)
+  => Tol r
+  -> a -- ^ actual
+  -> a -- ^ desired
+  -> Assertion
+assertAllClose tol a b =
+  case getAp (allCloseRaw tol a b) of
+    Nothing ->
+      assertString $ unlines $ header ++ ["x and y nan location mismatch:"] ++ footer
+    Just (Sum numMismatch, Sum numTotal, Max absDiff, Max relDiff)
+      | numMismatch == 0 -> return ()
+      | otherwise ->
+          assertString $ unlines $
+            header ++
+            [ printf "Mismatched elements: %d / %d (%f%%)" numMismatch numTotal (fromIntegral numMismatch * 100 / fromIntegral numTotal :: Double)
+            , " Max absolute difference: " ++ show absDiff
+            , " Max relative difference: " ++ show relDiff
+            ] ++ footer
+   where
+     header, footer :: [String]
+     header = [printf "Not equal to tolerance rtol=%s, atol=%s" (show (rtol tol)) (show (atol tol)), ""]
+     footer = [" x: " ++ show a, " y: " ++ show b]
+
+-- ------------------------------------------------------------------------
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,26 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE RankNTypes #-}
+import Test.Hspec
+
+import Numeric.Optimization.Backprop
+import Lens.Micro
+import IsClose
+
+
+main :: IO ()
+main = hspec $ do
+  describe "minimize" $ do
+    context "when given rosenbrock function" $
+      it "returns the global optimum" $ do
+        result <- minimize LBFGS def rosenbrock Nothing [] (-3,-4)
+        resultSuccess result `shouldBe` True
+        assertAllClose (def :: Tol Double) (resultSolution result) (1,1)
+
+
+-- https://en.wikipedia.org/wiki/Rosenbrock_function
+rosenbrock :: forall s. Reifies s W => BVar s (Double, Double) -> BVar s Double
+rosenbrock t = sq (1 - x) + 100 * sq (y - sq x)
+  where
+    sq x = x ** 2
+    x = t ^^. _1
+    y = t ^^. _2