diff --git a/Data/Eigen/Internal.hsc b/Data/Eigen/Internal.hsc
--- a/Data/Eigen/Internal.hsc
+++ b/Data/Eigen/Internal.hsc
@@ -1,4 +1,5 @@
 {-# OPTIONS_HADDOCK hide #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
 
 {-# LANGUAGE CPP #-} 
 {-# LANGUAGE EmptyDataDecls  #-}
@@ -22,11 +23,12 @@
 #endif
 import System.IO.Unsafe
 import Data.Complex
-import Data.IORef
 import Data.Bits
+import Data.Binary
+import Data.Binary.Get
+import Data.Binary.Put
 import qualified Data.Vector.Storable as VS
 import qualified Data.ByteString as BS
-import qualified Data.ByteString.Lazy as BSL
 import qualified Data.ByteString.Internal as BSI
 
 class (Num a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where
@@ -39,6 +41,18 @@
 class Cast a b where
     cast :: a -> b
 
+instance Storable a => Binary (VS.Vector a) where
+    put vs = put (BS.length bs) >> putByteString bs where
+        (fp,fs) = VS.unsafeToForeignPtr0 vs
+        es = sizeOf (VS.head vs)
+        bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)
+        
+    get = get >>= getByteString >>= \bs -> let
+        (fp,fo,fs) = BSI.toForeignPtr bs
+        es = sizeOf (VS.head vs)
+        vs = VS.unsafeFromForeignPtr0 (plusForeignPtr fp fo) (fs `div` es)
+        in return vs
+
 -- | Complex number for FFI with the same memory layout as std::complex\<T\>
 data CComplex a = CComplex !a !a
 
@@ -77,6 +91,9 @@
 instance Cast (CComplex CDouble) (Complex Double) where; cast (CComplex x y) = cast x :+ cast y
 instance Cast (Complex Double) (CComplex CDouble) where; cast (x :+ y) = CComplex (cast x) (cast y)
 
+instance Cast a b => Cast (CTriplet a) (Int, Int, b) where; cast (CTriplet x y z) = (cast x, cast y, cast z)
+instance Cast a b => Cast (Int, Int, a) (CTriplet b) where; cast (x,y,z) = CTriplet (cast x) (cast y) (cast z)
+
 intSize :: Int
 intSize = sizeOf (undefined :: CInt)
 
@@ -121,9 +138,15 @@
 instance Code (CComplex CFloat) where; code _ = 2
 instance Code (CComplex CDouble) where; code _ = 3
 
-magicCode :: Code a => a -> CInt
-magicCode x = code x `xor` 0x45696730
+newtype MagicCode = MagicCode CInt deriving Eq
 
+instance Binary MagicCode where
+    put (MagicCode code) = putWord32be $ fromIntegral code
+    get = MagicCode . fromIntegral <$> getWord32be
+
+magicCode :: Code a => a -> MagicCode
+magicCode x = MagicCode (code x `xor` 0x45696730)
+
 #let api1 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name
 
 #api1 random,        "Ptr b -> CInt -> CInt -> IO CString"
@@ -191,6 +214,7 @@
 #api2 sparse_setIdentity,   "CSparseMatrixPtr a b -> IO CString"
 #api2 sparse_reserve,       "CSparseMatrixPtr a b -> CInt -> IO CString"
 #api2 sparse_resize,        "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"
+
 #api2 sparse_conservativeResize,    "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"
 #api2 sparse_compressInplace,       "CSparseMatrixPtr a b -> IO CString"
 #api2 sparse_uncompressInplace,     "CSparseMatrixPtr a b -> IO CString"
@@ -223,22 +247,3 @@
 #api3 sparse_la_logAbsDeterminant,  "CSolverPtr a b -> Ptr b -> IO CString"
 #api3 sparse_la_absDeterminant,     "CSolverPtr a b -> Ptr b -> IO CString"
 #api3 sparse_la_signDeterminant,    "CSolverPtr a b -> Ptr b -> IO CString"
-
-
-openStream :: BSL.ByteString -> IO (IORef BSL.ByteString)
-openStream = newIORef
-
-readStream :: IORef BSL.ByteString -> Int -> IO BS.ByteString
-readStream ref size = readIORef ref >>= \a ->
-    let (b,c) = BSL.splitAt (fromIntegral size) a
-    in if BSL.length b /= fromIntegral size
-        then fail "readStream: stream exhausted"
-        else do
-            writeIORef ref c
-            return . BS.concat . BSL.toChunks $ b
-
-closeStream :: IORef BSL.ByteString -> IO ()
-closeStream ref = BSL.null <$> readIORef ref >>= (`unless` fail "closeStream: stream underrun")
-
-readInt :: IORef BSL.ByteString -> IO CInt
-readInt st = decodeInt <$> readStream st intSize
diff --git a/Data/Eigen/Matrix.hs b/Data/Eigen/Matrix.hs
--- a/Data/Eigen/Matrix.hs
+++ b/Data/Eigen/Matrix.hs
@@ -106,8 +106,9 @@
 import Prelude hiding (null, sum, all, any, map, filter)
 import Data.Tuple
 import Data.Complex hiding (conjugate)
+import Data.Binary hiding (encode, decode)
+import qualified Data.Binary as B
 import Foreign.Ptr
-import Foreign.ForeignPtr
 import Foreign.C.Types
 import Foreign.C.String
 import Foreign.Storable
@@ -125,7 +126,6 @@
 import qualified Data.Eigen.Internal as I
 import qualified Data.Eigen.Matrix.Mutable as M
 import qualified Data.ByteString.Lazy as BSL
-import qualified Data.ByteString.Internal as BSI
 
 -- | Matrix to be used in pure computations, uses column major memory layout, features copy-free FFI with C++ <http://eigen.tuxfamily.org Eigen> library.
 
@@ -158,6 +158,26 @@
     abs = map abs
     negate = map negate
 
+-- | Matrix binary serialization
+instance I.Elem a b => Binary (Matrix a b) where
+    put (Matrix rows cols vals) = do
+        put $ I.magicCode (undefined :: b)
+        put rows
+        put cols
+        put vals
+
+    get = do
+        get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))
+        Matrix <$> get <*> get <*> get
+
+-- | Encode the matrix as a lazy byte string
+encode :: I.Elem a b => Matrix a b -> BSL.ByteString
+encode = B.encode
+
+-- | Decode matrix from the lazy byte string
+decode :: I.Elem a b => BSL.ByteString -> Matrix a b
+decode = B.decode
+
 -- | Empty 0x0 matrix
 {-# INLINE empty #-}
 empty :: I.Elem a b => Matrix a b
@@ -556,31 +576,6 @@
 unsafeWith m@(Matrix rows cols vals) f
     | not (valid m) = fail "Matrix.unsafeWith: matrix layout is invalid"
     | otherwise = VS.unsafeWith vals $ \p -> f p (I.cast rows) (I.cast cols)
-
--- | Encode the matrix as a lazy byte string
-encode :: I.Elem a b => Matrix a b -> BSL.ByteString
-encode m@(Matrix rows cols vals)
-    | valid m = BSL.fromChunks [
-            I.encodeInt (I.magicCode $ VS.head vals),
-            I.encodeInt (I.cast rows),
-            I.encodeInt (I.cast cols),
-            let (fp, fs) = VS.unsafeToForeignPtr0 vals in BSI.PS (castForeignPtr fp) 0 (fs * sizeOf (VS.head vals))]
-    | otherwise = error "Matrix.encode: matrix layout is invalid"
-
--- | Decode matrix from the lazy byte string
-decode :: I.Elem a b => BSL.ByteString -> Matrix a b
-decode st = Matrix rows cols vals where
-    (rows, cols, vals) = I.performIO $ do
-        st <- I.openStream st
-        code <- I.readInt st
-        when (code /= I.magicCode (VS.head vals)) $
-            fail "Matrix.decode: wrong matrix type"
-        rows <- I.cast <$> I.readInt st
-        cols <- I.cast <$> I.readInt st
-        BSI.PS fp fo _ <- I.readStream st (rows * cols * sizeOf (VS.head vals))
-        I.closeStream st
-        return (rows, cols, VS.unsafeFromForeignPtr0 (I.plusForeignPtr fp fo) (rows * cols))
-
 
 {-# INLINE _prop #-}
 _prop :: I.Elem a b => (Ptr b -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> a
diff --git a/Data/Eigen/SparseMatrix.hs b/Data/Eigen/SparseMatrix.hs
--- a/Data/Eigen/SparseMatrix.hs
+++ b/Data/Eigen/SparseMatrix.hs
@@ -23,6 +23,8 @@
     -- * Matrix conversions
     fromList,
     toList,
+    fromVector,
+    toVector,
     fromDenseList,
     toDenseList,
     fromMatrix,
@@ -62,6 +64,8 @@
 import Prelude hiding (map)
 import qualified Data.List as L
 import Data.Complex
+import Data.Binary hiding (encode, decode)
+import qualified Data.Binary as B
 import Foreign.C.Types
 import Foreign.C.String
 import Foreign.Storable
@@ -80,9 +84,7 @@
 import qualified Data.Eigen.Internal as I
 import qualified Data.Vector.Storable as VS
 import qualified Data.Vector.Storable.Mutable as VSM
-import qualified Data.ByteString as BS
 import qualified Data.ByteString.Lazy as BSL
-import qualified Data.ByteString.Internal as BSI
 
 {-| A versatible sparse matrix representation.
 
@@ -158,17 +160,30 @@
     (+) = add
     (-) = sub
     fromInteger x = fromList 1 1 [(0,0,fromInteger x)]
-    signum = map signum
-    abs = map abs
-    negate = map negate
+    signum = _map signum
+    abs = _map abs
+    negate = _map negate
 
--- | Not exposed, For internal use donly
-map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b
-map f m = fromList (rows m) (cols m) . P.map (\(r,c,v) -> (r,c,f v)) . toList $ m
+instance I.Elem a b => Binary (SparseMatrix a b) where
+    put m = do
+        put $ I.magicCode (undefined :: b)
+        put $ rows m
+        put $ cols m
+        put $ toVector m
 
-mk :: I.Elem a b => Ptr (I.CSparseMatrix a b) -> IO (SparseMatrix a b)
-mk p = SparseMatrix <$> FC.newForeignPtr p (I.call $ I.sparse_free p)
+    get = do
+        get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))
+        fromVector <$> get <*> get <*> get
 
+-- | Encode the sparse matrix as a lazy byte string
+encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString
+encode = B.encode
+
+
+-- | Decode sparse matrix from the lazy byte string
+decode :: I.Elem a b => BSL.ByteString -> SparseMatrix a b
+decode = B.decode
+
 -- | Stores the coefficient values of the non-zeros.
 values :: I.Elem a b => SparseMatrix a b -> VS.Vector b
 values = _getvec I.sparse_values
@@ -224,7 +239,7 @@
 
 -- | Extract rectangular block from sparse matrix defined by startRow startCol blockRows blockCols
 block :: I.Elem a b => Int -> Int -> Int -> Int -> SparseMatrix a b -> SparseMatrix a b
-block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) mk
+block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) _mk
 
 -- | Number of non-zeros elements in the sparse matrix
 nonZeros :: I.Elem a b => SparseMatrix a b -> Int
@@ -232,11 +247,11 @@
 
 -- | The matrix in the compressed format
 compress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
-compress = _unop I.sparse_makeCompressed mk
+compress = _unop I.sparse_makeCompressed _mk
 
 -- | The matrix in the uncompressed mode
 uncompress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
-uncompress = _unop I.sparse_uncompress mk
+uncompress = _unop I.sparse_uncompress _mk
 
 -- | Is this in compressed form?
 compressed :: I.Elem a b => SparseMatrix a b -> Bool
@@ -250,51 +265,57 @@
 outerSize :: I.Elem a b => SparseMatrix a b -> Int
 outerSize = _unop I.sparse_outerSize (return . I.cast)
 
--- | Suppresses all nonzeros which are much smaller than reference under the tolerence epsilon
+-- | Suppresses all nonzeros which are much smaller than reference under the tolerence @epsilon@
 pruned :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b
-pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (I.cast r) >> I.sparse_prunedRef p pr pq) mk
+pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (I.cast r) >> I.sparse_prunedRef p pr pq) _mk
 
 -- | Multiply matrix on a given scalar
 scale :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b
