# Differentiating through a QP - order of constraints
**URL:** https://discourse.julialang.org/t/differentiating-through-a-qp-order-of-constraints/119275
**Category:** Optimization (Mathematical)
**Tags:** jump, diffopt
**Created:** [September 11, 2024, 9:48am UTC](https://discourse.julialang.org/t/differentiating-through-a-qp-order-of-constraints/119275 "2024-09-11T09:48:53Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@blob](https://discourse.julialang.org/u/blob)
#### Post date: [September 11, 2024, 9:48am UTC](https://discourse.julialang.org/t/differentiating-through-a-qp-order-of-constraints/119275/1 "2024-09-11T09:48:53Z")
This is a follow-up to [this topic](https://discourse.julialang.org/t/diffopt-differentiating-through-a-qp-different-from-manual/118141).
I have a QP that I want to differentiate through using DiffOpt. In contrast to the [example from the tutorial](https://jump.dev/DiffOpt.jl/dev/examples/matrix-inversion-manual/), I have two constraints:
```julia
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**
>
> ```julia
> 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**
>
> ```julia
> 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**
>
> ```julia
> #####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**
>
> ```julia
> 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
>
> ```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [September 12, 2024, 9:42pm UTC](https://discourse.julialang.org/t/differentiating-through-a-qp-order-of-constraints/119275/2 "2024-09-12T21:42:47Z")
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
```Julia
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`?
---
### Author:  [@blob](https://discourse.julialang.org/u/blob)
#### Post date: [September 13, 2024, 7:03am UTC](https://discourse.julialang.org/t/differentiating-through-a-qp-order-of-constraints/119275/3 "2024-09-13T07:03:46Z")
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:
```julia
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](https://jump.dev/DiffOpt.jl/dev/examples/sensitivity-analysis-svm/#Gradient-of-hyperplane-wrt-the-data-point-coordinates)) is
```julia
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)
```
---
### Author:  [@blegat](https://discourse.julialang.org/u/blegat)
#### Post date: [September 20, 2024, 6:04am UTC](https://discourse.julialang.org/t/differentiating-through-a-qp-order-of-constraints/119275/4 "2024-09-20T06:04:17Z")
Thanks, that’s indeed confusing, I have opened an issue: [Resetting differentiation input in-between differentiations · Issue #259 · jump-dev/DiffOpt.jl · GitHub](https://github.com/jump-dev/DiffOpt.jl/issues/259)
I’m tempted to say that it’s expected behavior that the previous gradient that was set stays but given that it’s not so easy to unset, maybe we should reset them in-between calls to `DiffOpt.forward_differentiate!(model)`