Building Rome with Convex Optimization

用凸优化建造罗马

Haoyu Han and Heng Yang Affiliation:  School of Engineering and Applied Sciences, Harvard University Affiliation:  https://computationalrobotics.seas.harvard.edu/XM

Abstract

摘要

Global bundle adjustment is made easy by depth prediction and convex optimization. We (i) propose a scaled bundle adjustment (SBA) formulation that lifts 2D keypoint measurements to 3D with learned depth, (ii) design an empirically tight convex semidefinite programming (SDP) relaxation that solves SBA to certifiable global optimality, (iii) solve the SDP relaxation at extreme scale with Burer-Monteiro factorization and a CUDA-based trust-region Riemannian optimizer (dubbed XM), (iv) build a structure from motion pipeline with XM as the optimization engine and show that XM-SfM compares favorably with existing pipelines in terms of reconstruction quality while being significantly faster, more scalable, and initialization-free.

全局束调整因深度预测与凸优化而变得容易。我们 (i) 提出缩放束调整（SBA）公式，用学习深度把二维关键点测量提升到三维；(ii) 设计经验上紧的凸半定规划（SDP）松弛，以可证明的全局最优性求解 SBA；(iii) 用 Burer-Monteiro 分解与基于 CUDA 的信赖域黎曼优化器（称为 XM）在极端规模上求解 SDP 松弛；(iv) 以 XM 为优化引擎构建运动恢复结构管线，并表明 XM-SfM 在重建质量上可与现有管线媲美，同时显著更快、更可扩展，且无需初始化。

(a) BAL Dataset. From left to right: $\mathchar 28721\mathchar 28720\mathchar 28721\mathchar 28725\mathchar 28725$, $\mathchar 28721\mathchar 28729\mathchar 28723\mathchar 28724$, and $\mathchar 28723\mathchar 28729\mathchar 28722$ camera frames.

(a) BAL 数据集。从左到右：$\mathchar 28721\mathchar 28720\mathchar 28721\mathchar 28725\mathchar 28725$、$\mathchar 28721\mathchar 28729\mathchar 28723\mathchar 28724$ 与 $\mathchar 28723\mathchar 28729\mathchar 28722$ 个相机帧。

(b) Replica Dataset. All three cases have $\mathchar 28722\mathchar 28720\mathchar 28720\mathchar 28720$ camera frames.

(b) Replica 数据集。三种情形均有 $\mathchar 28722\mathchar 28720\mathchar 28720\mathchar 28720$ 个相机帧。

(c) Mip-Nerf Dataset (left: reconstruction, right: novel view synthesis). $\mathchar 28721\mathchar 28728\mathchar 28725$ camera frames.

(c) Mip-Nerf 数据集（左：重建，右：新视角合成）。$\mathchar 28721\mathchar 28728\mathchar 28725$ 个相机帧。

(d) IMC Dataset. From left to right: $\mathchar 28723\mathchar 28727\mathchar 28726\mathchar 28725$, $\mathchar 28722\mathchar 28720\mathchar 28726\mathchar 28723$, and $\mathchar 28729\mathchar 28720\mathchar 28724$ camera frames.

(d) IMC 数据集。从左到右：$\mathchar 28723\mathchar 28727\mathchar 28726\mathchar 28725$、$\mathchar 28722\mathchar 28720\mathchar 28726\mathchar 28723$ 与 $\mathchar 28729\mathchar 28720\mathchar 28724$ 个相机帧。

Left (e) TUM Dataset ($\mathchar 28727\mathchar 28729\mathchar 28728$ and $\mathchar 28726\mathchar 28721\mathchar 28723$ camera frames); Right (f) C3VD Dataset ($\mathchar 28723\mathchar 28727\mathchar 28720$ and $\mathchar 28726\mathchar 28721\mathchar 28723$ camera frames).

左 (e) TUM 数据集（$\mathchar 28727\mathchar 28729\mathchar 28728$ 与 $\mathchar 28726\mathchar 28721\mathchar 28723$ 个相机帧）；右 (f) C3VD 数据集（$\mathchar 28723\mathchar 28727\mathchar 28720$ 与 $\mathchar 28726\mathchar 28721\mathchar 28723$ 个相机帧）。

Fig. 1: Faster, scalable, and initialization-free 3D reconstruction powered by conveX bundle adjustMent (XM).

图 1：由 conveX bundle adjustMent（XM）驱动的更快、可扩展且无需初始化的三维重建。

I Introduction

I 引言

At the heart of modern structure from motion (SfM) and simultaneous localization and mapping (SLAM) sits bundle adjustment (BA), the procedure of reconstructing camera poses and 3D landmarks from 2D image keypoints.

现代运动恢复结构（SfM）与同时定位与建图（SLAM）的核心是束调整（BA）：从二维图像关键点重建相机位姿与三维路标。

Classical BA formulation. Consider a so-called view graph illustrated in Fig. 2 with two types of nodes: 3D points $\mathchar 29040_{\mathchar 29035}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723}}\mathchar 24891\mathchar 29035\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29005$ and camera poses $\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785\mathchar 12850\mathrm{\mathchar 29011\mathchar 28997}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006$. The set of edges ${\cal\mathchar 28997}$ contains visibility information. An edge $\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}$ indicates the $\mathchar 29035$-th 3D point is visible to the $\mathchar 29033$-th camera, and a 2D keypoint measurement $\mathchar 29045_{\mathchar 29033\mathchar 29035}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28722}}$ has been obtained regarding the projection of the point onto the camera frame.[^1] The bundle adjustment problem consists of estimating points and poses (i.e., $\mathchar 29040_{\mathchar 29035}$’s and $\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785$’s on the nodes) using the 2D keypoint measurements (i.e., $\mathchar 29045_{\mathchar 29033\mathchar 29035}$’s on the edges). This is often formulated as an optimization problem:

经典 BA 公式。考虑图 2 所示的视角图，其中有两类节点：三维点 $\mathchar 29040_{\mathchar 29035}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723}}\mathchar 24891\mathchar 29035\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29005$ 与相机位姿 $\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785\mathchar 12850\mathrm{\mathchar 29011\mathchar 28997}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006$。边集 ${\cal\mathchar 28997}$ 包含可见性信息。边 $\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}$ 表示第 $\mathchar 29035$ 个三维点对第 $\mathchar 29033$ 个相机可见，并已获得关于该点投影到相机帧的二维关键点测量 $\mathchar 29045_{\mathchar 29033\mathchar 29035}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28722}}$。[^1] 束调整即用边上的二维关键点测量（即 $\mathchar 29045_{\mathchar 29033\mathchar 29035}$）估计节点上的点与位姿（即 $\mathchar 29040_{\mathchar 29035}$ 与 $\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785$）。这通常写成优化问题：

$$\displaystyle\min_{\begin{\mathchar29043\mathchar29045\mathchar29026\mathchar29025\mathchar29042\mathchar29042\mathchar29025\mathchar29049}{\mathchar29027}\mathchar 29040_{\mathchar 29035}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723}}\mathchar 24891\mathchar 29035\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29005\
\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785\mathchar 12850\mathrm{\mathchar 29011\mathchar 28997}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006\end{\mathchar29043\mathchar29045\mathchar29026\mathchar29025\mathchar29042\mathchar29042\mathchar29025\mathchar29049}}\mathchar 4944\displaylimits_{\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}}\left\delimiter 0\mathchar 29045_{\mathchar 29033\mathchar 29035}\mathchar 8704\mathchar 28953\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 29040_{\mathchar 29035}\mathchar 8235\mathchar 29044_{\mathchar 29033}\delimiter 84054785\right\delimiter 0^{\mathchar 28722}\mathchar 24891 \tag{1}$$

where $\mathchar 29010_{\mathchar 29033}\mathchar 29040_{\mathchar 29035}\mathchar 8235\mathchar 29044_{\mathchar 29033}$ transforms the point $\mathchar 29040_{\mathchar 29035}$ to the $\mathchar 29033$-th camera frame, $\mathchar 28953\delimiter 67273472\delimiter 69640972\delimiter 84054785\mathchar 12346\mathchar 12349\delimiter 67482370\delimiter 69640972_{\mathchar 28721}\mathchar 24635\delimiter 69640972_{\mathchar 28722}\delimiter 84267779\delimiter 68408078\delimiter 69640972_{\mathchar 28723}$ divides the first two coordinates by the depth and projects the 3D point to 2D, and the sum of squared errors evaluates how well the reprojected 2D points agree with the measurements $\mathchar 29045_{\mathchar 29033\mathchar 29035}$ for all $\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}$. Problem (1) assumes the availability of a view graph, which is often obtained through feature detection and matching in SfM and SLAM pipelines (to be detailed in §IV). An important observation is that problem (1) can only be solved up to scale. This is because $\mathchar 28953\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 29040_{\mathchar 29035}\mathchar 8235\mathchar 29044_{\mathchar 29033}\delimiter 84054785\mathchar 12817\mathchar 28953\delimiter 67273472\mathchar 29010_{\mathchar 29033}\delimiter 67273472\mathchar 28939\mathchar 29040_{\mathchar 29035}\delimiter 84054785\mathchar 8235\delimiter 67273472\mathchar 28939\mathchar 29044_{\mathchar 29033}\delimiter 84054785\delimiter 84054785$ for any scalar $\mathchar 28939\mathchar 12606\mathchar 28720$, i.e., scaling the points and translations by a factor of $\mathchar 28939$ does not change the objective value of problem (1).

