Bisect along a sphere

dune/tectonic/curvedfunction.hh

+#ifndef CURVED_FUNCTION_HH
+#define CURVED_FUNCTION_HH
+
+namespace Dune {
+template <class NonlinearityType> class CurvedFunction {
+  typedef typename NonlinearityType::VectorType VectorType;
+  typedef typename NonlinearityType::MatrixType MatrixType;
+
+public:
+  CurvedFunction(MatrixType const &A, VectorType const &b,
+                 NonlinearityType const &phi, VectorType const &x,
+                 VectorType const &dir)
+      : A(A), b(b), phi(phi), x(x), dir(dir) {}
+
+  double quadraticPart() const { return 0; }
+
+  double linearPart() const { return 0; }
+
+  void subDiff(double m, Interval<double> &D) const {
+    VectorType x;
+    evaluate(m, x);
+    VectorType tangent;
+    tangentialDirection(m, tangent);
+
+    phi.directionalSubDiff(x, tangent, D);
+
+    VectorType tmp;
+    A.mv(x, tmp);                      //  Ax
+    tmp -= b;                          //  Ax - b
+    double const dotp = tmp * tangent; // <Ax - b,t>
+    D[0] += dotp;
+    D[1] += dotp;
+  }
+
+  void domain(Interval<double> &domain) const {
+    // TODO
+    domain[0] = 0;
+    domain[1] = 1;
+  }
+
+  void evaluate(double m, VectorType &y) const {
+    y = 0;
+    y.axpy(1 - m, x);
+    y.axpy(m, dir);
+    y /= y.two_norm();
+  }
+
+private:
+  MatrixType const &A;
+  VectorType const &b;
+  NonlinearityType const &phi;
+  VectorType const &x;
+  VectorType const &dir;
+
+  /*
+    Tangential vector within the plane, with positive direction.
+
+    < (1-m)x + m*y,-m*|y|^2*x+(1-m)*|x|^2*y>
+    = -m(1-m)|y|^2<x,x> + m(1-m)|x|^2<y,y> (because <x,y> = 0)
+    = 0
+  */
+  void tangentialDirection(double m, VectorType &y) const {
+    y = 0;
+    y.axpy(-m * dir.two_norm2(), x);
+    y.axpy((1 - m) * x.two_norm2(), dir);
+  }
+};
+}
+
+#endif

dune/tectonic/samplefunctional.hh

@@ -13,6 +13,7 @@
 #include <dune/tnnmg/problem-classes/directionalconvexfunction.hh>
 
 #include "localnonlinearity.hh"
