From f5835dadd3532c9eb7c0355f915320f6d9fe4b32 Mon Sep 17 00:00:00 2001 From: Paul Schrimpf Date: Fri, 4 Feb 2022 09:45:16 -0800 Subject: [PATCH] Made IPNewton() compatible with alternate types (e.g. sparse) of hessians and jacobians --- .../solvers/constrained/ipnewton/interior.jl | 6 +- .../solvers/constrained/ipnewton/ipnewton.jl | 36 +++++--- .../solvers/constrained/ipnewton/interface.jl | 2 +- .../constrained/ipnewton/matrixtypes.jl | 92 +++++++++++++++++++ 4 files changed, 118 insertions(+), 18 deletions(-) create mode 100644 test/multivariate/solvers/constrained/ipnewton/matrixtypes.jl diff --git a/src/multivariate/solvers/constrained/ipnewton/interior.jl b/src/multivariate/solvers/constrained/ipnewton/interior.jl index 3406f3c2f..7c5877f97 100644 --- a/src/multivariate/solvers/constrained/ipnewton/interior.jl +++ b/src/multivariate/solvers/constrained/ipnewton/interior.jl @@ -149,13 +149,13 @@ Base.convert(::Type{BarrierLineSearch{T}}, bsl::BarrierLineSearch) where T = Parameters for interior-point line search methods that exploit the slope. """ -struct BarrierLineSearchGrad{T} +struct BarrierLineSearchGrad{T, TJ<:AbstractMatrix} c::Vector{T} # value of constraints-functions at trial point - J::Matrix{T} # constraints-Jacobian at trial point + J::TJ # constraints-Jacobian at trial point bstate::BarrierStateVars{T} # trial point for slack and λ variables bgrad::BarrierStateVars{T} # trial point's gradient end -Base.convert(::Type{BarrierLineSearchGrad{T}}, bsl::BarrierLineSearchGrad) where T = +Base.convert(::Type{BarrierLineSearchGrad{T,TJ}}, bsl::BarrierLineSearchGrad) where T where TJ = BarrierLineSearchGrad(convert(Vector{T}, bsl.c), convert(Matrix{T}, bsl.J), convert(BarrierStateVars{T}, bsl.bstate), diff --git a/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl b/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl index e56a8365b..3e26c6899 100644 --- a/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl +++ b/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl @@ -1,8 +1,9 @@ -struct IPNewton{F,Tμ<:Union{Symbol,Number}} <: IPOptimizer{F} +struct IPNewton{F,Tμ<:Union{Symbol,Number}, M} <: IPOptimizer{F} linesearch!::F μ0::Tμ # Initial value for the barrier penalty coefficient μ show_linesearch::Bool # TODO: μ0, and show_linesearch were originally in options + conJstorage::M end Base.summary(::IPNewton) = "Interior Point Newton" @@ -46,17 +47,18 @@ The algorithm was [originally written by Tim Holy](https://github.com/JuliaNLSol """ IPNewton(; linesearch::Function = backtrack_constrained_grad, μ0::Union{Symbol,Number} = :auto, - show_linesearch::Bool = false) = - IPNewton(linesearch, μ0, show_linesearch) + show_linesearch::Bool = false, + conJstorage::Union{Symbol,AbstractMatrix} = :auto) = + IPNewton(linesearch, μ0, show_linesearch, conJstorage) -mutable struct IPNewtonState{T,Tx} <: AbstractBarrierState +mutable struct IPNewtonState{T,Tx,TH, TJ, THtilde} <: AbstractBarrierState x::Tx f_x::T x_previous::Tx g::Tx f_x_previous::T - H::Matrix{T} # Hessian of the Lagrangian? - HP # TODO: remove HP? It's not used + H::TH # Hessian of the Lagrangian? + HP::Nothing # TODO: remove HP? It's not used Hd::Vector{Int8} # TODO: remove Hd? It's not used s::Tx # step for x # Barrier penalty fields @@ -68,12 +70,12 @@ mutable struct IPNewtonState{T,Tx} <: AbstractBarrierState bgrad::BarrierStateVars{T} # gradient of slack and λ variables at current "position" bstep::BarrierStateVars{T} # search direction for slack and λ constr_c::Vector{T} # value of the user-supplied constraints at x - constr_J::Matrix{T} # value of the user-supplied Jacobian at x + constr_J::TJ # value of the user-supplied Jacobian at x ev::T # equality violation, ∑_i λ_Ei (c*_i - c_i) Optim.@add_linesearch_fields() # x_ls and alpha b_ls::BarrierLineSearchGrad{T} gtilde::Tx - Htilde # Positive Cholesky factorization of H from PositiveFactorizations.jl + Htilde::THtilde # Positive Cholesky factorization of H from PositiveFactorizations.jl end # TODO: Do we need this convert thing? (It seems to be used with `show(IPNewtonState)`) @@ -123,14 +125,20 @@ function initial_state(method::IPNewton, options, d::TwiceDifferentiable, constr s = similar(initial_x) f_x_previous = NaN f_x, g_x = value_gradient!(d, initial_x) - g .= g_x # needs to be a separate copy of g_x - H = Matrix{T}(undef, n, n) + g .= g_x # needs to be a separate copy of g_x Hd = Vector{Int8}(undef, n) - hessian!(d, initial_x) + hd = hessian!(d, initial_x) + H = similar(hd) copyto!(H, hessian(d)) + Htilde = cholesky(Positive, H, Val{true}) + # More constraints - constr_J = Array{T}(undef, mc, n) + if (method.conJstorage == :auto) + constr_J = Array{T}(undef, mc, n) + else + constr_J = similar(method.conJstorage) + end gtilde = copy(g) constraints.jacobian!(constr_J, initial_x) μ = T(1) @@ -147,7 +155,7 @@ function initial_state(method::IPNewton, options, d::TwiceDifferentiable, constr g, # Store current gradient in state.g (TODO: includes Lagrangian calculation?) T(NaN), # Store previous f in state.f_x_previous H, - 0, # will be replaced + nothing, # will be replaced Hd, similar(initial_x), # Maintain current x-search direction in state.s μ, @@ -163,7 +171,7 @@ function initial_state(method::IPNewton, options, d::TwiceDifferentiable, constr Optim.@initial_linesearch()..., # Maintain a cache for line search results in state.lsr b_ls, gtilde, - 0) + Htilde) Hinfo = (state.H, hessianI(initial_x, constraints, 1 ./ bstate.slack_c, 1)) initialize_μ_λ!(state, constraints.bounds, Hinfo, method.μ0) diff --git a/test/multivariate/solvers/constrained/ipnewton/interface.jl b/test/multivariate/solvers/constrained/ipnewton/interface.jl index 69e4ebcdb..730f84fab 100644 --- a/test/multivariate/solvers/constrained/ipnewton/interface.jl +++ b/test/multivariate/solvers/constrained/ipnewton/interface.jl @@ -49,4 +49,4 @@ end optimize(exponential, exponential_gradient!, exponential_hessian!, lb, ub, initial_x, IPNewton(), Optim.Options()) optimize(TwiceDifferentiable(od, initial_x), lb, ub, initial_x) optimize(TwiceDifferentiable(od, initial_x), lb, ub, initial_x, Optim.Options()) -end \ No newline at end of file +end diff --git a/test/multivariate/solvers/constrained/ipnewton/matrixtypes.jl b/test/multivariate/solvers/constrained/ipnewton/matrixtypes.jl new file mode 100644 index 000000000..def889a9f --- /dev/null +++ b/test/multivariate/solvers/constrained/ipnewton/matrixtypes.jl @@ -0,0 +1,92 @@ +using Optim, Test +import SparseArrays: sparse +import LinearAlgebra: Tridiagonal + +let + # Why do the tests fail when this function is defined inside the @testset? + """ + Create TwiceDifferentiable objective and TwiceDifferentiableConstraints + representing + + min x'A*x s.t. x[i]^2 + x[i+1]^2 = 1 for i=1:2:(n-1) + + """ + function problem(A::AbstractMatrix) + n = size(A,1) + (mod(n,2) == 0) || error("size(A,1) must be even") + ncon = n ÷ 2 + h0 = A + g0 = zeros(eltype(A),n) + x0 = zeros(n) + function grad!(g,x) + for i in 1:n + g[i] = 2*x'*A[:,i] + end + end + function hess!(h,x) + copyto!(h,A) + end + function obj(x) + return(x'*A*x) + end + f = Optim.TwiceDifferentiable(obj, grad!, + hess!, x0, zero(eltype(A)), g0, h0 ) + function con_c!(c,x) + for i in 1:ncon + c[i] = x[2*i-1]^2 + x[2*i]^2 + end + c + end + Jstore = similar(A,n÷2,n) + Jstore .= zero(eltype(A)) + function con_j!(J,x) + for i in 1:ncon + J[i,2*i-1] = 2x[2*i-1] + J[i,2*i] = 2x[2*i] + end + end + function con_hl!(h,x,λ) + for i in 1:n + h[i,i] += 2λ[(i+1)÷2] + end + end + lx = fill(-Inf,n) + ux = fill(Inf,n) + lc = ones(n÷2) + uc = ones(n÷2) + fc = Optim.TwiceDifferentiableConstraints(con_c!, con_j!, con_hl!, lx, ux, lc, uc) + return(f, fc, Jstore) + end + + + @testset "hessian and jacobian alternate types" begin + n = 8 + A = zeros(n,n) + x0 = ones(n) + for i in 1:n + A[i,i] = 2.0 + if (i>1) + A[i-1,i] = 1.0 + A[i,i-1] = 1.0 + end + end + + f,con,J = problem(A) + dense_sol = Optim.optimize(f,con,x0,Optim.IPNewton(conJstorage=J)) + + @testset "sparse" begin + f,con,J = problem(sparse(A)) + sparse_sol = Optim.optimize(f,con,x0,Optim.IPNewton(conJstorage=J)) + @test dense_sol.minimum ≈ sparse_sol.minimum + @test dense_sol.minimizer ≈ sparse_sol.minimizer + end + + @testset "tridiagonal" begin + f,con,J = problem(Tridiagonal(A)) + tridiagonal_sol = Optim.optimize(f,con,x0,Optim.IPNewton(conJstorage=J)) + @test dense_sol.minimum ≈ tridiagonal_sol.minimum + @test dense_sol.minimizer ≈ tridiagonal_sol.minimizer + end + end + +end