Use PositiveFactorizations in Newton's method.

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.4
 Calculus
 DualNumbers 0.2
+PositiveFactorizations

src/Optim.jl

@@ -2,6 +2,7 @@ VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
 
 module Optim
     using Calculus
+    using PositiveFactorizations
     using Compat
 
     import Base.length,

src/newton.jl

@@ -95,9 +95,13 @@ function optimize{T}(d::TwiceDifferentiableFunction,
         # Increment the number of steps we've had to perform
         iteration += 1
 
-        # Search direction is always the negative gradient divided by H
-        # TODO: Do this calculation in place
-        @inbounds s[:] = -(H \ gr)
+        # Search direction is always the negative gradient divided by
+        # a matrix encoding the absolute values of the curvatures
+        # represented by H. It deviates from the usual "add a scaled
+        # identity matrix" version of the modified Newton method. More
+        # information can be found in the discussion at issue #153.
+        F, Hd = ldltfact!(Positive, H)
+        s[:] = -(F\gr)
 
         # Refresh the line search cache
         dphi0 = vecdot(gr, s)

test/newton.jl

@@ -55,3 +55,10 @@ for (name, prob) in Optim.UnconstrainedProblems.examples
 		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
 	end
 end
+
+let
+    prob=Optim.UnconstrainedProblems.examples["Himmelblau"]
+    ddf = TwiceDifferentiableFunction(prob.f, prob.g!, prob.h!)
+    res = optimize(ddf, [0., 0.], Newton())
+    @assert norm(res.minimum - prob.solutions) < 1e-10
+end