Initial PR2: BM solver for QCQP-Representable Factor Graphs.#2542
Initial PR2: BM solver for QCQP-Representable Factor Graphs.#2542zhexin1904 wants to merge 1 commit into
Conversation
There was a problem hiding this comment.
Pull request overview
Note
Copilot was unable to run its full agentic suite in this review.
Adds a Burer–Monteiro Riemannian Staircase optimizer for QCQP-representable factor graphs, with certificate verification (Spectra / dense eigen) and a rotation-averaging example.
Changes:
- Introduces
RiemannianStaircaseOptimizer(outer staircase + ALM inner loop) and certificate construction/verification utilities. - Adds a certifiable rotation averaging example executable (2D/3D) using the new optimizer.
- Adds a comprehensive CppUnitLite test suite covering layout, lifting, rounding, verification, and rejection cases.
Reviewed changes
Copilot reviewed 5 out of 5 changed files in this pull request and generated 7 comments.
Show a summary per file
| File | Description |
|---|---|
| gtsam/constrained/tests/testRiemannianStaircase.cpp | New unit tests for staircase behavior (certificate, lift, layout, rounding, extractors). |
| gtsam/constrained/RiemannianStaircaseOptimizer.h | Public API for staircase optimizer, params, layout, and result types. |
| gtsam/constrained/RiemannianStaircaseOptimizer.cpp | Implementation of layout stacking, certificate assembly, Spectra/dense verification, and staircase loop. |
| examples/CertifiableRotationAveraging.cpp | Example program demonstrating certifiable SO(2)/SO(3) synchronization via the staircase. |
| THANKS.md | Adds contributor acknowledgements for the new work. |
| Eigen::SparseMatrix<double> Id(n, n); | ||
| Id.setIdentity(); | ||
| Eigen::SparseMatrix<double> M = S + params_.eta * Id; | ||
| M.makeCompressed(); |
| EXPECT(assert_equal(Matrix(Yi), Matrix(YiNew.leftCols(Yi.cols())), 1e-15)); | ||
|
|
||
| const auto& slice = layout.sliceOf(key); | ||
| const Matrix expectedCol = alpha * v.segment(slice.offset, slice.rowDim); | ||
| EXPECT(assert_equal(expectedCol, Matrix(YiNew.rightCols(1)), 1e-15)); |
| // ============================================================================ | ||
|
|
||
| RiemannianStaircaseOptimizer::Layout | ||
| RiemannianStaircaseOptimizer::Layout::From(const Values& values) { |
There was a problem hiding this comment.
This Layout class is very similar I think to the existing (and powerful) VerticalBlockMatrix we have, except transposed? Pls evaluate whether we could just use that.
There was a problem hiding this comment.
Yeah I prefer to keep Layout. VerticalBlockMatrix is an int-indexed matrix container, but what we need here is a Key-indexed map, so most of Layout's methods would still be needed even if it sat on top of VerticalBlockMatrix.
| /// | ||
| /// Probes `QCQP(pMin)` once at construction so a non-QCQP-representable | ||
| /// graph fails here, not mid-iteration. | ||
| RiemannianStaircaseOptimizer(std::function<QcqpProblem(size_t)> QCQP, |
There was a problem hiding this comment.
This API places a big burden on the user: they have to know the internals of the solver, and know how to create the correct parameterized QCQP. I think maybe that’s inevitable if we start from the QCQP, but if we start from the NonlinearFactorGraph, we can supply the D parameter ourselves, right? Come to think of it, though - maybe not if variables are created using templating. We need a runtime D.
There was a problem hiding this comment.
Yeah, the reason we originally took a QCQP (more precisely, a rank-parameterized QCQP family) was to make the input requirement explicit at the type level: not every NonlinearFactorGraph can be passed to this solver, only one whose factors all override NonlinearFactor::qcqpFactors.
But I agree, using a NonlinearFactorGraph as the input would make the API simpler, so the new interface would look something like this?
- Riemannian Staircase takes a NonlinearFactorGraph directly.
- Internally rebuilds the rank-p QCQP at each staircase level via QcqpProblem(graph_, p).
- QCQP-representability is still checked at construction: we call QcqpProblem(graph_, pMin) once inside a try-catch; if any factor hits NonlinearFactor::qcqpFactors's base (which throws), we re-throw as invalid_argument with a clear "graph is not QCQP-representable" message.
The user does pick D at compile time, but only once. When building initialValues, they call InsertQcqpValue<RotT, D>( ). Most APIs used in the Riemannian staircase are already runtime?
QcqpProblem(graph, k) and qcqpFactors(K) are both runtime.
| /// `RotT` via `RotT::ClosestTo`. Slices of other sizes are skipped. | ||
| /// More constrained-variable extractors will be added as needed. | ||
| template <typename RotT> | ||
| static std::vector<std::pair<Key, RotT>> extractRotations( |
There was a problem hiding this comment.
Ow, that’s not great. We have type-specific machinery in what is supposed to be type agnostic by now. That rounding is something that will have to go in the conversion branch along with QcQpValue, no?
There was a problem hiding this comment.
OHH! Yes thanks! This is an important point.
There was a problem hiding this comment.
RiemannianStaircaseOptimizer is (should be) type-agnostic now. The QCQP traits may need small additions for more complex variable / factor types beyond rotation averaging (per-type conversion methods at the variable level), but the staircase itself should stay type-agnostic. Will see how that shakes out as I add other variables/factors.
|
Yep, I like that plan. By the way, about co-pilot comments, please evaluate each of them for merit, and resolve all of them either by fixing them or by saying not relevant. Please explicitly resolve so I know which ones you addressed and which one you didn't |
de9e45d to
228feb7
Compare
|
Hi Frank, the changes have been pushed and address most of the review comments. One point I am still thinking about is the output of the Riemannian Staircase solver. Currently, it returns the rounded rank-(d) solution as a stacked variable matrix(totally type-agnostic), but extracting the individual variable blocks should likely be type-specific, which we have not proposed yet. I think the appropriate design may become clearer once I add (SE(d)) and other variable types to this framework soon. |
|
Hi Jason, Great. I'll take a look tomorrow morning! |
|
Hi Frank, I just discussed this again with Dave. This is basically along the lines of what you discussed yesterday. Here are a few points we would like to add soon. Does this look good? 1. Add more usage examples / test cases
2. Provide different interfacesa. Optimizer interfaceWe can provide both interfaces: RiemannianStaircaseOptimizer(NonlinearFactorGraph)
RiemannianStaircaseOptimizer(QCQP)Internally, the RiemannianStaircaseOptimizer(NonlinearFactorGraph) {
// Convert the nonlinear factor graph to a QCQP
RiemannianStaircaseOptimizer(QCQP);
}b. Verification interfaceWe should also provide verification interfaces that do not require multipliers. ALM can return multipliers, but other solvers, including the one used in our paper, cannot. In those cases, we will compute the dual variables as part of the verification procedure. Therefore, we can provide multiple interfaces with different arguments, while keeping the output the same: the certificate matrix. verify(QCQP problem, Solution solution, Multipliers multipliers)
verify(QCQP problem, Solution solution)
// Additional verification interfaces as neededInternally, all of these interfaces should eventually construct the same certificate matrix and call the same eigenvalue-based verification routine: verify(Certificate certificate)This essentially enables a general workflow in which a solution produced by any local solver can be passed to the dual-certification procedure. |
|
Thanks ! Exciting progress ! About rotation-specific code in BM: that is a problem I think. Was your paper not generic on type once you had a QCQP? New plan looks good.
|
PR 2 — Burer-Monteiro Riemannian Staircase solver and RA example
Branch:
feature/QCQP_BM_solver. Depends on PR 1.Summary: add a Riemannian Staircase solver for the BM-factored SDP
relaxation of a QCQP, using ALM as the inner solver and Spectra for the
minimum-eigenvalue certificate check, with a usage example on 2D / 3D
rotation averaging.
Benchmark vs SE-Sync (random initialziation, initiaz rank = d)