Added slower and more accurate central differencing method

src/estimate_gradient.jl

@@ -1,38 +1,66 @@
-# Simple forward-differencing based on Nocedal and Wright and
-# suggestions from Nathaniel Daw.
-#
-# Will need to implement central-differencing as well. Wrap both.
+# Simple forward-differencing and central-differencing
+# based on Nocedal and Wright and
+# suggestions from Nathaniel Daw. 
 
 # Arguments:
 # f: A function.
 # x: A vector in the domain of f.
 
+# Forward-differencing
 function estimate_gradient(f::Function)
-  function g(x::Vector)
+  # Construct and return a closure.
+  function g(x::Vector)    
+    # How far do we move in each direction?
+    # See Nodedal and Wright for motivation.
+    epsilon = sqrt(10e-16)
+
+    # What is the dimension of x?
+    n = length(x)
+
+    # Compute forward differences.
+    diff = zeros(n)
+
     # Establish a baseline value of f(x).
     f_x = f(x)
 
+    # Iterate over each dimension of the gradient separately.
+    for j = 1:n
+      dx = zeros(n)
+      dx[j] = epsilon
+      diff[j] = f(x + dx) - f_x
+      diff[j] /= epsilon
+    end
+
+    # Return the estimated gradient.
+    diff
+  end
+  g
+end
+
+# Central-differencing.
+function estimate_gradient2(f::Function)
+  # Construct and return a closure.
+  function g(x::Vector)    
     # How far do we move in each direction?
     # See Nodedal and Wright for motivation.
-    alpha = sqrt(10e-16)
+    epsilon = sqrt(10e-16)
 
     # What is the dimension of x?
     n = length(x)
 
     # Compute forward differences.
     diff = zeros(n)
 
-    # Iterate over each dimension separately.
+    # Iterate over each dimension of the gradient separately.
     for j = 1:n
       dx = zeros(n)
-      dx[j] = alpha
-      diff[j] = f(x + dx)
+      dx[j] = epsilon
+      diff[j] = f(x + dx) - f(x - dx)
+      diff[j] /= 2 * epsilon
     end
 
-    # Estimate gradient by dividing out by forward movement.
-    (diff - f_x) ./ alpha
+    # Return the estimated gradient.
+    diff
   end
-
-  # Return a closure.
   g
 end

tests/estimate_gradient.jl

@@ -5,6 +5,10 @@ function f(x::Vector)
   (100.0 - x[1])^2 + (50.0 - x[2])^2
 end
 
+#
+# Forward differencing
+#
+
 # Check that gradient is approximately accurate.
 g = estimate_gradient(f)
 @assert norm(g([100.0, 50.0]) - [0.0, 0.0]) < 0.01
@@ -13,6 +17,20 @@ g = estimate_gradient(f)
 results = optimize(f, g, [0.0, 0.0])
 @assert norm(results.minimum - [100.0, 50.0]) < 0.01
 
-# Compare with running optimize() without using finite differencing.
+#
+# Central differencing
+#
+
+# Check that gradient is approximately accurate.
+g = estimate_gradient2(f)
+@assert norm(g([100.0, 50.0]) - [0.0, 0.0]) < 0.01
+
+# Run optimize() using approximate gradient.
+results = optimize(f, g, [0.0, 0.0])
+@assert norm(results.minimum - [100.0, 50.0]) < 0.01
+
+# Comparison
+
+# Compare both results with running optimize() without using any form of finite differencing.
 results = optimize(f, [0.0, 0.0])
 @assert norm(results.minimum - [100.0, 50.0]) < 0.01