其中 $\mathchar 29010_{\mathchar 29033}\mathchar 29040_{\mathchar 29035}\mathchar 8235\mathchar 29044_{\mathchar 29033}$ 把点 $\mathchar 29040_{\mathchar 29035}$ 变换到第 $\mathchar 29033$ 个相机坐标系，$\mathchar 28953\delimiter 67273472\delimiter 69640972\delimiter 84054785\mathchar 12346\mathchar 12349\delimiter 67482370\delimiter 69640972_{\mathchar 28721}\mathchar 24635\delimiter 69640972_{\mathchar 28722}\delimiter 84267779\delimiter 68408078\delimiter 69640972_{\mathchar 28723}$ 用深度去除前两维坐标并把三维点投影到二维，平方误差和则衡量重投影二维点与测量 $\mathchar 29045_{\mathchar 29033\mathchar 29035}$ 对所有 $\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}$。问题 (1) 假定视角图已知，通常由 SfM 与 SLAM 管线中的特征检测与匹配得到（详见 §IV）。一个重要观察是：问题 (1) 只能在尺度意义下求解。这是因为 $\mathchar 28953\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 29040_{\mathchar 29035}\mathchar 8235\mathchar 29044_{\mathchar 29033}\delimiter 84054785\mathchar 12817\mathchar 28953\delimiter 67273472\mathchar 29010_{\mathchar 29033}\delimiter 67273472\mathchar 28939\mathchar 29040_{\mathchar 29035}\delimiter 84054785\mathchar 8235\delimiter 67273472\mathchar 28939\mathchar 29044_{\mathchar 29033}\delimiter 84054785\delimiter 84054785$ 对任意标量 $\mathchar 28939\mathchar 12606\mathchar 28720$，亦即把点与平移同乘 $\mathchar 28939$ 不改变问题 (1) 的目标值。

Optimization challenges. Problem (1) is intuitive to formulate but extremely difficult to optimize. The difficulty comes from two challenges. First, problem (1) is highly nonconvex. The nonconvexity comes from both the nonconvex feasible set $\mathrm{\mathchar 29011\mathchar 28997}\delimiter 67273472\mathchar 28723\delimiter 84054785$ and the nonconvex objective function due to the 2D reprojection function $\mathchar 28953\delimiter 67273472\mathchar 8705\delimiter 84054785$. Second, problem (1) can have an extremely large scale. Both the number of camera poses $\mathchar 29006$ and the number of points $\mathchar 29005$ can range from hundreds to tens of thousands (cf. examples in Fig. 1). Due to these challenges, BA solvers such as Ceres [3], Gtsam [17], and accelerated variants [19, 36] require good initializations and can often get stuck in poor local minima (cf. §V-A). To address this, the popular SfM pipeline Colmap [39] employs an incremental strategy which starts by reconstructing only two views (i.e., (1) with $\mathchar 29006\mathchar 12349\mathchar 28722$) and then sequentially registers additional camera images and associated 3D structure. Incremental SfM ensures stable initialization but makes the pipeline slow and hard to scale to large datasets (cf. §V where Colmap requires several hours runtime). The recent global SfM pipeline Glomap [32] replaces the incremental process with rotation averaging and global positioning, but such initialization can still be time-consuming (see results in §V).

优化挑战。问题 (1) 表述直观但极难优化，困难来自两点。第一，问题 (1) 高度非凸。非凸性既来自非凸可行集 $\mathrm{\mathchar 29011\mathchar 28997}\delimiter 67273472\mathchar 28723\delimiter 84054785$，也来自二维重投影函数导致的非凸目标 $\mathchar 28953\delimiter 67273472\mathchar 8705\delimiter 84054785$。第二，问题 (1) 规模可能极大。相机位姿数 $\mathchar 29006$ 与点数 $\mathchar 29005$ 可从数百到数万（参见图 1 示例）。因此 Ceres [3]、Gtsam [17] 及加速变体 [19, 36] 等 BA 求解器需要良好初始化，并常陷入差局部极小（参见 §V-A）。为此，流行 SfM 管线 Colmap [39] 采用增量策略：先只重建两视图（即 (1) 取 $\mathchar 29006\mathchar 12349\mathchar 28722$）再依次注册更多相机图像与三维结构。增量 SfM 保证稳定初始化，但使管线缓慢且难扩展到大数据集（参见 §V，Colmap 需数小时）。近期全局 SfM 管线 Glomap [32] 用旋转平均与全局定位替换增量过程，但此类初始化仍可能耗时（见 §V 结果）。

Therefore, the motivating question of this paper is:

因此，本文的动机问题是：

Can we design an algorithm that solves the global bundle adjustment problem (1) without initialization and at scale?

我们能否设计一种算法，在无需初始化且可扩展的前提下求解全局束调整问题 (1)？

Fig. 2: A view graph for the bundle adjustment formulation (1). We propose a scaled bundle adjustment formulation (3) by lifting the 2D keypoints to 3D with learned depth.

图 2：束调整公式 (1) 的视角图。我们提出缩放束调整公式 (3)，用学习深度把二维关键点提升到三维。

BA with learned depth. We start by designing an approximate formulation to problem (1) that is simpler to optimize. The key insight is to lift the 2D keypoints $\mathchar 29045_{\mathchar 29033\mathchar 29035}$ as approximate (and noisy) 3D keypoints leveraging off-the-shelf large-scale pretrained monocular depth prediction models, such as [35, 49, 10]. Formally, let $\mathchar 29028_{\mathchar 29033\mathchar 29035}\mathchar 12606\mathchar 28720$ be the predicted depth of the 2D keypoint $\mathchar 29045_{\mathchar 29033\mathchar 29035}$, we generate a 3D keypoint measurement

带学习深度的 BA。我们先为问题 (1) 设计更易优化的近似公式。关键洞察是：借助现成大规模预训练单目深度预测模型（如 [35, 49, 10]），把二维关键点 $\mathchar 29045_{\mathchar 29033\mathchar 29035}$ 提升为近似（且有噪声）的三维关键点。形式化地，令 $\mathchar 29028_{\mathchar 29033\mathchar 29035}\mathchar 12606\mathchar 28720$ 为二维关键点 $\mathchar 29045_{\mathchar 29033\mathchar 29035}$ 的预测深度，我们生成三维关键点测量

$$\displaystyle\tilde{\mathchar 29045}{\mathchar 29033\mathchar 29035}\mathchar 12349\mathchar 29028{\mathchar 29033\mathchar 29035}\begin{\mathchar29026\mathchar29037\mathchar29025\mathchar29044\mathchar29042\mathchar29033\mathchar29048}\mathchar 29045_{\mathchar 29033\mathchar 29035}\
\mathchar 28721\end{\mathchar29026\mathchar29037\mathchar29025\mathchar29044\mathchar29042\mathchar29033\mathchar29048}\mathchar 24891\quad\mathchar 568\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}\mathchar 314 \tag{2}$$

In our work, we select a metric depth model, which is designed to produce depth estimates that match the true depth when the model is perfectly trained. Further, due to imperfect depth prediction caused by potentially out-of-distribution test images, we incorporate a scaling parameter that allows subsequent optimization to correct erroneous predictions, leading to the following scaled bundle adjustment (SBA) formulation

在本文中，我们选用度量深度模型：在完美训练时，其深度估计与真值深度一致。进一步，由于测试图像可能分布外导致深度预测不完美，我们引入缩放参数，使后续优化能纠正错误预测，从而得到如下缩放束调整（SBA）公式

$$\displaystyle\min_{\begin{\mathchar29043\mathchar29045\mathchar29026\mathchar29025\mathchar29042\mathchar29042\mathchar29025\mathchar29049}{\mathchar29027}\mathchar 29040_{\mathchar 29035}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723}}\mathchar 24891\mathchar 29035\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29005\
\mathchar 29043_{\mathchar 29033}\mathchar 12606\mathchar 28720\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006\
\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785\mathchar 12850\mathrm{\mathchar 29011\mathchar 28997}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006\end{\mathchar29043\mathchar29045\mathchar29026\mathchar29025\mathchar29042\mathchar29042\mathchar29025\mathchar29049}}\mathchar 4944\displaylimits_{\delimiter 67273472\mathchar 29033\mathchar 24891\mathchar 29035\delimiter 84054785\mathchar 12850{\cal\mathchar 28997}}\mathchar 29047_{\mathchar 29033\mathchar 29035}\left\delimiter 0\mathchar 29010_{\mathchar 29033}\delimiter 67273472\mathchar 29043_{\mathchar 29033}\tilde{\mathchar 29045}{\mathchar 29033\mathchar 29035}\delimiter 84054785\mathchar 8235\mathchar 29044{\mathchar 29033}\mathchar 8704\mathchar 29040_{\mathchar 29035}\right\delimiter 0^{\mathchar 28722}\mathchar 24891 \tag{3}$$

