Quick Start Guide

Installation

POPR can be installed by running from a terminal:

git clone --recurse-submodules git@github.com:duembgen/popcor
cd popcor
conda env create -f environment.yml

Problem Formulation

We start with polynomial optimization problems (POPs) of the form:

\[\begin{split}\begin{align} q^\star =&\min_{\theta} f(\theta) \\ \text{s.t. } &g(\theta) = 0 \\ &h(\theta) \geq 0 \end{align}\end{split}\]

where \(f,g,\text{ and } h\) are polynomial functions, and both \(g\) and \(h\) can be vector-valued. Many maximum-a-posteriori or maximum-likelihood estimation problems can be formulated as such, for example range-only localization and range-aided SLAM, (matrix-weighted) SLAM, and outlier-robust estimation. The same is true for many control and planning problems, for example the inverted pendulum and other classical dynamical systems, and even contact-rich problems such as slider-pusher planning problems.

Any POP can be equivalently written in the following QCQP form:

\[\begin{split}\begin{align} q^\star =&\min_{x} x^\top Q x \\ \text{s.t. } &(\forall i): x^\top A_i x = b_i \\ &(\forall j): x^\top B_j x \geq 0 \end{align}\end{split}\]

with cost matrix \(Q\), known constraint matrices \(A_i,B_j\). Note that

  • We always include the so-called homogenization variable, which enables to write linear and constant terms as quadratics. By convention, we set the first element of \(x\) to one, and we use \(b_0=1, A_0\) to encorce this constraint.

  • All inequality and some equality constraints correspond to the constraints from the original POP.

  • Some additional equality constraints correspond to new substitution variables that need to be added to formulate the problem as a quadratic.

Warning

Note that while inequality constraints can be added to the problem formulation, there is no implementation yet to add find and add redundant inequality constraints to the relaxation.

For the standard usage, the user first needs to define a custom Lifter class which essentially contains all elements related to the QCQP problem formulation. This class should inherit from StateLifter. A basic skeleton of such a Lifter class is provided in Example for AutoTight. The main purpose of this class is that it provides all basic operations related to the problem formulation, such as:

For a bit more advanced functionality (for example for the SDP Relaxation in the next section), you also need to define functions such as

Many example lifters are provided, you can find them under Examples.

Example: instantiating and using lifter

The following code snippet shows some basic operations (and useful sanity checks) for the example lifter class popcor.examples.Poly4Lifter. Note that this and all following examples can be found in the file ../../tests/test_quickstart.py.

from popcor.examples import Poly4Lifter

lifter = Poly4Lifter()

Q = lifter.get_Q()

# theta corresponds to the ground truth; in this case, the global minimum.
x = lifter.get_x(lifter.theta).flatten()
cost_optimum = float(x.T @ Q @ x)

# the cost at any other randomly sampled point has to be larger.
for i in range(10):
    theta_random = lifter.sample_theta()
    x_random = lifter.get_x(theta_random).flatten()
    assert float(x_random.T @ Q @ x_random) > cost_optimum

SDP Relaxation

It is straightforward to derive a convex relaxation of the original QCQP, using the reformulation \(x^\top Qx=\langle x, Qx\rangle = \langle Q, xx^\top \rangle\), where \(\langle \cdot, \cdot \rangle\) denotes the trace inner product. Then introducing \(X:=xx^\top\) and relaxing its rank, we obtain the following convex relaxation, in the form of an SDP:

\[\begin{split}\begin{align} p^\star = &\min_{X \succeq 0} \langle Q, X \rangle \\ \text{s.t. } &(\forall i): \langle A_i, X \rangle = b_i \\ &(\forall j): \langle B_j, X \rangle \geq 0 \end{align}\end{split}\]

Example: solving the QCQP using rank relaxation

The following code snippet shows how you can use the simple lifter from earlier to find the global optimum of the nonconvex polynomial problem, by solving an SDP.

from cert_tools.linalg_tools import rank_project
from cert_tools.sdp_solvers import solve_sdp

from popcor.examples import Poly4Lifter

lifter = Poly4Lifter()

# the cost matrix
Q = lifter.get_Q()

# the equality constraints
A_known = lifter.get_A_known()

# the homogenization constraint
A_0 = lifter.get_A0()

X, info = solve_sdp(Q, [(A_0, 1.0)] + [(A_i, 0.0) for A_i in A_known])
assert X is not None

# if X is rank one, the global optimum can be found in element X_10 of the matrix.
theta_pick = X[1, 0]
assert abs(theta_pick - lifter.theta) < 1e-5

