From c70ebe6544d3a03dd950007bbfbdfd145ff31ff0 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 7 Sep 2011 15:43:36 +0200
Subject: [PATCH] Directional derivatives only for unit vectors
---
src/bisection-example-flexible.cc | 28 +++++++++++++++++++++-------
src/bisection-example.cc | 25 +++++++++++++++++++------
2 files changed, 40 insertions(+), 13 deletions(-)
diff --git a/src/bisection-example-flexible.cc b/src/bisection-example-flexible.cc
index e8add65b..61398aee 100644
--- a/src/bisection-example-flexible.cc
+++ b/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;
diff --git a/src/bisection-example.cc b/src/bisection-example.cc
index 1ae3c95e..37c345ba 100644
--- a/src/bisection-example.cc
+++ b/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;
--
GitLab