where $\mathchar 29043_{\mathchar 29033}$ scales the predicted 3D keypoint $\tilde{\mathchar 29045}_{\mathchar 29033\mathchar 29035}$ in (2), $\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785$ transforms the scaled 3D keypoint to the global frame,[^2] and the squared errors in the objective computes 3D distances to the landmarks $\mathchar 29040_{\mathchar 29035}$ without the reprojection $\mathchar 28953\delimiter 67273472\mathchar 8705\delimiter 84054785$. In (3), we allow weighting the squared errors by confidence values ($\mathchar 29047_{\mathchar 29033\mathchar 29035}\mathchar 12606\mathchar 28720$) from the depth prediction model.

其中 $\mathchar 29043_{\mathchar 29033}$ 缩放 (2) 中的预测三维关键点 $\tilde{\mathchar 29045}_{\mathchar 29033\mathchar 29035}$，$\delimiter 67273472\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\delimiter 84054785$ 把缩放后的三维关键点变换到全局坐标系，[^2]目标中的平方误差计算到路标 $\mathchar 29040_{\mathchar 29035}$ 的三维距离，而不经重投影 $\mathchar 28953\delimiter 67273472\mathchar 8705\delimiter 84054785$。在 (3) 中，我们允许用深度预测模型给出的置信度（$\mathchar 29047_{\mathchar 29033\mathchar 29035}\mathchar 12606\mathchar 28720$）加权平方误差。

Remark 1(Connection to Other Perception Problems).

It is worth noting that (i) without the per-frame scaling, problem (3) recovers the multiple point cloud registration problem [15, 26]; (ii) when $\mathchar 29006\mathchar 12349\mathchar 28722$ (two frames), problem (3) reduces to (scaled) point cloud registration [25, 48].

注 1（与其他感知问题的联系）。

值得注意：(i) 去掉逐帧缩放时，问题 (3) 退化为多点云配准问题 [15, 26]；(ii) 当 $\mathchar 29006\mathchar 12349\mathchar 28722$（两帧）时，问题 (3) 化为（缩放）点云配准 [25, 48]。

Through depth prediction, we removed the reprojection function $\mathchar 28953\delimiter 67273472\mathchar 8705\delimiter 84054785$ from (1) and resolved one challenge. However, nonconvexity and large scale remain in problem (3).

通过深度预测，我们从 (1) 中去掉了重投影函数 $\mathchar 28953\delimiter 67273472\mathchar 8705\delimiter 84054785$，解决了其中一个挑战。然而问题 (3) 中非凸性与大规模仍然存在。

Contributions. First, we show the remaining nonconvexity in problem (3) is “benign”. Our strategy is to first rewrite (3) as a quadratically constrained quadratic program (QCQP), and then design a convex semidefinite program (SDP) relaxation, in the spirit of a growing family of SDP-enabled certifiable algorithms [47, 38, 52, 28, 6, 33, 44]. We show the SDP relaxation is empirically tight, i.e., globally optimal solutions of the nonconvex (3) can be computed from the convex SDP relaxation with optimality certificates. Second, we show the convex SDP relaxations can be solved at extreme scale and speed—faster and more scalable than even the best local solvers such as Ceres. The enabling technique is to exploit the low-rankness of (tight) SDP optimal solutions via Burer-Monteiro (BM) factorization [13] and solve the resulting Riemannian optimization using a trust-region algorithm [11]. For the first time, we implemented the trust-region Riemannian optimizer directly in C++/CUDA and show the GPU implementation is up to $\mathchar 28721\mathchar 28720\mathchar 28720$ times faster than the state-of-the-art CPU-based Manopt package [12] and can solve extreme-scale problems beyond the reach of Manopt (e.g., $\mathchar 29006\mathchar 12606\mathchar 28721\mathchar 28720\mathchar 24891\mathchar 28720\mathchar 28720\mathchar 28720$ camera frames, see Fig. 1). We name our GPU solver XM (conveX bundle adjustMent). Third, we build a full SfM pipeline with XM as the optimization engine and various techniques from prior work, such as feature matching from Colmap, view graph creation from Glomap, Ceres refinement (i.e., warmstart Ceres with XM’s solutions for solving (1)), and outlier-robust filtering and estimation schemes [4, 20, 31, 40]. We test XM-SfM across six popular datasets and demonstrate that XM-SfM compares favorably with existing SfM pipelines in terms of reconstruction quality while being significantly faster, more scalable, and initialization-free, thanks to convex optimization and GPU-based implementation.

贡献。首先，我们证明问题 (3) 中剩余的非凸性是「良性」的。策略是先把 (3) 改写为二次约束二次规划（QCQP），再设计凸半定规划（SDP）松弛，呼应日益增多的可证明 SDP 算法族 [47, 38, 52, 28, 6, 33, 44]。我们表明 SDP 松弛在经验上是紧的，亦即非凸 (3) 的全局最优解可从凸 SDP 松弛连同最优性证书一并算出。其次，我们证明凸 SDP 松弛可在极端规模与速度下求解——甚至比 Ceres 等最佳局部求解器更快、更可扩展。关键技术是借 Burer-Monteiro（BM）分解 [13] 利用（紧）SDP 最优解的低秩性，并用信赖域算法 [11] 求解所得黎曼优化。我们首次直接在 C++/CUDA 中实现信赖域黎曼优化器，并表明 GPU 实现相对当前最优的基于 CPU 的 Manopt 包 [12] 可快至 $\mathchar 28721\mathchar 28720\mathchar 28720$ 倍，并能求解 Manopt 力所不及的极端规模问题（例如 $\mathchar 29006\mathchar 12606\mathchar 28721\mathchar 28720\mathchar 24891\mathchar 28720\mathchar 28720\mathchar 28720$ 个相机帧，见图 1）。我们将 GPU 求解器命名为 XM（conveX bundle adjustMent）。第三，我们以 XM 为优化引擎，结合 Colmap 特征匹配、Glomap 视角图构建、Ceres 精修（即以 XM 解为 warmstart 求解 (1)）以及抗野值滤波与估计方案 [4, 20, 31, 40] 等先前技术，构建完整 SfM 管线。我们在六个流行数据集上测试 XM-SfM，表明其在重建质量上可与现有 SfM 管线媲美，同时因凸优化与 GPU 实现而显著更快、更可扩展且无需初始化。

In summary, our contributions are:

总结起来，我们的贡献是：

• • designing an empirically tight convex SDP relaxation for the scaled bundle adjustment problem (3);

• • 为缩放束调整问题 (3) 设计经验上紧的凸 SDP 松弛；

• • solving the convex SDP at extreme scales using BM factorization paired with a trust-region Riemannian optimizer directly implemented in C++/CUDA, i.e., XM;

• • 用 BM 分解配合直接以 C++/CUDA 实现的信赖域黎曼优化器在极端规模上求解凸 SDP，即 XM；

• • creating a full SfM pipeline called XM-SfM that “builds Rome with convex optimization”.

• • 创建完整 SfM 管线 XM-SfM，「用凸优化建造罗马」。

Paper organization. We derive the QCQP formulation for (3) and design the SDP relaxation in §II. We present BM factorization and the CUDA-based trust-region Riemannian optimizer in §III. We describe the XM-SfM pipeline in §IV and present experimental results in §V. We conclude in §VI.

论文组织。§II 推导 (3) 的 QCQP 公式并设计 SDP 松弛。§III 给出 BM 分解与基于 CUDA 的信赖域黎曼优化器。§IV 描述 XM-SfM 管线，§V 给出实验结果。§VI 作结。

II QCQP and Convex SDP Relaxation

II QCQP 与凸 SDP 松弛

In this section, let us focus on solving the SBA problem (3) to certifiable global optimality. We proceed in two steps. In §II-A, we simplify the original formulation as a quadratically constrained quadratic program (QCQP) through a sequence of mathematical manipulations. In §II-B, we apply Shor’s semidefinite relaxation to “convexify” the nonconvex QCQP.

本节聚焦把 SBA 问题 (3) 求解到可证明的全局最优。分两步：§II-A 经一系列数学处理把原公式化简为二次约束二次规划（QCQP）；§II-B 应用 Shor 半定松弛以「凸化」非凸 QCQP。

Before we get started, we remove the ambiguity of problem (3) through anchoring.

开始之前，我们通过锚定消除问题 (3) 的歧义。

Anchoring. Observe that one can choose $\mathchar 29043_{\mathchar 29033}\mathchar 12833\mathchar 28720\mathchar 24891\mathchar 568\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006$, $\mathchar 29044_{\mathchar 28721}\mathchar 12349\dots\mathchar 12349\mathchar 29044_{\mathchar 29006}\mathchar 12349\mathchar 29040_{\mathchar 28721}\mathchar 12349\dots\mathchar 12349\mathchar 29040_{\mathchar 29005}\mathchar 12349{\bf\mathchar 28720}$, and the objective value of (3) can be set arbitrarily close to zero. Additionally, multiplying an arbitary rotation matrix on $\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\mathchar 24891\mathchar 29040_{\mathchar 29033}$ does not change the objective of (3), leading to infinitely many solutions. To resolve these issues, we anchor the first frame and set $\mathchar 29010_{\mathchar 28721}\mathchar 12349{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\mathchar 29044_{\mathchar 28721}\mathchar 12349{\mathbf{\mathchar 28720}}\mathchar 24891\mathchar 29043_{\mathchar 28721}\mathchar 12349\mathchar 28721$.

锚定。观察到可选取 $\mathchar 29043_{\mathchar 29033}\mathchar 12833\mathchar 28720\mathchar 24891\mathchar 568\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006$, $\mathchar 29044_{\mathchar 28721}\mathchar 12349\dots\mathchar 12349\mathchar 29044_{\mathchar 29006}\mathchar 12349\mathchar 29040_{\mathchar 28721}\mathchar 12349\dots\mathchar 12349\mathchar 29040_{\mathchar 29005}\mathchar 12349{\bf\mathchar 28720}$，并使 (3) 的目标值任意接近零。此外，对 $\mathchar 29010_{\mathchar 29033}\mathchar 24891\mathchar 29044_{\mathchar 29033}\mathchar 24891\mathchar 29040_{\mathchar 29033}$ 乘任意旋转矩阵不改变 (3) 的目标，导致无穷多解。为消除这些问题，我们锚定第一帧并令 $\mathchar 29010_{\mathchar 28721}\mathchar 12349{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\mathchar 29044_{\mathchar 28721}\mathchar 12349{\mathbf{\mathchar 28720}}\mathchar 24891\mathchar 29043_{\mathchar 28721}\mathchar 12349\mathchar 28721$.

II-A QCQP Formulation

II-A QCQP 公式

We first show that the unconstrained variables in (3), namely the translations $\mathchar 29044_{\mathchar 29033}$ and the 3D landmarks $\mathchar 29040_{\mathchar 29035}$, can be “marginalized out”, leading to an optimization problem only concerning the scaling factors and 3D rotations.

我们先证明 (3) 中的无约束变量，即平移 $\mathchar 29044_{\mathchar 29033}$ 与三维路标 $\mathchar 29040_{\mathchar 29035}$，可被「边缘化掉」，从而得到只关乎缩放因子与三维旋转的优化问题。

Proposition 2(Scaled-Rotation-Only Formulation).

Problem(3)is equivalent to the following optimization

命题 2（仅缩放–旋转公式）。

问题 (3) 等价于如下优化

\displaystyle\mathchar 28954^{\mathchar 8511}\mathchar 12349\min \\ \displaystyle\ \ \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013\right\delimiter 84054785 \\ \displaystyle\ \ \mathchar 29013\mathchar 12349\left\delimiter 67482370\begin{array}[]{cccc}{\mathbf{\mathchar 29001}}_{\mathchar 28723}&\mathchar 29043_{\mathchar 28722}\mathchar 29010_{\mathchar 28722}&\mathinner{\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}}&\mathchar 29043_{\mathchar 29006}\mathchar 29010_{\mathchar 29006}\end{array}\right\delimiter 84267779 \\ \displaystyle\ \ \mathchar 29043_{\mathchar 29033}\mathchar 12606\mathchar 28720\mathchar 24891\mathchar 29010_{\mathchar 29033}\mathchar 12850\mathrm{\mathchar 29011\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 24891\quad\mathchar 29033\mathchar 12349\mathchar 28722\mathchar 24891\mathinner{\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}}\mathchar 24891\mathchar 29006 \tag{4a}

where $\mathchar 29009\mathchar 12850\mathbb{\mathchar 29011}^{\mathchar 28723\mathchar 29006}$ is a constant and symmetric “data matrix” whose expression is given in Appendix A-A.

其中 $\mathchar 29009\mathchar 12850\mathbb{\mathchar 29011}^{\mathchar 28723\mathchar 29006}$ 是常数对称「数据矩阵」，表达式见附录 A-A。

Let $\mathchar 29013^{\mathchar 8511}$ represent the optimal solution of(4)and let $\mathchar 29044\mathchar 12349\delimiter 67482370\mathchar 29044_{\mathchar 28721}\mathchar 24635\dots\mathchar 24635\mathchar 29044_{\mathchar 29006}\delimiter 84267779\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723\mathchar 29006}}$ be the concatenation of translations, $\mathchar 29040\mathchar 12349\delimiter 67482370\mathchar 29040_{\mathchar 28721}\mathchar 24635\dots\mathchar 24635\mathchar 29040_{\mathchar 29005}\delimiter 84267779\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723\mathchar 29005}}$ be the concatenation of landmark positions. The optimal translations $\mathchar 29044^{\mathchar 8511}$ and landmark positions $\mathchar 29040^{\mathchar 8511}$ of problem (3) can be recovered from $\mathchar 29013^{\mathchar 8511}$ as follows:

令 $\mathchar 29013^{\mathchar 8511}$ 表示 (4) 的最优解，令 $\mathchar 29044\mathchar 12349\delimiter 67482370\mathchar 29044_{\mathchar 28721}\mathchar 24635\dots\mathchar 24635\mathchar 29044_{\mathchar 29006}\delimiter 84267779\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723\mathchar 29006}}$ 为平移的拼接，$\mathchar 29040\mathchar 12349\delimiter 67482370\mathchar 29040_{\mathchar 28721}\mathchar 24635\dots\mathchar 24635\mathchar 29040_{\mathchar 29005}\delimiter 84267779\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723\mathchar 29005}}$ 为路标位置的拼接。问题 (3) 的最优平移 $\mathchar 29044^{\mathchar 8511}$ 与路标位置 $\mathchar 29040^{\mathchar 8511}$ 可由 $\mathchar 29013^{\mathchar 8511}$ 如下恢复：

$$\displaystyle\left\delimiter 67482370\begin{array}[]{c}\mathchar 29044^{\mathchar 8511}\
\mathchar 29040^{\mathchar 8511}\end{array}\right\delimiter 84267779\mathchar 12349\delimiter 67273472\mathchar 28993\mathchar 8714{\mathbf{\mathchar 29001}}_{\mathchar 28723}\delimiter 84054785\mathrm{\delimiter 69640972\mathchar 29029\mathchar 29027}\left\delimiter 67273472\mathchar 29013^{\mathchar 8511}\right\delimiter 84054785\mathchar 314$$

where the expression of $\mathchar 28993$ can be found in Appendix A-A.

其中 $\mathchar 28993$ 的表达式见附录 A-A。

Proof.

See Appendix A-A. ∎

证明。

见附录 A-A。∎

Proposition 2 reformulates the SBA problem (3) as a lower-dimensional problem (4). However, problem (4) is not a QCQP because the objective function is a degree-four polynomial in $\mathchar 29043_{\mathchar 29033}$ and $\mathchar 29010_{\mathchar 29033}$. In the next step, we show that it is possible to combine the “scaled rotation” $\mathchar 29043_{\mathchar 29033}\mathchar 29010_{\mathchar 29033}$ together as a new variable, effectively reducing the degree of the polynomial.

命题 2 把 SBA 问题 (3) 改写为低维问题 (4)。然而问题 (4) 不是 QCQP，因为目标是关于 $\mathchar 29043_{\mathchar 29033}$ 与 $\mathchar 29010_{\mathchar 29033}$ 的四次多项式。下一步我们证明可以把「缩放旋转」 $\mathchar 29043_{\mathchar 29033}\mathchar 29010_{\mathchar 29033}$ 合并为新变量，从而有效降低多项式次数。

Proposition 3(QCQP Formulation).

Define the set of scaled orthogonal group as

命题 3（QCQP 公式）。

定义缩放正交群集合为

$$\displaystyle\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 12349\{\bar{\mathchar 29010}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 28723\mathchar 8706\mathchar 28723}}\mathchar 12906\mathchar 569\mathchar 29043\mathchar 12606\mathchar 28720\mathchar 24891\mathchar 29010\mathchar 12850\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785~\text{s.t.}~\bar{\mathchar 29010}\mathchar 12349\mathchar 29043\mathchar 29010\}\mathchar 314 \tag{8}$$

The set $\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ can be described by quadratic constraints

集合 $\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ 可由二次约束刻画

$$\begin{split}\bar{\mathchar 29010}\mathchar 12349\left\delimiter 67482370\begin{array}[]{ccc}\mathchar 29027_{\mathchar 28721}&\mathchar 29027_{\mathchar 28722}&\mathchar 29027_{\mathchar 28723}\end{array}\right\delimiter 84267779\mathchar 12850\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathrel{\mathchar 552\mkern-3.0mu\mathchar 553}\
\begin{cases}\mathchar 29027_{\mathchar 28721}^{\mathsf{\mathchar 29012}}\mathchar 29027_{\mathchar 28721}\mathchar 12349\mathchar 29027_{\mathchar 28722}^{\mathsf{\mathchar 29012}}\mathchar 29027_{\mathchar 28722}\mathchar 12349\mathchar 29027_{\mathchar 28723}^{\mathsf{\mathchar 29012}}\mathchar 29027_{\mathchar 28723}\
\mathchar 29027_{\mathchar 28721}^{\mathsf{\mathchar 29012}}\mathchar 29027_{\mathchar 28722}\mathchar 12349\mathchar 29027_{\mathchar 28722}^{\mathsf{\mathchar 29012}}\mathchar 29027_{\mathchar 28723}\mathchar 12349\mathchar 29027_{\mathchar 28723}^{\mathsf{\mathchar 29012}}\mathchar 29027_{\mathchar 28721}\mathchar 12349\mathchar 28720\mathchar 314\end{cases}\end{split} \tag{9}$$

Consider the following QCQP

考虑如下 QCQP

\displaystyle\mathchar 28954_{{\sf\mathchar 29009\mathchar 28995\mathchar 29009\mathchar 29008}}^{\mathchar 8511}\mathchar 12349\min \\ \displaystyle\ \ \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013\right\delimiter 84054785 \\ \displaystyle\ \ \mathchar 29013\mathchar 12349\left\delimiter 67482370\begin{array}[]{cccc}{\mathbf{\mathchar 29001}}_{\mathchar 28723}&\bar{\mathchar 29010}_{\mathchar 28722}&\mathinner{\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}}&\bar{\mathchar 29010}_{\mathchar 29006}\end{array}\right\delimiter 84267779 \\ \displaystyle\ \ \bar{\mathchar 29010}_{\mathchar 29033}\mathchar 12850\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 24891~\mathchar 29033\mathchar 12349\mathchar 28722\mathchar 24891\dots\mathchar 24891\mathchar 29006\mathchar 314 \tag{10a}

and let $\mathchar 29013^{\mathchar 8511}\mathchar 12349\delimiter 67482370{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\bar{\mathchar 29010}_{\mathchar 28722}^{\mathchar 8511}\mathchar 24891\dots\mathchar 24891\bar{\mathchar 29010}_{\mathchar 29006}^{\mathchar 8511}\delimiter 84267779$ be a global optimizer. If

并令 $\mathchar 29013^{\mathchar 8511}\mathchar 12349\delimiter 67482370{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\bar{\mathchar 29010}_{\mathchar 28722}^{\mathchar 8511}\mathchar 24891\dots\mathchar 24891\bar{\mathchar 29010}_{\mathchar 29006}^{\mathchar 8511}\delimiter 84267779$ 为全局最优解。若

$$\displaystyle\det{\bar{\mathchar 29010}_{\mathchar 29033}^{\mathchar 8511}}\mathchar 12606\mathchar 28720\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28722\mathchar 24891\dots\mathchar 24891\mathchar 29006\mathchar 24891 \tag{11}$$

then $\mathchar 29013^{\mathchar 8511}$ is a global minimizer to problem(4).

则 $\mathchar 29013^{\mathchar 8511}$ 是问题 (4) 的全局极小点。

Proof.

See Appendix A-B. ∎

证明。

见附录 A-B。∎

It is clear that

显然有

$$\displaystyle\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29009\mathchar 28995\mathchar 29009\mathchar 29008}}\mathchar 12820\mathchar 28954^{\mathchar 8511} \tag{12}$$

because from (4) to (10) we have relaxed the $\mathrm{\mathchar 29011\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ constraint to $\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$. After solving (10), if there exists some index $\mathchar 29033$ such that $\det{\bar{\mathchar 29010}_{\mathchar 29033}^{\mathchar 8511}}\mathchar 12604\mathchar 28720$, we can extract a feasible solution to (4) by projecting onto $\mathrm{\mathchar 29011\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ after extracting the scaling from $\bar{\mathchar 29010}_{\mathchar 29033}^{\mathchar 8511}$.

因为从 (4) 到 (10) 我们把 $\mathrm{\mathchar 29011\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ 约束松弛为 $\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$。求解 (10) 后，若存在某个指标 $\mathchar 29033$ 使得 $\det{\bar{\mathchar 29010}_{\mathchar 29033}^{\mathchar 8511}}\mathchar 12604\mathchar 28720$，则可先从 $\mathrm{\mathchar 29011\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ 提取缩放，再投影到 $\bar{\mathchar 29010}_{\mathchar 29033}^{\mathchar 8511}$，从而得到 (4) 的可行解。

Remark 4(Bad Local Minimum).

Problem (10) is a nonconvex QCQP. In fact, it is a smooth Riemannian optimization by identifying $\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ as a product manifold of the positive manifold and the orthogonal group. Therefore, one can directly use, e.g., Manopt to solve (10) (and also (4)). However, as we will show in §V on the Mip-Nerf 360 dataset [7], directly solving (10) can get stuck in bad local minima.

注 4（坏局部极小）。

问题 (10) 是非凸 QCQP。实际上，通过把 $\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ 视为正流形与正交群的乘积流形，它是光滑黎曼优化。因此可直接用 Manopt 等求解 (10)（以及 (4)）。但如我们将在 §V 于 Mip-Nerf 360 数据集 [7] 上所示，直接求解 (10) 可能陷入坏局部极小。

This motivates and necessitates convex relaxation.

这促使并使得凸松弛成为必要。

II-B Convex SDP Relaxation

II-B 凸 SDP 松弛

For any QCQP, there exists a convex relaxation known as Shor’s semidefinite relaxation [46, Chapter 3]. The basic idea is fairly simple: by creating a matrix variable $\mathchar 29016\mathchar 12346\mathchar 12349\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013$ that is quadratic in the original variable $\mathchar 29013$, problem (10) becomes linear in $\mathchar 29016$. The convex relaxation proceeds by using convex positive semidefinite constraints and linear constraints to properly enforce the matrix variable $\mathchar 29016$.

对任意 QCQP，存在称为 Shor 半定松弛的凸松弛 [46, 第 3 章]。基本想法很简单：构造相对原变量 $\mathchar 29013$ 二次的矩阵变量 $\mathchar 29016\mathchar 12346\mathchar 12349\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013$，则问题 (10) 对 $\mathchar 29016$ 变为线性。凸松弛再用凸正半定约束与线性约束恰当强制矩阵变量 $\mathchar 29016$。

Proposition 5(SDP Relaxation).

The following semidefinite program (SDP)

命题 5（SDP 松弛）。

如下半定规划（SDP）

$$\displaystyle\mathchar 28954^{\mathchar 8511}{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}\mathchar 12349\min{\mathchar 29016\mathchar 12850\mathbb{\mathchar 29011}^{\mathchar 28723\mathchar 29006}} \ \displaystyle\ \ \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29016\right\delimiter 84054785 \ \displaystyle\ \ \mathchar 29016\mathchar 12349\left\delimiter 67482370\begin{array}[]{ccc}{\mathbf{\mathchar 29001}}{\mathchar 28723}&\mathinner{\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}}&\mathchar 8707\
\mathord{\vbox{\kern 6.0pt\hbox{.}\hbox{.}\hbox{.}}}&\mathinner{\mkern 1.0mu\raise 7.0pt\vbox{\kern 7.0pt\hbox{.}}\mkern 2.0mu\raise 4.0pt\hbox{.}\mkern 2.0mu\raise 1.0pt\hbox{.}\mkern 1.0mu}&\mathord{\vbox{\kern 6.0pt\hbox{.}\hbox{.}\hbox{.}}}\
&\mathinner{\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}}&\mathchar 28939{\mathchar 29006}{\mathbf{\mathchar 29001}}_{\mathchar 28723}\end{array}\right\delimiter 84267779\mathchar 12823\mathchar 28720 \tag{13a}$$

is a convex relaxation to(10), i.e.,

是对 (10) 的凸松弛，即

$$\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}\mathchar 12820\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29009\mathchar 28995\mathchar 29009\mathchar 29008}}\mathchar 314 \tag{14}$$

Let $\mathchar 29016^{\mathchar 8511}$ be a global minimizer of(13).

令 $\mathchar 29016^{\mathchar 8511}$ 为 (13) 的全局极小点。

If \mathrm{\mathchar 29042\mathchar 29025\mathchar 29038\mathchar 29035}\left\delimiter 67273472\mathchar 29016^{\mathchar 8511}\right\delimiter 84054785\mathchar 12349\mathchar 28723, then $\mathchar 29016^{\mathchar 8511}$ can be factorized as $\mathchar 29016^{\mathchar 8511}\mathchar 12349\delimiter 67273472\bar{\mathchar 29013}^{\mathchar 8511}\delimiter 84054785^{\mathsf{\mathchar 29012}}\bar{\mathchar 29013}^{\mathchar 8511}$, and it holds[^3]

若 \mathrm{\mathchar 29042\mathchar 29025\mathchar 29038\mathchar 29035}\left\delimiter 67273472\mathchar 29016^{\mathchar 8511}\right\delimiter 84054785\mathchar 12349\mathchar 28723，则 $\mathchar 29016^{\mathchar 8511}$ 可分解为 $\mathchar 29016^{\mathchar 8511}\mathchar 12349\delimiter 67273472\bar{\mathchar 29013}^{\mathchar 8511}\delimiter 84054785^{\mathsf{\mathchar 29012}}\bar{\mathchar 29013}^{\mathchar 8511}$，并且有[^3]

\displaystyle\bar{\mathchar 29013}^{\mathchar 8511}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cccc}\bar{\mathchar 29010}_{\mathchar 28721}^{\mathchar 8511}&\bar{\mathchar 29010}_{\mathchar 28722}^{\mathchar 8511}&\mathinner{\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}\mathpunct{\mathchar 513}}&\bar{\mathchar 29010}_{\mathchar 29006}^{\mathchar 8511}\end{array}\right\delimiter 84267779\mathchar 12850\mathchar 29007\delimiter 67273472\mathchar 28723\delimiter 84054785\mathchar 8706\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785^{\mathchar 29006\mathchar 8704\mathchar 28721}\mathchar 314

Define

定义

$$\displaystyle\mathchar 29013^{\mathchar 8511}\mathchar 12349\delimiter 67273472\bar{\mathchar 29010}_{\mathchar 28721}^{\mathchar 8511}\delimiter 84054785^{\mathsf{\mathchar 29012}}\bar{\mathchar 29013}^{\mathchar 8511}\mathchar 24891 \tag{17}$$

then $\mathchar 29013^{\mathchar 8511}$ is a global optimizer to(10). In this case, we say the relaxation is tight or exact.

则 $\mathchar 29013^{\mathchar 8511}$ 是 (10) 的全局最优解。此时称松弛是紧的或精确的。

Proof.

See Appendix A-C. ∎

证明。

见附录 A-C。∎

If \mathrm{\mathchar 29042\mathchar 29025\mathchar 29038\mathchar 29035}\left\delimiter 67273472\mathchar 29016^{\mathchar 8511}\right\delimiter 84054785\mathchar 12606\mathchar 28723, then we can extract feasible solutions to (10) by taking the top three eigenvectors of $\mathchar 29016^{\mathchar 8511}$ and project the corresponding entries to $\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$ and further to scaled rotations, a step that is typically called rounding. Denote $\hat{\mathchar 29013}$ as the rounded solution that is feasible for (4), we can evaluate the objective of (4) at $\hat{\mathchar 29013}$ and denote it $\hat{\mathchar 28954}$. Combining the chain of inequalities from (12) and (14), we get

若 \mathrm{\mathchar 29042\mathchar 29025\mathchar 29038\mathchar 29035}\left\delimiter 67273472\mathchar 29016^{\mathchar 8511}\right\delimiter 84054785\mathchar 12606\mathchar 28723，则可取 $\mathchar 29016^{\mathchar 8511}$ 的前三个特征向量，并把对应条目投影到 $\mathfrak{\mathchar 29043}\mathrm{\mathchar 29007}\delimiter 67273472\mathchar 28723\delimiter 84054785$，再进一步投影到缩放旋转，这一步通常称为取整。记 $\hat{\mathchar 29013}$ 为对 (4) 可行的取整解，在 $\hat{\mathchar 29013}$ 处计算 (4) 的目标并记为 $\hat{\mathchar 28954}$。合并 (12) 与 (14) 的不等式链，得

$$\displaystyle\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}\mathchar 12820\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29009\mathchar 28995\mathchar 29009\mathchar 29008}}\mathchar 12820\mathchar 28954^{\mathchar 8511}\mathchar 12820\hat{\mathchar 28954}\mathchar 24891 \tag{18}$$

where the last inequality follows from $\hat{\mathchar 29013}$ is a feasible solution for the minimization problem (4). Assuming we can solve the convex SDP, we compute $\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}$ and $\hat{\mathchar 28954}$ at both ends of the inequality (18), allowing us to evaluate a relative suboptimality:

其中最后一个不等式源于 $\hat{\mathchar 29013}$ 是极小化问题 (4) 的可行解。假定能求解凸 SDP，我们在不等式 (18) 两端计算 $\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}$ 与 $\hat{\mathchar 28954}$，从而可评估相对次优性：

$$\displaystyle\mathchar 28945\mathchar 12349{{\hat{\mathchar 28954}\mathchar 8704\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}\over\mathchar 28721\mathchar 8235\delimiter 69640972\hat{\mathchar 28954}\delimiter 69640972\mathchar 8235\delimiter 69640972\mathchar 28954^{\mathchar 8511}_{{\sf\mathchar 29011\mathchar 28996\mathchar 29008}}\delimiter 69640972}}\mathchar 314 \tag{19}$$

$\mathchar 28945\mathchar 12833\mathchar 28720$ certifies global optimality of the rounded solution $\hat{\mathchar 29013}$, and tightness of the SDP relaxation.

$\mathchar 28945\mathchar 12833\mathchar 28720$ 证明取整解 $\hat{\mathchar 29013}$ 的全局最优性，以及 SDP 松弛的紧性。

Remark 6(Connection to Prior Work).

For readers familiar with SDP relaxations, this section should not be surprising at all—that is the reason why we kept this section very brief and only focus on milestone results listed in the propositions. The closest two works to ours areSE-Sync*[38] andSIM-Sync[52], and the SDP relaxation technique traces back to at least the work by Carlone et al. [14]. The novelty of our formulation are twofold: (a) we estimate scaling factors with the motivation to correct learned depth whileSE-Syncestimates rotations and translations only; (b) we jointly estimate 3D landmarks and (scaled) camera poses whileSIM-Syncdoes not estimate 3D landmarks. These differences allow us to solve the long-standing bundle adjustment problem with convex optimization.*

注 6（与先前工作的联系）。

对熟悉 SDP 松弛的读者，本节并不意外——这也是我们写得很短、只聚焦命题中里程碑结果的原因。与我们最接近的两篇工作是 SE-Sync [38] 与 SIM-Sync [52]，SDP 松弛技术至少可追溯到 Carlone 等人 [14]。我们公式的新意有两点：(a) 我们估计缩放因子以纠正学习深度，而 SE-Sync 只估计旋转与平移；(b) 我们联合估计三维路标与（缩放）相机位姿，而 SIM-Sync 不估计三维路标。这些差异使我们能用凸优化求解长期存在的束调整问题。

We shall focus on how to solve the convex SDP (13).

我们将聚焦如何求解凸 SDP (13)。

III Burer-Monteiro Factorization and CUDA-based Riemannian Optimizer

III Burer-Monteiro 分解与基于 CUDA 的黎曼优化器

Let us first write the SDP (13) in standard primal form:

先把 SDP (13) 写成标准原形式：

\displaystyle\min_{\mathchar 29016\mathchar 12850\mathbb{\mathchar 29011}^{\mathchar 28723\mathchar 29006}} \\ \displaystyle\ \ \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29016\right\delimiter 84054785 \\ \displaystyle\ \ \delimiter 69632778\mathchar 28993_{\mathchar 29033}\mathchar 24891\mathchar 29016\delimiter 86414091\mathchar 12349\mathchar 29026_{\mathchar 29033}\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29037 \\ \displaystyle\ \ \mathchar 29016\mathchar 12823\mathchar 28720 \tag{20a}

where we have rewritten the constraint (13) as a positive semidefinite (PSD) constraint (20c) and $\mathchar 29037\mathchar 12349\mathchar 28725\mathchar 29006\mathchar 8235\mathchar 28721$ linear equality constraints (20b). To see why this reformulation is true, note that the diagonal $\mathchar 28723\mathchar 8706\mathchar 28723$ blocks of $\mathchar 29016$ are either ${\mathbf{\mathchar 29001}}_{\mathchar 28723}$ or scaled ${\mathbf{\mathchar 29001}}_{\mathchar 28723}$—the former can be enforced using 6 linear equalities and the latter can be enforced using 5, summing up to $\mathchar 28725\mathchar 29006\mathchar 8235\mathchar 28721$ linear equalities. Associated with the primal standard SDP (20) is the following dual standard SDP:

其中我们把约束 (13) 改写为正半定（PSD）约束 (20c) 以及 $\mathchar 29037\mathchar 12349\mathchar 28725\mathchar 29006\mathchar 8235\mathchar 28721$ 个线性等式约束 (20b)。该改写成立是因为 $\mathchar 28723\mathchar 8706\mathchar 28723$ 的对角 $\mathchar 29016$ 块要么是 ${\mathbf{\mathchar 29001}}_{\mathchar 28723}$，要么是缩放的 ${\mathbf{\mathchar 29001}}_{\mathchar 28723}$——前者可用 6 个线性等式强制，后者用 5 个，合计 $\mathchar 28725\mathchar 29006\mathchar 8235\mathchar 28721$ 个线性等式。与原标准 SDP (20) 相伴的是如下对偶标准 SDP：

$$\displaystyle\max_{\mathchar 29049\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 29037}}} \\ \displaystyle\ \ \mathchar 4944\displaylimits_{\mathchar 29033\mathchar 12349\mathchar 28721}^{\mathchar 29037}\mathchar 29026_{\mathchar 29033}\mathchar 29049_{\mathchar 29033} \\ \displaystyle\ \ \mathchar 29018\delimiter 67273472\mathchar 29049\delimiter 84054785\mathchar 12346\mathchar 12349\mathchar 29009\mathchar 8704\mathchar 4944\displaylimits_{\mathchar 29033\mathchar 12349\mathchar 28721}^{\mathchar 29037}\mathchar 29049_{\mathchar 29033}\mathchar 28993_{\mathchar 29033}\mathchar 12823\mathchar 28720\mathchar 314 \tag{21a}$$

Since Slater’s condition holds ($\mathchar 29016\mathchar 12349{\mathbf{\mathchar 29001}}_{\mathchar 28723\mathchar 29006}\mathchar 12831\mathchar 28720$ is feasible for the primal (20)), we know strong duality holds between primal (20) and dual (21), i.e., their optimal values are equal to each other [46, 45].

由于 Slater 条件成立（$\mathchar 29016\mathchar 12349{\mathbf{\mathchar 29001}}_{\mathchar 28723\mathchar 29006}\mathchar 12831\mathchar 28720$ 对原问题 (20) 可行），故原问题 (20) 与对偶 (21) 强对偶成立，即最优值相等 [46, 45]。

Scalability. Once the SDP (13) is formulated in standard forms, it can be solved by off-the-shelf SDP solvers such as MOSEK [5], provided that the SDP’s scale is not so large. Our SDP relaxation leads to a matrix variable with size $\mathchar 28723\mathchar 29006\mathchar 8706\mathchar 28723\mathchar 29006$. This means that when $\mathchar 29006$ is in the order of hundreds, MOSEK can solve the SDP without any problem. However, when $\mathchar 29006$ is in the order of thousands or tens of thousands, which is not uncommon in bundle adjustment problems, MOSEK will become very slow or even runs out of memory (e.g., MOSEK becomes unresponsive when $\mathchar 29006\mathchar 12606\mathchar 28722\mathchar 28720\mathchar 28720\mathchar 28720$). Therefore, we decided to customize a solver for our SDP relaxation.

可扩展性。一旦 SDP (13) 写成标准形式，只要规模不太大，就可用 MOSEK [5] 等现成 SDP 求解器求解。我们的 SDP 松弛导致矩阵变量规模为 $\mathchar 28723\mathchar 29006\mathchar 8706\mathchar 28723\mathchar 29006$。这意味着当 $\mathchar 29006$ 为数百量级时，MOSEK 可顺利求解。但当 $\mathchar 29006$ 为数千乃至数万（在束调整中并不少见）时，MOSEK 会非常慢甚至内存耗尽（例如当 $\mathchar 29006\mathchar 12606\mathchar 28722\mathchar 28720\mathchar 28720\mathchar 28720$ 时 MOSEK 无响应）。因此我们决定为该 SDP 松弛定制求解器。

III-A Burer-Monteiro Factorization

III-A Burer-Monteiro 分解

The key structure we will leverage is, as stated in Proposition 5, the optimal solution $\mathchar 29016^{\mathchar 8511}$ has its rank equal to three when the SDP relaxation is tight. In other words, the effective dimension of $\mathchar 29016^{\mathchar 8511}$ can be $\mathchar 28723\mathchar 8706\mathchar 28723\mathchar 29006$ instead of $\mathchar 28723\mathchar 29006\mathchar 8706\mathchar 28723\mathchar 29006$.

我们利用的关键结构是：如命题 5 所述，当 SDP 松弛紧时，最优解 $\mathchar 29016^{\mathchar 8511}$ 的秩等于三。换言之，$\mathchar 29016^{\mathchar 8511}$ 的有效维数可以是 $\mathchar 28723\mathchar 8706\mathchar 28723\mathchar 29006$ 而非 $\mathchar 28723\mathchar 29006\mathchar 8706\mathchar 28723\mathchar 29006$。

To exploit the low-rank structure, we will leverage the famous Burer-Monteiro (BM) factorization [13].

为利用低秩结构，我们将借助著名的 Burer-Monteiro（BM）分解 [13]。

Proposition 7(BM Factorization).

For a fixed rank $\mathchar 29042\mathchar 12821\mathchar 28723$, the Burer-Monteiro factorization of(13)and (20) reads:

命题 7（BM 分解）。

对固定秩 $\mathchar 29042\mathchar 12821\mathchar 28723$，(13) 与 (20) 的 Burer-Monteiro 分解为：

$$\displaystyle\mathchar 28954^{\mathchar 8511}{{\sf\mathchar 28994\mathchar 29005}\mathchar 24891\mathchar 29042}\mathchar 12349\min{\bar{\mathchar 29010}{\mathchar 29033}\mathchar 12850{{\mathbb{\mathchar 29010}}^{\mathchar 29042\mathchar 8706\mathchar 28723}}\mathchar 24891\mathchar 29033\mathchar 12349\mathchar 28721\mathchar 24891\dots\mathchar 24891\mathchar 29006} \ \displaystyle\ \ \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013\right\delimiter 84054785 \ \displaystyle\ \ \mathchar 29013\mathchar 12349\left\delimiter 67482370\begin{array}[]{ccc}\bar{\mathchar 29010}{\mathchar 28721}&\dots&\bar{\mathchar 29010}{\mathchar 29006}\end{array}\right\delimiter 84267779 \ \displaystyle\ \ \bar{\mathchar 29010}{\mathchar 29033}^{\mathsf{\mathchar 29012}}\bar{\mathchar 29010}{\mathchar 29033}\mathchar 12349\begin{cases}{\mathbf{\mathchar 29001}}{\mathchar 28723}&\mathchar 29033\mathchar 12349\mathchar 28721\
\mathchar 28939_{\mathchar 29033}{\mathbf{\mathchar 29001}}_{\mathchar 28723}&\mathchar 29033\mathchar 12821\mathchar 28722\end{cases}\mathchar 314 \tag{22a}$$

A few comments are in order. First, by factorizing $\mathchar 29016\mathchar 12349\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013$, $\mathchar 29016\mathchar 12823\mathchar 28720$ holds by construction. Second, note that when $\mathchar 29042\mathchar 12349\mathchar 28723$, problem (22) is exactly the same as the original QCQP (10) (up to the difference in $\bar{\mathchar 29010}_{\mathchar 28721}$), and thus is NOT convex. Third, as long as $\mathchar 29042\mathchar 12821\mathchar 29042^{\mathchar 8511}$ where $\mathchar 29042^{\mathchar 8511}$ is the minimum rank of all optimal solutions of the SDP (13), the nonconvex BM factorization has the same global minimum as the convex SDP (13) [46, 13]. Note that the factor $\mathchar 29013$ is a matrix of size $\mathchar 29042\mathchar 8706\mathchar 28723\mathchar 29006$, i.e., a flat matrix as shown in Fig. 3 bottom right.

有几点评注。第一，通过分解 $\mathchar 29016\mathchar 12349\mathchar 29013^{\mathsf{\mathchar 29012}}\mathchar 29013$，正半定 $\mathchar 29016\mathchar 12823\mathchar 28720$ 自然成立。第二，当 $\mathchar 29042\mathchar 12349\mathchar 28723$ 时，问题 (22) 与原 QCQP (10) 实质相同（仅 $\bar{\mathchar 29010}_{\mathchar 28721}$ 有差异），因而并非凸的。第三，只要 $\mathchar 29042\mathchar 12821\mathchar 29042^{\mathchar 8511}$（$\mathchar 29042^{\mathchar 8511}$ 为 SDP (13) 所有最优解的最小秩），非凸 BM 分解与凸 SDP (13) 有相同的全局最小值 [46, 13]。注意因子 $\mathchar 29013$ 是规模为 $\mathchar 29042\mathchar 8706\mathchar 28723\mathchar 29006$ 的扁平矩阵，如图 3 右下所示。

Counterintuitive? The reader might find this confusing. In §II, we applied a series of techniques to relax a nonconvex optimization problem into a convex SDP, enabling us to solve it to global optimality. Surprisingly, the BM factorization appears to “undo” this effort, pulling us back into nonconvex optimization. While the low-rank factorization offers a clear scalability advantage, it remains uncertain whether this benefit outweighs the drawbacks of reintroducing nonconvexity.

反直觉？ 读者或许困惑：§II 中我们用一系列技术把非凸问题松弛成凸 SDP，以便求到全局最优；BM 分解却似乎「撤销」了这些努力，把我们拉回非凸优化。虽然低秩分解有明确的可扩展优势，但这是否足以抵消重新引入非凸性的代价，仍不显然。

Rank “staircase”. The secret ingredient of BM factorization to tackle nonconvexity is that we will solve the factorized problem (22) at increasing ranks, like stepping up a “staircase” (cf. Fig. 3 bottom right). We will start with the lowest rank $\mathchar 29042\mathchar 12349\mathchar 28723$ and solve problem (22) using local optimization. One of two cases will happen. (a) The local optimizer is the same as the global optimizer of the SDP, in which case we can leverage the dual SDP (21) to certify global optimality and declare victory against the SDP. (b) The local optimizer is not the same as the global optimizer of the SDP, in which case we can again leverage the dual SDP (21) to “escape” the bad local minimum via increasing the rank, i.e., going up the staircase.

秩「阶梯」。BM 分解对付非凸性的秘诀是：我们按递增的秩求解分解问题 (22)，如同登上「阶梯」（参见图 3 右下）。从最低秩 $\mathchar 29042\mathchar 12349\mathchar 28723$ 出发，用局部优化求解 (22)。将出现两种情形之一。(a) 局部最优与 SDP 全局最优相同，则可借助对偶 SDP (21) 证明全局最优，宣告相对 SDP 获胜。(b) 局部最优不是 SDP 全局最优，则同样可借助对偶 SDP (21)，通过升秩（登上阶梯）「逃离」坏局部极小。

We formalize this in Algorithm 1 and prove its correctness.

我们把这形式化为算法 1 并证明其正确性。

1:  Input: the data matrix $\mathchar 29009$ 2:  Output: optimal solution $\mathchar 29013^{\mathchar 8511}$ 3:  # Initialization 4:  Set $\mathchar 29042\mathchar 12349\mathchar 28723$, $\mathchar 29013_{\mathchar 29042}^{\mathchar 28720}\mathchar 12349\delimiter 67482370{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\dots\mathchar 24891{\mathbf{\mathchar 29001}}_{\mathchar 28723}\delimiter 84267779$ 5:  while True do 6:   # Local optimization of (22) 7:   $\mathchar 29013_{\mathchar 29042}^{\mathchar 8511}\mathchar 12349\textsc{LocalOptimizer}\delimiter 67273472\mathchar 29013_{\mathchar 29042}^{\mathchar 28720}\delimiter 84054785$ 8:   # Compute dual certificate 9:   $\mathchar 29049_{\mathchar 29042}\mathchar 24891\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\mathchar 12832\textsc{Solve}~\eqref{eq:bm_first_order_optimality}$ 1 10:   # Certify global optimality 2 11:   if $\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\mathchar 12823\mathchar 28720$ then 12:    return $\mathchar 29013_{\mathchar 29042}$ 13:   # Escape local minimum 14:   $\delimiter 69640972\mathchar 12832\textsc{LeastEigenvector}\delimiter 67273472\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\delimiter 84054785$ 15:   \mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cc}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}&\mathchar 28939\delimiter 69640972\end{array}\right\delimiter 84267779^{\mathsf{\mathchar 29012}} with $\mathchar 28939\mathchar 12349\mathchar 28721$ 16:   # Line search 17:   while \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}^{\mathsf{\mathchar 29012}}\right\delimiter 84054785\mathchar 12821\mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013_{\mathchar 29042}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}\right\delimiter 84054785 do 18:    $\mathchar 28939\mathchar 12349\mathchar 28939\delimiter 68408078\mathchar 28722$ 19:    \mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cc}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}&\mathchar 28939\delimiter 69640972\end{array}\right\delimiter 84267779^{\mathsf{\mathchar 29012}} 20:   $\mathchar 29042\mathchar 12832\mathchar 29042\mathchar 8235\mathchar 28721$, $\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}^{\mathchar 28720}\mathchar 12832\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}$
1:  Input: the data matrix $\mathchar 29009$
2:  Output: optimal solution $\mathchar 29013^{\mathchar 8511}$
3:  # Initialization
4:  Set $\mathchar 29042\mathchar 12349\mathchar 28723$, $\mathchar 29013_{\mathchar 29042}^{\mathchar 28720}\mathchar 12349\delimiter 67482370{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\dots\mathchar 24891{\mathbf{\mathchar 29001}}_{\mathchar 28723}\delimiter 84267779$
5:  while True do
6:   # Local optimization of (22)
7:   $\mathchar 29013_{\mathchar 29042}^{\mathchar 8511}\mathchar 12349\textsc{LocalOptimizer}\delimiter 67273472\mathchar 29013_{\mathchar 29042}^{\mathchar 28720}\delimiter 84054785$
8:   # Compute dual certificate
9:   $\mathchar 29049_{\mathchar 29042}\mathchar 24891\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\mathchar 12832\textsc{Solve}~\eqref{eq:bm_first_order_optimality}$ 1

10:   # Certify global optimality 2

11:   if $\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\mathchar 12823\mathchar 28720$ then
12:    return $\mathchar 29013_{\mathchar 29042}$
13:   # Escape local minimum
14:   $\delimiter 69640972\mathchar 12832\textsc{LeastEigenvector}\delimiter 67273472\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\delimiter 84054785$
15:   \mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cc}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}&\mathchar 28939\delimiter 69640972\end{array}\right\delimiter 84267779^{\mathsf{\mathchar 29012}} with $\mathchar 28939\mathchar 12349\mathchar 28721$
16:   # Line search
17:   while \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}^{\mathsf{\mathchar 29012}}\right\delimiter 84054785\mathchar 12821\mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013_{\mathchar 29042}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}\right\delimiter 84054785 do
18:    $\mathchar 28939\mathchar 12349\mathchar 28939\delimiter 68408078\mathchar 28722$
19:    \mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cc}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}&\mathchar 28939\delimiter 69640972\end{array}\right\delimiter 84267779^{\mathsf{\mathchar 29012}}
20:   $\mathchar 29042\mathchar 12832\mathchar 29042\mathchar 8235\mathchar 28721$, $\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}^{\mathchar 28720}\mathchar 12832\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}$

1:  Input: the data matrix $\mathchar 29009$
2:  Output: optimal solution $\mathchar 29013^{\mathchar 8511}$
3:  # Initialization
4:  Set $\mathchar 29042\mathchar 12349\mathchar 28723$, $\mathchar 29013_{\mathchar 29042}^{\mathchar 28720}\mathchar 12349\delimiter 67482370{\mathbf{\mathchar 29001}}_{\mathchar 28723}\mathchar 24891\dots\mathchar 24891{\mathbf{\mathchar 29001}}_{\mathchar 28723}\delimiter 84267779$
5:  while True do
6:   # Local optimization of (22)
7:   $\mathchar 29013_{\mathchar 29042}^{\mathchar 8511}\mathchar 12349\textsc{LocalOptimizer}\delimiter 67273472\mathchar 29013_{\mathchar 29042}^{\mathchar 28720}\delimiter 84054785$
8:   # Compute dual certificate
9:   $\mathchar 29049_{\mathchar 29042}\mathchar 24891\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\mathchar 12832\textsc{Solve}~\eqref{eq:bm_first_order_optimality}$ 1

10:   # Certify global optimality 2

11:   if $\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\mathchar 12823\mathchar 28720$ then
12:    return $\mathchar 29013_{\mathchar 29042}$
13:   # Escape local minimum
14:   $\delimiter 69640972\mathchar 12832\textsc{LeastEigenvector}\delimiter 67273472\mathchar 29018\delimiter 67273472\mathchar 29049_{\mathchar 29042}\delimiter 84054785\delimiter 84054785$
15:   \mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cc}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}&\mathchar 28939\delimiter 69640972\end{array}\right\delimiter 84267779^{\mathsf{\mathchar 29012}} with $\mathchar 28939\mathchar 12349\mathchar 28721$
16:   # Line search
17:   while \mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}^{\mathsf{\mathchar 29012}}\right\delimiter 84054785\mathchar 12821\mathrm{\mathchar 29044\mathchar 29042}\left\delimiter 67273472\mathchar 29009\mathchar 29013_{\mathchar 29042}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}\right\delimiter 84054785 do
18:    $\mathchar 28939\mathchar 12349\mathchar 28939\delimiter 68408078\mathchar 28722$
19:    \mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}\mathchar 12349\left\delimiter 67482370\begin{array}[]{cc}\mathchar 29013_{\mathchar 29042}^{\mathsf{\mathchar 29012}}&\mathchar 28939\delimiter 69640972\end{array}\right\delimiter 84267779^{\mathsf{\mathchar 29012}}
20:   $\mathchar 29042\mathchar 12832\mathchar 29042\mathchar 8235\mathchar 28721$, $\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}^{\mathchar 28720}\mathchar 12832\mathchar 29013_{\mathchar 29042\mathchar 8235\mathchar 28721}$

Algorithm 1 Riemannian Staircase