# We can also first extract the rank-1 estimate (X=xx') and then extract theta.
x, info_rank = rank_project(X, p=1)
theta_round = x[1]

assert abs(theta_round - lifter.theta) < 1e-5

AutoTight Method

AutoTight is used to find all possible constraints matrices \(A_r\), which are also automatically satisfied by solutions of the QCQP. They are also called redundant constraints because they do not change the feasible set of the original problem, but when adding those constraints to the SDP (rank-)relaxation, they often improve tightness. Denoting by \(A_r\) the redundant constraints, we can solve the following SDP:

\[\begin{split}\begin{align} p_r^\star = &\min_{X \succeq 0} \langle Q, X \rangle \\ \text{s.t. } &(\forall i): \langle A_i, X \rangle = b_i \\ &(\forall r): \langle A_r, X \rangle = 0 \\ &(\forall j): \langle B_j, X \rangle \geq 0 \end{align}\end{split}\]

We use the term cost-tight to say that strong duality holds (\(p_r^\star = q^\star\)) while by rank-tight we denote the fact that the SDP solver returns a rank-one solution. If successful, the output is a set of constraints that leads to a tight SDP relaxation of the original problem, which can be used to solve the problem to global optimality (if we have rank tightness) or certify given solutions (if we have cost tightness).

More information on how to use AutoTight can be found here and a simple example is given next.

Example: tightening the SDP relaxation using AutoTight

from cert_tools.sdp_solvers import solve_sdp

from popcor.auto_tight import AutoTight
from popcor.examples.stereo1d_lifter import Stereo1DLifter

lifter = Stereo1DLifter(n_landmarks=5)

auto_tight = AutoTight()

# solve the SDP -- it is not cost tight!
Q = lifter.get_Q(noise=1e-1)
A_known = lifter.get_A_known()
A_0 = lifter.get_A0()

# solve locally starting at ground truth
x, info_local, cost = lifter.local_solver(lifter.theta, y=lifter.y_)
assert x is not None
assert info_local["success"]

constraints = lifter.get_A_b_list(A_list=A_known)
X, info_sdp = solve_sdp(Q, constraints)

gap = auto_tight.get_duality_gap(
    info_local["cost"], info_sdp["cost"], relative=True
)
assert gap > 0.1

# learn matrices and solve again
A_learned = auto_tight.get_A_learned(lifter)
constraints = lifter.get_A_b_list(A_list=A_learned)
X, info_sdp_learned = solve_sdp(Q, constraints)
gap = auto_tight.get_duality_gap(
    info_local["cost"], info_sdp_learned["cost"], relative=True
)

# note that the gap can be slightly negative because of mismatch in convergence tolerances etc.
assert abs(gap) < 1e-2

# problem: if we change landmarks, the constraints do not generalize!
new_lifter = Stereo1DLifter(n_landmarks=5)
assert not np.all(new_lifter.landmarks == lifter.landmarks)

AutoTemplate Method

AutoTemplate follows the same principle as AutoTight, but its output are templates rather than constraint matrices. These templates can be seen as “parametrized” versions of the constraint matrices, and can be applied to new problem instances of any size without having to learn the constraints again from scratch.

More information on how to use AutoTemplate can be found here and a simple example is given next.

Example: tightening a different problem using AutoTemplate



test_autotemplate():
"""Test the AutoTemplate example from the documentation."""
from cert_tools.sdp_solvers import solve_sdp

from popcor.auto_template import AutoTemplate
from popcor.auto_tight import AutoTight
from popcor.examples.stereo1d_lifter import Stereo1DLifter

# important: we need to use param_level="p", otherwise the parameters
# are not factored out and the constraints are not generalizable.
lifter = Stereo1DLifter(n_landmarks=3, param_level="p")

# learn the template matrices
auto_template = AutoTemplate(lifter)
data, success = auto_template.run(use_known=False, plot=True)
assert success

# apply the templates to a different lifter
new_lifter = Stereo1DLifter(n_landmarks=5, param_level="p")
A_learned = auto_template.apply(new_lifter, use_known=True)

Q = new_lifter.get_Q(noise=1e-1)

# adds homogenization constraint and all b_i terms
constraints = new_lifter.get_A_b_list(A_list=A_learned)
X, info_sdp = solve_sdp(Q, constraints)

# evaluate duality gap

References

[1] F. Dümbgen, C. Holmes, B. Agro and T. Barfoot, “Toward Globally Optimal State Estimation Using Automatically Tightened Semidefinite Relaxations,” in IEEE Transactions on Robotics, vol. 40, pp. 4338-4358, 2024, doi: 10.1109/TRO.2024.3454570.