Optimized hess coord

benchmark/figures_benchmark_no_tex.jl

@@ -13,7 +13,7 @@ using Random
 gr() #pgfplots()
 
 function figure()
-  names = ["2021-05-21__FPS_FPSFF_knitro_ipopt_45"]
+  names = ["2021-05-27_all__FPS_FPSFF_knitro_ipopt_105"]
 
   tod = string(today()) * ""
   dsolvers = [:ipopt, :knitro, :FPS, :FPSFF]
@@ -53,10 +53,10 @@ function figure()
       #p = performance_profile(stats, cost)
       p = performance_profile(
         stats,
-        cost, 
+        cost,
         title = perf_title(col),
         legend = :bottomright,
-        linestyles=styles,
+        linestyles = styles,
       )
       #savefig("$(tod)_$(list)_perf-$col.tex")
       png("$(tod)_$(list)_perf-$col")
@@ -70,7 +70,7 @@ function figure()
         ipairs = 0
         # combinations of solvers 2 by 2
         for i = 2:nsolvers
-          for j = 1:(i-1)
+          for j = 1:(i - 1)
             ipairs += 1
             pair = [solvers[i], solvers[j]]
             dfs = (stats[solver] for solver in pair)

src/model-Fletcherpenaltynlp.jl

@@ -15,6 +15,10 @@ import NLPModels:
 
 include("solve_two_systems_struct.jl")
 
