Skip to content

Commit

Permalink
Made IPNewton() compatible with alternate types (e.g. sparse) of hess…
Browse files Browse the repository at this point in the history
…ians and jacobians
  • Loading branch information
schrimpf committed Feb 4, 2022
1 parent cdebfca commit f5835da
Show file tree
Hide file tree
Showing 4 changed files with 118 additions and 18 deletions.
6 changes: 3 additions & 3 deletions src/multivariate/solvers/constrained/ipnewton/interior.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
36 changes: 22 additions & 14 deletions src/multivariate/solvers/constrained/ipnewton/ipnewton.jl
Original file line number Diff line number Diff line change
@@ -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"
Expand Down Expand Up @@ -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
Expand All @@ -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)`)
Expand Down Expand Up @@ -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)
Expand All @@ -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
μ,
Expand All @@ -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)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
end
92 changes: 92 additions & 0 deletions test/multivariate/solvers/constrained/ipnewton/matrixtypes.jl
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit f5835da

Please sign in to comment.