IPNewton() sparsity support

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),

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)

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
+end

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