Warm start using Dualization with JuMP and COSMO

Specific Domains Optimization (Mathematical)
1

I’m using the same example here as in my previous question about warm starting of HermitianPSDCone:

using JuMP, LinearAlgebra, COSMO, Dualization

function set_hermitian_start(x, start)
    for j in 1:size(x, 2), i in 1:size(x, 1)
        set_start_value(real(x[i, j]), real(start[i, j]))
        if i < j
            set_start_value(imag(x[i, j]), imag(start[i, j]))
        end
    end
    return
end

solver = optimizer_with_attributes(COSMO.Optimizer, "check_termination" => 1)
model = Model(dual_optimizer(solver))

ρs = [ [[1 0]; [0 0]],
       [[1 1]; [1 1]]/2 ]

@variable(model, E1[1:2, 1:2] in HermitianPSDCone())
@variable(model, E2[1:2, 1:2] in HermitianPSDCone())
@constraint(model, E1+E2 == LinearAlgebra.I)
@objective(model, Max, real(dot(ρs, [E1,E2])) / 2)

E1s = [[0.853553 -0.353553]; [-0.353553  0.146447]]
E2s = [[0.146447 0.353553]; [0.353553 0.853553]]
set_hermitian_start(E1, E1s)
set_hermitian_start(E2, E2s)

optimize!(model)

which COSMO then reports

------------------------------------------------------------------
          COSMO v0.8.8 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{16},
          constraints: A ∈ R^{20x16} (28 nnz),
          matrix size to factor: 36x36,
          Floating-point precision: Float64
Sets:     DensePsdConeTriangle of dim: 10 (4x4)
          DensePsdConeTriangle of dim: 10 (4x4)
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 1 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,	
          Safeguarded: true, tol: 2.0
Setup Time: 1.95ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-2.5007e+01	6.5518e+00	5.9999e-01	1.0000e-01
2	 3.1288e+01	7.7821e+00	3.5996e-01	1.0000e-01
3	-2.7261e+01	7.0431e+00	2.1595e-01	1.0000e-01
4	 2.2729e+01	5.6282e+00	1.2954e-01	1.0000e-01
5	 5.0000e-01	8.8388e-02	2.0329e-07	1.0000e-01
6	 5.0000e-01	8.8388e-02	1.0083e-12	1.0000e-01
7	 5.0000e-01	8.8388e-02	6.0496e-13	1.0000e-01
8	 5.0000e-01	8.8388e-02	3.6293e-13	1.0000e-01
9	 5.0000e-01	8.8388e-02	2.1794e-13	1.0000e-01
10	 5.0000e-01	8.8388e-02	1.3090e-13	1.0000e-01
11	 5.0000e-01	8.8388e-02	7.8382e-14	1.0000e-01
12	 5.0000e-01	8.8388e-02	4.7518e-14	1.0000e-01
13	 5.0000e-01	8.8388e-02	2.8422e-14	1.0000e-01
14	 5.0000e-01	8.8388e-02	1.7097e-14	1.0000e-01
15	 5.0000e-01	8.8388e-02	1.0214e-14	1.0000e-01
16	 5.0000e-01	8.8388e-02	6.4393e-15	1.0000e-01
17	 5.0000e-01	8.8388e-02	3.7748e-15	1.0000e-01
18	 5.0000e-01	8.8388e-02	2.5535e-15	1.0000e-01
19	 5.0000e-01	8.8388e-02	1.4433e-15	1.0000e-01
20	 5.0000e-01	8.8388e-02	8.8818e-16	1.0000e-01
21	 5.0000e-01	8.8388e-02	5.5511e-16	1.0000e-01
22	 5.0000e-01	8.8388e-02	4.1633e-16	1.0000e-01
23	 5.0000e-01	8.8388e-02	3.3307e-16	1.0000e-01
24	 5.0000e-01	8.8388e-02	2.2204e-16	1.0000e-01
25	 5.0000e-01	8.8388e-02	1.1102e-16	1.0000e-01
26	 5.0000e-01	8.8388e-02	4.4409e-16	1.0000e-01
27	 5.0000e-01	8.8388e-02	3.3307e-16	1.0000e-01
28	 5.0000e-01	8.8388e-02	1.1102e-16	1.0000e-01
29	 5.0000e-01	8.8388e-02	1.1102e-16	1.0000e-01
30	 5.0000e-01	8.8388e-02	1.4004e-16	1.0000e-01
31	 5.0000e-01	8.2924e-02	1.7293e-02	1.0000e-01
32	 1.3926e+00	1.1951e-01	7.1599e-03	1.0000e-01
33	 1.0481e+00	3.8818e-02	8.8949e-03	1.0000e-01
34	 7.2174e-01	2.6514e-02	8.4270e-03	1.0000e-01
35	 8.5355e-01	1.0243e-06	2.6155e-08	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 35
Optimal objective: 0.8536
Runtime: 0.014s (13.79ms)

Note I’m setting the start values to their optimal values.

I suspect that the start values are not being passed to COSMO via dualization since running without warm starting gives the same number iterations and COSMO reports the same objective value on the first iteration. If I instead just solve the primal solution then COSMO will report a different starting objective value and take a different number of iterations.

Curiously COSMO seems to take 14 iterations with warm starting as opposed to 12 without and I’m not sure why. Although COSMO reports the same dual residual with and without warm starting so perhaps it has to do with that?

1 Like
2

Correct. I don’t think Dualization.jl supports warm starts.

Edit: I opened an issue: Support starting values · Issue #165 · jump-dev/Dualization.jl · GitHub

1 Like
3

Thanks for reporting this.
The starting values are implemented in Add support for starting values by blegat · Pull Request #166 · jump-dev/Dualization.jl · GitHub
You can try it out with:

using Pkg
pkg"add Dualization#bl/start"

restart Julia for the changes to take effect.

1 Like
4

Thanks @blegat for implementing this so quickly!

I’m trying this example with your PR and it does not appear that COSMO is taking a different number of iterations or starting with a different objective value between being warm started or not.

Is there some other way I can check if the solver is actually using the warm start values or not?

5

@blegat’s PR is still a work in progress. It probably isn’t ready to try out just yet.

1 Like
6

Looking at the code of COSMO, it looks like it initializes the solution at zero and then replace these zero values by something once you set them.
At the end of the optimization, you see the value it converged to that has overwritten the starting value so you need to check this before optimize!.
You need to attach_optimizer the first CachingOptimizer and the second one deep down so that the model is copied to COSMO otherwise you will just see an empty model:

julia> using JuMP, COSMO

julia> model = Model(COSMO.Optimizer);

julia> @variable(model, x, start = 2);

julia> con_ref = @constraint(model, x == 1);

julia> set_dual_start_value(con_ref, 3.0)

julia> MOI.Utilities.attach_optimizer(model)

julia> MOI.Utilities.attach_optimizer(backend(model).optimizer.model)

julia> unsafe_backend(model).inner.vars
COSMO.Variables{Float64}([0.0, 0.0], [2.0, 0.0], [2.0], [0.0], [-3.0])

To avoid the attach_optimizer hack, you can also optimize with zero iterations but then the constraint dual might be reset:

julia> set_silent(model)

julia> set_attribute(model, "max_iter", 0)

julia> optimize!(model)

julia> unsafe_backend(model).inner.vars
COSMO.Variables{Float64}([0.4480000155199999, 0.01800001551999985], [2.0, 0.0], [2.0], [0.0], [0.0])