Added a modified Newton's algorithm to handle both convex and concave…

src/Optim.jl

@@ -50,6 +50,7 @@ module Optim
     include("newton.jl")
     include("bfgs.jl")
     include("l_bfgs.jl")
+    include("modified_newton.jl")
 
     # Constrained optimization
     include("fminbox.jl")

src/modified_newton.jl

@@ -0,0 +1,145 @@
+#Nocedal and Wright page 145
+function modified_newton{T}(d::TwiceDifferentiableFunction,
+                   initial_x::Vector{T};
+                   xtol::Real = 1e-32,
+                   ftol::Real = 1e-8,
+                   grtol::Real = 1e-8,
+                   iterations::Integer = 1_000,
+                   store_trace::Bool = false,
+                   show_trace::Bool = false,
+                   extended_trace::Bool = false,
+                   linesearch!::Function = hz_linesearch!)
+
+    # Maintain current state in x and previous state in x_previous
+    x, x_previous = copy(initial_x), copy(initial_x)
+
+    # Count the total number of iterations
+    iteration = 0
+
+    # Track calls to function and gradient
+    f_calls, g_calls = 0, 0
+
+    # Count number of parameters
+    n = length(x)
+
+    # Maintain current gradient in gr
+    gr = Array(T, n)
+
+    # The current search direction
+    # TODO: Try to avoid re-allocating s
+    s = Array(T, n)
+
+    # Buffers for use in line search
+    x_ls, gr_ls = Array(T, n), Array(T, n)
+
+    # Store f(x) in f_x
+    f_x_previous, f_x = NaN, d.fg!(x, gr)
+    f_calls, g_calls = f_calls + 1, g_calls + 1
+
+    # Store h(x) in H
+    H = Array(T, n, n)
+    d.h!(x, H)
+
+    # Keep track of step-sizes
+    alpha = alphainit(one(T), x, gr, f_x)
+
+    # TODO: How should this flag be set?
+    mayterminate = false
+
+    # Maintain a cache for line search results
+    lsr = LineSearchResults(T)
+
+    # Trace the history of states visited
+    tr = OptimizationTrace()
+    tracing = store_trace || show_trace || extended_trace
+    @newtontrace
+
+    # Assess multiple types of convergence
+    x_converged, f_converged, gr_converged = false, false, false
+
+    # Starting point for the Frobenius norm
+    succ = 0
+    β = 0.
+    diag_correct = eye(similar(H))
+    # Iterate until convergence
+    converged = false
+    while !converged && iteration < iterations
+        # Increment the number of steps we've had to perform
+        iteration += 1
+        if minimum(eigvals(H))<0
+        β = vecnorm(H)
+        τ = β/2
+          for i = 1:1:5000
+              succ = 1
+              try
+                  cholfact(H+τ*diag_correct)
+              catch
+                  τ = max(2τ, β/2)
+                  succ = 0
+              end
+              if succ == 1
+                  H=H+τ*diag_correct
+                  break
+              end
+          end
+        end
+        # Search direction is always the negative gradient divided by H
+        # TODO: Do this calculation in place
+        @inbounds s[:] = -(H \ gr)
+
+        # Refresh the line search cache
+        dphi0 = dot(gr, s)
+        clear!(lsr)
+        push!(lsr, zero(T), f_x, dphi0)
+
+        # Determine the distance of movement along the search line
+        alpha, f_update, g_update =
+          linesearch!(d, x, s, x_ls, gr_ls, lsr, alpha, mayterminate)
+        f_calls, g_calls = f_calls + f_update, g_calls + g_update
+
+        # Maintain a record of previous position
+        copy!(x_previous, x)
+
+        # Update current position
+        for i in 1:n
+            @inbounds x[i] = x[i] + alpha * s[i]
+        end
+
+        # Update the function value and gradient
+        f_x_previous, f_x = f_x, d.fg!(x, gr)
+        f_calls, g_calls = f_calls + 1, g_calls + 1
+
+        # Update the Hessian
+        d.h!(x, H)
+
+        x_converged,
+        f_converged,
+        gr_converged,
+        converged = assess_convergence(x,
+                                       x_previous,
+                                       f_x,
+                                       f_x_previous,
+                                       gr,
+                                       xtol,
+                                       ftol,
+                                       grtol)
+
+        @newtontrace
+    end
+
+    return MultivariateOptimizationResults("Modified Newton's Method",
+                                           initial_x,
+                                           x,
+                                           @compat Float64(f_x),
+                                           iteration,
+                                           iteration == iterations,
+                                           x_converged,
+                                           xtol,
+                                           f_converged,
+                                           ftol,
+                                           gr_converged,
+                                           grtol,
+                                           tr,
+                                           f_calls,
+                                           g_calls)
+end