Drake: Extract dual variables from MathematicalProgramResult

[RussTedrake]

I feel like this is a major missing feature in our MathematicalProgram offering. I would imagine the interface would be something like:

prog = MathematicalProgram()
x = prog.NewContinuousVariables(...)
c = prog.AddContraint(...)
...
result = Solve(prog)
xsol = result.GetSolution(x)
ysol = result.GetDualSolution(c)

(that is similar to what we have in CVX, Yalmip, and JuMP.

@hongkai-dai -- have you thought about this feature yet?

[hongkai-dai]

May I ask what is the use case for the dual variables? Each dual variable is associated with a constraint, should we return the dual variable with a given constraint? Namely the API could be

Eigen::VectorXd GetDualVariableForConstraint(const Binding<Constraint>& constraint);

Another question is that we need to reformulate some constraint for each solver. For example, in Mosek, the rotated Lorentz cone constraint is

2 * y1 * y2 >= (z1^2 + ... + zn^2)

while in Gurobi, the rotated Lorentz cone constraint doesn't have the factor 2 on the left hand-side

y1 * y2 >= (z1^2 + ... + zn^2)

So the meaning of the dual variable is different for the same RotatedLorentzConeConstraint in Mosek versus Gurobi.

The simplest solution is just to return the vector containing all the dual variables, but it would be hard for the user to parse this dual variable, as they don't know which constraint is associated with each dual variable.

[RussTedrake]

There are many use cases. Today I was trying to code up an example of using the dual solution of the LP relaxation as a heuristic for solving MIP. Can tell you more details soon.

I think that your proposed API is almost the same as what I proposed, with a different name? And you agree that we want the method inside MathematicalProgramResult, right?

I see your point about the subtlety of those more complicated constraints. I think that the user already has to deal with that subtlety if they are going to do anything with the constraint bindings that are returned from those methods? I don’t feel like adding the ability to get the dual solution for any one of those returned constraints actually adds any new complications to the API? (Perhaps it reveals an existing complexity)?

[RussTedrake]

Ah.. maybe our API’s were different. You are proposing to return the dual “variables”? And not the values? But as a VectorXd? I would have thought you just return the values as VectorXd?

[hongkai-dai]

Sorry I misread the comment, we proposed the same API, to return the value of the dual variables associated with a constraint.

This feature might take some time to implement. For each of the solver (Gurobi / Mosek / Snopt etc), it has its own ordering of the dual variables. We need to store the mapping between each Binding<Constraint> object with the indices of the dual variables inside each solver class GurobiSolver, MosekSolver, SnoptSolver, etc. I will add GetDualSolution for each solver. Which solver is more urgent for you?

[hongkai-dai]

BTW, as a quick work-around, if you know which solver is used, the result.get_solver_details() could contain the dual variable values (like for SNOPT). It packs all dual variable values for all constraints into a single vector, so it is not very interpretable. But Benoit and I have been using this dual solution since we don't care the associated constraint.

[RussTedrake]

I know that it's a lot of work, and I don't mean to ask you to stop and implement everything. But perhaps it's worth adding the plumbing to enable this feature, and implementing it for just a few cases (e.g. linear constraints in snopt or gurobi...) and throw "not implemented yet" in the others. That would let others contribute the missing solver bindings as needed?

[hongkai-dai]

I am adding the dual solution for each solver

• [x] OSQP
• [x] EqualityConstrainedQPSolver
• [x] Snopt
• [x] Ipopt
• [x] Gurobi
• [ ] Mosek
• [ ] SCS
• [ ] CSDP

Some of the solvers don't provide a dual solution (like MobyLCP).

[TobiaMarcucci]

Probably these are trivial observations... Just thought I could share my experience here.

I've been playing for a while with the dual solutions from Gurobi, and I remember I needed a lot of reverse engineering to understand the convention they use (e.g. sign of the multipliers). Nowhere in the docs they say how, e.g., the Lagrangian is defined. Moreover other solvers have other conventions. So this might require extra work than just exposing the user to the dual solutions from the solvers?

I think it'd be great to have somewhere in the Drake documentation a place where it is explained clearly how the signs are taken, and then reformat the outputs from the solvers in a coherent way?

It might be a ton of work though... I'd be happy to help in some way, despite my scarce coding skills.

[hongkai-dai]

Thanks Tobia!

I also observe that interpreting dual solution is hard, especially that for the same constraint, the dual variable can have different meanings in different solvers. I think maybe the easiest way to explain how to interpret the dual solution is through a jupyter notebook tutorial, where we can both explain it in text, and also the user can play with code blocks. What do you think?

[TobiaMarcucci]

Yes, I think a Jupyter notebook would be great! Probably that's more appropriate than the Drake docs for explaining that?

I think that a Jupyter notebook with a mathematical program with many different constraints (e.g. inequalities in all directions etc.) to show what are the signs of the resulting multipliers would be super useful, especially for the ones that need to do some algebra with the duals returned by MathematicalProgram.