Directional derivatives only for unit vectors

src/bisection-example-flexible.cc

@@ -34,12 +34,11 @@ class SampleFunctional {
     return y * v + func_(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
   }
 
-  double directionalDerivative(const SmallVector x,
-                               const SmallVector dir) const {
+  // |dir|-times the directional derivative wrt dir/|dir|.
+  double pseudoDirectionalDerivative(const SmallVector x,
+                                     const SmallVector dir) const {
     if (x == SmallVector(0.0))
-      // FIXME: Well, not in this way -- but can we compute them?
-      DUNE_THROW(Dune::Exception,
-                 "Directional derivatives cannot be computed at zero.");
+      return func_.rightDifferential(0) * dir.two_norm();
 
     if (x * dir > 0)
       return PlusGrad(x) * dir;
@@ -47,6 +46,20 @@ class SampleFunctional {
       return MinusGrad(x) * dir;
   }
 
+  double directionalDerivative(const SmallVector x,
+                               const SmallVector dir) const {
+    double norm = dir.two_norm();
+
+    if (norm == 0)
+      DUNE_THROW(Dune::Exception, "Directional derivatives cannot be computed "
+                                  "w.r.t. the zero-direction.");
+
+    SmallVector tmp = dir;
+    tmp /= norm;
+
+    return pseudoDirectionalDerivative(x, tmp);
+  }
+
   SmallVector minimise(const SmallVector x, unsigned int iterations) const {
     SmallVector descDir = ModifiedGradient(x);
     if (descDir == SmallVector(0.0))
@@ -160,6 +173,7 @@ class SampleFunctional {
     return ret;
   }
 
+  // No normalising is done!
   SmallVector project(const SmallVector z, const SmallVector x) const {
     SmallVector y = z;
     y.axpy(-(z * x) / x.two_norm2(), x);
@@ -183,7 +197,7 @@ void testSampleFunction() {
   SampleFunctional J(A, b);
 
   std::cout << J.directionalDerivative(b, b) << std::endl;
-  assert(J.directionalDerivative(b, b) == 10 + 2 * sqrt(5));
+  assert(J.pseudoDirectionalDerivative(b, b) == 10 + 2 * sqrt(5));
 
   SampleFunctional::SmallVector start = b;
   start *= 17;
@@ -209,7 +223,7 @@ void testTrivialFunction() {
   SampleFunctional J(A, b);
 
   std::cout << J.directionalDerivative(b, b) << std::endl;
-  assert(J.directionalDerivative(b, b) == 10);
+  assert(J.pseudoDirectionalDerivative(b, b) == 10);
 
   SampleFunctional::SmallVector start = b;
   start *= 17;

src/bisection-example.cc

@@ -32,12 +32,11 @@ template <int dimension> class SampleFunctional {
     return y * v + H(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
   }
 
-  double directionalDerivative(const SmallVector x,
-                               const SmallVector dir) const {
+  // |dir|-times the directional derivative wrt dir/|dir|.
+  double pseudoDirectionalDerivative(const SmallVector x,
+                                     const SmallVector dir) const {
     if (x == SmallVector(0.0))
-      // Well, not in this way -- but can we compute them?
-      DUNE_THROW(Dune::Exception,
-                 "Directional derivatives cannot be computed at zero.");
+      return HPrimePlus(0) * dir.two_norm();
 
     if (x * dir > 0)
       return PlusGrad(x) * dir;
@@ -45,6 +44,19 @@ template <int dimension> class SampleFunctional {
       return MinusGrad(x) * dir;
   }
 
+  double directionalDerivative(const SmallVector x,
+                               const SmallVector dir) const {
+    double norm = dir.two_norm();
+
+    if (norm == 0)
+      DUNE_THROW(Dune::Exception, "Directional derivatives cannot be computed "
+                                  "w.r.t. the zero-direction.");
+
+    SmallVector tmp = dir;
+    tmp /= norm;
+    return pseudoDirectionalDerivative(x, tmp);
+  }
+
   SmallVector minimise(const SmallVector x, unsigned int iterations) const {
     SmallVector descDir = ModifiedGradient(x);
 
@@ -153,6 +165,7 @@ template <int dimension> class SampleFunctional {
     return ret;
   }
 
+  // No normalising is done!
   SmallVector project(const SmallVector z, const SmallVector x) const {
     SmallVector y = z;
     y.axpy(-(z * x) / x.two_norm2(), x);
@@ -177,7 +190,7 @@ int main() {
     SampleFunctional J(A, b);
 
     std::cout << J.directionalDerivative(b, b) << std::endl;
-    assert(J.directionalDerivative(b, b) == 10 + 2 * sqrt(5));
+    assert(J.pseudoDirectionalDerivative(b, b) == 10 + 2 * sqrt(5));
 
     SampleFunctional::SmallVector start = b;
     start *= 17;