diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 36c60d719cbc7b0cca5b1bd0d5067c383c4462ec..7a9995d1228b19f99ee420d3f4cc61a24848d4b4 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -380,6 +380,10 @@ public:
     /** \brief Assemble the local stiffness matrix at the current position
 
     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
+
     */
     virtual void assembleHessian(const Entity& e,
                   const std::vector<TargetSpace>& localSolution);
@@ -454,7 +458,7 @@ public:
         }
 
     }
-
+    
     // assembled data
     Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > A_;
     
@@ -496,43 +500,65 @@ assembleHessian(const Entity& element,
 
     A_ = 0;
 
-    // first compute the Hessian in the embedding space
-    Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
 
-    std::vector<Dune::FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > embeddedGradient;
+    
+    
+    const double eps = 1e-4;
+        
+    std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
+    for (size_t i=0; i<B.size(); i++)
+        B[i] = localSolution[i].orthonormalFrame();
 
+        
+    // finite-difference approximation
     for (size_t i=0; i<localSolution.size(); i++) {
+        for (size_t i2=0; i2<blocksize; i2++) {
+            for (size_t j=0; j<localSolution.size(); j++) {
+                for (size_t j2=0; j2<blocksize; j2++) {
 
-        embeddedGradientOfEmbeddedGradient(element,localSolution, i, embeddedGradient);
-
-        for (size_t j=0; j<localSolution.size(); j++)
-            embeddedHessian[i][j] = embeddedGradient[j];
-
-    }
-
-    // transform to local tangent space bases
-    std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > orthonormalFrames(nDofs);
-    std::vector<Dune::FieldMatrix<double,embeddedBlocksize,blocksize> > orthonormalFramesTransposed(nDofs);
-
-    for (size_t i=0; i<nDofs; i++) {
-        orthonormalFrames[i] = localSolution[i].orthonormalFrame();
-
-        for (int j=0; j<embeddedBlocksize; j++)
-            for (int k=0; k<blocksize; k++)
-                orthonormalFramesTransposed[i][j][k] = orthonormalFrames[i][k][j];
-
-    }
-
-    for (size_t i=0; i<nDofs; i++)
-        for (size_t j=0; j<nDofs; j++) {
-
-            Dune::FieldMatrix<double,blocksize,embeddedBlocksize> tmp;
-            Dune::FMatrixHelp::multMatrix(orthonormalFrames[i],embeddedHessian[i][j],tmp);
-            A_[i][j] = tmp.rightmultiplyany(orthonormalFramesTransposed[j]);
+                    Dune::FieldVector<double,embeddedBlocksize> epsXi  = B[i][i2];    epsXi *= eps;
+                    Dune::FieldVector<double,embeddedBlocksize> epsEta = B[j][j2];   epsEta *= eps;
+            
+                    Dune::FieldVector<double,embeddedBlocksize> minusEpsXi  = epsXi;   minusEpsXi  *= -1;
+                    Dune::FieldVector<double,embeddedBlocksize> minusEpsEta = epsEta;  minusEpsEta *= -1;
+                        
+                    std::vector<TargetSpace> forwardSolutionXiEta  = localSolution;
+                    std::vector<TargetSpace> forwardSolutionXi     = localSolution;
+                    std::vector<TargetSpace> forwardSolutionEta    = localSolution;
+                    std::vector<TargetSpace> backwardSolutionXiEta  = localSolution;
+                    std::vector<TargetSpace> backwardSolutionXi     = localSolution;
+                    std::vector<TargetSpace> backwardSolutionEta    = localSolution;
             
+                    if (i==j)
+                        forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
+                    else {
+                        forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi);
+                        forwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],epsEta);
+                    }
+                    forwardSolutionXi[i]    = TargetSpace::exp(localSolution[i],epsXi);
+                    forwardSolutionEta[j]   = TargetSpace::exp(localSolution[j],epsEta);
+                        
+                    if (i==j)
+                        backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi+minusEpsEta);
+                    else {
+                        backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
+                        backwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],minusEpsEta);
+                    }
+                        
+                    backwardSolutionXi[i]    = TargetSpace::exp(localSolution[i],minusEpsXi);
+                    backwardSolutionEta[j]   = TargetSpace::exp(localSolution[j],minusEpsEta);
 
+                    double forwardValue  = energy(element, forwardSolutionXiEta) - energy(element, forwardSolutionXi) - energy(element, forwardSolutionEta);
+                    double centerValue   = -energy(element, localSolution);
+                    double backwardValue = energy(element, backwardSolutionXiEta) - energy(element, backwardSolutionXi) - energy(element, backwardSolutionEta);
+            
+                    A_[i][j][i2][j2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
+                        
+                }
+            }
         }
-    
+    }
+        
 }
 
 #endif