diff --git a/Numeric/AD.hs b/Numeric/AD.hs
--- a/Numeric/AD.hs
+++ b/Numeric/AD.hs
@@ -32,6 +32,10 @@
     -- * Jacobians (Forward Mode)
     , jacobianT, jacobianWithT
 
+    -- * Hessians (Forward-On-Reverse Mode)
+    , hessianProduct
+    , hessianProduct'
+
     -- * Derivatives (Forward Mode)
     , diff
     , diffF
@@ -83,6 +87,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 qualified Numeric.AD.Forward as Forward
 import qualified Numeric.AD.Reverse as Reverse
@@ -135,4 +140,21 @@
         size :: Foldable f => f a -> Int
         size = foldr' (\_ b -> 1 + b) 0
 {-# INLINE jacobianWith' #-}
+
+-- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:
+--
+-- > H v = (d/dr) grad_w (w + r v) | r = 0
+-- 
+-- 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 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.
+--
+-- > H v = (d/dr) grad_w (w + r v) | r = 0
+-- 
+-- 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))
 
diff --git a/Numeric/AD/Forward.hs b/Numeric/AD/Forward.hs
--- a/Numeric/AD/Forward.hs
+++ b/Numeric/AD/Forward.hs
@@ -27,6 +27,9 @@
     -- * Transposed Jacobian
     , jacobianT
     , jacobianWithT
+    -- * Hessian Product
+    , hessianProduct
+    , hessianProduct'
     -- * Derivatives
     , diff
     , diff'
@@ -50,6 +53,7 @@
 import Control.Monad (liftM)
 import Numeric.AD.Classes
 import Numeric.AD.Internal
+import Numeric.AD.Internal.Composition
 import Numeric.AD.Internal.Forward
 
 du :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> a
@@ -159,3 +163,11 @@
 gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
 gradWith' g f = bindWith' g (tangent . f)
 {-# INLINE gradWith' #-}
+
+-- | 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))
+
+-- | 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))
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 Name:         ad
-Version:      0.20
+Version:      0.21
 License:      BSD3
 License-File: LICENSE
 Copyright:    Edward Kmett 2010