From 7de078cac316ae7a76a1d4997c746deda4d25edd Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 31 Jul 2012 09:24:01 +0200
Subject: [PATCH] Modularise minimise()
---
dune/tectonic/samplefunctional.hh | 147 +++++++++++++++++-------------
1 file changed, 83 insertions(+), 64 deletions(-)
diff --git a/dune/tectonic/samplefunctional.hh b/dune/tectonic/samplefunctional.hh
index 750bba10..e88dd863 100644
--- a/dune/tectonic/samplefunctional.hh
+++ b/dune/tectonic/samplefunctional.hh
@@ -133,6 +133,87 @@ template <int dim> class SampleFunctional {
}
};
+template <class Functional>
+void descentMinimisation(Functional const &J,
+ typename Functional::SmallVector &x,
+ typename Functional::SmallVector const &descDir,
+ Bisection const &bisection) {
+ typedef typename Functional::SmallVector SmallVector;
+ typedef typename Functional::NonlinearityType LocalNonlinearityType;
+
+ // {{{ 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 = J.b; // b
+ J.A.mmv(x, tmp); // b-Au
+ double const JRestb = tmp * descDir; // <b-Au,v>
+
+ J.A.mv(descDir, tmp); // Av
+ double const JRestA = tmp * descDir; // <Av,v>
+
+ MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
+ JRestA, JRestb, *J.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.
+ */
+ }
+
+ int count;
+ double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
+ dverb << "Number of iterations in the bisection method: " << count
+ << std::endl;
+ ;
+
+ x.axpy(stepsize, descDir);
+}
+
+template <class Functional>
+void tangentialMinimisation(Functional const &J,
+ typename Functional::SmallVector &x,
+ typename Functional::SmallVector const &descDir,
+ Bisection const &bisection) {
+ typedef typename Functional::NonlinearityType LocalNonlinearityType;
+ typedef typename Functional::SmallVector SmallVector;
+
+ CircularConvexFunction<LocalNonlinearityType> const JRest(J.A, J.b, *J.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;
+ ;
+
+ // Since x is used in the computation of the rhs, do not write to it directly
+ SmallVector tmp;
+ JRest.cartesian(stepsize, tmp);
+ x = tmp;
+}
+
template <class Functional>
void minimise(Functional const &J, typename Functional::SmallVector &x,
size_t steps, Bisection const &bisection) {
@@ -145,75 +226,13 @@ void minimise(Functional const &J, typename Functional::SmallVector &x,
if (descDir == SmallVector(0.0))
return;
- typedef typename Functional::NonlinearityType LocalNonlinearityType;
if (linesearchp) {
- // {{{ 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 = J.b; // b
- J.A.mmv(x, tmp); // b-Au
- double const JRestb = tmp * descDir; // <b-Au,v>
-
- J.A.mv(descDir, tmp); // Av
- double const JRestA = tmp * descDir; // <Av,v>
-
- MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
- JRestA, JRestb, *J.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.
- */
- }
-
- int count;
- double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
- dverb << "Number of iterations in the bisection method: " << count
- << std::endl;
- ;
-
- x.axpy(stepsize, descDir);
+ descentMinimisation(J, x, descDir, bisection);
} else {
Bisection slowBisection(bisection);
slowBisection.setFastQuadratic(false);
- CircularConvexFunction<LocalNonlinearityType> const JRest(
- J.A, J.b, *J.phi, x, descDir);
-
- int count;
- double const stepsize = slowBisection.minimize(JRest, 0.0, 1.0, count);
- dverb << "Number of iterations in the bisection method: " << count
- << std::endl;
- ;
-
- // Since x is used in the computation of the rhs, do not write to it
- // directly
- SmallVector tmp;
- JRest.cartesian(stepsize, tmp);
- x = tmp;
+ tangentialMinimisation(J, x, descDir, slowBisection);
}
}
}
--
GitLab