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:
- Doing everything iteratively - differentiation through cons2 is wrong. I understand why, we modify the model for
cons1
, somodel
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
- Doing differentiation through a function. The values for
cons2
are incorrect. If I change the order and use firstcons2
thencons1
, then the results forcons1
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:
- 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