From ab9dedc3434b60b998fa8dd3d202946f296681ca Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 9 Jul 2019 11:09:57 +0200
Subject: [PATCH] Compare gradient and Hesse matrix with an FD approximation
Previously the code only tested whether the AverageDistanceAssembler
computed gradient vector and Hesse matrix. Now the test actually
tests whether these are correct.
---
test/averagedistanceassemblertest.cc | 111 ++++++++++++++++++++++++---
1 file changed, 101 insertions(+), 10 deletions(-)
diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc
index 863b8d0a..983b9707 100644
--- a/test/averagedistanceassemblertest.cc
+++ b/test/averagedistanceassemblertest.cc
@@ -17,6 +17,84 @@ const int dim = 2;
using namespace Dune;
+// Compute FD approximations to the gradient and the Hesse matrix
+template <class DistanceAssembler, class TargetSpace>
+void assembleGradientAndHessianApproximation(const DistanceAssembler& assembler,
+ const TargetSpace& argument,
+ typename TargetSpace::TangentVector& gradient,
+ FieldMatrix<double, TargetSpace::TangentVector::dimension, TargetSpace::TangentVector::dimension>& hesseMatrix)
+{
+ using field_type = typename TargetSpace::field_type;
+ constexpr auto blocksize = TargetSpace::TangentVector::dimension;
+ constexpr auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
+
+ const field_type eps = 1e-4;
+
+ FieldMatrix<double,blocksize,embeddedBlocksize> B = argument.orthonormalFrame();
+
+ // Precompute negative energy at the current configuration
+ // (negative because that is how we need it as part of the 2nd-order fd formula)
+ field_type centerValue = -assembler.value(argument);
+
+ // Precompute energy infinitesimal corrections in the directions of the local basis vectors
+ std::array<field_type,blocksize> forwardEnergy;
+ std::array<field_type,blocksize> backwardEnergy;
+
+ for (size_t i2=0; i2<blocksize; i2++)
+ {
+ typename TargetSpace::EmbeddedTangentVector epsXi = B[i2];
+ epsXi *= eps;
+ typename TargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
+ minusEpsXi *= -1;
+
+ TargetSpace forwardSolution = argument;
+ TargetSpace backwardSolution = argument;
+
+ forwardSolution = TargetSpace::exp(argument,epsXi);
+ backwardSolution = TargetSpace::exp(argument,minusEpsXi);
+
+ forwardEnergy[i2] = assembler.value(forwardSolution);
+ backwardEnergy[i2] = assembler.value(backwardSolution);
+ }
+
+ //////////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ //////////////////////////////////////////////////////////////
+
+ for (int j=0; j<blocksize; j++)
+ gradient[j] = (forwardEnergy[j] - backwardEnergy[j]) / (2*eps);
+
+ ///////////////////////////////////////////////////////////////////////////
+ // Compute Riemannian Hesse matrix by finite-difference approximation.
+ // We loop over the lower left triangular half of the matrix.
+ // The other half follows from symmetry.
+ ///////////////////////////////////////////////////////////////////////////
+ for (size_t i2=0; i2<blocksize; i2++)
+ {
+ for (size_t j2=0; j2<i2+1; j2++)
+ {
+ TargetSpace forwardSolutionXiEta = argument;
+ TargetSpace backwardSolutionXiEta = argument;
+
+ typename TargetSpace::EmbeddedTangentVector epsXi = B[i2]; epsXi *= eps;
+ typename TargetSpace::EmbeddedTangentVector epsEta = B[j2]; epsEta *= eps;
+
+ typename TargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1;
+ typename TargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1;
+
+ forwardSolutionXiEta = TargetSpace::exp(argument, epsXi+epsEta);
+ backwardSolutionXiEta = TargetSpace::exp(argument, minusEpsXi+minusEpsEta);
+
+ field_type forwardValue = assembler.value(forwardSolutionXiEta) - forwardEnergy[i2] - forwardEnergy[j2];
+ field_type backwardValue = assembler.value(backwardSolutionXiEta) - backwardEnergy[i2] - backwardEnergy[j2];
+
+ hesseMatrix[i2][j2] = hesseMatrix[j2][i2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
+ }
+ }
+}
+
+
+
/** \brief Test whether interpolation is invariant under permutation of the simplex vertices
*/
template <class TargetSpace>
@@ -35,26 +113,39 @@ void testPoint(const std::vector<TargetSpace>& corners,
// test the gradient
typename TargetSpace::TangentVector gradient;
assembler.assembleGradient(argument, gradient);
+ typename TargetSpace::TangentVector gradientApproximation;
// test the hessian
FieldMatrix<double, TargetSpace::TangentVector::dimension, TargetSpace::TangentVector::dimension> hessian;
FieldMatrix<double, TargetSpace::TangentVector::dimension, TargetSpace::TangentVector::dimension> hessianApproximation(0);
assembler.assembleHessian(argument, hessian);
- //assembler.assembleHessianApproximation(argument, hessianApproximation);
+ assembleGradientAndHessianApproximation(assembler, argument,
+ gradientApproximation, hessianApproximation);
+
+ // Check gradient
+ for (size_t i=0; i<gradient.size(); i++)
+ {
+ if (std::isnan(gradient[i]))
+ DUNE_THROW(Dune::Exception, "Gradient contains NaN");
+ if (std::isnan(gradientApproximation[i]))
+ DUNE_THROW(Dune::Exception, "Gradient approximation contains NaN");
+ if (std::abs(gradient[i] - gradientApproximation[i]) > 1e-6)
+ DUNE_THROW(Dune::Exception, "Gradient and its approximation do not match");
+ }
+ // Check Hesse matrix
for (size_t i=0; i<hessian.N(); i++)
- for (size_t j=0; j<hessian.M(); j++) {
- assert(!std::isnan(hessian[i][j]));
- assert(!std::isnan(hessianApproximation[i][j]));
+ for (size_t j=0; j<hessian.M(); j++)
+ {
+ if (std::isnan(hessian[i][j]))
+ DUNE_THROW(Dune::Exception, "Hesse matrix contains NaN");
+ if (std::isnan(hessianApproximation[i][j]))
+ DUNE_THROW(Dune::Exception, "Hesse matrix approximation contains NaN");
+ if (std::abs(hessian[i][j] - hessianApproximation[i][j]) > 1e-6)
+ DUNE_THROW(Dune::Exception, "Hesse matrix and its approximation do not match");
}
- std::cout << "WARNING: no approximation of the Hessian available, not testing" << std::endl;
- exit(1);
-
- FieldMatrix<double, TargetSpace::TangentVector::dimension, TargetSpace::TangentVector::dimension> diff = hessian;
- diff -= hessianApproximation;
- std::cout << "Matrix inf diff: " << diff.infinity_norm() << std::endl;
}
--
GitLab