From 4dd7b29eba815116ad8d2ae074f3bdc6f3ccbba9 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 11 Nov 2011 15:53:15 +0100
Subject: [PATCH] Bisect along a sphere
---
dune/tectonic/curvedfunction.hh | 70 +++++++++++++++
dune/tectonic/samplefunctional.hh | 137 ++++++++++++++++++------------
src/test-gradient-method.cc | 18 ++--
3 files changed, 161 insertions(+), 64 deletions(-)
create mode 100644 dune/tectonic/curvedfunction.hh
diff --git a/dune/tectonic/curvedfunction.hh b/dune/tectonic/curvedfunction.hh
new file mode 100644
index 00000000..9d5fb281
--- /dev/null
+++ b/dune/tectonic/curvedfunction.hh
@@ -0,0 +1,70 @@
+#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 φ
+ 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
diff --git a/dune/tectonic/samplefunctional.hh b/dune/tectonic/samplefunctional.hh
index a19d40a5..50d8f785 100644
--- a/dune/tectonic/samplefunctional.hh
+++ b/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;
}
}
}
diff --git a/src/test-gradient-method.cc b/src/test-gradient-method.cc
index fa3b5827..ffb7e9ca 100644
--- a/src/test-gradient-method.cc
+++ b/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);
}
--
GitLab