Make dual program #21565
Make dual program #21565
Conversation
There was a problem hiding this comment.
+@hongkai-dai for feature review! Thanks!
A lot of this code should look familiar from clarabel_solver and scs_clarabel_common
Reviewable status: LGTM missing from assignee hongkai-dai, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits), missing label for release notes (waiting on @AlexandreAmice)
There was a problem hiding this comment.
Checkpoint.
Reviewed 1 of 17 files at r1.
Reviewable status: 9 unresolved discussions, LGTM missing from assignee hongkai-dai, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits), missing label for release notes (waiting on @AlexandreAmice)
solvers/aggregate_costs_constraints.h line 106 at r1 (raw file):
/** * Aggregate all convex (conic) constraints into conic standard form i.e. a * single constraint Ax + b ∈ K where K is the product of cones.
Should also explain what Aeq and beq are.
solvers/aggregate_costs_constraints.h line 108 at r1 (raw file):
* single constraint Ax + b ∈ K where K is the product of cones. */ void AggregateConvexConstraints(const MathematicalProgram& prog,
I think it is better to call this AggregateConvexConicConstraints. Not every convex constraint is a conic constraint.
I don't see the implementation of this constraint. Do you want to implement it?
solvers/aggregate_costs_constraints.h line 116 at r1 (raw file):
namespace internal { struct ConvexConstraintAggregationOptions {
Document each of these options?
solvers/aggregate_costs_constraints.h line 129 at r1 (raw file):
// needed by SCS and Clarabel to transcribe the problem to the solvers and // extract dual solutions. struct ConvexConstraintAggregationInfo {
Why supporting the form Ax + s = b; s in K here? I think only SCS or Clarabel use this form, so we can leave it in scs_clarabel_common.h
solvers/dual_convex_program.h line 13 at r1 (raw file):
/** Compute a dual of the convex program. In particular, given the convex * program * min c^T x + d s.t.
Use unicode for transpose here?
solvers/dual_convex_program.h line 30 at r1 (raw file):
* 2) The lorentz cone which is self-dual * 3) The rotated lorentz cone. Drake's rotated lorentz cone is not self-dual, * so λ ∈ (2, 2, I) K where K is the rotated lorentz cone.
I don't think the meaning of (2, 2, I) K is very clear.
If you prefer self-dual cone, then I think it is better not to use Drake's rotated Lorentz cone, but the version of
2 * x1 * x2 >= sqrt(x3^2 + ... + xn^2)
as the SelfDualRotatedLorentzCone.
solvers/dual_convex_program.h line 33 at r1 (raw file):
* 4) The positive semidefinite cone which is self-dual. Note that for λ ∈ the * positive semidefinite cone, we take the convention that λ corresponds to the * lower triangular elements of the symmetric matrix Λ. To obtain the matrix Λ
Also mention that the symmetric matrix it is in the column major?
solvers/dual_convex_program.h line 52 at r1 (raw file):
* The only except is if binding is a linear equality constraint in which case * primal_result.GetDualSolution(binding) = * -dual_result.GetSolution(constraint_to_dual_variable_map.at(binding)).
I suppose if you write the dual as
max beqᵀy + d - bᵀλ s.t.
c + Aeqᵀy - Aᵀλ = 0
λ ∈ K*
then we do not need to take the negation here, and we always have primal_result.GetDualSolution(binding) = dual_result.GetSolution(constraint_to_dual_variable_map.at(binding)).
solvers/dual_convex_program.h line 55 at r1 (raw file):
* This is because Drake's GetDualSolution adopts the shadow price * interpretation for linear equality constraints, which is the negative of the * dual formulation provided by this method.s
extra s in the end.
There was a problem hiding this comment.
Address checkpoint review.
Reviewable status: 9 unresolved discussions, LGTM missing from assignee hongkai-dai, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits), missing label for release notes (waiting on @AlexandreAmice)
solvers/aggregate_costs_constraints.h line 106 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
Should also explain what Aeq and beq are.
Done. Removed method.
solvers/aggregate_costs_constraints.h line 108 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
I think it is better to call this AggregateConvexConicConstraints. Not every convex constraint is a conic constraint.
I don't see the implementation of this constraint. Do you want to implement it?
No I don't actually. Good catch.
solvers/aggregate_costs_constraints.h line 116 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
Document each of these options?
Done.
solvers/aggregate_costs_constraints.h line 129 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
Why supporting the form
Ax + s = b; s in Khere? I think only SCS or Clarabel use this form, so we can leave it in scs_clarabel_common.h
I plan on using this to support both constructing dual programs and primal standard form. I have it in this form to avoid even more changes in this PR.
solvers/dual_convex_program.h line 13 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
Use unicode for transpose here?
Done.
solvers/dual_convex_program.h line 30 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
I don't think the meaning of
(2, 2, I) Kis very clear.If you prefer self-dual cone, then I think it is better not to use Drake's rotated Lorentz cone, but the version of
2 * x1 * x2 >= sqrt(x3^2 + ... + xn^2)as the SelfDualRotatedLorentzCone.
Is the new spelling better? It is also more accurate.
solvers/dual_convex_program.h line 33 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
Also mention that the symmetric matrix it is in the column major?
Done.
solvers/dual_convex_program.h line 52 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
I suppose if you write the dual as
max beqᵀy + d - bᵀλ s.t. c + Aeqᵀy - Aᵀλ = 0 λ ∈ K*then we do not need to take the negation here, and we always have primal_result.GetDualSolution(binding) = dual_result.GetSolution(constraint_to_dual_variable_map.at(binding)).
This would break from the standard conic form convention. I am not sure which is best.
solvers/dual_convex_program.h line 55 at r1 (raw file):
Previously, hongkai-dai (Hongkai Dai) wrote…
extra
sin the end.
Done.
There was a problem hiding this comment.
Checkpoint
Reviewed 1 of 17 files at r1, 3 of 9 files at r2.
Reviewable status: 12 unresolved discussions, LGTM missing from assignee hongkai-dai, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits), missing label for release notes (waiting on @AlexandreAmice)
solvers/aggregate_costs_constraints.h line 114 at r2 (raw file):
// triangular part. If we wish tr(XY) = xᵀy, then we need to scale the off // diagonal entries of X by sqrt(2). If this option is true, then we perform // this scaling. See ParsePositiveSemidefiniteConstraints for more details.s
Extra s at the end of the line.
solvers/aggregate_costs_constraints.h line 135 at r2 (raw file):
std::vector<double> b_std; // The number of rows in the matrix A. int A_row_count{0};
It is redundant to keep A_row_count here as that is the same as b_std.size().
If you do need A_row_count, then I suggest to have a function
int A_row_count() const {
return b_std.size();
}solvers/aggregate_costs_constraints.h line 136 at r2 (raw file):
// The number of rows in the matrix A. int A_row_count{0}; // The total number of variables with bounding box constraints such that
BTW, I think you would also need to store the number of columns in A?
solvers/aggregate_costs_constraints.h line 156 at r2 (raw file):
std::vector<std::vector<std::pair<int, int>>> linear_constraint_dual_indices; // The lengths of each second order cone s in the Cartesian product K. See // // ParseSecondOrderConeConstraints for more details.
Extra //.
solvers/aggregate_costs_constraints.h line 191 at r2 (raw file):
// (0, 1) // This convention is compatible with both the SCS and Clarabel solvers. void DoAggregateConicConstraints(
Why call it DoAggregateConicConstraints instead of simpler AggregateConicConstraints?
solvers/aggregate_costs_constraints.h line 454 at r2 (raw file):
// triangular part of the symmetric matrix. SCS uses the lower triangular part, // and Clarabel uses the upper triangular part. // @param[in] upper_triangular Whether to scale the off-diagonal entries of the
Combine this line with the line above, they are both for the parameter upper_triangular.
solvers/aggregate_costs_constraints.cc line 276 at r2 (raw file):
// Parse the Bounding Box constraints. The bounding box constraints which are // equalities are immediately added to A and b. The bounding box constraints
Better to be more explicit, "Aeq and beq" instead of "A and b"
solvers/aggregate_costs_constraints.cc line 288 at r2 (raw file):
prog, // We can directly add the bounding box equality constraints // to A and b.
Ditto
solvers/aggregate_costs_constraints.cc line 912 at r2 (raw file):
scale_factor = sqrt2; } else if (i != j) { scale_factor = 2;
Sorry, why use scale_factor = 2 instead of scale_factor =1? I thought that when preserve_psd_inner_product_vectorization = False, we then only store the lower/upper triangular part directly.
solvers/dual_convex_program.h line 56 at r2 (raw file):
* interpretation for linear equality constraints, which is the negative of the * dual formulation provided by this method. */
Should this function throw an exception if prog contains non-conic constraint? If yes, better to mention it in the documentation.
solvers/dual_convex_program.h line 59 at r2 (raw file):
std::unique_ptr<MathematicalProgram> CreateDualConvexProgram( const MathematicalProgram& prog, std::unordered_map<Binding<Constraint>, MatrixX<symbolic::Expression>>*
It is surprising to use Matrix<symbolic::Expression> instead of Matrix<symbolic::Variable>. Could you explain a bit?
solvers/dual_convex_program.cc line 25 at r2 (raw file):
std::initializer_list<ProgramAttribute>{ // Supported Constraints. ProgramAttribute::kLinearEqualityConstraint,
Why not to support bounding box constraint as well?
There was a problem hiding this comment.
Reviewable status: 15 unresolved discussions, LGTM missing from assignee hongkai-dai, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits), missing label for release notes (waiting on @AlexandreAmice)
solvers/aggregate_costs_constraints.h line 125 at r2 (raw file):
// the form. // A x + s = b // s in K
Is there any special ordering on the cones in K? From the code in CreateDualConvexProgram, it seems that you assume the linear equality constraints come first, namely K start with the zero cones. If that is the case, we should be very explicit in the documentation.
solvers/dual_convex_program.cc line 59 at r2 (raw file):
auto A_triplets_end_of_equalities_iterator = info.A_triplets.begin(); std::advance(A_triplets_end_of_equalities_iterator, info.num_linear_equality_constraint_rows);
I am not sure about this code. Although the first info.num_linear_equality_constraint_rows correspond to the equality constraints, that doesn't mean the first info.num_linear_equality_constraint_rows entries in A_triplets correspond to equality constraint. It is possible that one equality constraint have multiple non-zero entries in A.
I think this part of the code is brittle, it basically makes these assumptions
- The first num_linear_equality_constraint_rows in
A*x+s=bcorrespond tos=0. - The first chunk of entries in A_triplets are obtained from the first chunk of rows in
A*x+s=b, but this doesn't have to be true.
I would suggest to change the code in ConicConstraintAggregationInfo to be more specific, namely have a separate Aeq, beq fields in this struct to store the equality constraints.
solvers/dual_convex_program.cc line 62 at r2 (raw file):
// The first info.num_linear_equality_constraint_rows correspond to equality // constraints.
This assumption is never explicitly stated in ConicConstraintAggregationInfo. The code is brittle by making this assumption but never explicitly state it.
|
This seems to have stalled out. Is there a plan to continue work? If not, I will plan to close it in a few weeks, to conserve CI resources. |
jwnimmer-tri
added
the
status: defer ci
This change is