Tutorial
Input
RipQP uses the package QuadraticModels.jl to model convex quadratic problems.
Here is a basic example:
using QuadraticModels, LinearAlgebra, SparseMatricesCOO
Q = [6. 2. 1.
2. 5. 2.
1. 2. 4.]
c = [-8.; -3; -3]
A = [1. 0. 1.
0. 2. 1.]
b = [0.; 3]
l = [-1.0;0;0]
u = [Inf; Inf; Inf]
QM = QuadraticModel(
c,
SparseMatrixCOO(tril(Q)),
A=SparseMatrixCOO(A),
lcon=b,
ucon=b,
lvar=l,
uvar=u,
c0=0.,
name="QM"
)QuadraticModels.QuadraticModel{Float64, Vector{Float64}, SparseMatricesCOO.SparseMatrixCOO{Float64, Int64}, SparseMatricesCOO.SparseMatrixCOO{Float64, Int64}}
Problem name: QM
All variables: ████████████████████ 3 All constraints: ████████████████████ 2
free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ████████████████████ 3 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ████████████████████ 2
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: ( 0.00% sparsity) 6 linear: ████████████████████ 2
nonlinear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzj: ( 33.33% sparsity) 4
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
Once your QuadraticModel is loaded, you can simply solve it RipQP:
using RipQP
stats = ripqp(QM)
println(stats)[ Info: Solving in Float64 using K2LDL with LDLFactorizationData
[ Info: iter obj rgap ‖rb‖ ‖rc‖ α_pri α_du μ ρ δ
[ Info: 0 3.6e+00 9.9e-01 1.4e+00 3.8e+00 0.0e+00 0.0e+00 5.9e+00 1.5e-03 1.5e-03
[ Info: 1 1.6e+00 2.7e-01 4.1e-03 1.4e-03 1.0e+00 1.0e+00 2.4e-01 1.5e-03 1.5e-03
[ Info: 2 1.2e+00 2.1e-02 5.5e-05 4.5e-05 1.0e+00 1.0e+00 1.5e-02 1.5e-04 1.5e-04
[ Info: 3 1.1e+00 1.0e-04 6.9e-07 1.9e-05 1.0e+00 1.0e+00 6.1e-05 1.5e-05 1.5e-05
[ Info: 4 1.1e+00 1.0e-07 1.0e-09 1.9e-08 1.0e+00 1.0e+00 6.1e-08 1.5e-06 1.5e-06
[ Info: 5 1.1e+00 1.0e-10 1.1e-12 1.8e-11 1.0e+00 1.0e+00 6.1e-11 1.5e-07 1.5e-07
Generic Execution stats
status: first-order stationary
objective value: 1.1250000001819647
primal feasibility: 5.253818735984856e-13
dual feasibility: 4.5240388878182354e-11
solution: [-6.602556960437742e-11 1.499999999967075 6.550018773077893e-11]
multipliers: [-5.000000000443094 2.249999999907261]
multipliers_L: [1.3498983056074131e-12 9.350004367856938e-13 2.7500000006802963]
multipliers_U: [0.0 0.0 0.0]
iterations: 5
elapsed time: 13.493134021759033
solver specific:
iters_sp: 5
pdd: 1.002449365261327e-10
psoperations: ∅
iters_sp3: 0
iters_sp2: 0
relative_iter_cnt: 20The stats output is a GenericExecutionStats.
It is also possible to use the package QPSReader.jl in order to read convex quadratic problems in MPS or SIF formats: (download QAFIRO)
using QPSReader, QuadraticModels
QM = QuadraticModel(readqps("assets/QAFIRO.SIF"))QuadraticModels.QuadraticModel{Float64, Vector{Float64}, SparseMatricesCOO.SparseMatrixCOO{Float64, Int64}, SparseMatricesCOO.SparseMatrixCOO{Float64, Int64}}
Problem name: Generic
All variables: ████████████████████ 32 All constraints: ████████████████████ 27
free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ████████████████████ 32 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ███████████████⋅⋅⋅⋅⋅ 19
low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ██████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 8
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: ( 98.86% sparsity) 6 linear: ████████████████████ 27
nonlinear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzj: ( 90.39% sparsity) 83
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
Logging
RipQP displays some logs at each iterate.
stats = ripqp(QM)
println(stats)[ Info: Solving in Float64 using K2LDL with LDLFactorizationData
[ Info: iter obj rgap ‖rb‖ ‖rc‖ α_pri α_du μ ρ δ
[ Info: 0 5.1e+03 5.2e+00 5.3e+02 6.2e+01 0.0e+00 0.0e+00 8.2e+02 1.5e-03 1.5e-03
[ Info: 1 2.8e+02 4.7e+01 7.7e+01 8.7e+00 8.6e-01 2.9e-01 3.4e+02 1.5e-03 1.5e-03
[ Info: 2 2.8e+02 3.7e+01 4.3e+01 7.3e+00 4.4e-01 2.8e-01 2.6e+02 1.5e-04 1.5e-04
[ Info: 3 2.2e+02 3.0e+00 2.5e+00 5.4e-01 9.4e-01 9.9e-01 1.5e+01 1.5e-05 1.5e-05
[ Info: 4 3.6e+01 2.1e+00 8.2e-03 6.4e-02 1.0e+00 9.9e-01 1.6e+00 1.5e-06 1.5e-06
[ Info: 5 4.7e+00 1.9e+00 3.5e-06 1.4e-06 1.0e+00 1.0e+00 2.2e-01 1.5e-07 1.5e-07
[ Info: 6 -4.3e-01 7.7e-01 1.1e-07 3.1e-01 1.0e+00 7.7e-01 3.5e-02 1.5e-08 1.5e-08
[ Info: 7 -1.5e+00 5.3e-02 1.1e-08 1.1e-02 9.1e-01 9.8e-01 3.1e-03 1.5e-09 1.5e-09
[ Info: 8 -1.6e+00 6.1e-05 1.0e-10 1.3e-05 1.0e+00 1.0e+00 3.7e-06 1.5e-10 1.5e-09
[ Info: 9 -1.6e+00 6.1e-08 2.6e-13 1.3e-08 1.0e+00 1.0e+00 3.7e-09 1.5e-11 1.5e-09
[ Info: 10 -1.6e+00 6.1e-11 7.1e-15 1.3e-11 1.0e+00 1.0e+00 3.7e-12 1.5e-12 1.5e-09
Generic Execution stats
status: first-order stationary
objective value: -1.590781793787462
primal feasibility: 7.105427357601002e-15
dual feasibility: 4.0220152292944574e-11
solution: [0.38028479685706695 8.806189977065276e-13 0.3802847968561864 0.40310188466849095 ⋯ 1.6759073619030825e-13]
multipliers: [-4.466932765449747 -2.2377603513257897e-12 -1.980617639140097e-14 -3.642178610385085 ⋯ -6.781571768245936e-15]
multipliers_L: [3.0441783478982717e-12 1.185323748792536 2.340230871153934e-12 2.242879722664324e-12 ⋯ 10.000000000000178]
multipliers_U: [0.0 0.0 0.0 0.0 ⋯ 0.0]
iterations: 10
elapsed time: 0.09997987747192383
solver specific:
iters_sp: 10
pdd: 6.100201586626379e-11
psoperations: [QuadraticModels.SingletonRow{Float64, Vector{Float64}}(3, 1, 1.0, false, true) QuadraticModels.SingletonRow{Float64, Vector{Float64}}(13, 16, 1.0, false, true)]
iters_sp3: 0
iters_sp2: 0
relative_iter_cnt: 40You can deactivate logging with
stats = ripqp(QM, display = false)
println(stats)Generic Execution stats
status: first-order stationary
objective value: -1.590781793787462
primal feasibility: 7.105427357601002e-15
dual feasibility: 4.0220152292944574e-11
solution: [0.38028479685706695 8.806189977065276e-13 0.3802847968561864 0.40310188466849095 ⋯ 1.6759073619030825e-13]
multipliers: [-4.466932765449747 -2.2377603513257897e-12 -1.980617639140097e-14 -3.642178610385085 ⋯ -6.781571768245936e-15]
multipliers_L: [3.0441783478982717e-12 1.185323748792536 2.340230871153934e-12 2.242879722664324e-12 ⋯ 10.000000000000178]
multipliers_U: [0.0 0.0 0.0 0.0 ⋯ 0.0]
iterations: 10
elapsed time: 0.0003490447998046875
solver specific:
iters_sp: 10
pdd: 6.100201586626379e-11
psoperations: [QuadraticModels.SingletonRow{Float64, Vector{Float64}}(3, 1, 1.0, false, true) QuadraticModels.SingletonRow{Float64, Vector{Float64}}(13, 16, 1.0, false, true)]
iters_sp3: 0
iters_sp2: 0
relative_iter_cnt: 40It is also possible to get a history of several quantities such as the primal and dual residuals and the relative primal-dual gap. These quantites are available in the dictionary solver_specific of the stats.
stats = ripqp(QM, display = false, history = true)
pddH = stats.solver_specific[:pddH]11-element Vector{Float64}:
5.202001221940632
47.09291657152063
37.04715154649735
2.9705628772130837
2.0982950800102707
1.949512880288423
0.7749630402562204
0.05312630164594008
6.100127510228456e-5
6.100081602192295e-8
6.100201586626379e-11Change configuration and tolerances
You can use RipQP without presolve and scaling with:
stats = ripqp(QM, ps=false, scaling = false)You can also change the RipQP.InputTol type to change the tolerances for the stopping criteria:
stats = ripqp(QM, itol = InputTol(max_iter = 5))
println(stats)[ Info: Solving in Float64 using K2LDL with LDLFactorizationData
[ Info: iter obj rgap ‖rb‖ ‖rc‖ α_pri α_du μ ρ δ
[ Info: 0 5.1e+03 5.2e+00 5.3e+02 6.2e+01 0.0e+00 0.0e+00 8.2e+02 1.5e-03 1.5e-03
[ Info: 1 2.8e+02 4.7e+01 7.7e+01 8.7e+00 8.6e-01 2.9e-01 3.4e+02 1.5e-03 1.5e-03
[ Info: 2 2.8e+02 3.7e+01 4.3e+01 7.3e+00 4.4e-01 2.8e-01 2.6e+02 1.5e-04 1.5e-04
[ Info: 3 2.2e+02 3.0e+00 2.5e+00 5.4e-01 9.4e-01 9.9e-01 1.5e+01 1.5e-05 1.5e-05
[ Info: 4 3.6e+01 2.1e+00 8.2e-03 6.4e-02 1.0e+00 9.9e-01 1.6e+00 1.5e-06 1.5e-06
[ Info: 5 4.7e+00 1.9e+00 3.5e-06 1.4e-06 1.0e+00 1.0e+00 2.2e-01 1.5e-07 1.5e-07
Generic Execution stats
status: maximum iteration
objective value: 4.658954498519296
primal feasibility: 3.4981491245367202e-6
dual feasibility: 1.6824715908247967e-6
solution: [0.743878138254691 0.004674605650346285 0.7392033550578112 0.78851071023687 ⋯ 0.002243926732604866]
multipliers: [-7.826615009307647 -0.44355741689044687 -0.0002826585049596877 -7.445204284037653 ⋯ -9.745279751756744e-5]
multipliers_L: [0.1452930232869297 1.5112382757139278 0.6014690171102718 0.44363007256944875 ⋯ 10.002732536283734]
multipliers_U: [0.0 0.0 0.0 0.0 ⋯ 0.0]
iterations: 5
elapsed time: 0.007070064544677734
solver specific:
iters_sp: 5
pdd: 1.949512880288423
psoperations: [QuadraticModels.SingletonRow{Float64, Vector{Float64}}(3, 1, 1.0, false, true) QuadraticModels.SingletonRow{Float64, Vector{Float64}}(13, 16, 1.0, false, true)]
iters_sp3: 0
iters_sp2: 0
relative_iter_cnt: 20Save the Interior-Point system
At every iteration, RipQP solves two linear systems with the default Predictor-Corrector method (the affine system and the corrector-centering system), or one linear system with the Infeasible Path-Following method.
To save these systems, you can use:
w = SystemWrite(write = true, name="test_", kfirst = 4, kgap=3)
stats = ripqp(QM, w = w)This will save one matrix and the associated two right hand sides of the PC method every three iterations starting at iteration four. Then, you can read the saved files with:
using DelimitedFiles, MatrixMarket
K = MatrixMarket.mmread("test_K_iter4.mtx")
rhs_aff = readdlm("test_rhs_iter4_aff.rhs", Float64)[:]
rhs_cc = readdlm("test_rhs_iter4_cc.rhs", Float64)[:];Timers
You can see the elapsed time with:
stats.elapsed_timeFor more advance timers you can use TimerOutputs.jl:
using TimerOutputs
TimerOutputs.enable_debug_timings(RipQP)
reset_timer!(RipQP.to)
stats = ripqp(QM)
TimerOutputs.complement!(RipQP.to) # print complement of timed sections
show(RipQP.to, sortby = :firstexec)