-scale x = _unop (\p pq -> alloca $ \px -> poke px (I.cast x) >> I.sparse_scale p px pq) mk
+scale x = _unop (\p pq -> alloca $ \px -> poke px (I.cast x) >> I.sparse_scale p px pq) _mk
 
 -- | Transpose of the sparse matrix
 transpose :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
-transpose = _unop I.sparse_transpose mk
+transpose = _unop I.sparse_transpose _mk
 
 -- | Adjoint of the sparse matrix
 adjoint :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
-adjoint = _unop I.sparse_adjoint mk
+adjoint = _unop I.sparse_adjoint _mk
 
 -- | Adding two sparse matrices by adding the corresponding entries together. You can use @(+)@ function as well.
 add :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b
-add = _binop I.sparse_add mk
+add = _binop I.sparse_add _mk
 
 -- | Subtracting two sparse matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.
 sub :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b
-sub = _binop I.sparse_sub mk
+sub = _binop I.sparse_sub _mk
 
 -- | Matrix multiplication. You can use @(*)@ function as well.
 mul :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b
-mul = _binop I.sparse_mul mk
+mul = _binop I.sparse_mul _mk
 
 -- | Construct sparse matrix of given size from the list of triplets (row, col, val)
 fromList :: I.Elem a b => Int -> Int -> [(Int, Int, a)] -> SparseMatrix a b