+#include "curvedfunction.hh"
 
 namespace Dune {
 template <int dim> class SampleFunctional {
@@ -34,7 +35,7 @@ template <int dim> class SampleFunctional {
     return y * v + phi(v); // <1/2 Av - b,v> + H(|v|)
   }
 
-  void descentDirection(SmallVector const x, SmallVector &ret) const {
+  bool descentDirection(SmallVector const x, SmallVector &ret) const {
     // Check the squared norm rather than each component because
     // complementaryProjection() divides by it
     if (x.two_norm2() == 0.0) {
@@ -54,7 +55,7 @@ template <int dim> class SampleFunctional {
       } else {
         ret = 0;
       }
-      return;
+      return true;
     }
 
     SmallVector pg;
@@ -68,11 +69,15 @@ template <int dim> class SampleFunctional {
       dverb << "## Directional derivative (as per scalar product w/ "
                "semigradient): " << -(ret * mg)
             << " (coordinates of the restriction)" << std::endl;
+      ret *= -1;
+      return true;
     } else if (pgx <= 0 && mgx <= 0) {
       ret = mg;
       dverb << "## Directional derivative (as per scalar product w/ "
                "semigradient): " << -(ret * pg)
             << " (coordinates of the restriction)" << std::endl;
+      ret *= -1;
+      return true;
     } else {
       // Includes the case that pg points in direction x and mg
       // points in direction -x. The projection will then be zero.
@@ -82,8 +87,9 @@ template <int dim> class SampleFunctional {
       dverb << "## Directional derivative (as per scalar product w/ "
                "semigradient): " << -(ret * ret)
             << " (coordinates of the restriction)" << std::endl;
+      ret *= -1;
+      return false;
     }
-    ret *= -1;
   }
 
   SmallMatrix const &A;
@@ -136,7 +142,7 @@ void minimise(Functional const J, typename Functional::SmallVector &x,
                                  // become smaller than this number
                             1.0, // acceptFactor: ?
                             1e-12, // requiredResidual: ?
-                            true,  // fastQuadratic
+                            false, // fastQuadratic
                             0))    // safety: acceptance factor for inexact
                                    // minimization
 {
@@ -144,68 +150,89 @@ void minimise(Functional const J, typename Functional::SmallVector &x,
 
   for (size_t step = 0; step < steps; ++step) {
     SmallVector descDir;
-    J.descentDirection(x, descDir);
+    double linesearchp = J.descentDirection(x, descDir);
 
     if (descDir == SmallVector(0.0))
       return;
 
-    // {{{ Construct a restriction of J to the line x + t * descDir
-
-    /* We have
-
-       1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
-
-       since A is symmetric.
-    */
-    SmallVector tmp;
+    if (linesearchp) {
+      // {{{ Construct a restriction of J to the line x + t * descDir
 
-    J.A.mv(descDir, tmp);                //  Av
-    double const JRestA = tmp * descDir; // <Av,v>
+      /* We have
 
-    tmp = J.b;                           //  b
-    J.A.mmv(x, tmp);                     //  b-Au
-    double const JRestb = tmp * descDir; // <b-Au,v>
+         1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2
+         Au-b,u>
 
-    typedef typename Functional::NonlinearityType LocalNonlinearityType;
-    LocalNonlinearityType phi = J.phi;
-    typedef DirectionalConvexFunction<LocalNonlinearityType>
-    MyDirectionalConvexFunctionType;
-    // FIXME: We cannot pass J.phi directly because the constructor
-    // does not allow for constant arguments
-    MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir);
-    // }}}
-
-    { // Debug
-      Interval<double> D;
-      JRest.subDiff(0, D);
+         since A is symmetric.
+      */
+      SmallVector tmp;
+
+      J.A.mv(descDir, tmp);                //  Av
+      double const JRestA = tmp * descDir; // <Av,v>
+
+      tmp = J.b;                           //  b
+      J.A.mmv(x, tmp);                     //  b-Au
+      double const JRestb = tmp * descDir; // <b-Au,v>
+
+      typedef typename Functional::NonlinearityType LocalNonlinearityType;
+      LocalNonlinearityType phi = J.phi;
+      typedef DirectionalConvexFunction<LocalNonlinearityType>
+      MyDirectionalConvexFunctionType;
+      // FIXME: We cannot pass J.phi directly because the constructor
+      // does not allow for constant arguments
+      MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir);
+      // }}}
+
+      { // Debug
+        Interval<double> D;
+        JRest.subDiff(0, D);
+
+        dverb << "## Directional derivative (as per subdifferential of "
+                 "restriction): " << D[1] << " (coordinates of the restriction)"
+              << std::endl;
+        /*
+          It is possible that this differs quite a lot from the
+          directional derivative computed in the descentDirection()
+          method:
+
+          If phi is nonsmooth at x, so that the directional
+          derivatives jump, and |x| is computed to be too small or too
+          large globally or locally, the locally computed
+          subdifferential and the globally computed subdifferential
+          will no longer coincide!
+
+          The assertion D[1] <= 0 may thus fail.
+        */
+      }
 
-      dverb << "## Directional derivative (as per subdifferential of "
-               "restriction): " << D[1] << " (coordinates of the restriction)"
+      int count;
+      // FIXME: The value of x_old should not matter if the factor is 1.0,
+      // correct?
+      double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
+      dverb << "Number of iterations in the bisection method: " << count
             << std::endl;
-      /*
-        It is possible that this differs quite a lot from the
-        directional derivative computed in the descentDirection()
-        method:
-
-        If phi is nonsmooth at x, so that the directional
-        derivatives jump, and |x| is computed to be too small or too
-        large globally or locally, the locally computed
-        subdifferential and the globally computed subdifferential
-        will no longer coincide!
-
-        The assertion D[1] <= 0 may thus fail.
-      */
-    }
+      ;
 
-    int count;
-    // FIXME: The value of x_old should not matter if the factor is 1.0,
-    // correct?
-    double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
-    dverb << "Number of iterations in the bisection method: " << count
-          << std::endl;
-    ;
+      x.axpy(stepsize, descDir);
+    } else {
+      typedef typename Functional::NonlinearityType LocalNonlinearityType;
+      LocalNonlinearityType phi = J.phi;
+      typedef Dune::CurvedFunction<LocalNonlinearityType> MyCurvedFunctionType;
+      MyCurvedFunctionType JRest(J.A, J.b, phi, x, descDir);
+
+      int count;
+      double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
+      dverb << "Number of iterations in the bisection method: " << count
+            << std::endl;
+      ;
 
-    x.axpy(stepsize, descDir);
+      // Since x is used in the computation of the rhs, do not write to it
+      // directly
+      SmallVector tmp;
+      JRest.evaluate(stepsize, tmp);
+      x = tmp;
+    }
+    // std::cout << "Resulting vector: " << x << std::endl;
   }
 }
 }

src/test-gradient-method.cc

@@ -379,19 +379,19 @@ void testSampleFunctionSteep1() {
   start[0] = 0;
   start[1] = 1;
 
-  double const ret1 = functionTester(J, start, 7);
+  double const ret1 = functionTester(J, start, 3);
 
   // Something random
   start[0] = 279;
   start[1] = -96;
 
-  double const ret2 = functionTester(J, start, 7);
+  double const ret2 = functionTester(J, start, 3);
   assert(std::abs(ret1 - ret2) < 1e-5);
 
   start[0] = 0;
   start[1] = 0;
 
-  double const ret3 = functionTester(J, start, 2);
+  double const ret3 = functionTester(J, start, 1);
   assert(std::abs(ret1 - ret3) < 1e-5);
 }
 
@@ -416,19 +416,19 @@ void testSampleFunctionSteep2() {
   start[0] = 0;
   start[1] = 1;
 
-  double const ret1 = functionTester(J, start, 7);
+  double const ret1 = functionTester(J, start, 1);
 
   // Something random
   start[0] = 279;
   start[1] = -96;
 
-  double const ret2 = functionTester(J, start, 7);
+  double const ret2 = functionTester(J, start, 3);
   assert(std::abs(ret1 - ret2) < 1e-5);
 
   start[0] = 0;
   start[1] = 0;
 
-  double const ret3 = functionTester(J, start, 3);
+  double const ret3 = functionTester(J, start, 1);
   assert(std::abs(ret1 - ret3) < 1e-5);
 }
 
@@ -453,19 +453,19 @@ void testSteepFunction() {
   start[0] = 0;
   start[1] = 1;
 
-  double const ret1 = functionTester(J, start, 1000);
+  double const ret1 = functionTester(J, start, 1);
 
   // Something random
   start[0] = 279;
   start[1] = -96;
 
-  double const ret2 = functionTester(J, start, 1000);
+  double const ret2 = functionTester(J, start, 4);
   assert(std::abs(ret1 - ret2) < 1e-8);
 
   start[0] = 0;
   start[1] = 0;
 
-  double const ret3 = functionTester(J, start, 3);
+  double const ret3 = functionTester(J, start, 1);
   assert(std::abs(ret1 - ret3) < 1e-5);
 }