Fix Krylov and remove stuff from Project.toml.

Project.toml

@@ -4,11 +4,8 @@ authors = ["Patrick Kofod Mogensen <[email protected]>"]
 version = "0.1.0"
 
 [deps]
-IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
@@ -19,6 +16,7 @@ ArbNumerics = "7e558dbc-694d-5a72-987c-6f4ebed21442"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -28,4 +26,4 @@ StaticOptim = "814864de-89d5-11e8-1c72-cfea42f978d7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ArbNumerics", "DoubleFloats", "ForwardDiff", "GeometryTypes", "Printf", "StaticArrays", "SparseDiffTools", "SparseArrays", "Random"]
+test = ["Test", "ArbNumerics", "DoubleFloats", "ForwardDiff", "GeometryTypes", "IterativeSolvers", "Printf", "StaticArrays", "SparseDiffTools", "SparseArrays", "Random"]

src/NLSolvers.jl

@@ -34,9 +34,6 @@ using LinearAlgebra:  dot, I, norm, # used everywhere in updates, convergence, e
                       axpy! # for Anderson
 import LinearAlgebra: mul!, dot # need to extend for preconditioners
 
-# For better random number generators and rand!
-using RandomNumbers
-
 using Printf
 
 function solve end
@@ -56,7 +53,6 @@ end
 export objective_return
 isallfinite(x) = mapreduce(isfinite, *, x)
 
-using StaticArrays
 abstract type MutateStyle end
 
 abstract type AbstractProblem end

src/nlsolve/linesearch/krylov.jl

@@ -5,7 +5,6 @@
 # https://pdfs.semanticscholar.org/b321/3084f663260076dcb92f2fa6031b362dc5bc.pdf
 # https://www.sciencedirect.com/science/article/abs/pii/0098135483800027
 # 
-using IterativeSolvers
 abstract type ForcingSequence end
 
 struct FixedForceTerm{T} <: ForcingSequence