-fromList rows cols tris = I.performIO $ VS.unsafeWith vs $ \p -> alloca $ \pq -> do
-    I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length vs) pq
-    peek pq >>= mk
-    where vs = VS.fromList $ P.map (\(row,col,val) -> I.CTriplet (I.cast row) (I.cast col) (I.cast val)) tris
+fromList rows cols = fromVector rows cols . VS.fromList . P.map I.cast
 
+-- | Construct sparse matrix of given size from the storable vector of triplets (row, col, val)
+fromVector :: I.Elem a b => Int -> Int -> VS.Vector (I.CTriplet b) -> SparseMatrix a b
+fromVector rows cols tris = I.performIO $ VS.unsafeWith tris $ \p -> alloca $ \pq -> do
+    I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length tris) pq
+    peek pq >>= _mk
+
 -- | Convert sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included
 toList :: I.Elem a b => SparseMatrix a b -> [(Int, Int, a)]
-toList m@(SparseMatrix fp) = I.performIO $ do
+toList = P.map I.cast . VS.toList . toVector
+
+-- | Convert sparse matrix to the storable vector of triplets (row, col, val). Compressed elements will not be included
+toVector :: I.Elem a b => SparseMatrix a b -> VS.Vector (I.CTriplet b)
+toVector m@(SparseMatrix fp) = I.performIO $ do
     let size = nonZeros m
     tris <- VSM.new size
     withForeignPtr fp $ \p ->
         VSM.unsafeWith tris $ \q ->
             I.call $ I.sparse_toList p q (I.cast size)