+#=
+We need to create a workspace structure Val1/Val2/matrix-free Val2
+=#
+
 """
 We consider here the implementation of Fletcher's exact penalty method for
 the minimization problem:
@@ -77,6 +81,17 @@ mutable struct FletcherPenaltyNLP{
   Jcρ::T
   Jv::T
   Ss::Array{S, 2} # only when Val(1)
+  Arows
+  Acols
+  Avals
+  Hrows
+  Hcols
+  Hvals
+  Ax # jacobian matrix (remove for matrix free implementation) (ncon x nvar)
+  F
+  Pt # dense projector (nvar x nvar)
+  invAtA # dense pseudo-inverse of AA' (ncon x ncon)
+  Hs # = Symmetric(hess(nlp.nlp, x, -ys), :L) (nvar x nvar)
 
   # Problem parameter
   σ::P
@@ -100,7 +115,7 @@ function FletcherPenaltyNLP(
   x0::AbstractVector{S};
   qds = LDLtSolver(nlp, S(0)),
 ) where {S}
-  nvar = nlp.meta.nvar
+  nvar, ncon = nlp.meta.nvar, nlp.meta.ncon
 
   meta = NLPModelMeta(
     nvar,
@@ -113,6 +128,19 @@ function FletcherPenaltyNLP(
     name = "Fletcher penalization of $(nlp.meta.name)",
   )
 
+  Arows, Acols = jac_structure(nlp)
+  Avals = rand(S, nlp.meta.nnzj)
+  Hrows, Hcols = hess_structure(nlp)
+  Hvals = Vector{S}(undef, nlp.meta.nnzh)
+
+  tempA = sparse(Arows, Acols, Avals, ncon, nvar)
+  invAtA = sparse(tempA * tempA' + 1e-14 * Matrix(I, ncon, ncon))
+  F = if nlp.meta.ncon > 0
+    lu(invAtA)
+  else
+    nothing
+  end
+
   return FletcherPenaltyNLP(
     meta,
     Counters(),
@@ -134,6 +162,17 @@ function FletcherPenaltyNLP(
     Vector{S}(undef, nlp.meta.nvar),
     Vector{S}(undef, nlp.meta.ncon),
     Array{S, 2}(undef, nlp.meta.ncon, nlp.meta.nvar),
+    Arows,
+    Acols,
+    Avals,
+    Hrows,
+    Hcols,
+    Hvals,
+    Array{S, 2}(undef, nlp.meta.ncon, nlp.meta.nvar),
+    F,
+    Array{S, 2}(undef, nlp.meta.nvar, nlp.meta.nvar),
+    invAtA,
+    Array{S, 2}(undef, nlp.meta.nvar, nlp.meta.nvar),
     σ,
     zero(typeof(σ)),
     zero(typeof(σ)),
@@ -156,7 +195,7 @@ function FletcherPenaltyNLP(
   x0::AbstractVector{S};
   qds = LDLtSolver(nlp, S(0)), #IterativeSolver(nlp, S(NaN)),
 ) where {S}
-  nvar = nlp.meta.nvar
+  nvar, ncon = nlp.meta.nvar, nlp.meta.ncon
 
   meta = NLPModelMeta(
     nvar,
@@ -169,6 +208,20 @@ function FletcherPenaltyNLP(
     name = "Fletcher penalization of $(nlp.meta.name)",
   )
   counters = Counters()
+
+  Arows, Acols = jac_structure(nlp)
+  Avals = rand(S, nlp.meta.nnzj)
+  Hrows, Hcols = hess_structure(nlp)
+  Hvals = Vector{S}(undef, nlp.meta.nnzh)
+
+  tempA = sparse(Arows, Acols, Avals, ncon, nvar)
+  invAtA = sparse(tempA * tempA' + 1e-14 * Matrix(I, ncon, ncon))
+  F = if nlp.meta.ncon > 0
+    lu(invAtA)
+  else
+    nothing
+  end
+
   return FletcherPenaltyNLP(
     meta,
     counters,
@@ -190,6 +243,17 @@ function FletcherPenaltyNLP(
     Vector{S}(undef, nlp.meta.nvar),
     Vector{S}(undef, nlp.meta.ncon),
     Array{S, 2}(undef, nlp.meta.ncon, nlp.meta.nvar),
+    Arows,
+    Acols,
+    Avals,
+    Hrows,
+    Hcols,
+    Hvals,
+    Array{S, 2}(undef, nlp.meta.ncon, nlp.meta.nvar),
+    F,
+    Array{S, 2}(undef, nlp.meta.nvar, nlp.meta.nvar),
+    invAtA,
+    Array{S, 2}(undef, nlp.meta.nvar, nlp.meta.nvar),
     σ,
     ρ,
     δ,
@@ -364,11 +428,11 @@ function hess_structure!(
 end
 
 function hess_coord!(
-  nlp::FletcherPenaltyNLP,
+  nlp::FletcherPenaltyNLP{S, Tt, Val{1}, P, QDS},
   x::AbstractVector{T},
   vals::AbstractVector{T};
   obj_weight::Real = one(T),
-) where {T}
+) where {T, S, Tt, P, QDS}
   @lencheck nlp.meta.nvar x
   @lencheck nlp.meta.nnzh vals
   increment!(nlp, :neval_hess)
@@ -381,9 +445,13 @@ function hess_coord!(
   g = nlp.gx
   c = nlp.cx
   A = jac(nlp.nlp, x) # If used, this should be allocated probably
+  # nlp.Ax .= jac(nlp.nlp, x) # do in-place ?
+  # A = nlp.Ax
   σ, ρ, δ = nlp.σ, nlp.ρ, nlp.δ
 
   Hs = Symmetric(hess(nlp.nlp, x, -ys), :L)
+  # nlp.Hs .= Symmetric(hess(nlp.nlp, x, -ys), :L) # do in-place ?
+  # Hs = nlp.Hs
   Im = Matrix(I, ncon, ncon)
   τ = T(max(nlp.δ, 1e-14))
   invAtA = pinv(Matrix(A * A') + τ * Im) #inv(Matrix(A*A') + τ * Im) # Euh... wait !
@@ -397,11 +465,136 @@ function hess_coord!(
     Hx += Hc + T(ρ) * A' * A
   end
 
-  if nlp.hessian_approx == Val(1)
-    for j = 1:ncon
-      nlp.Ss[j, :] = gs' * Symmetric(jth_hess(nlp.nlp, x, j), :L)
+  # for 1st hessian approximation only
+  for j = 1:ncon
+    nlp.Ss[j, :] = gs' * Symmetric(jth_hess(nlp.nlp, x, j), :L) # could we use hprod ?
+  end
+  Hx += -AinvAtA * nlp.Ss - nlp.Ss' * invAtA * A
+
+  #=
+  if nlp.η > 0.0
+    In = Matrix(I, nvar, nvar)
+    Hx += T(nlp.η) * In
+  end
+  =#
+
+  k = 1
+  for j = 1:nvar
+    for i = j:nvar
+      vals[k] = obj_weight * Hx[i, j]
+      if i == j && nlp.η > 0.0
+        vals[k] += obj_weight * T(nlp.η)
+      end
+      k += 1
+    end
+  end
+
+  return vals
+end
+
+function _compute_pseudo_inverse_dense!(nlp, x::AbstractVector{T}) where {T}
+  nvar = nlp.meta.nvar
+  ncon = nlp.nlp.meta.ncon
+  A = jac(nlp.nlp, x)
+  # put in a separate function
+  # Im = Matrix(I, ncon, ncon)
+  τ = T(max(nlp.δ, 1e-14))
+  # invAtA = pinv(Matrix(A * A') + τ * Im) #inv(Matrix(A*A') + τ * Im) # Euh... wait !
+  nlp.invAtA .= Matrix(I, ncon, ncon) # spdiagm(ones(ncon)) #  
+  mul!(nlp.invAtA, A, A', 1, τ)
+  # lu!(nlp.F, nlp.invAtA)
+  nlp.invAtA = pinv(Matrix(nlp.invAtA)) # in-place
+  # AinvAtA = A' * invAtA
+  # Pt = AinvAtA * A
+  mul!(nlp.Ax, nlp.invAtA, A)
+  # ldiv!(nlp.Ax, nlp.F, A)
+  mul!(nlp.Pt, A', nlp.Ax)
+  return nlp.Pt
+end
+
+function _compute_pseudo_inverse_sparse!(nlp, x::AbstractVector{T}) where {T}
+  nvar = nlp.meta.nvar
+  ncon = nlp.nlp.meta.ncon
+  jac_coord!(nlp.nlp, x, nlp.Avals)
+  A = sparse(nlp.Arows, nlp.Acols, nlp.Avals, ncon, nvar)
+  # nlp.Ax .= jac(nlp.nlp, x) # do in-place ?
+  # A = nlp.Ax
+  # put in a separate function
+  # Im = Matrix(I, ncon, ncon)
+  τ = T(max(nlp.δ, 1e-14))
+  # invAtA = pinv(Matrix(A * A') + τ * Im) #inv(Matrix(A*A') + τ * Im) # Euh... wait !
+  nlp.invAtA .= Matrix(I, ncon, ncon) # spdiagm(ones(ncon)) #  
+  mul!(nlp.invAtA, A, A', 1, τ)
+  lu!(nlp.F, nlp.invAtA)
+  # AinvAtA = A' * invAtA
+  # Pt = AinvAtA * A
+  ldiv!(nlp.Ax, nlp.F, A)
+  mul!(nlp.Pt, A', nlp.Ax)
+  return A, nlp.Pt
+end
+
+function hess_coord!(
+  nlp::FletcherPenaltyNLP{S, Tt, Val{2}, P, QDS},
+  x::AbstractVector{T},
+  vals::AbstractVector{T};
+  obj_weight::Real = one(T),
+) where {T, S, Tt, P, QDS}
+  @lencheck nlp.meta.nvar x
+  @lencheck nlp.meta.nnzh vals
+  increment!(nlp, :neval_hess)
+
+  nvar = nlp.meta.nvar
+  ncon = nlp.nlp.meta.ncon
+
+  gs, ys, _, _ = _compute_ys_gs!(nlp, x)
+  f = nlp.fx
+  g = nlp.gx
+  c = nlp.cx
+  σ, ρ, δ = nlp.σ, nlp.ρ, nlp.δ
+
+  Hs = Symmetric(hess(nlp.nlp, x, -ys), :L) # do in-place ?
+  # Hs = nlp.Hs
+
+  if ncon > 0
+    A = jac(nlp.nlp, x) # If used, this should be allocated probably
+    # jac_coord!(nlp.nlp, x, nlp.Avals)
+    # A = sparse(nlp.Arows, nlp.Acols, nlp.Avals, ncon, nvar)
+    # nlp.Ax .= jac(nlp.nlp, x) # do in-place ?
+    # A = nlp.Ax
+    # put in a separate function
+    # Im = Matrix(I, ncon, ncon)
+    τ = T(max(nlp.δ, 1e-14))
+    # invAtA = pinv(Matrix(A * A') + τ * Im) #inv(Matrix(A*A') + τ * Im) # Euh... wait !
+    nlp.invAtA .= Matrix(I, ncon, ncon) # spdiagm(ones(ncon)) #  
+    mul!(nlp.invAtA, A, A', 1, τ)
+    # lu!(nlp.F, nlp.invAtA)
+    nlp.invAtA = pinv(Matrix(nlp.invAtA)) # in-place
+    # AinvAtA = A' * invAtA
+    # Pt = AinvAtA * A
+    mul!(nlp.Ax, nlp.invAtA, A)
+    # ldiv!(nlp.Ax, nlp.F, A)
+    mul!(nlp.Pt, A', nlp.Ax)
+    # end of separate function
+    # A, nlp.Pt = _compute_pseudo_inverse_dense!(nlp, x) # doesn't work...
+    # A, Pt = _compute_pseudo_inverse_sparse!(nlp, x)
+
+    # mul!(C, A, B, α, β) -> C = A * B * α + C * β
+    # Hx = Hs - nlp.Pt * Hs - Hs * nlp.Pt + 2 * T(σ) * nlp.Pt
+    mul!(nlp.Hs, nlp.Pt, Hs, -1, 0) # - nlp.Pt * Hs
+    mul!(nlp.Hs, Hs, nlp.Pt, -1, 1) # nlp.Hs -= nlp.Pt * Hs
+    @. nlp.Hs += Hs + 2 * T(σ) * nlp.Pt
+    Hx = nlp.Hs
+    #regularization term
+    if ρ > 0.0
+      Hc = hess(nlp.nlp, x, c * T(ρ), obj_weight = zero(T))
+      Hx += Hc + T(ρ) * A' * A
+    end
+  else
+    Hx = Hs
+    if ρ > 0.0
+      Hc = hess(nlp.nlp, x, c * T(ρ), obj_weight = zero(T))
+      Hx += Hc
     end
-    Hx += -AinvAtA * nlp.Ss - nlp.Ss' * invAtA * A
   end
 
   #=

src/solve_linear_system.jl

@@ -115,7 +115,7 @@ function solve_two_mixed(
   if !stats1.solved
     @warn "Failed solving 1st linear system lsqr in mixed."
   end
-  
+
   if nlp.δ != 0.0
     (p2, q2, stats2) = craig!(
       nlp.qdsolver.solver_struct_least_norm,
@@ -131,8 +131,8 @@ function solve_two_mixed(
     )
   else
     (p2, q2, stats2) = craig!(
-      nlp.qdsolver.solver_struct_least_norm, 
-      nlp.Aop, 
+      nlp.qdsolver.solver_struct_least_norm,
+      nlp.Aop,
       -rhs2,
       atol = nlp.qdsolver.ln_atol,
       rtol = nlp.qdsolver.ln_rtol,

src/solve_two_systems_struct.jl

@@ -90,8 +90,8 @@ function IterativeSolver(
   ln_btol::T = √eps(T),
   ln_conlim::T = 1 / √eps(T),
   ln_itmax::Integer = 5 * (nlp.meta.ncon + nlp.meta.nvar),
-  ne_atol::T = √eps(T)/100,
-  ne_rtol::T = √eps(T)/100,
+  ne_atol::T = √eps(T) / 100,
+  ne_rtol::T = √eps(T) / 100,
   ne_ratol::T = zero(T),
   ne_rrtol::T = zero(T),
   ne_etol::T = √eps(T),

test/problems/controlinvestment.jl

@@ -0,0 +1,18 @@
+function controlinvestment_autodiff(args...; n::Int=200, type::Val{T}=Val(Float64), kwargs...) where T
+  N = div(n, 2)
+  h = 1/N
+  x0 = 1.0
+  gamma = 3
+  function f(y)
+    x, u = y[1:N], y[N+1:end]
+    return 0.5 * h * sum((u[k] - 1) * x[k] + (u[k+1] - 1) * x[k+1] for k = 1 : N-1)
+  end
+  function c(y)
+    x, u = y[1:N], y[N+1:end]
+    return [x[k+1] - x[k] - 0.5 * h * gamma * (u[k] * x[k] + u[k+1] * x[k+1]) for k = 1 : N-1]
+  end
+  lvar = vcat(x0, -Inf*ones(T, N-1), zeros(T, N))
+  uvar = vcat(x0, Inf*ones(T, N-1), ones(T, N))
+  xi = vcat(ones(T, N), zeros(T, N))
+  return ADNLPModel(f, xi, lvar, uvar, c, zeros(T, N-1), zeros(T, N-1), name="controlinvestment_autodiff")
+end