diff --git a/CHANGES.md b/CHANGES.md
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,6 +1,10 @@
 Revision history for Haskell package learning-hmm
 ===
 
+## Version 0.2.1.0
+- Performance improvements with employing 'hmatrix'
+- `withEmission` now does iterative re-estimation based on the Viterbi path
+
 ## Version 0.2.0.0
 - Remove dependency on the 'logfloat' package
 - Performance improvements
diff --git a/learning-hmm.cabal b/learning-hmm.cabal
--- a/learning-hmm.cabal
+++ b/learning-hmm.cabal
@@ -1,5 +1,5 @@
 name:                learning-hmm
-version:             0.2.0.0
+version:             0.2.1.0
 stability:           experimental
 
 synopsis:            Yet another library for hidden Markov models
@@ -28,12 +28,12 @@
   other-modules:     Data.Random.Distribution.Categorical.Util
                    , Data.Random.Distribution.Simplex
                    , Data.Vector.Generic.Util
-                   , Data.Vector.Generic.Util.LinearAlgebra
                    , Learning.HMM.Internal
   -- other-extensions:  
   build-depends:     base >=4.7 && <4.8
                    , containers
                    , deepseq
+                   , hmatrix < 0.16.1
                    , random-fu
                    , random-source
                    , vector
diff --git a/src/Data/Vector/Generic/Util/LinearAlgebra.hs b/src/Data/Vector/Generic/Util/LinearAlgebra.hs
deleted file mode 100644
--- a/src/Data/Vector/Generic/Util/LinearAlgebra.hs
+++ /dev/null
@@ -1,141 +0,0 @@
--- | Operators commonly used in the basic linear algebra. Note that all the
---   functions defined here do not check the dimension/length of
---   vectors/matrices.
-module Data.Vector.Generic.Util.LinearAlgebra (
-  -- * Pairwise operators
-    (>+>)
-  -- , (>->)
-  , (>.>)
-  , (>/>)
-  , (#+#)
-  -- , (#-#)
-  -- , (#.#)
-  -- , (#/#)
-
-  -- * Scalar-vector/vector-scalar operators
-  -- , (+>)
-  -- , (->)
-  , (.>)
-  -- , (/>)
-  -- , (>+)
-  -- , (>-)
-  -- , (>.)
-  , (>/)
-
-  -- * Scalar-matrix/matrix-scalar operators
-  -- , (+#)
-  -- , (-#)
-  -- , (.#)
-  -- , (/#)
-  -- , (#+)
-  -- , (#-)
-  -- , (#.)
-  , (#/)
-
-  -- * Dot and matrix-vector/vector-matrix products
-  , (<.>)
-  , (#.>)
-  , (<.#)
-
-  -- * Unary operators
-  , transpose
-  ) where
-
-import Prelude hiding (any, head, map, null, sum, tail, zipWith)
-import Data.Vector.Generic (
-    Vector, any, cons, convert, empty, head, map, null, sum, tail, zipWith
-  )
-
--- $setup
--- >>> :module + Data.Vector
-
--- | Pairwise addition between two vectors
---
--- >>> fromList [1, 2] >+> fromList [3, 4 :: Int]
--- fromList [4,6]
-(>+>) :: (Num a, Vector v a) => v a -> v a -> v a
-{-# INLINE (>+>) #-}
-u >+> v = zipWith (+) u v
-
--- | Pairwise product between two vectors
---
--- >>> fromList [1, 2] >.> fromList [3, 4 :: Double]
--- fromList [3.0,8.0]
-(>.>) :: (Num a, Vector v a) => v a -> v a -> v a
-{-# INLINE (>.>) #-}
-u >.> v = zipWith (*) u v
-
--- | Pairwise division between two vectors
---
--- >>> fromList [1, 2] >/> fromList [3, 4 :: Double]
--- fromList [0.3333333333333333,0.5]
-(>/>) :: (Fractional a, Vector v a) => v a -> v a -> v a
-{-# INLINE (>/>) #-}
-u >/> v = zipWith (/) u v
-
--- | Pairwise addition between two matrices
---
--- >>> fromList [fromList [1, 2], fromList [3, 4]] #+# fromList [fromList [5, 6], fromList [7, 8 :: Int]]
--- fromList [fromList [6,8],fromList [10,12]]
-(#+#) :: (Num a, Vector v a, Vector w (v a)) => w (v a) -> w (v a) -> w (v a)
-{-# INLINE (#+#) #-}
-m #+# n = zipWith (>+>) m n
-
--- | Scalar-vector product
---
--- >>> 2 .> fromList [1, 2 :: Integer]
--- fromList [2,4]
-(.>) :: (Num a, Vector v a) => a -> v a -> v a
-{-# INLINE (.>) #-}
-s .> v = map (s *) v
-
--- | Vector-scalar division
---
--- >>> fromList [1, 2 :: Double] >/ 2
--- fromList [0.5,1.0]
-(>/) :: (Fractional a, Vector v a) => v a -> a -> v a
-{-# INLINE (>/) #-}
-v >/ s = map (/ s) v
-
--- | Matrix-scalar division
---
--- >>> fromList [fromList [1, 2], fromList [3, 4 :: Double]] #/ 2
--- fromList [fromList [0.5,1.0],fromList [1.5,2.0]]
-(#/) :: (Fractional a, Vector v a, Vector w (v a)) => w (v a) -> a -> w (v a)
-{-# INLINE (#/) #-}
-m #/ s = map (>/ s) m
-
--- | Dot product
---
--- >>> fromList [1, 2] <.> fromList [3, 4 :: Int]
--- 11
-(<.>) :: (Num a, Vector v a) => v a -> v a -> a
-{-# INLINE (<.>) #-}
-u <.> v = sum $ u >.> v
-
--- | Matrix-vector product
---
--- >>> fromList [fromList [1, 2], fromList [3, 4]] #.> fromList [1, 2 :: Double]
--- fromList [5.0,11.0]
-(#.>) :: (Num a, Vector v a, Vector w (v a), Vector w a) => w (v a) -> v a -> v a
-{-# INLINE (#.>) #-}
-m #.> v = convert $ map (<.> v) m
-
--- | Vector-matrix product
---
--- >>> fromList [1, 2 :: Double] <.# fromList [fromList [1, 2], fromList [3, 4]]
--- fromList [7.0,10.0]
-(<.#) :: (Num a, Vector v a, Vector w (v a), Vector w a) => v a -> w (v a) -> v a
-{-# INLINE (<.#) #-}
-v <.# m | any null m = empty
-        | otherwise  = (v <.> convert (map head m)) `cons` (v <.# map tail m)
-
--- | Matrix transpose
---
--- >>> transpose $ fromList [fromList "ab", fromList "cd"]
--- fromList [fromList "ac",fromList "bd"]
-transpose :: (Vector v a, Vector w (v a), Vector w a) => w (v a) -> w (v a)
-{-# INLINE transpose #-}
-transpose m
-  | any null m = empty
-  | otherwise  = convert (map head m) `cons` transpose (map tail m)
diff --git a/src/Learning/HMM.hs b/src/Learning/HMM.hs
--- a/src/Learning/HMM.hs
+++ b/src/Learning/HMM.hs
@@ -22,10 +22,15 @@
 import Data.Random.Distribution.Categorical.Util ()
 import Data.Random.RVar (RVar)
 import Data.Random.Sample (sample)
-import qualified Data.Vector as V ((!), elemIndex, fromList, map, toList)
+import qualified Data.Vector as V (
+    elemIndex, fromList, map, toList, unsafeIndex
+  )
 import qualified Data.Vector.Generic as G (convert)
-import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)
-import qualified Data.Vector.Unboxed as U (fromList, toList)
+import qualified Data.Vector.Unboxed as U (fromList)
+import qualified Numeric.LinearAlgebra.Data as H (
+    (!), fromList, fromLists, toList
+  )
+import qualified Numeric.LinearAlgebra.HMatrix as H (tr)
 import Learning.HMM.Internal
 
 -- | Parameter set of the hidden Markov model. Direct use of the
@@ -33,13 +38,13 @@
 data HMM s o = HMM { states  :: [s] -- ^ Hidden states
                    , outputs :: [o] -- ^ Outputs
                    , initialStateDist :: Categorical Double s
-                     -- ^ Categorical distribution of initial states
+                     -- ^ Categorical distribution of initial state
                    , transitionDist :: s -> Categorical Double s
-                     -- ^ Categorical distribution of next states
-                     --   conditioned by the previous states
+                     -- ^ Categorical distribution of next state
+                     --   conditioned by the previous state
                    , emissionDist :: s -> Categorical Double o
-                     -- ^ Categorical distribution of outputs conditioned
-                     --   by the hidden states
+                     -- ^ Categorical distribution of output conditioned
+                     --   by the hidden state
                    }
 
 instance (Show s, Show o) => Show (HMM s o) where
@@ -92,10 +97,10 @@
 init ss os = fromHMM' ss os <$> init' (length ss) (length os)
 
 -- | @model \`withEmission\` xs@ returns a model in which the
---   'emissionDist' is updated by using the observed outputs @xs@. The
---   'emissionDist' is set to be normalized histograms each of which is
---   calculated from a partial set of @xs@ for each state. The partition is
---   based on the most likely state path obtained by the Viterbi algorithm.
+--   'emissionDist' is updated by re-estimations using the observed outputs
+--   @xs@. The 'emissionDist' is set to be normalized histograms each of
+--   which is calculated from segumentations of @xs@ based on the Viterbi
+--   state path.
 withEmission :: (Eq s, Eq o) => HMM s o -> [o] -> HMM s o
 withEmission model xs = fromHMM' ss os $ withEmission' model' xs'
   where
@@ -112,12 +117,13 @@
 viterbi model xs =
   checkModelIn "viterbi" model `seq`
   checkDataIn "viterbi" model xs `seq`
-  first (V.toList . V.map (ss V.!) . G.convert) $ viterbi' model' xs'
+  first toStates $ viterbi' model' xs'
   where
-    ss     = V.fromList $ states model
+    ss'    = V.fromList $ states model
     os'    = V.fromList $ outputs model
     model' = toHMM' model
     xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs
+    toStates = V.toList . V.map (V.unsafeIndex ss') . G.convert
 
 -- | @baumWelch model xs@ iteratively performs the Baum-Welch algorithm
 --   using the observed outputs @xs@, and returns a list of updated models
@@ -188,29 +194,29 @@
   where
     pi0 = initialStateDist' hmm'
     w   = transitionDist' hmm'
-    phi = G.transpose $ emissionDistT' hmm'
-    pi0'   = zip (U.toList pi0) ss
-    w' i   = zip (U.toList $ w V.! i) ss
-    phi' i = zip (U.toList $ phi V.! i) os
+    phi = H.tr $ emissionDistT' hmm'
+    pi0'   = zip (H.toList pi0) ss
+    w' i   = zip (H.toList $ w H.! i) ss
+    phi' i = zip (H.toList $ phi H.! i) os
 
 -- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transitionDist'',
 --   and 'emissionDistT'' are normalized.
 toHMM' :: (Eq s, Eq o) => HMM s o -> HMM'
 toHMM' hmm = HMM' { nStates'          = length ss
                   , nOutputs'         = length os
-                  , initialStateDist' = U.fromList pi0'
-                  , transitionDist'   = V.fromList w'
-                  , emissionDistT'    = V.fromList phi'
+                  , initialStateDist' = pi0
+                  , transitionDist'   = w
+                  , emissionDistT'    = phi'
                   }
   where
-    ss  = states hmm
-    os  = outputs hmm
-    pi0 = C.normalizeCategoricalPs $ initialStateDist hmm
-    w   = C.normalizeCategoricalPs . transitionDist hmm
-    phi = C.normalizeCategoricalPs . emissionDist hmm
-    pi0' = [pdf pi0 s | s <- ss]
-    w'   = [U.fromList [pdf (w s) s' | s' <- ss] | s <- ss]
-    phi' = [U.fromList [pdf (phi s) o | s <- ss] | o <- os]
+    ss   = states hmm
+    os   = outputs hmm
+    pi0_ = C.normalizeCategoricalPs $ initialStateDist hmm
+    w_   = C.normalizeCategoricalPs . transitionDist hmm
+    phi_ = C.normalizeCategoricalPs . emissionDist hmm
+    pi0  = H.fromList [pdf pi0_ s | s <- ss]
+    w    = H.fromLists [[pdf (w_ s) s' | s' <- ss] | s <- ss]
+    phi' = H.fromLists [[pdf (phi_ s) o | s <- ss] | o <- os]
 
 errorIn :: String -> String -> a
 errorIn fun msg = error $ "Learning.HMM." ++ fun ++ ": " ++ msg
diff --git a/src/Learning/HMM/Internal.hs b/src/Learning/HMM/Internal.hs
--- a/src/Learning/HMM/Internal.hs
+++ b/src/Learning/HMM/Internal.hs
@@ -18,21 +18,28 @@
 import Data.Random.RVar (RVar)
 import Data.Random.Distribution.Simplex (stdSimplex)
 import qualified Data.Vector as V (
-    Vector, (!), filter, foldl1', freeze, fromList, generate, map , tail
-  , zip, zipWith, zipWith3
+    Vector, filter, foldl1', map, unsafeFreeze, unsafeIndex, unsafeTail
+  , zip, zipWith3
   )
 import qualified Data.Vector.Generic as G (convert)
 import qualified Data.Vector.Generic.Util as G (frequencies)
-import Data.Vector.Generic.Util.LinearAlgebra (
-    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#.>), (<.#)
+import qualified Data.Vector.Mutable as MV (
+    unsafeNew, unsafeRead, unsafeWrite
   )
-import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)
-import qualified Data.Vector.Mutable as MV (new, read, write)
 import qualified Data.Vector.Unboxed as U (
-    Vector, (!), freeze, fromList, generate, length, map, maxIndex, maximum
-  , replicate, sum, tail, zip
+    Vector, fromList, length, map, sum, unsafeFreeze, unsafeIndex
+  , unsafeTail, zip
   )
-import qualified Data.Vector.Unboxed.Mutable as MU (new, read, write)
+import qualified Data.Vector.Unboxed.Mutable as MU (
+    unsafeNew, unsafeRead, unsafeWrite
+  )
+import qualified Numeric.LinearAlgebra.Data as H (
+    (!), Matrix, Vector, diag, fromColumns, fromList, fromLists
+  , fromRows, konst, maxElement, maxIndex, toColumns, tr
+  )
+import qualified Numeric.LinearAlgebra.HMatrix as H (
+    (<>), (#>), sumElements
+  )
 
 type LogLikelihood = Double
 
@@ -43,43 +50,69 @@
 --   in order to simplify the calculation.
 data HMM' = HMM' { nStates'          :: Int -- ^ Number of states
                  , nOutputs'         :: Int -- ^ Number of outputs
-                 , initialStateDist' :: U.Vector Double
-                 , transitionDist'   :: V.Vector (U.Vector Double)
-                 , emissionDistT'    :: V.Vector (U.Vector Double)
+                 , initialStateDist' :: H.Vector Double
+                 , transitionDist'   :: H.Matrix Double
+                 , emissionDistT'    :: H.Matrix Double
                  }
 
 instance NFData HMM' where
-  rnf hmm' = rnf n `seq` rnf m `seq` rnf pi0 `seq` rnf w `seq` rnf phi'
+  rnf hmm' = rnf k `seq` rnf l `seq` rnf pi0 `seq` rnf w `seq` rnf phi'
     where
-      n    = nStates' hmm'
-      m    = nOutputs' hmm'
+      k    = nStates' hmm'
+      l    = nOutputs' hmm'
       pi0  = initialStateDist' hmm'
       w    = transitionDist' hmm'
       phi' = emissionDistT' hmm'
 
 init' :: Int -> Int -> RVar HMM'
-init' n m = do
-  pi0 <- U.fromList <$> stdSimplex (n-1)
-  w   <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (n-1))
-  phi <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (m-1))
-  return HMM' { nStates'          = n
-              , nOutputs'         = m
+init' k l = do
+  pi0 <- H.fromList <$> stdSimplex (k-1)
+  w   <- H.fromLists <$> replicateM k (stdSimplex (k-1))
+  phi <- H.fromLists <$> replicateM k (stdSimplex (l-1))
+  return HMM' { nStates'          = k
+              , nOutputs'         = l
               , initialStateDist' = pi0
               , transitionDist'   = w
-              , emissionDistT'    = G.transpose phi
+              , emissionDistT'    = H.tr phi
               }
 
 withEmission' :: HMM' -> U.Vector Int -> HMM'
-withEmission' model xs = model { emissionDistT' = phi' }
+withEmission' model xs = model'
   where
-    ss   = V.generate (nStates' model) id
-    os   = U.generate (nOutputs' model) id
-    phi' = let (path, _) = viterbi' model xs
-               freqs     = G.frequencies $ U.zip path xs
-               hists     = V.map (\s -> U.map (\o ->
-                                 M.findWithDefault 0 (s, o) freqs) os) ss
-           in V.map (\f -> f >/ U.sum f) hists
+    n = U.length xs
+    k = nStates' model
+    l = nOutputs' model
+    ss = [0..(k-1)]
+    os = [0..(l-1)]
 
+    step m = fst $ baumWelch1' (m { emissionDistT' = H.tr phi }) n xs
+      where
+        phi :: H.Matrix Double
+        phi = let zs  = fst $ viterbi' m xs
+                  fs  = G.frequencies $ U.zip zs xs
+                  hs  = H.fromLists $ map (\s -> map (\o ->
+                          M.findWithDefault 0 (s, o) fs) os) ss
+                  -- hs' is needed to not yield NaN vectors
+                  hs' = hs + H.konst 1e-9 (k, l)
+                  ns  = hs' H.#> H.konst 1 k
+              in hs' / H.fromColumns (replicate l ns)
+
+    ms  = iterate step model
+    ms' = tail ms
+    ds  = zipWith euclideanDistance ms ms'
+
+    model' = fst $ head $ dropWhile ((> 1e-9) . snd) $ zip ms' ds
+
+-- | Return the Euclidean distance between two models.
+euclideanDistance :: HMM' -> HMM' -> Double
+euclideanDistance model model' =
+  sqrt $ (H.sumElements $ (w - w') ** 2) + (H.sumElements $ (phi - phi') ** 2)
+  where
+    w    = transitionDist' model
+    w'   = transitionDist' model'
+    phi  = emissionDistT' model
+    phi' = emissionDistT' model'
+
 viterbi' :: HMM' -> U.Vector Int -> (U.Vector Int, LogLikelihood)
 viterbi' model xs = (path, logL)
   where
@@ -87,36 +120,39 @@
 
     -- First, we calculate the value function and the state maximizing it
     -- for each time.
-    deltas :: V.Vector (U.Vector Double)
+    deltas :: V.Vector (H.Vector Double)
     psis   :: V.Vector (U.Vector Int)
     (deltas, psis) = runST $ do
-      ds <- MV.new n
-      ps <- MV.new n
-      MV.write ds 0 $ U.map log (phi' V.! (xs U.! 0)) >+> U.map log pi0
-      MV.write ps 0 $ U.replicate k 0
+      ds <- MV.unsafeNew n
+      ps <- MV.unsafeNew n
+      let x0 = U.unsafeIndex xs 0
+      MV.unsafeWrite ds 0 $ log (phi' H.! x0) + log pi0
       forM_ [1..(n-1)] $ \t -> do
-        d <- MV.read ds (t-1)
-        let dws = V.map (\wj -> d >+> U.map log wj) w'
-        MV.write ds t $ U.map log (phi' V.! (xs U.! t)) >+> G.convert (V.map U.maximum dws)
-        MV.write ps t $ G.convert (V.map U.maxIndex dws)
-      ds' <- V.freeze ds
-      ps' <- V.freeze ps
+        d <- MV.unsafeRead ds (t-1)
+        let x   = U.unsafeIndex xs t
+            dws = map (\wj -> d + log wj) w'
+        MV.unsafeWrite ds t $ log (phi' H.! x) + H.fromList (map H.maxElement dws)
+        MV.unsafeWrite ps t $ U.fromList (map H.maxIndex dws)
+      ds' <- V.unsafeFreeze ds
+      ps' <- V.unsafeFreeze ps
       return (ds', ps')
       where
-        k    = nStates' model
         pi0  = initialStateDist' model
-        w'   = G.transpose $ transitionDist' model
+        w'   = H.toColumns $ transitionDist' model
         phi' = emissionDistT' model
 
-    -- The most likely path and corresponding log likelihood is as follows.
+    deltaE = V.unsafeIndex deltas (n-1)
+
+    -- The most likely path and corresponding log likelihood are as follows.
     path = runST $ do
-      ix <- MU.new n
-      MU.write ix (n-1) $ U.maxIndex (deltas V.! (n-1))
-      forM_ (reverse [0..(n-2)]) $ \t -> do
-        i <- MU.read ix (t+1)
-        MU.write ix t $ psis V.! (t+1) U.! i
-      U.freeze ix
-    logL = U.maximum $ deltas V.! (n-1)
+      ix <- MU.unsafeNew n
+      MU.unsafeWrite ix (n-1) $ H.maxIndex deltaE
+      forM_ [n-l | l <- [1..(n-1)]] $ \t -> do
+        i <- MU.unsafeRead ix t
+        let psi = V.unsafeIndex psis t
+        MU.unsafeWrite ix (t-1) $ U.unsafeIndex psi i
+      U.unsafeFreeze ix
+    logL = H.maxElement deltaE
 
 baumWelch' :: HMM' -> U.Vector Int -> [(HMM', LogLikelihood)]
 baumWelch' model xs = zip models (tail logLs)
@@ -130,6 +166,9 @@
 baumWelch1' :: HMM' -> Int -> U.Vector Int -> (HMM', LogLikelihood)
 baumWelch1' model n xs = force (model', logL)
   where
+    k = nStates' model
+    l = nOutputs' model
+
     -- First, we calculate the alpha, beta, and scaling values using the
     -- forward-backward algorithm.
     (alphas, cs) = forward' model n xs
@@ -142,16 +181,14 @@
     -- Using the gamma and xi values, we obtain the optimal initial state
     -- probability vector, transition probability matrix, and emission
     -- probability matrix.
-    pi0  = gammas V.! 0
-    w    = let ds = V.foldl1' (#+#) xis -- denominators
-               ns = V.map U.sum ds      -- numerators
-           in V.zipWith (>/) ds ns
+    pi0  = V.unsafeIndex gammas 0
+    w    = let ds = V.foldl1' (+) xis   -- denominators
+               ns = ds H.#> H.konst 1 k -- numerators
+           in H.diag (H.konst 1 k / ns) H.<> ds
     phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip (G.convert xs) gammas
-               ds    = V.foldl1' (>+>) . gs'  -- denominators
-               ns    = V.foldl1' (>+>) gammas -- numerators
-           in V.map (\o -> ds o >/> ns) os
-      where
-        os = V.generate (nOutputs' model) id
+               ds    = V.foldl1' (+) . gs'  -- denominators
+               ns    = V.foldl1' (+) gammas -- numerators
+           in H.fromRows $ map (\o -> ds o / ns) [0..(l-1)]
 
     -- We finally obtain the new model and the likelihood for the old model.
     model' = model { initialStateDist' = pi0
@@ -161,56 +198,61 @@
     logL = - (U.sum $ U.map log cs)
 
 -- | Return alphas and scaling variables.
-forward' :: HMM' -> Int -> U.Vector Int -> (V.Vector (U.Vector Double), U.Vector Double)
+forward' :: HMM' -> Int -> U.Vector Int -> (V.Vector (H.Vector Double), U.Vector Double)
 {-# INLINE forward' #-}
 forward' model n xs = runST $ do
-  as <- MV.new n
-  cs <- MU.new n
-  let a0 = (phi' V.! (xs U.! 0)) >.> pi0
-      c0 = 1 / U.sum a0
-  MV.write as 0 (c0 .> a0)
-  MU.write cs 0 c0
+  as <- MV.unsafeNew n
+  cs <- MU.unsafeNew n
+  let x0 = U.unsafeIndex xs 0
+      a0 = (phi' H.! x0) * pi0
+      c0 = 1 / H.sumElements a0
+  MV.unsafeWrite as 0 (H.konst c0 k * a0)
+  MU.unsafeWrite cs 0 c0
   forM_ [1..(n-1)] $ \t -> do
-    a <- MV.read as (t-1)
-    let a' = (phi' V.! (xs U.! t)) >.> (a <.# w)
-        c' = 1 / U.sum a'
-    MV.write as t (c' .> a')
-    MU.write cs t c'
-  as' <- V.freeze as
-  cs' <- U.freeze cs
+    a <- MV.unsafeRead as (t-1)
+    let x  = U.unsafeIndex xs t
+        a' = (phi' H.! x) * (w' H.#> a)
+        c' = 1 / H.sumElements a'
+    MV.unsafeWrite as t (H.konst c' k * a')
+    MU.unsafeWrite cs t c'
+  as' <- V.unsafeFreeze as
+  cs' <- U.unsafeFreeze cs
   return (as', cs')
   where
+    k    = nStates' model
     pi0  = initialStateDist' model
-    w    = transitionDist' model
+    w'   = H.tr $ transitionDist' model
     phi' = emissionDistT' model
 
 -- | Return betas using scaling variables.
-backward' :: HMM' -> Int -> U.Vector Int -> U.Vector Double -> V.Vector (U.Vector Double)
+backward' :: HMM' -> Int -> U.Vector Int -> U.Vector Double -> V.Vector (H.Vector Double)
 {-# INLINE backward' #-}
 backward' model n xs cs = runST $ do
-  bs <- MV.new n
-  let bE = U.replicate k 1
-      cE = cs U.! (n-1)
-  MV.write bs (n-1) $ cE .> bE
-  forM_ (reverse [0..(n-2)]) $ \t -> do
-    b <- MV.read bs (t+1)
-    let b' = w #.> ((phi' V.! (xs U.! (t+1))) >.> b)
-        c' = cs U.! t
-    MV.write bs t $ c' .> b'
-  V.freeze bs
+  bs <- MV.unsafeNew n
+  let bE = H.konst 1 k
+      cE = U.unsafeIndex cs (n-1)
+  MV.unsafeWrite bs (n-1) (H.konst cE k * bE)
+  forM_ [n-l | l <- [1..(n-1)]] $ \t -> do
+    b <- MV.unsafeRead bs t
+    let x  = U.unsafeIndex xs t
+        b' = w H.#> ((phi' H.! x) * b)
+        c' = U.unsafeIndex cs (t-1)
+    MV.unsafeWrite bs (t-1) (H.konst c' k * b')
+  V.unsafeFreeze bs
   where
     k    = nStates' model
     w    = transitionDist' model
     phi' = emissionDistT' model
 
 -- | Return the posterior distribution.
-posterior' :: HMM' -> Int -> U.Vector Int -> V.Vector (U.Vector Double) -> V.Vector (U.Vector Double) -> U.Vector Double -> (V.Vector (U.Vector Double), V.Vector (V.Vector (U.Vector Double)))
+posterior' :: HMM' -> Int -> U.Vector Int -> V.Vector (H.Vector Double) -> V.Vector (H.Vector Double) -> U.Vector Double -> (V.Vector (H.Vector Double), V.Vector (H.Matrix Double))
 {-# INLINE posterior' #-}
 posterior' model _ xs alphas betas cs = (gammas, xis)
   where
-    gammas = V.zipWith3 (\a b c -> a >.> b >/ c) alphas betas (G.convert cs)
-    xis    = V.zipWith3 (\a b x -> let w' = V.zipWith (.>) (G.convert a) w
-                                   in V.map ((phi' V.! x) >.> b >.>) w')
-               alphas (V.tail betas) (G.convert $ U.tail xs)
+    gammas = V.zipWith3 (\a b c -> a * b / H.konst c k)
+               alphas betas (G.convert cs)
+    xis    = V.zipWith3 (\a b x -> H.diag a H.<> w H.<> H.diag (b * (phi' H.! x)))
+               alphas (V.unsafeTail betas) (G.convert $ U.unsafeTail xs)
+    k    = nStates' model
     w    = transitionDist' model
     phi' = emissionDistT' model
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -4,6 +4,5 @@
 
 main :: IO ()
 main = doctest [ "-isrc"
-               , "src/Data/Vector/Generic/Util/LinearAlgebra.hs"
                , "src/Learning/HMM.hs"
                ]