Differentiating through a QP - order of constraints

Specific Domains Optimization (Mathematical)
,
1

This is a follow-up to this topic.

I have a QP that I want to differentiate through using DiffOpt. In contrast to the example from the tutorial, I have two constraints:

n = 2 # variable dimension
m = 2; # no of inequality constraints

Q = [4.0 1.0; 1.0 2.0]
q = [1.0; 1.0]
G = [1.0 1.0]
h = [-1.0]   # initial values set
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)
@variable(model, x[1:2])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);

@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)

result = value.(x)

The constraints are the same just for the sake of computations. I want to differentiate both constraints wrt to the right-hand side, h.

Question: how do I do that?

What does not work:

  1. Doing everything iteratively - differentiation through cons2 is wrong. I understand why, we modify the model for cons1, so model changes and we get a different result for cons2.
Code 
MOI.set(
    model,
    DiffOpt.ForwardConstraintFunction(),
    cons1[1],
    0.0 * index(x[1]) - 1.0,  # to indicate the direction vector to get directional derivatives
)

DiffOpt.forward_differentiate!(model)
dxh1 = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), x)

MOI.set(
    model,
    DiffOpt.ForwardConstraintFunction(),
    cons2[1],
    0.0 * index(x[1]) - 1.0,  # to indicate the direction vector to get directional derivatives
)

DiffOpt.forward_differentiate!(model)
dxh2 = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), x)

@show dxh1, dxh2
  1. Doing differentiation through a function. The values for cons2 are incorrect. If I change the order and use first cons2 then cons1, then the results for cons1 are incorrect. This is what puzzles me, I thought that whatever was done inside the function was not available outside.
Code 
function directionbFcn(model,constr)
    dwdbOut = []
    for cs in constr
        MOI.set(
            model,
            DiffOpt.ForwardConstraintFunction(),
            cs,
            0.0 * index(model[:x][1]) - 1.0,  # to indicate the direction vector to get directional derivatives
        )
        DiffOpt.forward_differentiate!(model)
        dwdb = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), model[:x])
        push!(dwdbOut, dwdb)
    end
    return dwdbOut
end
dxh1Fcn = directionbFcn(model,cons1)
dxh2Fcn = directionbFcn(model,cons2) ###Wrong?
@show dwh, dxh1Fcn, dxh2Fcn

What does work:

  1. Redoing the model for every constraint. I see why but I want to avoid solving the problem multiple times.
Code 
#####Model 1
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)
@variable(model, x[1:2])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);

@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)

dxh1Fcn = directionbFcn(model,cons1)

#####Model 2
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)
@variable(model, x[1:2])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);

@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)

dxh2Fcn = directionbFcn(model,cons2)

The entire MWE, together with a comparison with a different computation method:

Code 
using JuMP
import DiffOpt
import Ipopt
##########
function directionbFcn(model,constr)
    dwdbOut = []
    for cs in constr
        MOI.set(
            model,
            DiffOpt.ForwardConstraintFunction(),
            cs,
            0.0 * index(model[:x][1]) - 1.0,  # to indicate the direction vector to get directional derivatives
        )
        DiffOpt.forward_differentiate!(model)
        dwdb = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), model[:x])
        push!(dwdbOut, dwdb)
    end
    return dwdbOut
end
###############

n = 2 # variable dimension
m = 2; # no of inequality constraints

Q = [4.0 1.0; 1.0 2.0]
q = [1.0; 1.0]
G = [1.0 1.0]
h = [-1.0]   # initial values set
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)
@variable(model, x[1:2])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);

@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)

result = value.(x)

######This uses duals from JuMP which can be negative. 
######However, Amos and Kolter (2017) assume λ≥0, so
######we need to adjust the signs
######Duals from JuMP are non-positive for ≤ constraints
q = q
G = G
h = h
Q = Q

LHM = [
    Q -transpose(vcat(G,G))
    -Diagonal(vcat(dual.(cons1), dual.(cons2)))*vcat(G,G) -Diagonal(vcat(G*result-h,G*result-h)) 
]
LHMinv = inv(LHM)
k1 = length(LHM[:,1])
k2 = length(Diagonal(vcat(dual.(cons1), dual.(cons2)))[1,:])
RHMb = [zeros(k1-k2,k2)
-Diagonal(vcat(dual.(cons1), dual.(cons2)))] #Matrix(1.0I, size(h)[1])  
RHMc = [-Matrix(1.0I, 2,2)
zeros(k2,2)]