-    let f (I.CTriplet row col val) = (I.cast row, I.cast col, I.cast val)
-    P.map f . VS.toList <$> VS.unsafeFreeze tris
+    VS.unsafeFreeze tris
 
 -- | Construct sparse matrix of two-dimensional list of values. Matrix dimensions will be detected automatically. Zero values will be compressed.
 fromDenseList :: (I.Elem a b, Eq a) => [[a]] -> SparseMatrix a b
@@ -316,7 +337,7 @@
 fromMatrix m1 = I.performIO $ alloca $ \pm0 ->
     M.unsafeWith m1 $ \vals rows cols -> do
         I.call $ I.sparse_fromMatrix vals rows cols pm0
-        peek pm0 >>= mk
+        peek pm0 >>= _mk
 
 -- | Construct dense matrix from sparse matrix
 toMatrix :: I.Elem a b => SparseMatrix a b -> M.Matrix a b
@@ -327,58 +348,6 @@
             I.call $ I.sparse_toMatrix pm1 vals rows cols
     M.unsafeFreeze m0
 
--- | Encode the sparse matrix as a lazy byte string
-encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString
-encode m@(SparseMatrix fp) = I.performIO $ do
-    let size = nonZeros m
-    tris <- VSM.new size
-    withForeignPtr fp $ \p ->
-        VSM.unsafeWith tris $ \q ->
-            I.call $ I.sparse_toList p q (I.cast size)
-    tris <- VS.unsafeFreeze tris
-    let
-        tri@(I.CTriplet _ _ val) = VS.head tris
-
-    return $ BSL.fromChunks [
-        encodeInt (I.magicCode val),
-        encodeInt (I.cast $ rows m),
-        encodeInt (I.cast $ cols m),
-        encodeInt (I.cast $ size),
-        let (fp, fs) = VS.unsafeToForeignPtr0 tris in BSI.PS (castForeignPtr fp) 0 (fs * sizeOf tri)]
-
-    where
-        encodeInt :: CInt -> BS.ByteString
-        encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr
-
-
--- | Decode sparse matrix from the lazy byte string
-decode :: forall a b . I.Elem a b => BSL.ByteString -> SparseMatrix a b
-decode st = I.performIO $ do
-    let
-        val = undefined :: b
-        tri = undefined :: I.CTriplet b
-        triSize = sizeOf tri
-
-    st <- I.openStream st
-
-    code <- I.decodeInt <$> I.readStream st I.intSize
-    when (code /= I.magicCode val) $
-        fail "SparseMatrix.decode: wrong matrix type"
-
-    rows <- I.readInt st
-    cols <- I.readInt st
-    size <- I.readInt st
-
-    BSI.PS fp fo _ <- I.readStream st (I.cast size * triSize)
-
-    I.closeStream st
-
-    let tris = VS.unsafeFromForeignPtr0 (I.plusForeignPtr fp fo) (I.cast size)
-
-    VS.unsafeWith tris $ \p -> alloca $ \pq -> do
-        I.call $ I.sparse_fromList rows cols p size pq
-        peek pq >>= mk
-
 -- | Yield an immutable copy of the mutable matrix
 freeze :: I.Elem a b => SMM.IOSparseMatrix a b -> IO (SparseMatrix a b)
 freeze (SMM.IOSparseMatrix fp) = SparseMatrix <$> _clone fp
