diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index b91bcc2081d8ec4eb61b1c8dabd8325961b2a397..7e332e90723ca6e318917f1d9457cf99bdca51b2 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -146,7 +146,12 @@ public:
     This default implementation used finite-difference approximations to compute the second derivatives
 
     The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
-    'Optimization algorithms on matrix manifolds', page 107
+    'Optimization algorithms on matrix manifolds', page 107.  There it says that
+    \f[
+        \langle Hess f(x)[\xi], \eta \rangle
+            = \frac 12 \frac{d^2}{dt^2} \Big(f(\exp_x(t(\xi + \eta))) - f(\exp_x(t\xi)) - f(\exp_x(t\eta))\Big)\Big|_{t=0}.
+    \f]
+    We compute that using a finite difference approximation.
 
     */
     virtual void assembleHessian(const Entity& e,