RHMc = [-Matrix(1.0I, 2,2)
zeros(k2,2)]

dallb = LHMinv*RHMb    ####both w and λ

dwh = dallb[1:2,:] ####picking only w

########Trying out DiffOpt
########Doesn't work, but that's understandable
# MOI.set(
#     model,
#     DiffOpt.ForwardConstraintFunction(),
#     cons1[1],
#     0.0 * index(x[1]) - 1.0,  # to indicate the direction vector to get directional derivatives
# )

# DiffOpt.forward_differentiate!(model)
# dxh1 = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), x)

# MOI.set(
#     model,
#     DiffOpt.ForwardConstraintFunction(),
#     cons2[1],
#     0.0 * index(x[1]) - 1.0,  # to indicate the direction vector to get directional derivatives
# )

# DiffOpt.forward_differentiate!(model)
# dxh2 = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), x)

# @show dwh, dxh1, dxh2
########These return wrong results - why?
dxh1Fcn = directionbFcn(model,cons1)
dxh2Fcn = directionbFcn(model,cons2)
@show dwh, dxh1Fcn, dxh2Fcn

#####REDOING THE MODEL FOR EVERY CONSTRAINT

#####Model 1
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)
@variable(model, x[1:2])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);

@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)

dxh1Fcn = directionbFcn(model,cons1)

#####Model 2
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
set_silent(model)
@variable(model, x[1:2])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h[j]);

@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)

dxh2Fcn = directionbFcn(model,cons2)
@show dwh,dxh1,dxh2, dxh1Fcn, dxh2Fcn
2

I think this is a question for @mbesancon @joaquimg or @blegat. I’ll happily admit that I don’t know much about DiffOpt.

If I have to guess though, you can’t differentiate through h in two places at the same time, because JuMP doesn’t know that they are the same h. It just sees the value -1.0.

Perhaps you could do

using JuMP, Ipopt
n = 2 # variable dimension
m = 2; # no of inequality constraints
Q = [4.0 1.0; 1.0 2.0]
q = [1.0; 1.0]
G = [1.0 1.0]
h = [-1.0]
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, x[1:2])
@variable(model, h_var[j in 1:1])
@constraint(model, consh[j in 1:1], h_var[j] == h[j])
@constraint(model, cons1[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h_var[j]);
@constraint(model, cons2[j in 1:1], sum(G[j, i] * x[i] for i in 1:2) <= h_var[j]);
@objective(
    model,
    Min,
    1 / 2 * sum(Q[j, i] * x[i] * x[j] for i in 1:2, j in 1:2) +
    sum(q[i] * x[i] for i in 1:2)
)
optimize!(model)
result = value.(x)
h_result = reduced_cost.(h_var)

And then differentiate with respect to the right-hand side of consh?

3

Thanks for the reply. I am not sure this is the problem, I don’t want to differentiate simultaneously. What I think is happening is when I do this:

    for cs in constr
        MOI.set(
            model,
            DiffOpt.ForwardConstraintFunction(),
            cs,
            0.0 * index(model[:x][1]) - 1.0,  # to indicate the direction vector to get directional derivatives
        )
        DiffOpt.forward_differentiate!(model)
        dwdb = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), model[:x])
        push!(dwdbOut, dwdb)
    end

in every iteration DiffOpt wants to differentiate through the current cs and that works in the first iteration. In the second iteration I would like to differentiate through the second constraints in a given direction, but DiffOpt keeps the direction for the first constraint that I set in the first iteration. This would explain why changing the order of constraints in the differentiation changes the results.

My current workaround (inspired by this example) is

MOI.set.(
    model,
    DiffOpt.ForwardConstraintFunction(),
    [cons1;cons2],
    [0.0 * index(x[1]) - 1.0;0.0 * index(x[1]) ],  # to indicate the direction vector to get directional derivatives
)

DiffOpt.forward_differentiate!(model)
dxh1 = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), x)

MOI.set.(
    model,
    DiffOpt.ForwardConstraintFunction(),
    [cons1;cons2],
    [0.0 * index(x[1]) ;0.0 * index(x[1]) - 1.0 ],  # to indicate the direction vector to get directional derivatives
)

DiffOpt.forward_differentiate!(model)
dxh3 = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), x)