Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix start value #138

Merged
merged 5 commits into from
Mar 4, 2019
Merged

Fix start value #138

merged 5 commits into from
Mar 4, 2019

Conversation

blegat
Copy link
Member

@blegat blegat commented Feb 19, 2019

Closes #136

@kalmarek
Copy link
Collaborator

maybe I'm using it wrong but here is what I tried:

function example_problem(m = Model())
    @variable(m, x >= 0)
    @variable(m, y <= 2)
    @variable(m, P[1:2,1:2] >= 0)
    @objective(m, Max, x+y + P[1,1] - P[2,2])
    @constraint(m, x <= y+4)
    @constraint(m, P[1,1] + P[2,2] <= 4)
    @constraint(m, P[2,1] - P[1,1] <= 4)
    
    return m
end

m = example_problem()
optimize!(m, with_optimizer(SCS.Optimizer, eps=1e-10))

solvedx = value(m[:x])
solvedy = value(m[:y])
solvedP = value.(m[:P])

m2 = example_problem()
set_start_value(m2[:x], solvedx)
set_start_value(m2[:y], solvedy)
set_start_value.(m2[:P], solvedP)

optimize!(m2, with_optimizer(SCS.Optimizer, eps=1e-10, max_iters=3, warm_start=true))
@assert termination_status(m2) == termination_status(m) # fails

I can however overwrite SCS.MOISolution:

primal = backend(m).optimizer.model.optimizer.sol.primal
dual = backend(m).optimizer.model.optimizer.sol.dual
slack = backend(m).optimizer.model.optimizer.sol.slack

m3 = example_problem(Model(with_optimizer(SCS.Optimizer, eps=1e-10, max_iters=3, warm_start=true)))

m3.moi_backend.optimizer.model.optimizer.sol.primal = primal
m3.moi_backend.optimizer.model.optimizer.sol.dual = dual
m3.moi_backend.optimizer.model.optimizer.sol.slack = slack

optimize!(m3)
@assert termination_status(m3) == termination_status(m) # this works

@blegat
Copy link
Member Author

blegat commented Feb 19, 2019

What is the termination status for m2 ?

@kalmarek
Copy link
Collaborator

m2 = example_problem()
set_start_value(m2[:x], solvedx)
set_start_value(m2[:y], solvedy)
set_start_value.(m2[:P], solvedP)

with_SCS_3iters = with_optimizer(SCS.Optimizer, eps=1e-10, max_iters=3, warm_start=true)
optimize!(m2, with_SCS_3iters)
termination_status(m2)
----------------------------------------------------------------------------
	SCS v2.0.2 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012-2017
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 12, CG tol ~ 1/iter^(2.00)
eps = 1.00e-10, alpha = 1.50, max_iters = 3, normalize = 1, scale = 1.00
acceleration_lookback = 20, rho_x = 1.00e-03
Variables n = 6, constraints m = 9
Cones:	linear vars: 9
Setup time: 3.03e-05s
SCS using variable warm-starting
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 8.03e+00  1.75e+00  7.88e-01 -7.37e+01 -8.29e+00  0.00e+00  1.89e-05 
     3| 7.00e-01  5.44e-01  9.19e-01  2.76e+00 -8.51e+00  8.83e-16  3.28e-05 
----------------------------------------------------------------------------
Status: Solved/Inaccurate
Hit max_iters, solution may be inaccurate
Timing: Solve time: 3.39e-05s
	Lin-sys: avg # CG iterations: 2.50, avg solve time: 9.99e-07s
	Cones: avg projection time: 1.84e-07s
	Acceleration: avg step time: 1.94e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 7.6374e-01, dist(y, K*) = 0.0000e+00, s'y/|s||y| = -1.3466e-01
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 6.9958e-01
dual res:   |A'y + c|_2 / (1 + |c|_2) = 5.4372e-01
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 9.1850e-01
----------------------------------------------------------------------------
c'x = 2.7591, -b'y = -8.5115
============================================================================

ALMOST_OPTIMAL::TerminationStatusCode = 7

compare with the same problem where start values are not set:

m2 = example_problem()
# set_start_value(m2[:x], solvedx)
# set_start_value(m2[:y], solvedy)
# set_start_value.(m2[:P], solvedP)

with_SCS_3iters = with_optimizer(SCS.Optimizer, eps=1e-10, max_iters=3, warm_start=true)
optimize!(m2, with_SCS_3iters)
termination_status(m2)
----------------------------------------------------------------------------
	SCS v2.0.2 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012-2017
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 12, CG tol ~ 1/iter^(2.00)
eps = 1.00e-10, alpha = 1.50, max_iters = 3, normalize = 1, scale = 1.00
acceleration_lookback = 20, rho_x = 1.00e-03
Variables n = 6, constraints m = 9
Cones:	linear vars: 9
Setup time: 3.53e-05s
SCS using variable warm-starting
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 6.15e+00  1.28e+00  7.53e-01 -5.71e+01 -7.61e+00  0.00e+00  2.28e-05 
     3| 4.22e-01  2.59e-01  2.60e-01 -6.17e+00 -1.09e+01  4.17e-16  4.02e-05 
----------------------------------------------------------------------------
Status: Solved/Inaccurate
Hit max_iters, solution may be inaccurate
Timing: Solve time: 4.18e-05s
	Lin-sys: avg # CG iterations: 2.25, avg solve time: 1.28e-06s
	Cones: avg projection time: 1.44e-07s
	Acceleration: avg step time: 2.35e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 5.9819e-16, dist(y, K*) = 0.0000e+00, s'y/|s||y| = 7.1798e-02
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 4.2229e-01
dual res:   |A'y + c|_2 / (1 + |c|_2) = 2.5927e-01
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 2.6037e-01
----------------------------------------------------------------------------
c'x = -6.1744, -b'y = -10.8735
============================================================================

ALMOST_OPTIMAL::TerminationStatusCode = 7

@blegat
Copy link
Member Author

blegat commented Feb 19, 2019

Isn't the difference between m2 and m3 the fact that in m3 you also set the start for slack and dual while in m2 you only set the start for primal ?

@kalmarek
Copy link
Collaborator

SCS will ignore primal when dual and slack are not set:
https://github.com/JuliaOpt/SCS.jl/blob/master/src/c_wrapper.jl#L65
and this is what I believe happening, based on Your comment.

So maybe the question should be: is it possible to set dual and slack from the level of MOI?

@kalmarek
Copy link
Collaborator

sorry for dumb questions, but the documentation is rather sparse at the moment:

I tried copying solution from m to m2:

with_SCS = with_optimizer(SCS.Optimizer,
    eps=1e-10,
#     max_iters=3,
    warm_start=true,
    verbose=false)

with_SCS_3iters = with_optimizer(SCS.Optimizer,
    eps=1e-10,
    max_iters=3,
    warm_start=true,
    verbose=true)

m = example_problem(Model(with_SCS))

optimize!(m)

m2 = example_problem(Model(with_SCS_3iters));
MOI.set(m2, MOI.VariablePrimalStart(), m2[:x], value(m[:x]))
MOI.set(m2, MOI.VariablePrimalStart(), m2[:y], value(m[:y]))
MOI.set.(m2, MOI.VariablePrimalStart(), m2[:P], value.(m[:P]))

# is this how you iterate over constraints?
for (F,S) in MOI.get(m.moi_backend, MOI.ListOfConstraints())
    cnstrs = MOI.get(m.moi_backend, MOI.ListOfConstraintIndices{F,S}())
    for c in cnstrs
        dual_val = MOI.get(m.moi_backend, MOI.ConstraintDual(), c)
        slack_val = MOI.get(m.moi_backend, MOI.ConstraintPrimal(), c)
        MOI.set(m2.moi_backend, MOI.ConstraintDualStart(), c, dual_val)
        MOI.set(m2.moi_backend, MOI.ConstraintPrimalStart(), c, slack_val)
    end
end

optimize!(m2)

however I still get warnings:

┌ Warning: MathOptInterface.ConstraintPrimalStart() is not supported by MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{SCS.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}. This 
│   information will be discarded. = information will be discarded.
└ @ MathOptInterface.Utilities /home/kalmar/.julia/packages/MathOptInterface/QVlW6/src/Utilities/copy.jl:133
┌ Warning: MathOptInterface.ConstraintDualStart() is not supported by MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{SCS.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}. This 
│   information will be discarded. = information will be discarded.
└ @ MathOptInterface.Utilities /home/kalmar/.julia/packages/MathOptInterface/QVlW6/src/Utilities/copy.jl:133
...

@blegat
Copy link
Member Author

blegat commented Feb 20, 2019

Only VariablePrimalStart is implemented for the MOI SCS wrapper.

@kalmarek
Copy link
Collaborator

as noted above implementing only VariablePrimalStart is pointless for SCS; what needs to be done for ConstraintDualStart and ConstraintPrimalStart?

(I see that load functions are already defined, with typos though)

@blegat
Copy link
Member Author

blegat commented Feb 20, 2019

(I see that load functions are already defined, with typos though)

Indeed, it seems we only need to add supports methods
https://github.com/JuliaOpt/SCS.jl/blob/980ef0f1032eb45be43e8a3d7a701f87257f9434/src/MOI_wrapper.jl#L82

@kalmarek
Copy link
Collaborator

You're missing comma after Union;

unfortunately I've already tried this and I get the same warnings (and results). note:

F = MOI.SingleVariable
S = MOI.EqualTo{Float64}

MOI.supports(Model(with_SCS).moi_backend, MOI.ConstraintDualStart(), MOI.ConstraintIndex{F,S})

still returns false;

the defined supports functions are for
MOI.supports(Model(with_SCS).moi_backend.optimizer.model.optimizer, MOI.ConstraintDualStart(), MOI.ConstraintIndex{F,S})) and the chain gets broken somewhere

@blegat
Copy link
Member Author

blegat commented Feb 21, 2019

Good point, we need to implement starting values for bridges

@blegat
Copy link
Member Author

blegat commented Feb 21, 2019

Let's first check if it works if all the constraints are natively supported, i.e.
https://github.com/JuliaOpt/SCS.jl/blob/980ef0f1032eb45be43e8a3d7a701f87257f9434/src/MOI_wrapper.jl#L8-L9
then we will add support for bridges in MOI.

kalmarek added a commit to kalmarek/PropertyT.jl that referenced this pull request Feb 21, 2019
@kalmarek
Copy link
Collaborator

@blegat what do you mean by natively supported?

reading jump-dev/MathOptInterface.jl#528 I understand that we need to wait to remove VectorOfVariables support (as this will be transformed by Bridge) is that correct?

@blegat
Copy link
Member Author

blegat commented Feb 22, 2019

Indeed, I was ambiguous, I meant natively supported by the MOI wrapper (i.e. no bridge are needed). So it should work if for a model like

model = Model(with_optimizer(SCS.Optimizer))
@variable(model, x[1:2])
@constraint(model, [x[1] - x[2]] in MOI.Nonnegatives(1))
@constraint(model, [1; x] in MOI.SecondOrderCone(2))
@objective(model, Max, sum(x))

@blegat blegat merged commit 71482c5 into master Mar 4, 2019
@blegat
Copy link
Member Author

blegat commented Mar 4, 2019

I have opened an issue for start values in the bridges: jump-dev/MathOptInterface.jl#684

@odow odow deleted the bl/fixstart branch October 8, 2020 01:09
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

Successfully merging this pull request may close these issues.

2 participants