@@ -426,3 +395,11 @@
     I.call $ I.sparse_clone p pq
     q <- peek pq
     FC.newForeignPtr q $ I.call $ I.sparse_free q
+
+_map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b
+_map f m = fromVector (rows m) (cols m) . VS.map g . toVector $ m where
+    g (I.CTriplet r c v) = I.CTriplet r c $ I.cast $ f $ I.cast v
+
+_mk :: I.Elem a b => Ptr (I.CSparseMatrix a b) -> IO (SparseMatrix a b)
+_mk p = SparseMatrix <$> FC.newForeignPtr p (I.call $ I.sparse_free p)
+
diff --git a/eigen.cabal b/eigen.cabal
--- a/eigen.cabal
+++ b/eigen.cabal
@@ -1,5 +1,5 @@
 name:           eigen
-version:        2.1.4
+version:        2.1.5
 homepage:       https://github.com/osidorkin/haskell-eigen
 synopsis:       Eigen C++ library (linear algebra: matrices, sparse matrices, vectors, numerical solvers).
 description:    This module provides Haskell binding for <http://eigen.tuxfamily.org/ Eigen C++ library>.
@@ -1308,6 +1308,7 @@
     build-depends:      base >= 3 && < 5,
                         vector >= 0.5 && < 0.12,
                         primitive >= 0.1 && < 0.7,
+                        binary,
                         bytestring,
                         transformers >= 0.3
 
@@ -1322,20 +1323,20 @@
 Test-Suite test-solve
     type:               exitcode-stdio-1.0
     main-is:            test/solve.hs
-    build-depends:      base, primitive, vector, eigen
+    build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen
 
 Test-Suite test-solve-sparse
     type:               exitcode-stdio-1.0
     main-is:            test/solve-sparse.hs
-    build-depends:      base, primitive, vector, eigen
+    build-depends:      base, primitive, vector, bytestring, transformers, mtl, binary, eigen
 
 Test-Suite test-rank
     type:               exitcode-stdio-1.0
     main-is:            test/rank.hs
-    build-depends:      base, primitive, vector, eigen
+    build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen
 
 Test-Suite test-regression
     type:               exitcode-stdio-1.0
     main-is:            test/regression.hs
-    build-depends:      base, primitive, vector, eigen
+    build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen