In Perfect Shape: Certifiably Optimal 3D Shape Reconstruction from 2D Landmarks

Heng Yang, Luca Carlone

Introduction

3D object detection and pose estimation from a single image is a fundamental problem in computer vision. Despite the progress in semantic segmentation , depth estimation , and pose estimation , reconstructing the 3D shape and pose of an object from a single image remains a challenging task .

A typical approach for 3D shape reconstruction is to first detect 2D landmarks in a single image, and then solve a model-based optimization to lift the 2D landmarks to form a 3D model . For the optimization to be well-posed, the unknown shape is assumed to be a 3D deformable model, composed by a linear combination of basis shapes, handcrafted or learned from a large corpus of training data . The optimization then seeks to jointly optimize the coefficients of the linear combination (shape parameters) and the camera pose to minimize the reprojection errors between the 3D model and the 2D landmarks. This model-based paradigm has been successful in several applications such as face recognition , car model fitting , and human pose estimation .

Despite its long history and broad range of applications, there is still no globally optimal solver for the non-convex optimization problem arising in 3D shape reconstruction. Therefore, most existing solutions adopt a local optimization strategy, which alternates between solving for the camera pose and the shape parameters. These techniques, as shown in prior works , require an initial guess for the solution and often get stuck in local minima. In addition, 2D landmark detectors are prone to produce outliers, causing existing methods to be brittle . Therefore, the motivation for this paper is two-fold: (i) to develop a certifiably optimal shape reconstruction solver, and (ii) to develop a robust reconstruction algorithm that is insensitive to a large amount of outlier 2D measurements (e.g., 70%70\%).

Contributions. Our first contribution is to formulate the shape reconstruction problem as a polynomial optimization problem and apply Lasserre’s hierarchy of Sums-of-Squares (SOS) relaxations to relax the non-convex polynomial optimization into a convex semidefinite program (SDP). We show the SOS relaxation of minimum order 2 empirically solves the non-convex shape reconstruction problem exactly and provides a global optimality certificate. The second contribution is to apply basis reduction, a technique that exploits the sparse structure of the polynomial in the objective function, to reduce the size of the resulting SDP. We show that basis reduction significantly improves the efficiency of the SOS relaxation without compromising global optimality. To the best of our knowledge, this is the first certifiably optimal solver for shape reconstruction, and we name it Shape⋆. Our third contribution is to robustify Shape⋆ by adopting a truncated least squares (TLS) robust cost function and solving the resulting robust estimation problem using graduated non-convexity . The resulting algorithm, named Shape#{\#}, is robust against 70%70\% outliers and does not require an initial guess.

The rest of this paper is organized as follows. Section 2 reviews related work. Section 3 introduces notation and preliminaries on SOS relaxations. Section 4 introduces the shape reconstruction problem. Section 5 develops our SOS solver (Shape⋆). Section 6 presents an algorithm (Shape#{\#}) to robustify the SOS relaxation against outliers. Section 7 provides experimental results in both simulations and real datasets, while Section 8 concludes the paper.

Related Work

We limit our review to optimization-based approaches for 3D shape reconstruction from 2D landmarks. The interested reader can find a review of end-to-end shape and pose reconstruction using deep learning in .

Notation and Preliminaries

We now give a brief summary of SOS relaxations for polynomial optimization. Our review is based on . We first introduce the notion of SOS polynomial.

and Q{\bm{Q}} is called the Gram matrix of p(x)p(\bm{x}).

Now consider the following polynomial optimization:

the ideal and the 2β2\beta-th truncated ideal of h\bm{h}, where deg(⋅)\text{deg}(\cdot) is the degree of a polynomial. The ideal is simply a summation of polynomials with polynomial coefficients, a construct that will simplify the notation later on. We call

the quadratic module and the β\beta-th truncated quadratic module generated from g\bm{g}. Note that the quadratic module is similar to the ideal, except now we require the polynomial coefficients to be SOS. Apparently, if p(x)∈⟨h⟩+Q(g)p(\bm{x})\in\langle\bm{h}\rangle+Q(\bm{g}), then p(x)p(\bm{x}) is nonnegative on X{\cal X}If p∈⟨h⟩+Q(g)p\in\langle\bm{h}\rangle+Q(\bm{g}), then p=h+gp=h+g, with h∈⟨h⟩h\in\langle\bm{h}\rangle and g∈Q(g)g\in Q(\bm{g}). For any x∈X\bm{x}\in{\cal X}, since hi(x)=0h_{i}(\bm{x})=0, so h(x)=∑λihi=0h(\bm{x})=\sum\lambda_{i}h_{i}=0; since gk(x)≥0g_{k}(\bm{x})\geq 0 and sk(x)≥0s_{k}(\bm{x})\geq 0, so g=∑skgk≥0g=\sum s_{k}g_{k}\geq 0. Therefore, p=h+g≥0p=h+g\geq 0. Putinar’s Positivstellensatz describes when the reverse is also true.

Based on Putinar’s Positivstellensatz, Lasserre derived a sequence of SOS relaxations that approximates the global minimum of problem (3) with increasing accuracy. The key insight behind Lasserre’s hierarchy is twofold. The first insight is that problem (3), which we can write succinctly as min⁡x∈Xf(x)\min_{\bm{x}\in{\cal X}}f(x), can be equivalently written as max⁡x,γγ,s.t.f(x)−γ≥0 on X\displaystyle\max_{\bm{x},\gamma}\gamma,s.t.f(x)-\gamma\geq 0\text{ on }{\cal X} (intuition: the latter pushes the lower bound γ\gamma to reach the global minimum of f(x)f(\bm{x})). The second intuition is that we can rewrite the condition f(x)−γ≥0 on Xf(x)-\gamma\geq 0\text{ on }{\cal X}, using Putinar’s Positivstellensatz (Theorem 2), leading to the following hierarchy of Sums-of-Squares relaxations.

Lasserre’s hierarchy of order β\beta is the following SOS program:

which can be written as a standard SDP. Moreover, let f⋆f^{\star} be the global minimum of (3) and fβ⋆f^{\star}_{\beta} be the optimal value of (8), then fβ⋆f^{\star}_{\beta} monotonically increases and fβ⋆→f⋆f^{\star}_{\beta}\rightarrow f^{\star} when β→∞\beta\rightarrow\infty. More recently, Nie proved that under Archimedeanness, Lasserre’s hierarchy has finite convergence generically (i.e., fβ⋆=f⋆f^{\star}_{\beta}=f^{\star} for some finite β\beta).

In computer vision, Lasserre’s hierarchy was first used by Kahl and Henrion to minimize rational functions arising in geometric reconstruction problems, and more recently by Probst et al. as a framework to solve a set of 3D vision problems. In this paper we will show that the SOS relaxation as written in eq. (8) allows using basis reduction to exploit the sparsity pattern of polynomials and leads to significantly smaller semidefinite programs.

Problem Statement: Shape Reconstruction

Certifiably Optimal Shape Reconstruction

This section shows how to develop a certifiably optimal solver for problem (13). Our first step is to algebraically eliminate the translation t\bm{t} and obtain a translation-free shape reconstruction problem, as shown below.

The shape reconstruction problem (13) is equivalent to the following translation-free optimization:

Further, let R⋆{\bm{R}}^{\star} and ck⋆,k=1,…,K,c_{k}^{\star},k=1,\dots,K, be the global minimizer of the above translation-free optimization (14), then the optimal translation t⋆\bm{t}^{\star} can be recovered as:

A formal proof of Theorem 4 can be found in the Supplementary Material. The intuition behind Theorem 4 is that if we express the landmark coordinates and 3D basis shapes with respect to their (weighted) centroids zˉw\bar{\bm{z}}^{w} and Bˉkw,k=1,…,K\bar{{\bm{B}}}^{w}_{k},k=1,\dots,K, we can remove the dependence on the translation t\bm{t}. This strategy is inspired by Horn’s method for point cloud registration , and generalizes to the weighted and non-centered case.

This section applies Lasserre’s hierarchy as described in Theorem 3 to solve the translation-free shape reconstruction problem (14). We do this in two steps: we first show problem (14) can be formulated as a polynomial optimization in the form (3); and then we add valid constraints to make the feasible set Archimedean.

then it becomes clear that qi(x)q_{i}(\bm{x}) is a polynomial function of x\bm{x} with degree 4. Because the Lasso regularization is linear in c\bm{c}, the objective function f(x)f(\bm{x}) is a degree-4 polynomial.

In summary, the translation-free problem (14) is equivalent to a polynomial optimization with a degree-4 objective f(x)f(\bm{x}), constrained by 15 quadratic equalities hi(x)h_{i}(\bm{x}) (eq. (19)) and KK linear inequalities gk(x)=ckg_{k}(\bm{x})=c_{k}.

Now we can certify the Archimedeanness of ⟨h⟩+Q(g)\langle\bm{h}\rangle+Q(\bm{g}):

with M=K+3M=K+3 and β=1\beta=1 (cf. Theorem 2).

Apply Lasserre’s Hierarchy. With Archimedeanness, we can now apply Lasserre’s hierarchy of SOS relaxations.

The SOS relaxation of order β\beta (β≥2\beta\geq 2)The minimum relaxation order is 2 because f(x)f(\bm{x}) has degree 4. for the translation-free shape reconstruction problem (14) is the following convex semidefinite program:

where f(x)f(\bm{x}) is the objective function defined in (14), gk(x),k=1,…,2Kg_{k}(\bm{x}),k=1,\dots,2K are the inequality constraints ck,1−ck2c_{k},1-c_{k}^{2}, hi(x),i=1,…,15h_{i}(\bm{x}),i=1,\dots,15 are the equality constraints defined in (19), and N0:=(K+9+ββ)N_{0}:=\left(\begin{subarray}{c}K+9+\beta\\ \beta\end{subarray}\right), Ns:=(K+8+ββ−1)N_{s}:=\left(\begin{subarray}{c}K+8+\beta\\ \beta-1\end{subarray}\right), Nλ:=(K+7+2β2β−2)N_{\lambda}:=\left(\begin{subarray}{c}K+7+2\beta\\ 2\beta-2\end{subarray}\right) are the sizes of matrices and vectors.

While a formal proof of Proposition 6 is given in the Supplementary Material, we observe that (22) immediately results from the application of Lasserre’s hierarchy to (8), by parametrizing Qβ(g)Q_{\beta}(\bm{g}) with monomial bases [x]β−1[\bm{x}]_{\beta-1}, [x]β[\bm{x}]_{\beta} and PSD matrices S0{\bm{S}}_{0}, Sk,k=1,…,2K{\bm{S}}_{k},k=1,\dots,2K (one for each gkg_{k}), and by parametrizing ⟨h⟩2β\langle\bm{h}\rangle_{2\beta} with monomial basis [x]2β−2[\bm{x}]_{2\beta-2} and coefficient vectors λi,i=1,…,15\bm{\lambda}_{i},i=1,\dots,15 (one for each hih_{i}). Problem (22) can be written as an SDP and solved globally using standard convex solvers (e.g. YALMIP ). We call the SDP written in (22) the primal SDP. The dual SDP of (22) can be derived using moment relaxation , which is readily available in GloptiPoly 3 .

Extract Solutions from SDP. After solving the SDP (22), we can extract solutions to the original non-convex problem (14), a procedure we call rounding.

Let fβ⋆=γ⋆f^{\star}_{\beta}=\gamma^{\star} and S0β⋆,Skβ⋆,λiβ⋆{\bm{S}}_{0}^{\beta\star},{\bm{S}}_{k}^{\beta\star},\bm{\lambda}_{i}^{\beta\star} be the optimal solutions to the SDP (22) at order β\beta; compute vβ⋆\bm{v}^{\beta\star} as the eigenvector corresponding to the minimum eigenvalue of S0β⋆{\bm{S}}_{0}^{\beta\star}, and then normalize vβ⋆\bm{v}^{\beta\star} such that the first entry is equal to 1. Then an approximate solution to problem (14) can be obtained as:

where f⋆f^{\star} is the true (unknown) global minimum of problem (14). We define the relative duality gap ηβ\eta_{\beta} as:

which quantifies the quality of the SOS relaxation.

Certifiable Global Optimality. Besides extracting solutions to the original problem, we can also verify when the SOS relaxation solves the original problem exactly.

Let fβ⋆=γ⋆f^{\star}_{\beta}=\gamma^{\star} and S0β⋆{\bm{S}}_{0}^{\beta\star} be the optimal solutions to the SDP (22) at order β\beta. If corank(S0β⋆)=1\text{corank}({\bm{S}}_{0}^{\beta\star})=1 (the corank is the dimension of the null space of S0β⋆{\bm{S}}_{0}^{\beta\star}), then fβ⋆f^{\star}_{\beta} is the global minimum of problem (14), and the relaxation is said to be tight at order β\beta. Moreover, the relative duality gap ηβ=0\eta_{\beta}=0 and the solution x^β\hat{\bm{x}}^{\beta} extracted using Proposition 7 is the unique global minimizer of problem (14).

The proof of Theorem 8 is given in the Supplementary Material. Empirically (Section 7), we observed that the relaxation is always tight at the minimum relaxation order β=2\beta=2. Note that even when the relaxation is not tight, one can still obtain an approximate solution using Proposition 7 and quantify how suboptimal the approximate solution is using the relative duality gap ηβ\eta_{\beta}.

2 Basis Reduction

Despite the theoretical soundness and finite convergence at order β=2\beta=2, the size of the SDP (22) is N0=(K+9+ββ)N_{0}=\left(\begin{subarray}{c}K+9+\beta\\ \beta\end{subarray}\right), which for β=2\beta=2 becomes (K+112)\left(\begin{subarray}{c}K+11\\ 2\end{subarray}\right), implying that the size of the SDP grows quadratically in the number of bases KK. Although there have been promising advances in improving the scalability of SDP solvers (see for a thorough review), such as exploiting sparsity and low-rankness , in this section we demonstrate a simple yet effective approach, called basis reduction, that exploits the structure of the objective function to significantly reduce the size of the SDP in (22).

In a nutshell, basis reduction methods seek to find a smaller, but still expressive enough, subset of the full vector of monomials [x]β[\bm{x}]_{\beta} on the right-hand side (RHS) of eq. (23), to explain the objective function f(x)f(\bm{x}) on the left-hand side (LHS). There exist standard approximation algorithms for basis reduction, discussed in and implemented in YALMIP . However, in practice we found the basis selection method in YALMIP failed to find any reduction for the SDP (22). Therefore, here we propose a problem-specific reduction, which follows from the examination of which monomials appear on the LHS of (23).

The SOS relaxation of order β=2\beta=2 with basis reduction for the translation-free shape reconstruction problem (14) is the following convex semidefinite program:

Comparing the SDP (27) and (22), the most significant change is replacing the full monomial basis [x]β[\bm{x}]_{\beta} in (22) with a much smaller monomial basis m2(x)m_{2}(\bm{x}) that excludes degree-2 monomials purely supported in c\bm{c} and r\bm{r}. This reduction is motivated by analyzing the monomial terms in f(x)f(\bm{x}). Although a formal proof of the equivalence between (22) and (27) remains open, we provide an intuitive explanation in the Supplementary Material. After basis reduction, the size of the SDP (27) is N0′=10K+10N_{0}^{\prime}=10K+10, which is linear in KK and much smaller than the size of the original SDP (22) Nk=(K+112)N_{k}=\left(\begin{subarray}{c}K+11\\ 2\end{subarray}\right)For K=5,10,20K=5,10,20, N0=120,210,465N_{0}=120,210,465, while N0′=60,110,210N_{0}^{\prime}=60,110,210.. Section 7 numerically shows that the SDP after basis reduction gives the same (tight) solution as the original SDP.

3 Shape⋆: Algorithm Summary

To summarize the derivation in this section, our solver for the shape reconstruction problem (13), named Shape⋆, works as follows. It first solves the SDP (27) and applies the rounding described in Proposition 7 to compute an estimate of the shape parameters ckc_{k} and rotation R{\bm{R}} and possibly certify its optimality. Then, Shape⋆ uses the closed-form expression (17) to retrieve the translation estimate t\bm{t}.

Robust Outlier Rejection

Section 5 proposed a certifiably optimal solver for problem (13). However, the least squares formulation (13) tends to be sensitive to outliers: the pixel measurements Z{\bm{Z}} in eq. (9) are typically produced by learning-based or handcrafted detectors , which might produce largely incorrect measurements (e.g. due to wrong data association zi↔Bki\bm{z}_{i}\leftrightarrow{\bm{B}}_{ki}), which in turn leads to poor shape reconstruction results. This section shows how to regain robustness by iteratively solving the weighted least squares problem (13) and adjusting the weights wiw_{i} to reject outliers.

The key insight is to substitute the least square penalty in (13) with a robust cost function, namely the truncated least squares (TLS) cost . Hence, we propose the following TLS shape reconstruction formulation:

where ri(ck,R,t):=∥zi ⁣− ⁣ΠR(∑k=1KckBki) ⁣− ⁣t∥r_{i}(c_{k},{\bm{R}},\bm{t}):=\left\|\bm{z}_{i}\!-\!\Pi{\bm{R}}\left(\sum_{k=1}^{K}c_{k}{\bm{B}}_{ki}\right)\!-\!\bm{t}\right\| (introduced for notational convenience), and ρcˉ(r)=min⁡(r2,cˉ2)\rho_{\bar{c}}(r)=\min(r^{2},\bar{c}^{2}) implements a truncated least squares cost, which is quadratic for small residuals and saturates to a constant value for residuals larger than a maximum error cˉ\bar{c}.

Our second insight is that ρcˉ(r)\rho_{\bar{c}}(r) can be written as ρcˉ(r)=min⁡w∈{0,1}wr2+(1−w)cˉ2\rho_{\bar{c}}(r)=\min_{w\in\{0,1\}}wr^{2}+(1-w)\bar{c}^{2}, by introducing extra slack binary variables w∈{0,1}w\in\{0,1\}. Therefore, we can write problem (29) equivalently as:

The final insight is that now we can minimize (30) by iteratively minimizing (i) over ck,R,tc_{k},{\bm{R}},\bm{t} (with fixed weights wiw_{i}), and (ii) over the weights wiw_{i} (with fixed ck,R,tc_{k},{\bm{R}},\bm{t}). The rationale for this approach is that step (i) can be implemented using Shape⋆ (since in this case the weights are fixed), and step (ii) can be implemented in closed-form. To improve convergence of this iterative algorithm, we adopt graduated non-convexity , which starts with a convex approximation of problem (30) and uses a control parameter μ\mu to gradually increase the amount of non-convexity, till (for large μ\mu) one solves (30). The resulting algorithm named Shape#{\#} is given in Algorithm 1. We refer the reader to the Supplementary Material and for a complete derivation of Algorithm 1 and for the closed-form expression of the weight update in line 1 of the algorithm.

Shape#{\#} is deterministic and does not require an initial guess. We remark that the graduated non-convexity scheme in Shape#{\#} (contrarily to Shape⋆) is not guaranteed to converge to an optimal solution of (30), but we show in the next section that it is empirically robust to 70%70\% outliers.

Experiments

Implementation details. Both Shape⋆ and Shape#{\#} are implemented in Matlab, with both SOS relaxations (22) and (27) implemented using YALMIP and the resulting SDPs solved using MOSEK .

2 Shape⋆ for Outlier-Free Reconstruction

3 Shape##{\#} for Robust Reconstruction

Conclusions

We presented Shape⋆, the first certifiably optimal solver for 3D shape reconstruction from 2D landmarks in a single image. Shape⋆ is developed by applying Lasserre’s hierarchy of SOS relaxations combined with basis reduction to improve efficiency. Experimental results show that the SOS relaxation of order 2 always achieves global optimality. To handle outlying measurements, we also proposed Shape#{\#}, which solves a truncated least squares robust estimation problem by iteratively running Shape⋆ without the need for an initial guess. We show that Shape#{\#} achieves robustness against 70%70\% outliers on the FG3DCar dataset and outperforms state-of-the-art solvers.

Proof of Theorem 4

Here we prove Theorem 4 in the main document. Recall the weighted least squares optimization for shape reconstruction in eq. (13) and denote its objective function as f(c,R,t)f(\bm{c},{\bm{R}},\bm{t}), with c=[c1,…,ck]T\bm{c}=[c_{1},\dots,c_{k}]^{\mathsf{T}}:

In order to marginalize out the translation t\bm{t}, we compute the derivative of f(c,R,t)f(\bm{c},{\bm{R}},\bm{t}) w.r.t. t\bm{t}:

and set it to 0\bm{0}, which allows us to write t⋆\bm{t}^{\star} in closed form using R⋆{\bm{R}}^{\star} and c⋆\bm{c}^{\star}:

with zˉw\bar{\bm{z}}^{w} and Bˉkw,k=1,…,K,\bar{{\bm{B}}}^{w}_{k},k=1,\dots,K, being the weighted centers of the 2D landmarks Z{\bm{Z}} and the 3D basis shapes Bk{\bm{B}}_{k}:

Then we can substitute the expression of t⋆\bm{t}^{\star} in (A3) back into the objective function in (A1) and obtain an objective function without translation:

we can see the equivalence between the objective function in eq. (A5) and the objective function in eq. (14) of Theorem 4. The constraints remain unchanged because we only marginalize out the unconstrained variable t\bm{t}. Therefore, the shape reconstruction problem (13) is equivalent to the translation-free problem (14), and the optimal translation can be recovered using eq. (A3). ∎

Proof of Proposition 6

Here we prove the SOS relaxation of order β\beta (β≥2\beta\geq 2) for the translation-free shape reconstruction problem (14) is the semidefinite program in (22). First, let us rewrite the general form of Lasserre’s hierarchy of order β\beta in eq. (8) in Theorem 3 as the following:

In words, the constraints of (A8) ask the polynomial f(x)−γf(\bm{x})-\gamma to be written as a sum of two polynomials hh and gg, with hh in the 2β2\beta-th truncated ideal of h\bm{h}, and gg in the β\beta-th truncated quadratic module of g\bm{g}.

Next, we use the definition of the 2β2\beta-th truncated ideal and the β\beta-th truncated quadratic module to explicitly represent hh and gg. First recall the definition of the 2β2\beta-th truncated ideal in eq. (5), which states that hh must be written as a sum of polynomial products between the equality constraints hih_{i}’s and the polynomial multipliers λi\lambda_{i}’s:

with λi\bm{\lambda}_{i} being the vector of unknown coefficients associated with the monomial basis [x]2β−2[\bm{x}]_{2\beta-2}. The size of λi\bm{\lambda}_{i} is equal to the length of [x]2β−2[\bm{x}]_{2\beta-2}, which can be computed by (n+dd)\left(\begin{subarray}{c}n+d\\ d\end{subarray}\right), with n=K+9n=K+9 being the number of variables in x\bm{x}, and d=2β−2d=2\beta-2 being the maximum degree of the monomial basis. Similarly for gg, we recall the definition of the β\beta-th truncated quadratic module in eq. (7), which states that gg must be written as a sum of polynomial products between the inequality constraints gkg_{k}’s and the SOS polynomial multipliers sks_{k}’s:

and the degree of each polynomial product skgks_{k}g_{k} must be no greater than 2β2\beta, i.e., deg⁡(skgk)≤2β\deg(s_{k}g_{k})\leq 2\beta. For our specific shape reconstruction problem, we have g0:=1g_{0}:=1, gk=ck,k=1,…,Kg_{k}=c_{k},k=1,\dots,K, and gK+k=1−ck2,k=1,…,Kg_{K+k}=1-c_{k}^{2},k=1,\dots,K. Since g0g_{0} has degree 0, s0s_{0} can have degree up to 2β2\beta. All gk,k=1,…,Kg_{k},k=1,\dots,K, have degree 1, so sk,k=1,…,Ks_{k},k=1,\dots,K, can have degree up to 2β−12\beta-1. However, because SOS polynomials can only have even degree, sk,k=1,…,Ks_{k},k=1,\dots,K can only have degree up to 2β−22\beta-2. For gK+k,k=1,…,Kg_{K+k},k=1,\dots,K, they have degree 2, so their corresponding SOS polynomial multipliers sK+k,k=1,…,Ks_{K+k},k=1,\dots,K can have degree up to 2β−22\beta-2. Now for each SOS polynomial sk,k=0,…,2Ks_{k},k=0,\dots,2K, from the Gram matrix representation in eq. (2), we can associate a PSD matrix Sk{\bm{S}}_{k} with it using corresponding monomial bases:

Finally, by inserting the expressions of sks_{k} in (A12) back to the expression of gg in (A11), and inserting the expression of λi\lambda_{i} in (A10) back to the expression of hh in (A9), we can convert the SOS relaxation of general form (A8) to the semidefinite program (22). ∎

Proof of Theorem 8

According to , the dual SDP of (22) is the following SDP:

and v⋆\bm{v}^{\star} is in the null-space of S0β⋆{\bm{S}}_{0}^{\beta\star}. Therefore, the solution extracted using Proposition 7 is also the unique global minimizer of problem (14). ∎

Derivation of Proposition 9

Here we show the intuition for using the basis reduction in Proposition 9. In the original SOS relaxation (22), the parametrization of the SOS polynomial multipliers sk,k=0,…,2Ks_{k},k=0,\dots,2K, and the polynomial multipliers λi,i=1,…,15\lambda_{i},i=1,\dots,15, uses the vector of all monomials up to their corresponding degrees (cf. (A10) and (A12)), which leads to an SDP of size N0=(K+9+ββ)N_{0}=\left(\begin{subarray}{c}K+9+\beta\\ \beta\end{subarray}\right) that grows quadratically with the number of basis shapes KK. In basis reduction, we do not limit ourselves to the vector of full monomials, but rather parametrize s0s_{0}, sks_{k} and λi\lambda_{i} with unknown monomials bases v0[x]v_{0}[\bm{x}], vs[x]v_{s}[\bm{x}] and vλ[x]v_{\lambda}[\bm{x}], which allows us to rewrite (23) as:

with the hope that v0[x]⊆[x]2v_{0}[\bm{x}]\subseteq[\bm{x}]_{2}, vs[x]⊆[x]1v_{s}[\bm{x}]\subseteq[\bm{x}]_{1} and vλ[x]⊆[x]2v_{\lambda}[\bm{x}]\subseteq[\bm{x}]_{2} have much smaller sizes (we limit ourselves to the case of β=2\beta=2, at which level the relaxation is empirically tight).

As described, one can see that the problem of finding smaller v0[x]v_{0}[\bm{x}], vs[x]v_{s}[\bm{x}] and vλ[x]v_{\lambda}[\bm{x}], while keeping the relaxation empirically tight, is highly combinatorial in general. Therefore, our strategy is to only consider the following case:

Expressive: choose v0[x]v_{0}[\bm{x}] such that s0s_{0} contains all the monomials in f(x)−γf(\bm{x})-\gamma,

Balanced: choose vs[x]v_{s}[\bm{x}] and vλ[x]v_{\lambda}[\bm{x}] such that the sum s0+∑skgk+∑λihis_{0}+\sum s_{k}g_{k}+\sum\lambda_{i}h_{i} can only have monomials from f(x)−γf(\bm{x})-\gamma.

In words, condition (i) ensures that the right-hand side (RHS) of (A20) contains all the monomials of the left-hand side (LHS). Condition (ii) asks the three terms of the RHS, i.e., s0s_{0}, ∑skgk\sum s_{k}g_{k} and ∑λihi\sum\lambda_{i}h_{i}, to be self-balanced in the types of monomials. For example, if s0s_{0} contains extra monomials that are not in the LHS, then those extra monomials better appear also in ∑skgk\sum s_{k}g_{k} and/or ∑λihi\sum\lambda_{i}h_{i} so that they could be canceled by summation. Under these two conditions, it is possible to have equation (A20) holdWhether or not these are sufficient or necessary conditions remains open. However, leveraging Theorem 8 we can still check optimality a posteriori..

The choices in both conditions depend on analyzing the monomials in f(x)−γf(\bm{x})-\gamma. Recall the expression of f(x)f(\bm{x}) in (14) and the expression qi(x)q_{i}(\bm{x}) in (18) for each term inside the summation, it can be seen that f(x)f(\bm{x}) only contains the following types of monomials:

and the key observation is that f(x)f(\bm{x}) does not contain degree-4 monomials purely in c\bm{c} or r\bm{r}, i.e., ck1ck2ck3ck4c_{k_{1}}c_{k_{2}}c_{k_{3}}c_{k_{4}} and rj1rj2rj3rj4r_{j_{1}}r_{j_{2}}r_{j_{3}}r_{j_{4}}, or any degree-3 monomials in c\bm{c} and r\bm{r}. Therefore, when choosing v0[x]v_{0}[\bm{x}], we can exclude degree-2 monomials purely in c\bm{c} and r\bm{r} from [x]2[\bm{x}]_{2}, and set v0[x]=m2(x)=[1,cT,rT,cT⊗rT]Tv_{0}[\bm{x}]=m_{2}(\bm{x})=[1,\bm{c}^{\mathsf{T}},\bm{r}^{\mathsf{T}},\bm{c}^{\mathsf{T}}\otimes\bm{r}^{\mathsf{T}}]^{\mathsf{T}}A more rigorous analysis should follow the rules of Newton Polytope , but the intuition is the same as what we describe here. as stated in Proposition 9. This will satisfy the expressive condition (i), because s0=m2[x]TS0m2[x]s_{0}=m_{2}[\bm{x}]^{\mathsf{T}}{\bm{S}}_{0}m_{2}[\bm{x}] can have the following monomials:

and those in (A25) cover the monomials in f(x)f(\bm{x}). Replacing [x]2[\bm{x}]_{2} with v0[x]=m2(x)v_{0}[\bm{x}]=m_{2}(\bm{x}) is the key step in reducing the size of the SDP, because it reduces the size of the SDP from (K+112)\left(\begin{subarray}{c}K+11\\ 2\end{subarray}\right) to 10K+1010K+10, i.e., from quadratic to linear in KK.

In order to satisfy condition (ii), when choosing vs[x]v_{s}[\bm{x}] and vλ[x]v_{\lambda}[\bm{x}], the goal is to have the product between sks_{k}, λi\lambda_{i} and gkg_{k}, hih_{i} result in monomials that appear in f(x)−γf(x)-\gamma, and ensure that monomials that do not appear in the latter can simplify our in the summation. For example, as stated in Proposition 9, we choose vs[x]=[r]1=[1,rT]Tv_{s}[\bm{x}]=[\bm{r}]_{1}=[1,\bm{r}^{\mathsf{T}}]^{\mathsf{T}} and sks_{k} will contain monomials 1, rjr_{j} and rj1rj2r_{j_{1}}r_{j_{2}}. Because gkg_{k}’s have monomials 11, ckc_{k} and ck2c_{k}^{2}, we can see that ∑skgk\sum s_{k}g_{k} will contain the following monomials:

This still satisfies the balanced condition, because monomials of ∑skgk\sum s_{k}g_{k} in (A27) balance with monomials of s0s_{0} in (A25), and monomials of ∑skgk\sum s_{k}g_{k} in (A27) balance with monomials of s0s_{0} in (A26). Similarly, choosing vλ[x]=[c]2v_{\lambda}[\bm{x}]=[\bm{c}]_{2} makes λi\lambda_{i} have monomials 11, ckc_{k} and ck1ck2c_{k_{1}}c_{k_{2}}, and because hih_{i}’s have monomials 11, rjr_{j} and rj1rj2r_{j_{1}}r_{j_{2}}, we see that ∑λihi\sum\lambda_{i}h_{i} contains the following monomials:

which balance with monomials in s0s_{0} from (A25) and (A26).

We remark that we cannot guarantee that the SOS relaxation resulting from basis reduction can achieve the same performance as the original SOS relaxation and we cannot guarantee our choice of basis is “optimal” in any sense. Therefore, in practice, one needs to check the solution and compute corank(S02⋆)\text{corank}({\bm{S}}_{0}^{2\star}) and η2\eta_{2} to check the optimality of the solution produced by (27). Moreover, it remains an open problem to find a better set of monomials bases to achieve better reduction (e.g., knowing more about the algebraic geometry of gkg_{k} and hih_{i} could possibly enable using the standard monomials as a set of bases ).

Derivation of Algorithm 1

For a complete discussion of graduated non-convexity and its applications for robust spatial perception, please see .

In the main document, for robust shape reconstruction, we adopt the TLS shape reconstruction formulation:

where ri(ck,R,t):=∥zi ⁣− ⁣ΠR(∑k=1KckBki) ⁣− ⁣t∥r_{i}(c_{k},{\bm{R}},\bm{t}):=\left\|\bm{z}_{i}\!-\!\Pi{\bm{R}}\left(\sum_{k=1}^{K}c_{k}{\bm{B}}_{ki}\right)\!-\!\bm{t}\right\| is called the residual, and ρcˉ(r)=min⁡(r2,cˉ2)\rho_{\bar{c}}(r)=\min(r^{2},\bar{c}^{2}) implements a truncated least squares cost. Recalling that ρcˉ(r)=min⁡(r2,cˉ2)=min⁡w∈{0,1}wr2+(1−w)cˉ2\rho_{\bar{c}}(r)=\min(r^{2},\bar{c}^{2})=\min_{w\in\{0,1\}}wr^{2}+(1-w)\bar{c}^{2}, we can rewrite the TLS shape reconstruction as a joint optimization of (c,R,t)(\bm{c},{\bm{R}},\bm{t}) and the binary variables wiw_{i}’s, as in eq. (30) in the main document. However, as hinted in the main document, due to the non-convexity of the TLS cost, directly solving the joint problem or alternating between solving for (c,R,t)(\bm{c},{\bm{R}},\bm{t}) and binary variables wiw_{i}’s would require an initial guess and is prone to bad local optima.

The idea of graduated non-convexity (GNC) is to introduce a surrogate function ρcˉμ(r)\rho^{\mu}_{\bar{c}}(r), governed by a control parameter μ\mu, such that changing μ\mu allows ρcˉμ(r)\rho^{\mu}_{\bar{c}}(r) to start from a convex proxy of ρcˉ(r)\rho_{\bar{c}}(r), and gradually increase the amount of non-convexity till the original TLS function ρcˉ(r)\rho_{\bar{c}}(r) is recovered. The surrogate function for TLS is stated below.

The truncated least squares function is defined as:

where cˉ\bar{c} is a given truncation threshold. The GNC surrogate function with control parameter μ\mu is:

By inspection, one can verify ρcˉμ(r)\rho^{\mu}_{\bar{c}}(r) is convex for μ\mu approaching zero ((ρcˉμ(r))′′=−2μ→0\rho^{\mu}_{\bar{c}}(r))^{\prime\prime}=-2\mu\rightarrow 0) and retrieves ρcˉ(r)\rho_{\bar{c}}(r) in (A32) for μ→+∞\mu\rightarrow+\infty. An illustration of ρcˉμ(r)\rho^{\mu}_{\bar{c}}(r) is given in Fig. A1.

The nice property of the GNC surrogate function is that when μ\mu is close to zero, ρcˉμ\rho^{\mu}_{\bar{c}} is convex, which means the only non-convexity of problem (A31) comes from the constraints and can be relaxed using the SOS relaxations.

For the GNC surrogate function ρcˉμ\rho^{\mu}_{\bar{c}}, the simple trick of introducing binary variables (ρcˉ(r)=min⁡w∈{0,1}wr2+(1−w)cˉ2\rho_{\bar{c}}(r)=\min_{w\in\{0,1\}}wr^{2}+(1-w)\bar{c}^{2}) would not work. However, Black and Rangarajan showed that this idea of introducing an outlier variableww can be thought of an outlier variable: when w=1w=1, the measurement is an inlier, when w=1w=1, the measurement is an outlier. can be generalized to many robust cost functions. In particular, for the GNC surrogate function, we have the following.

The GNC surrogate TLS shape reconstruction:

with ρcˉμ(r)\rho^{\mu}_{\bar{c}}(r) defined in (A33), is equivalent to the following optimization with outlier variables wiw_{i}’s:

where Φcˉμ(wi)\Phi^{\mu}_{\bar{c}}(w_{i}) is the following outlier process:

The derivation of Φcˉμ(wi)\Phi^{\mu}_{\bar{c}}(w_{i}) in (A36) follows the Black-Rangarajan procedure in Fig. 10 of . ∎

In words, the Black-Rangarajan duality allows us to rewrite the non-convex shape reconstruction problem as a joint optimization in (c,R,t)(\bm{c},{\bm{R}},\bm{t}) and outlier variables wiw_{i}’s. The interested readers can find closed-form outlier processes for many other robust cost functions in the original paper .

Leveraging the Black-Rangarajan duality, for any given choice of the control parameter μ\mu, we can solve problem (A35) in two steps: first we solve (c,R,t)(\bm{c},{\bm{R}},\bm{t}) using Shape⋆ with fixed weights wiw_{i}’s, and then we update the weights with fixed (c,R,t)(\bm{c},{\bm{R}},\bm{t}). In particular, at each iteration τ\tau (corresponding to a given control parameter μ\mu), we perform the following:

Variable update: minimize (A35) with respect to (c,R,t)(\bm{c},{\bm{R}},\bm{t}), with fixed weights wi(τ−1)w_{i}^{(\tau-1)}:

where we have dropped the term ∑i=1NΦcˉμ(wi)\sum_{i=1}^{N}\Phi^{\mu}_{\bar{c}}(w_{i}) because it is independent from (c,R,t)(\bm{c},{\bm{R}},\bm{t}). This problem is exactly the weighted least squares problem (13) and can be solved using Shape⋆ (cf. line 1 in Algorithm 1). Using the solutions (ck(τ),R(τ),t(τ))(c_{k}^{(\tau)},{\bm{R}}^{(\tau)},\bm{t}^{(\tau)}), we can compute the residuals ri(τ)r_{i}^{(\tau)} (cf. line 1 in Algorithm 1).

Weight update: minimize (A35) with respect to wiw_{i}, with fixed residuals ri(τ)r_{i}^{(\tau)}:

where we have dropped ∑k=1Kck(τ)\sum_{k=1}^{K}c_{k}^{(\tau)} because it is a constant for the optimization. This optimization, fortunately, can be solved in closed-form. We take the gradient of the objective function with respect to wiw_{i}:

and observe that ∇wi=(ri(τ))2−μ+1μcˉ2\nabla_{w_{i}}=(r_{i}^{(\tau)})^{2}-\frac{\mu+1}{\mu}\bar{c}^{2} when wi=0w_{i}=0, and ∇wi=(ri(τ))2−μμ+1cˉ2\nabla_{w_{i}}=(r_{i}^{(\tau)})^{2}-\frac{\mu}{\mu+1}\bar{c}^{2} when wi=1w_{i}=1. Therefore, the global minimizer wi⋆:=wi(τ)w_{i}^{\star}:=w_{i}^{(\tau)} is:

and this is the weight update rule in line 1 of Algorithm 1.

After both the variables and weights are updated using the shaperobust approach described above, we increase the control parameter μ\mu to increase the non-convexity of the surrogate function ρcˉμ\rho^{\mu}_{\bar{c}} (cf. line 1 of Algorithm 1). At the next iteration τ+1\tau+1, the updated weights are used to perform the variable update. The iterations terminate when the change in the objective function becomes negligible (cf. line 1 of Algorithm 1) or after a maximum number of iterations (cf. line 1 of Algorithm 1). Note that all weights are initialized to 11 (cf. line 1 in Algorithm 1), which means that initially all measurements are tentatively accepted as inliers, therefore no prior information about inlier/outlier is required.

FG3DCar Qualitative Results

Fig. 14 shows 9 full qualitative results comparing the performances of Altern+Robust , Convex+Robust and Shape#{\#} on the FG3DCar dataset under 10%10\% to 70%70\% outlier rates. One can further see that the performance of Shape#{\#} is sensitive to 70%70\% outliers, while the performances of Altern+Robust and Convex+Robust gradually degrade and fail at 50%50\% to 60%60\% outliers.

& Convex+Robust Shape#{\#} Honda CRV Altern+Robust Convex+Robust Shape#{\#} Mercedes-Benz GL450

& Convex+Robust Shape#{\#} Volvo V70 Altern+Robust Convex+Robust Shape#{\#} Saab 93