@@ -83,19 +82,21 @@ end
 Constructs a method type for the Inexact Newton's method with Linesearch.
 
 """
-struct InexactNewton{ForcingType<:ForcingSequence, Tη}
+struct InexactNewton{LinSolve, ForcingType<:ForcingSequence, Tη}
+    linsolve::LinSolve
     force_seq::ForcingType
     η₀::Tη
     maxiter::Int
 end
 #InexactNewton(; force_seq=FixedForceTinexacerm(1e-4), eta0 = 1e-4, maxiter=300)=InexactNewton(force_seq, eta0, maxiter)
-InexactNewton(; force_seq=DemboSteihaug(), eta0 = 1e-4, maxiter=300)=InexactNewton(force_seq, eta0, maxiter)
+
+InexactNewton(; linsolve, force_seq=DemboSteihaug(), eta0 = 1e-4, maxiter=300)=InexactNewton(linsolve, force_seq, eta0, maxiter)
 # map from method to forcing sequence
 η(fft::InexactNewton, info) = η(fft.force_seq, info)
 
 
 function solve(problem::NEqProblem, x, method::InexactNewton, options::NEqOptions)
-    if !(mstyle(prob) === InPlace())
+    if !(mstyle(problem) === InPlace())
         throw(ErrorException("solve() not defined for OutOfPlace() with InexactNewton"))
     end
     t0 = time()
@@ -105,7 +106,7 @@ function solve(problem::NEqProblem, x, method::InexactNewton, options::NEqOption
 
     Tx = eltype(x)
     xp, Fx = copy(x), copy(x)
-    Fx = F(x, Fx)
+    Fx = F(Fx, x)
     JvOp = JvGen(x)
 
     ρ2F0 = norm(Fx, 2)
@@ -126,21 +127,13 @@ function solve(problem::NEqProblem, x, method::InexactNewton, options::NEqOption
         if iter == 1 && !isa(method.force_seq, FixedForceTerm) 
             ηₖ = method.η₀
         else 
-           ηₖ = η(method, force_info)
+            ηₖ = η(method, force_info)
         end
-
         JvOp = JvGen(x)
 
         xp .= 0
-        krylov_iter = IterativeSolvers.gmres_iterable!(xp, JvOp, Fx; maxiter=50)
-        res = copy(Fx)
-        rhs = ηₖ*norm(Fx, 2)
-        for item in krylov_iter
-            res = krylov_iter.residual.current
-            if res <= rhs
-                break
-            end
-        end
+        xp, res = method.linsolve(xp, JvOp, Fx, ηₖ)
+
         ρs = norm(xp)
         t = 1e-4
         # use line searc functions here
@@ -152,10 +145,9 @@ function solve(problem::NEqProblem, x, method::InexactNewton, options::NEqOption
         it = 0
         while !btk_conv
             it += 1
-
-            z = retract(problem, z, x, xp)
-
-            Fx = problem.R.F(z, Fx)
+            z = retract(problem, z, x, xp, -1)
+
+            Fx = problem.R.F(Fx, z)
             btk_conv = norm(Fx, 2) ≤ (1-t*(1-ηₖ))*ρFx || it > 20
         end
         if norm(Fx, 2) < stoptol

src/optimize/randomsearch/simulatedannealing.jl

@@ -32,7 +32,7 @@ log_temperature(t) = 1 / log(t)^2
 function default_neighbor(x_best)
   T = eltype(x_best)
   n = length(x_best)
-  return x_best .+ T.(RandomNumbers.randn(n))
+  return x_best .+ T.(randn(n))
 end
 
 function solve(prob::OptimizationProblem, x0, method::SimulatedAnnealing, options::OptimizationOptions)
@@ -71,7 +71,7 @@ function solve(prob::OptimizationProblem, x0, method::SimulatedAnnealing, option
     else
       # If proposal is inferior, we move to it with probability p
       p = exp(-(f_candidate - f_now) / temperature)
-      if RandomNumbers.rand() <= p
+      if rand() <= p
         x_now = copy(x_candidate)
         f_now = f_candidate
       end

test/nlsolve/interface.jl

@@ -291,4 +291,78 @@ end
 # #       colorbar=true)
 
 
+end
+
+@testset "Krylov" begin
+    # Start with a stupid example
+    n = 10
+    A = sprand(10, 10, 0.1)
+    A = A*A' + I
+    F(Fx, x) = mul!(Fx, A, x)
+
+    x0, = rand(10)
+    xp = copy(x0)
+    Fx = copy(xp)
+
+    function theta(x)
+        if x[1] > 0
+            return atan(x[2] / x[1]) / (2.0 * pi)
+        else
+            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
+        end
+    end
+
+    function F_powell!(Fx, x)
+            if ( x[1]^2 + x[2]^2 == 0 )
+                dtdx1 = 0;
+                dtdx2 = 0;
+            else
+                dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
+                dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
+            end
+            Fx[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
+                200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
+            Fx[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
+                200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
+            Fx[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
+        Fx
+    end
+
+    function F_jacobian_powell!(Fx, Jx, x)
+        ForwardDiff.jacobian!(Jx, F_powell!, Fx, x)
+        Fx, Jx
+    end
+    x0 = [-1.0, 0.0, 0.0]
+
+    Fc, Jc = zeros(3), zeros(3,3)
+    F_jacobian_powell!(Fc, Jc, x0)
+
+    jv = JacVec(F_powell!, rand(3); autodiff=false)
+    function jvop(x)
+        jv.u .= x
+        jv
+    end
+    prob_obj = NLSolvers.VectorObjective(F_powell!, nothing, F_jacobian_powell!, jvop)
+
+
+    prob = NEqProblem(prob_obj)
+    res = solve(prob, copy(x0), LineSearch(Newton(), Backtracking()))
+    @test norm(res.info.best_residual) < 1e-15
+
+    using IterativeSolvers
+    function inexact_linsolve(x0, JvOp, Fx, ηₖ)
+        krylov_iter = IterativeSolvers.gmres_iterable!(x0, JvOp, Fx; maxiter=50)
+        res = copy(Fx)
+        rhs = ηₖ*norm(Fx, 2)
+        for item in krylov_iter
+            res = krylov_iter.residual.current
+            if res <= rhs
+                break
+            end
+        end
+        return x0, res
+    end
+    res = solve(prob, copy(x0), InexactNewton(inexact_linsolve, FixedForceTerm(0.001), 1e-4, 300), NEqOptions(maxiter=1000))
+    @test norm(res.info.best_residual) < 1e-10
+
 end

test/nlsolve/krylov.jl

@@ -24,7 +24,7 @@ function theta(x)
    end
 end
 
-function F_powell!(x, Fx)
+function F_powell!(Fx, x)
         if ( x[1]^2 + x[2]^2 == 0 )
             dtdx1 = 0;
             dtdx2 = 0;
@@ -40,20 +40,20 @@ function F_powell!(x, Fx)
     Fx
 end
 
-function F_jacobian_powell!(x, Fx, Jx)
-    ForwardDiff.jacobian!(Jx, (y,x)->F_powell!(x,y), Fx, x)
+function F_jacobian_powell!(Fx, Jx, x)
+    ForwardDiff.jacobian!(Jx, F_powell!, Fx, x)
     Fx, Jx
 end
 Fc, Jc = zeros(3), zeros(3,3)
-F_jacobian_powell!(x0, Fc, Jc)
+F_jacobian_powell!(Fc, Jc, x0)
 
 import NLSolvers: OnceDiffedJv
 function OnceDiffedJv(F; seed, autodiff=false)
     JacOp = JacVec(F, seed; autodiff=false)
     OnceDiffedJv(F, JacOp)
 end
 
-jv = JacVec((y,x)->F_powell!(x,y), rand(3); autodiff=false)
+jv = JacVec(F_powell!, rand(3); autodiff=false)
 function jvop(x)
     jv.u .= x
     jv
@@ -63,8 +63,24 @@ prob_obj = NLSolvers.VectorObjective(F_powell!, nothing, F_jacobian_powell!, jvo
 
 prob = NEqProblem(prob_obj)
 x0 = [-1.0, 0.0, 0.0]
-solve(prob, x0, LineSearch(Newton(), Backtracking()))
-solve(prob, x0, InexactNewton(FixedForceTerm(0.4), 1e-12, 300), NEqOptions(maxiter=1))
+res = solve(prob, copy(x0), LineSearch(Newton(), Backtracking()))
+@test norm(res.info.best_residual) < 1e-15
+
+using IterativeSolvers
+function inexact_linsolve(x0, JvOp, Fx, ηₖ)
+    krylov_iter = IterativeSolvers.gmres_iterable!(x0, JvOp, Fx; maxiter=50)
+    res = copy(Fx)
+    rhs = ηₖ*norm(Fx, 2)
+    for item in krylov_iter
+        res = krylov_iter.residual.current
+        if res <= rhs
+            break
+        end
+    end
+    return x0, res
+end
+res = solve(prob, copy(x0), InexactNewton(inexact_linsolve, FixedForceTerm(0.001), 1e-4, 300), NEqOptions(maxiter=1000))
+@test norm(res.info.best_residual) < 1e-10
 
 function f_2by2!(F, x)
     F[1] = (x[1]+3)*(x[2]^3-7)+18

test/nlsolve/runtests.jl

@@ -1,2 +1,3 @@
 include("interface.jl")
+include("krylov.jl")
 include("MGH.jl")