diff --git a/Numeric/AD.hs b/Numeric/AD.hs
--- a/Numeric/AD.hs
+++ b/Numeric/AD.hs
@@ -16,12 +16,16 @@
 module Numeric.AD
     (
     -- * Gradients (Reverse Mode)
-      grad, grad'
-    , gradWith, gradWith'
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
 
     -- * Jacobians (Mixed Mode)
-    , jacobian, jacobian'
-    , jacobianWith, jacobianWith'
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
 
     -- * Jacobians (Reverse Mode)
     , gradF
@@ -30,9 +34,16 @@
     , gradWithF'
 
     -- * Jacobians (Forward Mode)
-    , jacobianT, jacobianWithT
+    , jacobianT
+    , jacobianWithT
 
-    -- * Hessians (Forward-On-Reverse Mode)
+    -- * Hessian (Forward-On-Reverse)
+    , hessian
+
+    -- * Hessian Tensors (Forward-On-Mixed)
+    , hessianTensor
+
+    -- * Hessian Vector Products (Forward-On-Reverse)
     , hessianProduct
     , hessianProduct'
 
@@ -87,7 +98,7 @@
 import Numeric.AD.Forward  (diff, diff', diffF, diffF', du, du', duF, duF', diffM, diffM', jacobianT, jacobianWithT) 
 import Numeric.AD.Tower    (diffsF, diffs0F , diffs, diffs0, taylor, taylor0, maclaurin, maclaurin0)
 import Numeric.AD.Reverse  (grad, grad', gradWith, gradWith', gradM, gradM', gradWithM, gradWithM', gradF, gradF', gradWithF, gradWithF')
-import Numeric.AD.Internal.Composition (compose, decompose)
+import Numeric.AD.Internal.Composition
 
 import qualified Numeric.AD.Forward as Forward
 import qualified Numeric.AD.Reverse as Reverse
@@ -147,7 +158,7 @@
 -- 
 -- Or in other words, we take the directional derivative of the gradient.
 hessianProduct :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f a
-hessianProduct f = duF (grad (decompose . f . fmap compose))
+hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))
 
 -- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' <$> wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\". The outputs are returned wrapped in the same functor.
 --
@@ -156,5 +167,14 @@
 -- Or in other words, we take the directional derivative of the gradient.
 -- 
 hessianProduct' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f (a, a)
-hessianProduct' f = duF' (grad (decompose . f . fmap compose))
+hessianProduct' f = duF' (grad (decomposeMode . f . fmap composeMode))
 
+-- hessianProductWith' :: (Traversable f, Num a) => (a -> a -> a -> a -> b) -> (forall s. Mode s. f (AD s a) -> AD s a) -> f (a, a) -> f b
+
+-- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in forward mode.
+hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
+hessian f = Forward.jacobian (grad (decomposeMode . f . fmap composeMode))
+
+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the forward-mode Jacobian of the mixed-mode Jacobian of the function.
+hessianTensor :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
+hessianTensor f = decomposeFunctor . Forward.jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
diff --git a/Numeric/AD/Forward.hs b/Numeric/AD/Forward.hs
--- a/Numeric/AD/Forward.hs
+++ b/Numeric/AD/Forward.hs
@@ -166,8 +166,8 @@
 
 -- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD. 
 hessianProduct :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f a
-hessianProduct f = duF (grad (decompose . f . fmap compose))
+hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))
 
 -- | Compute the gradient and hessian product using forward-on-forward-mode AD. 
 hessianProduct' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f (a, a)
-hessianProduct' f = duF' (grad (decompose . f . fmap compose))
+hessianProduct' f = duF' (grad (decomposeMode . f . fmap composeMode))
diff --git a/Numeric/AD/Internal/Composition.hs b/Numeric/AD/Internal/Composition.hs
--- a/Numeric/AD/Internal/Composition.hs
+++ b/Numeric/AD/Internal/Composition.hs
@@ -12,98 +12,111 @@
 -----------------------------------------------------------------------------
 
 module Numeric.AD.Internal.Composition
