PartiallySeparableSolvers.jl Tutorial
Solvers exploiting partial separability
PartiallySeparableSolvers.jl defines a solver exploiting 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,\]
where $f$ is the sum of element functions $f_i$. This solver exploits the partitioned structure of the gradient
\[\nabla f(x) = \sum_{i=1}^N U_i^\top \nabla f_i (U_i x),\]
and the Hessian
\[\nabla^2 f(x) = \sum_{i=1}^N U_i^\top \nabla^2 f_i (U_i x) U_i,\]
to maintain a partitioned a quasi-Newton approximation
\[B = \sum_{i=1}^N U_i^\top B_{i} U_i,\]
where each $B_{i} \approx \nabla^2 f_i$ has size $n_i \times n_i$.
The usual quasi-Newton updates, such as SR1, BFGS, or their limited-memory variants LSR1 or LBFGS, approximate $\nabla^2 f(x)\approx B$ by low rank updates, which inevitably result in dense approximations. By relying on element Hessian approximations $B_i$ updated at each iteration, the partitioned updates respect the Hessian sparsity structure and perform updates of rank $\min(\lambda N,n), \; \lambda = 1,2$ depending on the update applied to $B_i$.
See PartitionedStructures.jl tutorial for more details about partitioned quasi-Newton approximations and how they compare with standard updates. Some of these partitioned quasi-Newton methods are detailed in the reference below, and some are new (see PartitionedStructures.jl).
Reference
- A. Griewank and Ph. L. Toint, Partitioned variable metric updates for large structured optimization problems, Numerische Mathematik volume, 39, pp. 119–137, 1982.
Running a partitioned quasi-Newton solver
For now, PartiallySeparableSolvers.jl supports ADNLPModels, which are defined by pure julia code, or MathOptNLPModels, which are based on JuMP models. Regardless of which model is used, the solver returns a GenericExecutionStats.
An ADNLPModel
Let's first define the function that we seek to minimize and wrap it into an ADNLPModel:
using ADNLPModels
function f(x)
n = length(x)
n < 2 && @error("length of x must be >= 2")
return sum((x[j] + x[j+1])^2 for j=1:n-1)
end
model = ADNLPModel(f, ones(10))ADNLPModel - Model with automatic differentiation backend ADNLPModels.ForwardDiffAD{ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}, Float64, 10}}}}(55, 0, ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}, Float64, 10}}}((Partials(1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Partials(0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Partials(0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Partials(0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0), Partials(0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0), Partials(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0), Partials(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0), Partials(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0), Partials(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0)), ForwardDiff.Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(1.0e-323,5.0e-324,6.92740947939756e-310,2.0e-323,1.0e-323,6.9274095445984e-310,2.5e-323,1.5e-323,6.9274095445984e-310,3.5e-323,2.0e-323), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(6.92740947939756e-310,1.1e-322,3.0e-323,6.92732565929036e-310,4.0e-323,3.0e-323,6.92740947939756e-310,7.0e-323,3.5e-323,2.33419537006e-312,4.4e-323), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(3.5e-323,2.6312747808e-312,6.4e-323,4.0e-323,2.121995790965e-312,5.0e-323,4.0e-323,2.10077583306e-312,5.4e-323,4.4e-323,3.182993686453e-312), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(6.4e-323,5.0e-323,3.628612802556e-312,7.0e-323,6.0e-323,3.96813212911e-312,1.0e-322,6.4e-323,4.095451876563e-312,7.4e-323,6.4e-323), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(4.392531287303e-312,9.4e-323,7.0e-323,3.883252297466e-312,8.0e-323,7.0e-323,3.86203233956e-312,8.4e-323,7.4e-323,4.944250192954e-312,9.4e-323), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(8.0e-323,5.389869309057e-312,1.0e-322,9.0e-323,5.49596909987e-312,1.1e-322,9.4e-323,2.33419538276e-313,1.04e-322,9.4e-323,5.920368258063e-312), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(4.0e-323,1.0e-322,5.92036825806e-312,5.92036825817e-312,5.81426846863e-312,6.30232750059e-312,6.408427290197e-312,6.472087163946e-312,6.15378779531e-312,6.15378779523e-312,6.15378779523e-312), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(1.433e-321,6.87526636432e-312,7.34210543838e-312,7.554305017556e-312,6.13256783766e-312,7.702844722983e-312,8.190903754984e-312,8.275783586667e-312,8.38188337626e-312,7.766504596736e-312,7.74528463883e-312), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(1.803e-321,8.530423081637e-312,8.63652287123e-312,8.678962787064e-312,8.74262266081e-312,8.763842618683e-312,8.82750249246e-312,9.082141987436e-312,9.12458190327e-312,9.188241777014e-312,9.20946173489e-312), Dual{ForwardDiff.Tag{typeof(Main.ex-PSSolver.f), Float64}}(9.25190165076e-312,8.99726215585e-312,7.68162476546e-312,3.8195924455e-313,9.0e-323,9.442881269894e-312,9.63386089323e-312,9.464101229943e-312,2.1219960172e-314,5.0e-324,1.3e-322)]))
Problem name: Generic
All variables: ████████████████████ 10 All constraints: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
free: ████████████████████ 10 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) 55 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzj: (------% sparsity)
Counters:
obj: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
Running the partitioned quasi-Newton trust-region method
You minimize your model by calling the partitioned-update solver PTRUNK:
using PartiallySeparableSolvers
stats = PTRUNK(model; verbose=1)
print(stats)Start: Trust-region plse [ Info: iter f(x) π Δ ratio cgstatus [ Info: 0 3.6e+01 2.3e+01 1.0e+00 -1 initial point [ Info: 1 1.6e+01 1.6e+01 2.0e+00 8.7e-01 on trust-region boundary [ Info: 2 3.2e-01 1.5e+00 4.0e+00 1.0e+00 on trust-region boundary [ Info: 3 8.6e-02 5.8e-01 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 4 3.2e-02 3.4e-01 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 5 1.5e-02 2.0e-01 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 6 8.4e-03 1.6e-01 4.0e+00 1.5e+00 solution good enough given atol and rtol [ Info: 7 3.8e-03 1.2e-01 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 8 1.4e-03 9.0e-02 4.0e+00 1.4e+00 solution good enough given atol and rtol [ Info: 9 1.8e-04 3.0e-02 4.0e+00 1.2e+00 solution good enough given atol and rtol [ Info: 10 2.5e-05 1.3e-02 4.0e+00 1.1e+00 solution good enough given atol and rtol [ Info: 11 2.6e-06 3.1e-03 4.0e+00 1.2e+00 solution good enough given atol and rtol [ Info: 12 5.9e-08 7.2e-04 4.0e+00 1.0e+00 solution good enough given atol and rtol [ Info: 13 2.6e-09 1.7e-04 4.0e+00 1.2e+00 solution good enough given atol and rtol [ Info: 14 7.8e-11 1.6e-05 4.0e+00 1.1e+00 solution good enough given atol and rtol [ Info: 15 2.3e-13 1.4e-06 4.0e+00 1.0e+00 solution good enough given atol and rtol [ Info: 16 3.3e-15 1.9e-07 4.0e+00 1.1e+00 solution good enough given atol and rtol Generic Execution stats status: first-order stationary objective value: 3.3153079341003517e-15 primal feasibility: 0.0 dual feasibility: 1.94972452963349e-7 solution: [-2.4992321185261527e-8 2.531928630003402e-8 -5.041820581345534e-9 6.34669864372573e-9 ⋯ -2.5893626388969547e-8] iterations: 16 elapsed time: 2.0059099197387695
The partitioned quasi-Newton update performed by PTRUNK is chosen by the optional argument name. Allowed values include:
name=:pbfgseach $B_i$ is updated by BFGS;name=:psr1each $B_i$ is updated by SR1;name=:pseeach $B_i$ may be updated by both BFGS and SR1;name=:plbfgseach $B_i$ is an LBFGS operator;name=:plsr1each $B_i$ is an LSR1 operator;name=:plse(by default) each $B_i$ may be an LBFGS or LSR1 operator.
Thus for example,
stats = PTRUNK(model; name=:pbfgs, verbose=1)
print(stats)Start: Trust-region pbfgs [ Info: iter f(x) π Δ ratio cgstatus [ Info: 0 3.6e+01 2.3e+01 1.0e+00 -1 initial point [ Info: 1 1.6e+01 1.6e+01 2.0e+00 8.7e-01 on trust-region boundary [ Info: 2 3.2e-01 1.5e+00 4.0e+00 1.0e+00 on trust-region boundary [ Info: 3 2.9e-02 2.7e-01 4.0e+00 1.2e+00 solution good enough given atol and rtol [ Info: 4 9.6e-03 1.6e-01 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 5 3.6e-03 8.3e-02 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 6 7.6e-04 7.1e-02 4.0e+00 1.2e+00 solution good enough given atol and rtol [ Info: 7 2.6e-04 4.7e-02 4.0e+00 1.4e+00 solution good enough given atol and rtol [ Info: 8 5.3e-06 6.4e-03 4.0e+00 1.1e+00 solution good enough given atol and rtol [ Info: 9 2.7e-07 4.5e-04 4.0e+00 1.0e+00 solution good enough given atol and rtol [ Info: 10 4.1e-08 3.1e-04 4.0e+00 1.3e+00 solution good enough given atol and rtol [ Info: 11 3.2e-11 1.7e-05 4.0e+00 1.0e+00 solution good enough given atol and rtol [ Info: 12 5.9e-12 6.7e-06 4.0e+00 1.4e+00 solution good enough given atol and rtol [ Info: 13 1.4e-13 2.5e-07 4.0e+00 1.0e+00 solution good enough given atol and rtol Generic Execution stats status: first-order stationary objective value: 1.398142667201933e-13 primal feasibility: 0.0 dual feasibility: 2.515472254209989e-7 solution: [5.36024642406355e-7 -4.779964750225963e-7 3.692278582302615e-7 -2.2396869251400482e-7 ⋯ 5.360246403698534e-7] iterations: 13 elapsed time: 0.2174239158630371
updates each element Hessian approximation with BFGS.
See PartitionedStructures.jl tutorial for more details about partitioned quasi-Newton approximations.