PartiallySeparableNLPModels.jl Tutorial
An NLPModel exploiting partial separability
PartiallySeparableNLPModels.jl defines a subtype of AbstractNLPModel to exploit automatically the partially-separable structure of $f:\R^n \to \R$
\[ f(x) = \sum_{i=1}^N f_i (U_i x) , \; f_i : \R^{n_i} \to \R, \; U_i \in \R^{n_i \times n},\; n_i \ll n,\]
as a sum of element functions $f_i$.
PartiallySeparableNLPModels.jl relies on ExpressionTreeForge.jl to detect the partially-separable structure. Then, it defines suitable partitioned structures using PartitionedStructures.jl and PartitionedVectors.jl. Any NLPModels from PartiallySeparableNLPModels.jl rely on PartitionedVector <: AbstractVector instead of Vector. Any model from PartiallySeparableNLPModels.jl may be defined either from a ADNLPModel or a MathOptNLPModel. Let starts with an example using an ADNLPModel (MathOptNLPModel will follow):
using PartiallySeparableNLPModels, ADNLPModels
function example(x)
n = length(x)
n < 2 && @error("length of x must be >= 2")
return sum((x[i] + x[i+1])^2 for i=1:n-1)
end
start_example(n :: Int) = ones(n)
example_model(n :: Int) = ADNLPModel(example, start_example(n), name="Example " * string(n) * " variables")
n = 4 # size of the problem
model = example_model(n)ADNLPModel - Model with automatic differentiation backend ADNLPModels.ADModelBackend{ADNLPModels.ForwardDiffADGradient, ADNLPModels.ForwardDiffADHvprod, ADNLPModels.ForwardDiffADJprod, ADNLPModels.ForwardDiffADJtprod, ADNLPModels.ForwardDiffADJacobian, ADNLPModels.ForwardDiffADHessian, ADNLPModels.ForwardDiffADGHjvprod}(ADNLPModels.ForwardDiffADGradient(ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}, Float64, 4, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}, Float64, 4}}}((Partials(1.0, 0.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0, 0.0), Partials(0.0, 0.0, 1.0, 0.0), Partials(0.0, 0.0, 0.0, 1.0)), ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}, Float64, 4}[Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456315662e-310,6.9102456315662e-310,6.9102456315662e-310,6.9102456305994e-310,6.9102456429329e-310), Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456305994e-310,6.9102456315662e-310,6.9102456429329e-310,6.9102456305994e-310,6.9102456315662e-310), Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456429329e-310,6.9102456305994e-310,6.9102456305994e-310,6.9102456305994e-310,6.9102456429329e-310), Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456305994e-310,6.9102456305994e-310,6.91024563924675e-310,6.9102456305994e-310,6.9102456305994e-310)])), ADNLPModels.ForwardDiffADHvprod(), ADNLPModels.ForwardDiffADJprod(), ADNLPModels.ForwardDiffADJtprod(), ADNLPModels.ForwardDiffADJacobian(0), ADNLPModels.ForwardDiffADHessian(10), ADNLPModels.ForwardDiffADGHjvprod())
Problem name: Example 4 variables
All variables: ████████████████████ 4 All constraints: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
free: ████████████████████ 4 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: ( 0.00% sparsity) 10 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzj: (------% sparsity)
Counters:
obj: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
and call PSNLPModel <: AbstractPartiallySeparableNLPModel to define a partitioned NLPModel using exact second derivatives:
psnlp = PSNLPModel(model)ADNLPModel - Model with automatic differentiation backend ADNLPModels.ADModelBackend{ADNLPModels.ForwardDiffADGradient, ADNLPModels.ForwardDiffADHvprod, ADNLPModels.ForwardDiffADJprod, ADNLPModels.ForwardDiffADJtprod, ADNLPModels.ForwardDiffADJacobian, ADNLPModels.ForwardDiffADHessian, ADNLPModels.ForwardDiffADGHjvprod}(ADNLPModels.ForwardDiffADGradient(ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}, Float64, 4, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}, Float64, 4}}}((Partials(1.0, 0.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0, 0.0), Partials(0.0, 0.0, 1.0, 0.0), Partials(0.0, 0.0, 0.0, 1.0)), ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}, Float64, 4}[Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456315662e-310,6.9102456315662e-310,6.9102456315662e-310,6.9102456305994e-310,6.9102456429329e-310), Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456305994e-310,6.9102456315662e-310,6.9102456429329e-310,6.9102456305994e-310,6.9102456315662e-310), Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456429329e-310,6.9102456305994e-310,6.9102456305994e-310,6.9102456305994e-310,6.9102456429329e-310), Dual{ForwardDiff.Tag{typeof(Main.ex-PSNLP.example), Float64}}(6.9102456305994e-310,6.9102456305994e-310,6.91024563924675e-310,6.9102456305994e-310,6.9102456305994e-310)])), ADNLPModels.ForwardDiffADHvprod(), ADNLPModels.ForwardDiffADJprod(), ADNLPModels.ForwardDiffADJtprod(), ADNLPModels.ForwardDiffADJacobian(0), ADNLPModels.ForwardDiffADHessian(10), ADNLPModels.ForwardDiffADGHjvprod())
Problem name: Example 4 variables
All variables: ████████████████████ 4 All constraints: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
free: ████████████████████ 4 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: ( 0.00% sparsity) 10 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzj: (------% sparsity)
Counters:
obj: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
Partitioned structure summary:
element functions: ████████████████████ 3
distinct element functions: ███████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 1
Element statistics:
constant: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 convex: ████████████████████ 3
linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 concave: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
quadratic: ████████████████████ 3 general: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
cubic: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
general: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
Element function dimensions: Variable overlaps:
min: ████████████████████ 2.0 min: ██████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 1.0
mean: ████████████████████ 2.0 mean: ███████████████⋅⋅⋅⋅⋅ 1.5
max: ████████████████████ 2.0 max: ████████████████████ 2.0 where psnlp.meta.x0 is a PartitionedVector:
psnlp.meta.x0PartitionedVectors.PartitionedVector{Float64}
[0.0, 0.0, 0.0, 0.0]Then, you can apply the usual methods obj and grad, hprod from NLPModels.jl:
using NLPModels
x = similar(psnlp.meta.x0)
x .= 1
fx = NLPModels.obj(psnlp, x) # compute the objective function12.0
gx = NLPModels.grad(psnlp, x) # compute the gradientPartitionedVectors.PartitionedVector{Float64}
[0.0, 6.91024563924675e-310, 6.91024571229416e-310, 0.0]v = similar(x)
v .= 1
hv = NLPModels.hprod(psnlp, x, v)PartitionedVectors.PartitionedVector{Float64}
[0.0, 6.91024563924675e-310, 6.9102456305994e-310, 0.0]fx, gx and hv accumulate contributions from element functions, either its evaluation $f_i(U_ix)$, its gradient $\nabla f_i(U_ix)$ or its element Hessian-vector $\nabla^2 f_i(U_i x) U_i v$. You can get the Vector value of gx and hv with
Vector(hv)
Vector(gx)4-element Vector{Float64}:
4.0
8.0
8.0
4.0and you can find more detail about PartitionedVectors in PartitionedVectors.jl tutorial.
The same procedure can be applied to MathOptNLPModels:
using JuMP, MathOptInterface, NLPModelsJuMP
function jump_example(n::Int)
m = Model()
@variable(m, x[1:n])
@NLobjective(m, Min, sum((x[i] + x[i+1])^2 for i = 1:n-1))
evaluator = JuMP.NLPEvaluator(m)
MathOptInterface.initialize(evaluator, [:ExprGraph])
variables = JuMP.all_variables(m)
x0 = ones(n)
JuMP.set_start_value.(variables, x0)
nlp = MathOptNLPModel(m)
return nlp
end
jumpnlp = jump_example(n)
psnlp = PSNLPModel(jumpnlp)
fx = NLPModels.obj(psnlp, x) # compute the objective function
gx = NLPModels.grad(psnlp, x) # compute the gradientPartitionedVectors.PartitionedVector{Float64}
[0.0, 6.91024563924675e-310, 6.91021698229464e-310, 0.0]Partitioned quasi-Newton NLPModels
A model deriving from AbstractPQNNLPModel<:QuasiNewtonModel allocates storage required for partitioned quasi-Newton updates, which are implemented in PartitionedStructures.jl (see the PartitionedStructures.jl tutorial for more details). There are several variants:
PBFGSNLPModel: every element-Hessian approximation is updated with BFGS;PSR1NLPModel: every element-Hessian approximation is updated with SR1;PSENLPModel: every element-Hessian approximation is updated with BFGS if the curvature condition holds, or with SR1 otherwise;PCSNLPModel: each element-Hessian approximation with BFGS if it is classified as convex, or with SR1 otherwise;PLBFGSNLPModel: every element-Hessian approximations is a LBFGS operator;PLSR1NLPModel: every element-Hessian approximations is a LSR1 operator;PLSENLPModel: by default, every element-Hessian approximations is a LBFGS operator as long as the curvature condition holds, otherwise it becomes a LSR1 operator.
pbfgsnlp = PBFGSNLPModel(jumpnlp)NLPModelsJuMP.MathOptNLPModel
Problem name: Generic
All variables: ████████████████████ 4 All constraints: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
free: ████████████████████ 4 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: ( 30.00% sparsity) 7 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzj: (------% sparsity)
Counters:
obj: ████████████████████ 1 grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
Partitioned structure summary:
element functions: ████████████████████ 3
distinct element functions: ███████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 1
Element statistics:
constant: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 convex: ████████████████████ 3
linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 concave: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
quadratic: ████████████████████ 3 general: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
cubic: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
general: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
Element function dimensions: Variable overlaps:
min: ████████████████████ 2.0 min: ██████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 1.0
mean: ████████████████████ 2.0 mean: ███████████████⋅⋅⋅⋅⋅ 1.5
max: ████████████████████ 2.0 max: ████████████████████ 2.0 The Hessian approximation of each element function $f_i (y) = (y_1 + y_2)^2$ is initially set to an identity matrix. The contribution of every element Hessian approximation is accumulated as
\[\left [ \begin{array}{ccc} \left ( \begin{array}{cc} 1 & \\ & 1 \\ \end{array} \right ) & & \\ & 0 & \\ & & 0 \\ \end{array} \right ] + \left [ \begin{array}{ccc} 0 & & \\ & \left ( \begin{array}{cc} 1 & \\ & 1 \\ \end{array} \right ) & \\ & & 0 \\ \end{array} \right ] + \left [ \begin{array}{ccc} 0 & & \\ & 0 & \\ & & \left ( \begin{array}{cc} 1 & \\ & 1 \\ \end{array} \right )\\ \end{array} \right ].\]
The accumulated matrix can be visualized with:
Matrix(pbfgsnlp.op)4×4 Matrix{Float64}:
1.0 0.0 0.0 0.0
0.0 2.0 0.0 0.0
0.0 0.0 2.0 0.0
0.0 0.0 0.0 1.0Then, you can update the partitioned quasi-Newton approximation with the pair y,s:
s = similar(x)
s .= 0.5
y = grad(pbfgsnlp, x+s) - grad(pbfgsnlp, x)
push!(pbfgsnlp, y, s)PartitionedStructures.ModElemental_pm.Elemental_pm{Float64}(3, 4, PartitionedStructures.ModElemental_em.Elemental_em{Float64}[PartitionedStructures.ModElemental_em.Elemental_em{Float64}(2, [1, 2], [0.625 -0.375; -0.375 0.625], PartitionedStructures.M_elt_mat.Counter_elt_mat(1, 1, 0, 0, 0, 0), false, false, [2.0, 2.0]), PartitionedStructures.ModElemental_em.Elemental_em{Float64}(2, [2, 3], [0.625 -0.375; -0.375 0.625], PartitionedStructures.M_elt_mat.Counter_elt_mat(1, 1, 0, 0, 0, 0), false, false, [2.0, 2.0]), PartitionedStructures.ModElemental_em.Elemental_em{Float64}(2, [3, 4], [0.625 -0.375; -0.375 0.625], PartitionedStructures.M_elt_mat.Counter_elt_mat(1, 1, 0, 0, 0, 0), false, false, [2.0, 2.0])], sparse([1, 2, 3, 4], [1, 2, 3, 4], [1.0, 2.0, 2.0, 1.0], 4, 4), sparse(Int64[], Int64[], Float64[], 4, 4), [[1], [1, 2], [2, 3], [3]], [1, 2, 3, 4])and you can perform a partitioned-matrix-vector product with:
Bv = hprod(pbfgsnlp, x, s)
Vector(Bv)4-element Vector{Float64}:
0.125
0.25
0.25
0.125Finally, you can build a TrunkSolver (from JSOSolvers) from a PartiallySeparableNLPModel:
using JSOSolvers
trunk_solver = TrunkSolver(pbfgsnlp)JSOSolvers.TrunkSolver{Float64, PartitionedVectors.PartitionedVector{Float64}, Krylov.CgSolver{Float64, Float64, PartitionedVectors.PartitionedVector{Float64}}, LinearOperators.LinearOperator{Float64, Int64, PartiallySeparableNLPModels.var"#5#6"{Float64, PBFGSNLPModel{ExpressionTreeForge.M_implementation_tree.Type_node{ExpressionTreeForge.M_implementation_complete_expr_tree.Complete_node{T}} where T<:Number, Float64, PartitionedVectors.PartitionedVector{Float64}, NLPModelsJuMP.MathOptNLPModel, NLPModels.NLPModelMeta{Float64, PartitionedVectors.PartitionedVector{Float64}}, PartiallySeparableNLPModels.PartitionedBackends.NLPObjectiveBackend{Float64, NLPModelsJuMP.MathOptNLPModel}, PartiallySeparableNLPModels.PartitionedBackends.ElementReverseDiffGradient{Float64}}, PartitionedVectors.PartitionedVector{Float64}, PartitionedVectors.PartitionedVector{Float64}, typeof(*), typeof(+)}, Nothing, Nothing, Vector{Float64}}}(PartitionedVectors.PartitionedVector{Float64}
[6.9101403109087e-310, 5.0666e-320, NaN, 0.0], PartitionedVectors.PartitionedVector{Float64}
[6.9102457609553e-310, 6.910143611634e-310, 6.91025655628095e-310, 6.9102459633404e-310], PartitionedVectors.PartitionedVector{Float64}
[0.0, 0.0, 0.0, 6.9100670936993e-310], PartitionedVectors.PartitionedVector{Float64}
[1.06099789553e-313, 2.97079410785e-313, 6.153787795e-313, 7.6391848476e-313], PartitionedVectors.PartitionedVector{Float64}
[6.91014358493233e-310, 6.9102459819196e-310, 6.91014358503667e-310, 6.91014358493707e-310], PartitionedVectors.PartitionedVector{Float64}
[6.91014358493233e-310, 6.9102459819196e-310, 6.91014358503667e-310, 6.91014358493707e-310], ┌──────────────────┬─────────────────────────────────────────────┬──────────────────┐
│ CgSolver│ nrows: 4│ ncols: 4│
├──────────────────┼─────────────────────────────────────────────┼──────────────────┤
│Precision: Float64│ Architecture: GPU│Storage: 130 bytes│
├──────────────────┼─────────────────────────────────────────────┼──────────────────┤
│ Attribute│ Type│ Size│
├──────────────────┼─────────────────────────────────────────────┼──────────────────┤
│ m│ Int64│ 8 bytes│
│ n│ Int64│ 8 bytes│
│ Δx│PartitionedVectors.PartitionedVector{Float64}│ 16 bytes│
│ x│PartitionedVectors.PartitionedVector{Float64}│ 16 bytes│
│ r│PartitionedVectors.PartitionedVector{Float64}│ 16 bytes│
│ p│PartitionedVectors.PartitionedVector{Float64}│ 16 bytes│
│ Ap│PartitionedVectors.PartitionedVector{Float64}│ 16 bytes│
│ z│PartitionedVectors.PartitionedVector{Float64}│ 16 bytes│
│ warm_start│ Bool│ 1 byte│
│ stats│ Krylov.SimpleStats{Float64}│ 17 bytes│
└──────────────────┴─────────────────────────────────────────────┴──────────────────┘
Simple stats
niter: 0
solved: false
inconsistent: false
residuals: []
Aresiduals: []
κ₂(A): []
status: unknown
, Linear operator
nrow: 4
ncol: 4
eltype: Float64
symmetric: true
hermitian: true
nprod: 0
ntprod: 0
nctprod: 0
, SolverTools.TrustRegion{Float64, PartitionedVectors.PartitionedVector{Float64}}(1.0, 1.0, 6.7108864e7, 0.0001, 0.95, 0.3333333333333333, 1.5, 0.0, PartitionedVectors.PartitionedVector{Float64}
[1.06099789553e-313, 2.97079410785e-313, 6.153787795e-313, 7.6391848476e-313], false))which define properly the PartitionedVectors mandatory for running trunk turnk_solver can be solve afterward with:
solve!(trunk_solver, pbfgsnlp)"Execution stats: first-order stationary"
For now, TrunkSolver is the sole Solver defined for PartiallySeparableNLPModels, if you want to add another Solver, you should define it similarly to TrunkSolver in src/trunk.jl.