-    ( (:.)(..)
-    , On
-    , compose
-    , decompose
+    ( ComposeFunctor(..)
+    , ComposeMode(..)
+    , composeMode
+    , decomposeMode
     ) where
 
+import Data.Traversable
+import Control.Applicative
+import Data.Foldable
 import Numeric.AD.Classes
 import Numeric.AD.Internal
 
-newtype (f :. g) a = O { runO :: f (AD g a) }
+-- * Functor composition
 
+newtype ComposeFunctor f g a = ComposeFunctor { decomposeFunctor :: f (g a) }
 
-type On f g = g :. f
+instance (Functor f, Functor g) => Functor (ComposeFunctor f g) where
+    fmap f (ComposeFunctor a) = ComposeFunctor (fmap (fmap f) a)
+    
+instance (Foldable f, Foldable g) => Foldable (ComposeFunctor f g) where
+    foldMap f (ComposeFunctor a) = foldMap (foldMap f) a
+    
+instance (Traversable f, Traversable g) => Traversable (ComposeFunctor f g) where
+    traverse f (ComposeFunctor a) = ComposeFunctor <$> traverse (traverse f) a
 
-compose :: AD f (AD g a) -> AD (f :. g) a
-compose (AD a) = AD (O a)
+newtype ComposeMode f g a = ComposeMode { runComposeMode :: f (AD g a) }
 
-decompose :: AD (f :. g) a -> AD f (AD g a)
-decompose (AD (O a)) = AD a
+composeMode :: AD f (AD g a) -> AD (ComposeMode f g) a
+composeMode (AD a) = AD (ComposeMode a)
 
-instance (Primal f, Mode g, Primal g) => Primal (f :. g) where
-    primal = primal . primal . runO
+decomposeMode :: AD (ComposeMode f g) a -> AD f (AD g a)
+decomposeMode (AD (ComposeMode a)) = AD a
 
-instance (Mode f, Mode g) => Mode (f :. g) where
-    lift = O . lift . lift
-    O a <+> O b = O (a <+> b) 
-    a *^ O b = O (lift a *^ b) 
-    O a ^* b = O (a ^* lift b)
-    O a ^/ b = O (a ^/ lift b)
+instance (Primal f, Mode g, Primal g) => Primal (ComposeMode f g) where
+    primal = primal . primal . runComposeMode
 
-instance (Mode f, Mode g) => Lifted (f :. g) where
-    showsPrec1 n (O a) = showsPrec1 n a
-    O a ==! O b  = a ==! b
-    compare1 (O a) (O b) = compare1 a b
-    fromInteger1 = O . lift . fromInteger1
-    O a +! O b = O (a +! b)
-    O a -! O b = O (a -! b)
-    O a *! O b = O (a *! b)
-    negate1 (O a) = O (negate1 a)
-    abs1 (O a) = O (abs1 a)
-    signum1 (O a) = O (signum1 a)
-    O a /! O b = O (a /! b) 
-    recip1 (O a) = O (recip1 a)
-    fromRational1 = O . lift . fromRational1
-    toRational1 (O a) = toRational1 a
-    pi1 = O pi1
-    exp1 (O a) = O (exp1 a)
-    log1 (O a) = O (log1 a) 
-    sqrt1 (O a) = O (sqrt1 a)
-    O a **! O b = O (a **! b)
-    logBase1 (O a) (O b) = O (logBase1 a b)
-    sin1 (O a) = O (sin1 a)
-    cos1 (O a) = O (cos1 a)
-    tan1 (O a) = O (tan1 a)
-    asin1 (O a) = O (asin1 a)
-    acos1 (O a) = O (acos1 a)
-    atan1 (O a) = O (atan1 a)
-    sinh1 (O a) = O (sinh1 a)
-    cosh1 (O a) = O (cosh1 a)
-    tanh1 (O a) = O (tanh1 a)
-    asinh1 (O a) = O (asinh1 a)
-    acosh1 (O a) = O (acosh1 a)
-    atanh1 (O a) = O (atanh1 a)
-    properFraction1 (O a) = (b, O c) where
+instance (Mode f, Mode g) => Mode (ComposeMode f g) where
+    lift = ComposeMode . lift . lift
+    ComposeMode a <+> ComposeMode b = ComposeMode (a <+> b) 
+    a *^ ComposeMode b = ComposeMode (lift a *^ b) 
+    ComposeMode a ^* b = ComposeMode (a ^* lift b)
+    ComposeMode a ^/ b = ComposeMode (a ^/ lift b)
+
+instance (Mode f, Mode g) => Lifted (ComposeMode f g) where
+    showsPrec1 n (ComposeMode a) = showsPrec1 n a
+    ComposeMode a ==! ComposeMode b  = a ==! b
+    compare1 (ComposeMode a) (ComposeMode b) = compare1 a b
+    fromInteger1 = ComposeMode . lift . fromInteger1
+    ComposeMode a +! ComposeMode b = ComposeMode (a +! b)
+    ComposeMode a -! ComposeMode b = ComposeMode (a -! b)
+    ComposeMode a *! ComposeMode b = ComposeMode (a *! b)
+    negate1 (ComposeMode a) = ComposeMode (negate1 a)
+    abs1 (ComposeMode a) = ComposeMode (abs1 a)
+    signum1 (ComposeMode a) = ComposeMode (signum1 a)
+    ComposeMode a /! ComposeMode b = ComposeMode (a /! b) 
+    recip1 (ComposeMode a) = ComposeMode (recip1 a)
+    fromRational1 = ComposeMode . lift . fromRational1
+    toRational1 (ComposeMode a) = toRational1 a
+    pi1 = ComposeMode pi1
+    exp1 (ComposeMode a) = ComposeMode (exp1 a)
+    log1 (ComposeMode a) = ComposeMode (log1 a) 
+    sqrt1 (ComposeMode a) = ComposeMode (sqrt1 a)
+    ComposeMode a **! ComposeMode b = ComposeMode (a **! b)
+    logBase1 (ComposeMode a) (ComposeMode b) = ComposeMode (logBase1 a b)
+    sin1 (ComposeMode a) = ComposeMode (sin1 a)
+    cos1 (ComposeMode a) = ComposeMode (cos1 a)
+    tan1 (ComposeMode a) = ComposeMode (tan1 a)
+    asin1 (ComposeMode a) = ComposeMode (asin1 a)
+    acos1 (ComposeMode a) = ComposeMode (acos1 a)
+    atan1 (ComposeMode a) = ComposeMode (atan1 a)
+    sinh1 (ComposeMode a) = ComposeMode (sinh1 a)
+    cosh1 (ComposeMode a) = ComposeMode (cosh1 a)
+    tanh1 (ComposeMode a) = ComposeMode (tanh1 a)
+    asinh1 (ComposeMode a) = ComposeMode (asinh1 a)
+    acosh1 (ComposeMode a) = ComposeMode (acosh1 a)
+    atanh1 (ComposeMode a) = ComposeMode (atanh1 a)
+    properFraction1 (ComposeMode a) = (b, ComposeMode c) where
         (b, c) = properFraction1 a
-    truncate1 (O a) = truncate1 a
-    round1 (O a) = round1 a
-    ceiling1 (O a) = ceiling1 a
-    floor1 (O a) = floor1 a
-    floatRadix1 (O a) = floatRadix1 a
-    floatDigits1 (O a) = floatDigits1 a
-    floatRange1 (O a) = floatRange1 a
-    decodeFloat1 (O a) = decodeFloat1 a
-    encodeFloat1 m e = O (encodeFloat1 m e)
-    exponent1 (O a) = exponent1 a
-    significand1 (O a) = O (significand1 a)
-    scaleFloat1 n (O a) = O (scaleFloat1 n a)
-    isNaN1 (O a) = isNaN1 a 
-    isInfinite1 (O a) = isInfinite1 a
-    isDenormalized1 (O a) = isDenormalized1 a
-    isNegativeZero1 (O a) = isNegativeZero1 a
-    isIEEE1 (O a) = isIEEE1 a
-    atan21 (O a) (O b) = O (atan21 a b)
-    succ1 (O a) = O (succ1 a)
-    pred1 (O a) = O (pred1 a)
-    toEnum1 n = O (toEnum1 n)
-    fromEnum1 (O a) = fromEnum1 a
-    enumFrom1 (O a) = map O $ enumFrom1 a
-    enumFromThen1 (O a) (O b) = map O $ enumFromThen1 a b
-    enumFromTo1 (O a) (O b) = map O $ enumFromTo1 a b
-    enumFromThenTo1 (O a) (O b) (O c) = map O $ enumFromThenTo1 a b c
-    minBound1 = O minBound1
-    maxBound1 = O maxBound1
+    truncate1 (ComposeMode a) = truncate1 a
+    round1 (ComposeMode a) = round1 a
+    ceiling1 (ComposeMode a) = ceiling1 a
+    floor1 (ComposeMode a) = floor1 a
+    floatRadix1 (ComposeMode a) = floatRadix1 a
+    floatDigits1 (ComposeMode a) = floatDigits1 a
+    floatRange1 (ComposeMode a) = floatRange1 a
+    decodeFloat1 (ComposeMode a) = decodeFloat1 a
+    encodeFloat1 m e = ComposeMode (encodeFloat1 m e)
+    exponent1 (ComposeMode a) = exponent1 a
+    significand1 (ComposeMode a) = ComposeMode (significand1 a)
+    scaleFloat1 n (ComposeMode a) = ComposeMode (scaleFloat1 n a)
+    isNaN1 (ComposeMode a) = isNaN1 a 
+    isInfinite1 (ComposeMode a) = isInfinite1 a
+    isDenormalized1 (ComposeMode a) = isDenormalized1 a
+    isNegativeZero1 (ComposeMode a) = isNegativeZero1 a
+    isIEEE1 (ComposeMode a) = isIEEE1 a
+    atan21 (ComposeMode a) (ComposeMode b) = ComposeMode (atan21 a b)
+    succ1 (ComposeMode a) = ComposeMode (succ1 a)
+    pred1 (ComposeMode a) = ComposeMode (pred1 a)
+    toEnum1 n = ComposeMode (toEnum1 n)
+    fromEnum1 (ComposeMode a) = fromEnum1 a
+    enumFrom1 (ComposeMode a) = map ComposeMode $ enumFrom1 a
+    enumFromThen1 (ComposeMode a) (ComposeMode b) = map ComposeMode $ enumFromThen1 a b
+    enumFromTo1 (ComposeMode a) (ComposeMode b) = map ComposeMode $ enumFromTo1 a b
+    enumFromThenTo1 (ComposeMode a) (ComposeMode b) (ComposeMode c) = map ComposeMode $ enumFromThenTo1 a b c
+    minBound1 = ComposeMode minBound1
+    maxBound1 = ComposeMode maxBound1
 
 -- deriveNumeric (conT `appT` varT (mkName "f") `appT` varT (mkName "g")) 
diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs
--- a/Numeric/AD/Newton.hs
+++ b/Numeric/AD/Newton.hs
@@ -106,11 +106,11 @@
 --
 -- > take 10 $ extremum cos 1 -- convert to 0 
 extremum :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
-extremum f = findZero (diff (decompose . f . compose))
+extremum f = findZero (diff (decomposeMode . f . composeMode))
 {-# INLINE extremum #-}
 
 extremumM :: (Monad m, Fractional a) => (forall s. Mode s => AD s a -> m (AD s a)) -> a -> MList m a
-extremumM f = findZeroM (diffM (liftM decompose . f . compose))
+extremumM f = findZeroM (diffM (liftM decomposeMode . f . composeMode))
 {-# INLINE extremumM #-}
 
 -- | The 'gradientDescent' function performs a multivariate
diff --git a/Numeric/AD/Reverse.hs b/Numeric/AD/Reverse.hs
--- a/Numeric/AD/Reverse.hs
+++ b/Numeric/AD/Reverse.hs
@@ -29,6 +29,11 @@
     , jacobian'
     , jacobianWith
     , jacobianWith'
+    -- * Hessian
+    , hessian
+    , hessianM
+    , hessianTensor
+    
     -- * Derivatives
     , diff
     , diff'
@@ -57,6 +62,7 @@
 
 import Numeric.AD.Classes
 import Numeric.AD.Internal
+import Numeric.AD.Internal.Composition
 import Numeric.AD.Internal.Reverse
 
 -- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
@@ -193,3 +199,20 @@
 gradWithM' :: (Traversable f, Monad m, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> m (AD s a)) -> f a -> m (a, f b)
 gradWithM' g f = unwrapMonad . jacobianWith' g (WrapMonad . f)
 
+-- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.
+--
+-- However, since the @'grad f :: f a -> f a'@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian' in "Numeric.AD".
+hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
+hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))
+
+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the forward-mode Jacobian of the mixed-mode Jacobian of the function.
+--
+-- While this is less efficient than 'Numeric.AD.hessianTensor' from "Numeric.AD" or 'Numeric.AD.Forward.hessianTensor' from "Numeric.AD.Forward", the type signature is more permissive with regards to the output non-scalar, and it may be more efficient if only a few coefficients of the result are consumed.
+hessianTensor :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
+hessianTensor f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
+
+-- | Compute the hessian via the reverse-mode jacobian of the reverse-mode gradient of a non-scalar-to-scalar monadic action. 
+--
+-- While this is less efficient than 'Numeric.AD.hessianTensor' from "Numeric.AD" or 'Numeric.AD.Forward.hessianTensor' from "Numeric.AD.Forward", the type signature is more permissive with regards to the output non-scalar, and it may be more efficient if only a few coefficients of the result are consumed.
+hessianM :: (Traversable f, Monad m, Num a) => (forall s. Mode s => f (AD s a) -> m (AD s a)) -> f a -> m (f (f a))
+hessianM f = unwrapMonad . hessianTensor (WrapMonad . f)
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 Name:         ad
-Version:      0.21
+Version:      0.22
 License:      BSD3
 License-File: LICENSE
 Copyright:    Edward Kmett 2010