From f0f9d0c54fc4bd390566e86012f0be0561a38578 Mon Sep 17 00:00:00 2001 From: Patrick Kofod Mogensen Date: Sat, 17 Oct 2020 23:42:01 +0200 Subject: [PATCH] Fix Krylov and remove stuff from Project.toml. --- Project.toml | 6 +- src/NLSolvers.jl | 4 - src/nlsolve/linesearch/krylov.jl | 32 +++----- .../randomsearch/simulatedannealing.jl | 4 +- test/nlsolve/interface.jl | 74 +++++++++++++++++++ test/nlsolve/krylov.jl | 30 ++++++-- test/nlsolve/runtests.jl | 1 + 7 files changed, 114 insertions(+), 37 deletions(-) diff --git a/Project.toml b/Project.toml index be46d2b..5523d75 100644 --- a/Project.toml +++ b/Project.toml @@ -4,11 +4,8 @@ authors = ["Patrick Kofod Mogensen "] 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"] diff --git a/src/NLSolvers.jl b/src/NLSolvers.jl index c4158a6..430599f 100644 --- a/src/NLSolvers.jl +++ b/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 diff --git a/src/nlsolve/linesearch/krylov.jl b/src/nlsolve/linesearch/krylov.jl index 87d90c8..06bb7d1 100644 --- a/src/nlsolve/linesearch/krylov.jl +++ b/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 diff --git a/src/optimize/randomsearch/simulatedannealing.jl b/src/optimize/randomsearch/simulatedannealing.jl index 26bef17..c9b319a 100644 --- a/src/optimize/randomsearch/simulatedannealing.jl +++ b/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 diff --git a/test/nlsolve/interface.jl b/test/nlsolve/interface.jl index c76347b..bc99194 100644 --- a/test/nlsolve/interface.jl +++ b/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 \ No newline at end of file diff --git a/test/nlsolve/krylov.jl b/test/nlsolve/krylov.jl index 981a4ba..7862cf5 100644 --- a/test/nlsolve/krylov.jl +++ b/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,12 +40,12 @@ 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) @@ -53,7 +53,7 @@ function OnceDiffedJv(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 diff --git a/test/nlsolve/runtests.jl b/test/nlsolve/runtests.jl index 3b71870..8096f73 100644 --- a/test/nlsolve/runtests.jl +++ b/test/nlsolve/runtests.jl @@ -1,2 +1,3 @@ include("interface.jl") +include("krylov.jl") include("MGH.jl") \ No newline at end of file