11institutetext: Megvii Inc., China 22institutetext: University of Edinburgh,
United Kingdom

# MegBA: A GPU-Based Distributed Library for Large-Scale Bundle Adjustment

Jie Ren∗ 1122 0000-0003-3332-2092 Wenteng Liang∗ 11 0000-0002-3750-7695
Ran Yan† 11 0000-0001-6705-1644 Luo Mai 22 0000-0002-3594-1092
Shiwen Liu 11 0000-0003-1779-409X Xiao Liu 11 0000-0002-6073-030X

###### Abstract

Large-scale Bundle Adjustment (BA) requires massive memory and computation
resources which are difficult to be fulfilled by existing BA libraries. In
this paper, we propose MegBA, a GPU-based distributed BA library. MegBA can
provide massive aggregated memory by automatically partitioning large BA
problems, and assigning the solvers of sub-problems to parallel nodes. The
parallel solvers adopt distributed Precondition Conjugate Gradient and
distributed Schur Elimination, so that an effective solution, which can match
the precision of those computed by a single node, can be efficiently computed.
To accelerate BA computation, we implement end-to-end BA computation using
high-performance primitives available on commodity GPUs. MegBA exposes easy-
to-use APIs that are compatible with existing popular BA libraries.
Experiments show that MegBA can significantly outperform state-of-the-art BA
libraries: Ceres (41.45$\times$), RootBA (64.576$\times$) and DeepLM
(6.769$\times$) in several large-scale BA benchmarks. The code of MegBA is
available at: https://github.com/MegviiRobot/MegBA.

††∗ Equal contribution, work was done during their internship in Megvii Inc.
† Corresponding author.

##  1 Introduction

Bundle Adjustment (BA) is the foundation for many real-world 3D vision
applications [34, 20], including structure-from-motion and simultaneous-
localization-and-mapping. A BA problem minimises the re-projection error
between camera poses and map points. The error is a non-linear square
function, and it is minimised through iterative methods, such as Gauss-Newton
(GN) [37], Leverberg-Marquardt (LM) [26] and Dog-Leg [30]. In each iteration,
a BA library differentiates the errors with respect to solving states and
constructs a linear system which is solved by optimisation algorithms, such as
Cholesky decomposition [4] and Precondition Conjugate Gradient (PCG) [9].

Large-scale BA is increasingly important given the recent rise of city-level
high-definition maps for autonomous driving [3, 21, 24, 22] and indoor maps
for augmentation reality [41, 36, 29]. A structure-from-motion application,
for example, can produce massive images [15, 2], resulting in billions of
points and observations to be adjusted. Such a BA problem is orders of
magnitude larger than those in conventional vision applications [42, 35].

Existing BA libraries (e.g., g2o [12] and Ceres [1]) however provide
insufficient support for large-scale BA. We observe several reasons: (i)
Existing libraries focus on single-node execution, and they lack algorithms to
distribute computation. They thus cannot provide massive aggregated memory
that is the key for large-scale BA. Even though there are algorithms, such as
RPBA [27], DPBA [6] and STBA [42], which explore distributed BA. These
algorithms adopt _approximation_ which can adversely affect the precision of
found solutions. (ii) Existing BA libraries are designed for CPU
architectures, and they under-utilise GPUs which is particularly useful for
large-scale BA. Even though there are systems, such as PBA [39], to accelerate
BA with GPUs. They leave key BA operations un-accelerated (e.g., error
differentiation and linear system construction). DeepLM [16] offloads error
differentiation into GPUs through PyTorch, but the performance is often sub-
optimised.

In this paper, we propose MegBA, a novel GPU-optimised distributed library for
large-scale BA. The design of MegBA makes several contributions:

(1) Distributed BA algorithms. MegBA provides a large amount of aggregated
memory by distributing BA computation to multiple nodes. To this end, we
propose a generic BA problem partitioning method. This method leverages a key
observation in BA problems: BA problems are often expressed as graphs where
nodes represent points/cameras, and edges represent the associations between
cameras and points. MegBA can thus automatically partition the graphs based on
edges, and ensure each sub-graph has an equal number of edges (with an aim of
achieving load balancing). MegBA further assigns sub-graphs to distributed
nodes and merges the local solutions to sub-graphs. To ensure that distributed
BA can offer the precision as those computed by single-node BA libraries, we
propose the distributed PCG algorithm and the distributed Schur elimination
algorithm. These two algorithms synchronise the states of solvers on parallel
nodes, and the synchronisation is realised using NCCL.

(2) GPU-Optimised BA computation. MegBA thoroughly optimise BA computation for
GPUs, thus providing massive computation power for large-scale BA.
Computation-intensive operators (e.g., inverse, inner project, etc) are
implemented as Single-Instruction-Multiple-Data (SIMD) operators. MegBA store
data in _JetVector_ , a data structure that stores BA data in SIMD-friendly
vectors, and JetVector minimises data serialisation cost between CPUs and
GPUs. To minimise data movement cost which could block GPU execution, MegBA
has algorithms that can predict the GPU memory usage of BA, thus pre-fetching
BA data if possible. It exposes easy-to-use APIs that are compatible with g2o
and Ceres. Ceres and g2o applications can be thus easily ported to MegBA.

We evaluate the performance of MegBA on servers and each server has 8 NVIDIA
V100 GPUs. Experiments with public large BA datasets (i.e., Final-13682 [2])
show that MegBA can out-perform Ceres by up to 41.45$\times$, and RootBA [7]
by up to 64.576$\times$, indicating the benefits of optimising BA computation
for GPUs. We further compare MegBA with DeepLM [16], a GPU-based BA library.
MegBA out-performs DeepLM by 5.213$\times$ on 4 GPUs. With 8 GPUs, MegBA out-
performs DeepLM by 6.769$\times$, making MegBA the state-of-the-art BA library
on GPUs.

To evaluate the scalability of MegBA, we construct an extremely large
synthetic BA dataset which is modelled after by the BA problems we have in
real-world applications. This dataset contains 80 million observations,
2.76$\times$ larger than Final-13682. DeepLM and RootBA incur out-of-memory
error and cannot handle such a dataset. On the contrary, MegBA can solve this
BA problem in 216.26 seconds by distributing BA computation to 8 GPUs, which
is 23.54$\times$ faster than Ceres.

##  2 Related Work

This section describes the related work of MegBA. g2o [12] and Ceres [1] are
exact BA libraries that can compute high-accuracy solutions to BA problems.
These libraries are designed for using parallel CPU cores, and they cannot use
GPUs. These libraries also fail to provide distributed execution, which makes
them suffer from out-of-memory issues in solving large-scale BA problems.

Approximated BA algorithms can substantially speed up BA, though often come
with a compromise in the quality of BA solutions. PBA [39] is limited to run
on a single device. $\sqrt{BA}$ [8] replaced Schur Complement with a memory-
efficient nullspace projection of Jacobian, thus improving its performance
with single-precision float numbers. iSAM[18] and iSAM2[17] exploit states
ordering; while ICE-BA [23] exploits the states in temporal orders. Though
fast in speed, approximated BA algorithms modify the original BA problems,
which adversely affect the quality of BA solutions. As a result, commercial 3D
vision software, such as PIX4D1 usually avoid any form of approximation and
adopt exact BA libraries if possible. ††1https://www.pix4d.com/

Distributed BA libraries have been recently designed for large-scale BA.
Anders Eriksson et al. [10] present consensus-based optimisation which
leverages proximal splitting. Runze Zhang et al. [40] purpose an Alternating
Direction Method of Multipliers to distribute the optimisation problem. Later
RPBA [27], DPBA [6], STBA [42] partition the BA problems based on ADMM. These
ADMM-based systems incur massive redundant computation on distributed devices,
making them sometimes under-perform single-node libraries. Further, their
users must manually partition BA problems, resulting in sub-optimal
distributed performance. BA-Net [33] and DeepLM [16] leverages GPUs to speed
up BA. They however rely on PyTorch to use GPU, which incurs non-trivial
performance overheads when using GPU and extra memory copies. Decentralised
SLAM libraries, such as DEDV-SLAM [5], often solve approximated BA problems on
distributed robots, then they merge local solutions. However, the merged
solution is not equivalent to the original global BA problem.

Custom hardware and algorithms are useful in accelerating BA [14]. GBP [28]
uses a neural processing unit (i.e., GraphCore IPU) to speed up BA; but the
limited availability of IPU makes GBP difficult to be used as a general
solution. Practitioners also propose an approximated BA solver tailored for
facial capture [11], and this solver cannot be used for arbitrary BA problems
such as structure-from-motion.

##  3 Preliminaries

Figure 1: MegBA overview. 1 MegBA partitions a BA problem based on edges. BA
sub-problems are in the same size, and they are dispatched to distributed
GPUs. Each GPU computes Jacobians 2 , constructs a linear system 3 , and solve
the linear system using the the distributed PCG algorithm 4 . The
communication involved in linear system construction and distributed PCG is
implemented using allreduce operations. Step 2 , 3 , and 4 are executed
iteratively until 5 convergence criteria has been met.

This section introduces the preliminaries of MegBA. A BA problem can be
expressed as a graph, and its solving is realised an iterative process which
minimises a non-linear square error objective function:

$\boldsymbol{x}^{*}=\arg\min_{\boldsymbol{x}}\sum\boldsymbol{e}_{i,j}^{\top}\boldsymbol{\Sigma}_{i,j}\boldsymbol{e}_{i,j},$
(1)

where $\boldsymbol{e}_{i,j}$ is the constraint (i.e. error or graph edge)
between state (i.e. parameters or graph nodes) $\boldsymbol{x}_{i}$ and
$\boldsymbol{x}_{j}$, $\boldsymbol{\Sigma}_{i,j}$ is an information matrix.

Solving Equation 1 is equivalent to iteratively updating the incremental
amount $\boldsymbol{\Delta x}$, given by the linear system
$\boldsymbol{H}\Delta\boldsymbol{x}=\boldsymbol{g}$, upon the current state
$\boldsymbol{x}$ until convergence. The Hessian matrix
$\boldsymbol{H}=\boldsymbol{J}^{T}\Sigma\boldsymbol{J}$ for GN method and
$\boldsymbol{H}=\boldsymbol{J}^{T}\Sigma\boldsymbol{J}+\lambda\boldsymbol{I}$
for LM method, the residual vector $\boldsymbol{g}$ equals to
$-\boldsymbol{J}^{T}\Sigma\boldsymbol{r}$, $\boldsymbol{J}$ is the Jacobian of
the error $\boldsymbol{e}$ with respect to current state $\boldsymbol{x}$.

To solve BA problems, BA libraries can use Schur Complement (SC):

$\begin{bmatrix}\boldsymbol{B}&\boldsymbol{E}\\\
\boldsymbol{E}^{T}&\boldsymbol{C}\end{bmatrix}\begin{bmatrix}\Delta\boldsymbol{x}_{c}\\\
\Delta\boldsymbol{x}_{p}\end{bmatrix}=\begin{bmatrix}\boldsymbol{v}\\\
\boldsymbol{w}\end{bmatrix}$ (2)

where $\boldsymbol{B}$ and $\boldsymbol{C}$ are block diagonal and they are
related to camera-camera and point-point edges, respectively. $\boldsymbol{E}$
refers to edges between camera and point. $\boldsymbol{v}$ and
$\boldsymbol{w}$ refer to the residual vectors for camera and point states.

Solving $\boldsymbol{H}\Delta\boldsymbol{x}=\boldsymbol{g}$ is equivalent to
compute the incremental update for states related to cameras
$\Delta\boldsymbol{x}_{c}$ by solving an alternative linear system, called
Reduced Camera System (RCS)

$[\boldsymbol{B}-\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{E}^{T}]\Delta\boldsymbol{x}_{c}=\boldsymbol{v}-\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{w},$
(3)

and followed by a substitution $\Delta\boldsymbol{x}_{c}$ into

$\Delta\boldsymbol{x}_{p}=\boldsymbol{C}^{-1}\left(\boldsymbol{w}-\boldsymbol{E}^{T}\Delta\boldsymbol{x}_{c}\right),$
(4)

to get the update for 3D map points.

BA libraries solve linear systems using either direct methods or iterative
methods. Direct methods, such as Gaussian-Elimination, LU, QR, and Cholesky
Decomposition, return optimised solution of $\boldsymbol{x}$ in one pass. They
however suffer from $O(n^{3})$ time and $O(n^{2})$ space complexity, making
them only suitable for small-scale BA problems. On the contrary, iterative
methods, such as PCG [34], are suitable for large-scale BA problems.
Specifically, PCG replaces the explicit computation of
$\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{E}^{T}$ with multiple iterative
sparse matrix-vector operations. It reduces the space complexity to $O(n)$,
thus saving memory.

##  4 MegBA Design

This section introduces the design of MegBA. A key design goal of MegBA is to
transparently distribute the solving of a given BA problem to multiple nodes,
thus addressing the memory wall of a single node.

Figure 1 presents an overview of the distributed execution of MegBA. A MegBA
user declares a BA problem as a graph. MegBA can automatically partition the
BA problem based on edges with an aim of each BA sub-problem to have an even
number of edges 1 . Specifically, each GPU first 2 computes the Jacobian (i.e.
differentiation of the edge for the node), and then 3 construct the linear
system, and finally, 4 apply PCG to compute the update for adjusting the
current BA sub-problem. The PCG intermediate state is synchronised so that
MegBA can eventually solve the shared global BA problem. The BA update step is
iteratively performed until a user-defined convergence criterion is met 5 .
Notably, the BA computation on each GPU is implemented as SIMD operations
which can best utilise GPUs (The details of SIMD-friendly data structure and
the memory-efficient runtime are given in Section 5)

###  4.1 Edge-Based Partitioning Method

We focus on partitioning the Hessian matrix produced in BA. For example, in a
BA dataset with 80M edges, a Hessian matrix $\boldsymbol{H}$ consume over 50G
memory, leading to over 99.9% storage to be allocated to $\boldsymbol{E}$ and
$\boldsymbol{E}^{T}$. We want to have a generic method that partitions the
Hessian matrix in a BA problem. This method needs to assign each parallel
device with a part of the matrix $\boldsymbol{E}$ and $\boldsymbol{E}^{T}$,
preferably in equal sizes. The partitioning needs to guarantee that the global
solution merged from local solutions is _equivalent_ to the one computed using
a single node. This equivalence property is the key to ensuring a high-
precision solution found by MegBA.

At the high level, MegBA achieves distributed BA using two major components:
(i) a method that can divide a BA problem into sub-problems, and (ii) an
algorithm that can coordinate distributed PCG algorithms to solve sub-problems
in parallel. Our partitioning method is based on a key observation in BA
problems: the non-zero blocks in $\boldsymbol{E}$ and $\boldsymbol{E}^{T}$ are
corresponding with edges, i.e., the $i$-th row $j$-th column non-zero block in
$\boldsymbol{E}$ is computed by $e_{i,j}$, we can partition edges based on the
number of available GPUs, and each GPU only store part of these non-zero
blocks (We provide an example to illustrate the partitioning process in the
supplementary materials).

Assume there are $K$ GPUs, given a BA problem, we tile edges to a vector
$\boldsymbol{e}=\begin{array}[]{ccc}[\dots&e_{i,j}&\dots]^{T}\end{array}$,
then we partition it to several blocks
$\boldsymbol{e}=\begin{array}[]{cccc}[\boldsymbol{e}_{1}^{T}&\boldsymbol{e}_{2}^{T}&\dots&\boldsymbol{e}_{K}^{T}]^{T}\end{array}$.
The Jacobian $\boldsymbol{J}$ could be partitioned into several blocks:

$\boldsymbol{J}=\frac{d\boldsymbol{e}}{d\boldsymbol{x}}=\left[\begin{array}[]{cccc}\frac{d\boldsymbol{e}_{1}}{d\boldsymbol{x}}^{T}&\frac{d\boldsymbol{e}_{2}}{d\boldsymbol{x}}^{T}&\dots&\frac{d\boldsymbol{e}_{K}}{d\boldsymbol{x}}^{T}\end{array}\right]^{T}=\left[\begin{array}[]{cccc}\boldsymbol{J}_{1}^{T}&\boldsymbol{J}_{2}^{T}&\dots&\boldsymbol{J}_{K}^{T}\end{array}\right]^{T}.$
(5)

Assuming identity information matrix is given here, Hessian $\boldsymbol{H}$
can be partitioned:

$\displaystyle\boldsymbol{H}$ $\displaystyle=\boldsymbol{J}^{T}\boldsymbol{J}$
$\displaystyle=\sum_{k=1}^{K}\boldsymbol{J}_{k}^{T}\boldsymbol{J}_{k}=\sum_{k=1}^{K}\boldsymbol{H}_{k}.$
(6)

To perform Schur elimination, we represented $\boldsymbol{H}_{k}$ as sub-
blocks:

$\boldsymbol{H}_{k}=\left[\begin{array}[]{cc}\boldsymbol{B}_{k}&\boldsymbol{E}_{k}\\\
\boldsymbol{E}_{k}^{T}&\boldsymbol{C}_{k}\end{array}\right].$ (7)

Matrix blocks in Equation 2 have the following equivalent forms in the edge-
based partition setting: $\boldsymbol{B}=\sum_{k=1}^{K}\boldsymbol{B}_{k}$,
$\boldsymbol{E}=\sum_{k=1}^{K}\boldsymbol{E}_{k}$,
$\boldsymbol{E}^{T}=\sum_{k=1}^{K}\boldsymbol{E}_{k}^{T}$, and
$\boldsymbol{C}=\sum_{k=1}^{K}\boldsymbol{C}_{k}$. The number of non-zero
parameter blocks in $\boldsymbol{E}$ or $\boldsymbol{E}^{T}$ equals the number
of edges. Notably, the sub-matrices $\boldsymbol{B}_{k}$ and
$\boldsymbol{C}_{k}$ have the same number of non-zero elements as
$\boldsymbol{B}$ and $\boldsymbol{C}$, respectively. Since we store matrices
in the Compressed Sparse Row (CSR) format, each GPU only stores $\frac{1}{K}$
non-zero blocks in $\boldsymbol{E}$ and $\boldsymbol{E}^{T}$. The blocking
strategy greatly alleviates the problem that $\boldsymbol{E}$ and
$\boldsymbol{E}^{T}$ are too large to be stored on a single device.

By applying the above partition method for Equation 2, an equivalent
distributed version can be formulated as follow:

$\boldsymbol{g}=-\boldsymbol{J}^{T}\boldsymbol{r}=-[\begin{array}[]{cccc}\boldsymbol{J}_{1}^{T}&\boldsymbol{J}_{2}^{T}&\dots\boldsymbol{J}_{K}^{T}\end{array}][\begin{array}[]{cccc}\boldsymbol{r}_{1}^{T}&\boldsymbol{r}_{2}^{T}&\dots&\boldsymbol{r}_{K}^{T}\end{array}]^{T}=-\sum_{k=1}^{K}\boldsymbol{J}_{k}^{T}\boldsymbol{r}.$
(8)

###  4.2 Distributed BA Algorithm

Algorithm 1 Distributed BA

0: BA initial state
$\boldsymbol{x}=[\begin{array}[]{cc}\boldsymbol{x}_{c}^{T}&\boldsymbol{x}_{p}^{T}\end{array}]^{T}$,
vector of edges $\boldsymbol{e}_{k}$, and local GPU rank $k$

0: Optimised state $\boldsymbol{x}$

1: while BA Convergence Criteria not satisfied do

2:
$\boldsymbol{r}_{k}=\boldsymbol{e}_{k}(\boldsymbol{x}),\boldsymbol{J}_{k}=d\boldsymbol{e}_{k}(\boldsymbol{x})/d\boldsymbol{x}$
/* Residual and Jacobian */

3: $\left[\begin{array}[]{cc}\boldsymbol{B}_{k}&\boldsymbol{E}_{k}\\\
\boldsymbol{E}_{k}^{T}&\boldsymbol{C}_{k}\end{array}\right]=\boldsymbol{J}_{k}^{T}\boldsymbol{J}_{k},[\begin{array}[]{cc}\boldsymbol{v}_{k}&\boldsymbol{w}_{k}\end{array}]=-\boldsymbol{J}_{k}^{T}\boldsymbol{r}_{k}$/*
Hessian and Constant vector */

4: $\boldsymbol{B}=allreduce(\boldsymbol{B}_{k})$,
$\boldsymbol{C}=allreduce(\boldsymbol{C}_{k})$,$\boldsymbol{v}=allreduce(\boldsymbol{v}_{k})$,
$\boldsymbol{w}=allreduce(\boldsymbol{w}_{k})$/*
$\boldsymbol{B}=\sum_{i=1}^{K}\boldsymbol{B}_{i}$,
$\boldsymbol{C}=\sum_{i=1}^{K}\boldsymbol{C}_{i}$,
$\boldsymbol{v}=\sum_{i=1}^{K}\boldsymbol{v}_{i}$,
$\boldsymbol{w}=\sum_{i=1}^{K}\boldsymbol{w}_{i}$ */

5:
$\boldsymbol{\alpha}_{k}=\boldsymbol{E}_{k}\boldsymbol{C}^{-1}\boldsymbol{w}$

6: $\boldsymbol{\alpha}=allreduce(\boldsymbol{\alpha}_{k})$ /*
$\sum_{i=1}^{K}\boldsymbol{\alpha}_{i}$ */

7: $\boldsymbol{g}=\boldsymbol{v}-\boldsymbol{\alpha}$/* Constant vector in
Equation 3 */

8:
$\Delta\boldsymbol{x}_{c}=\texttt{DPCG}(\boldsymbol{0},\boldsymbol{B},\boldsymbol{E}_{k},\boldsymbol{E}_{k}^{T},\boldsymbol{C},\boldsymbol{g},k)$/*
Update $\boldsymbol{x}_{c}$ using Algorithm 2 */

9: $\boldsymbol{\beta}_{k}=\boldsymbol{E}_{k}^{T}\Delta\boldsymbol{x}_{c}$

10: $\boldsymbol{\beta}=allreduce(\boldsymbol{\beta}_{k})$ /*
$\sum_{i=1}^{K}\boldsymbol{\beta}_{i}$ */

11:
$\Delta\boldsymbol{x}_{p}=\boldsymbol{C}^{-1}(\boldsymbol{w}-\boldsymbol{\beta})$
/* Increment of $\boldsymbol{x}_{p}$ */

12: $\boldsymbol{x}_{c}=\boldsymbol{x}_{c}+\Delta\boldsymbol{x}_{c}$,
$\boldsymbol{x}_{p}=\boldsymbol{x}_{p}+\Delta\boldsymbol{x}_{p}$/* Update
state */

13: end while

14: return
$\boldsymbol{x}=[\begin{array}[]{cc}\boldsymbol{x}_{c}^{T}&\boldsymbol{x}_{p}^{T}\end{array}]^{T}$

By far we have partitioned a BA problem and assigned sub-problems to all GPUs.
In the following, we will discuss how does MegBA coordinates the solving of
sub-problems in a distributed manner.

Algorithm 1 introduces the overall distributed BA algorithm in MegBA. The
distributed BA algorithm takes as initial state and partitioned edges as
described in Section 8. We use _JecVector_ to compute the Jacobian and
residual (Line 2). _JetVector_ is a novel data structure to represent BA data
in a SIMD format, it can make full use of the hardware characteristics of GPU
(e.g. coalesced memory loading) to do auto-differentiation over millions of
edges in parallel. We give a more detailed illustration in Section 5.1. Then
we build a linear system (Line 3). We perform allreduce [25] on diagonal-
blocks and constant vector (Line 4) before solving the linear system because
the size of diagonal-blocks and constant vector is small and they would be
used several times in the following procedures.

We then compute constant vector in Equation 3 (Line 5-7) and solve the linear
system by using a Distributed PCG (DPCG) algorithm (Line 8). Notably, we do
necessary allreduce in the DPCG algorithm to guarantee DPCG output the same
result as non-distributed PCG solver does in solving Equation 3, further
implementation details will be shown in Section 4.3. After solving the linear
system in Equation 3, we compute the increment of $\boldsymbol{x}_{p}$
following Equation 4 (Line 9-11). Once we have computed the increment
$\Delta\boldsymbol{x}_{c}$ and $\Delta\boldsymbol{x}_{p}$, we update the state
$\boldsymbol{x}_{c}$ and $\boldsymbol{x}_{p}$ (Line 12). If it doesn’t satisfy
the convergence criteria we will start another loop; otherwise, we will return
the optimised state $\boldsymbol{x}$.

Algorithm 2 Distributed PCG (DPCG)

0: Initial state $\boldsymbol{x}^{0}$, matrix block $\boldsymbol{B}$,
$\boldsymbol{E}_{k}$, $\boldsymbol{E}_{k}^{T}$, $\boldsymbol{C}$ of
$\boldsymbol{H}_{k}$, constant vector $\boldsymbol{b}$, and local GPU rank $k$

0: Solution $\boldsymbol{x}$ for linear system
[$\boldsymbol{B}-\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{E}^{T}]\boldsymbol{x}=\boldsymbol{b}$,
where $\boldsymbol{E}=\sum_{i=1}^{K}\boldsymbol{E}_{i}$ and
$\boldsymbol{E}^{T}=\sum_{i=1}^{K}\boldsymbol{E}_{i}$

1:
$\boldsymbol{r}^{0}=\boldsymbol{b}-\texttt{DSE}(\boldsymbol{x}^{0},\boldsymbol{B},\boldsymbol{E}_{k},\boldsymbol{E}_{k}^{T},\boldsymbol{C}^{-1},k)$
/* Algorithm 3 */

2: $n=0$

3: while Convergence Criteria not satisfied do

4: $\boldsymbol{z}^{n}=\boldsymbol{B}^{-1}\boldsymbol{r}^{n}$

5: $\rho^{n}={\boldsymbol{r}^{n}}^{T}\boldsymbol{z}^{n}$

6: if $n>1$ then

7: $\beta^{n}=\rho^{n}/\rho^{n-1}$

8: $\boldsymbol{p}^{n}=\boldsymbol{z}^{n}+\beta^{n}\boldsymbol{p}^{n}$

9: else

10: $\boldsymbol{p}^{n}=\boldsymbol{z}^{n}$

11: end if

12:
$\boldsymbol{q}^{n}=\texttt{DSE}(\boldsymbol{p}^{n},\boldsymbol{B},\boldsymbol{E}_{k},\boldsymbol{E}_{k}^{T},\boldsymbol{C}^{-1},k)$
/* Algorithm 3 */

13: $\alpha^{n}=\rho^{n}/{\boldsymbol{p}^{n}}^{T}\boldsymbol{q}^{n}$

14: $\boldsymbol{x}^{n+1}=\boldsymbol{x}^{n}+\alpha^{n}\boldsymbol{p}^{n}$

15: $\boldsymbol{r}^{n+1}=\boldsymbol{r}^{n}-\alpha^{n}\boldsymbol{q}^{n-1}$

16: $n=n+1$

17: end while

18: return $\boldsymbol{x}^{n}$

Algorithm 3 Distributed Schur Elimination (DSE)

0: Vector $\boldsymbol{x}$, matrix
$\boldsymbol{B},\boldsymbol{E}_{k},\boldsymbol{E}_{k}^{T},\boldsymbol{C}^{-1}$,
and local GPU rank $k$

0: Schur elimination result
$[\boldsymbol{B}-\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{E}]\boldsymbol{x}$,
where
$\boldsymbol{E}=\sum_{i=1}^{K}\boldsymbol{E}_{i},\boldsymbol{E}^{T}=\sum_{i=1}^{K}\boldsymbol{E}_{i}^{T}$

1: $\boldsymbol{a}_{k}=\boldsymbol{E}_{k}^{T}\boldsymbol{x}$

2: $\boldsymbol{a}=allreduce(\boldsymbol{a}_{k})$ /*
$\sum_{i=1}^{K}\boldsymbol{a}_{i}$ */

3: $\boldsymbol{b}=\boldsymbol{C}^{-1}\boldsymbol{a}$

4: $\boldsymbol{c}_{k}=\boldsymbol{E}_{k}\boldsymbol{b}$

5: $\boldsymbol{c}=allreduce(\boldsymbol{c}_{k})$ /*
$\sum_{i=1}^{K}\boldsymbol{c}_{i}$ */

6: $\boldsymbol{d}=\boldsymbol{B}\boldsymbol{x}$

7: return $\boldsymbol{d}-\boldsymbol{c}$

###  4.3 Distributed PCG

We then discuss how to distribute the PCG algorithm in BA, shown in Algorithm
1. This algorithm first constructs a linear system defined in Equation 2. It
then solves Equation 3 and computes increment following Equation 4. It finally
uses the increments update state $\boldsymbol{x}$, and tested if a convergence
criterion has been met. To guarantee that the distributed BA algorithm
achieves the convergence performance, we make Algorithm 1, named DPCG, return
the same result as the non-distributed linear solver.

In the following, we describe the execution of DPCG. DPCG takes BA initial
state $\boldsymbol{x}^{0}$, matrix block $\boldsymbol{B}$,
$\boldsymbol{E}_{k}$, $\boldsymbol{E}_{k}^{T}$, $\boldsymbol{C}$ of
$\boldsymbol{H}_{k}$, constant vector $\boldsymbol{b}$ as input and output
solution $\boldsymbol{x}$ for linear system
[$\boldsymbol{B}-\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{E}^{T}]\boldsymbol{x}=\boldsymbol{b}$,
where $\boldsymbol{E}=\sum_{k=1}^{K}\boldsymbol{E}_{k}$ and
$\boldsymbol{E}^{T}=\sum_{k=1}^{K}\boldsymbol{E}_{k}$. The procedures of DPCG
using Schur elimination is similar to single-node PCG. Notably, the
coefficient matrix of the linear system to be solved is Schur complement. The
matrix-vector multiplication operations (Line 1, 12 in Algorithm 2) is thus
the multiplication between Schur complement and vector. The difference between
distributed compared with non-distributed setting is that DPCG only assign
$\boldsymbol{E}_{k}$ and $\boldsymbol{E}_{k}^{T}$ rather than the complete
matrices $\boldsymbol{E}$ and $\boldsymbol{E}^{T}$ to GPU $k$, so we need to
guarantee operations that use $\boldsymbol{E}_{k}$ and
$\boldsymbol{E}_{k}^{T}$ have the same output compared with using
$\boldsymbol{E}$ and $\boldsymbol{E}^{T}$, these operations happen when doing
Schur elimination (Line 1, 12).

Our key idea of computing Schur elimination in a distributed manner is that:
the summation of matrix-vector multiplication is the same as the the result
matrix summation multiplies vector, i.e.,
$\sum_{k=1}^{K}(\boldsymbol{E}_{k}\boldsymbol{v})=\sum_{k=1}^{K}(\boldsymbol{E}_{k})\boldsymbol{v}$.
We compute an intermediate vector (Line 1) and reduce it (Line 2), then we
compute intermediate vectors sequentially (Line 3, 4). We perform all-reduce
operation over the intermediate vector (Line 5) and compute another
intermediate vector (Line 6. After those procedures, we do subtraction to the
last two intermediate vectors and output the final result. The result would be
the same as computing the complete Schur complement
$[\boldsymbol{B}-\boldsymbol{E}\boldsymbol{C}^{-1}\boldsymbol{E}]$ then
multiplying it with vector $\boldsymbol{x}$.

###  4.4 Complexity Analysis

In the end, we present the complexity analysis of MegBA. Assume that MegBA is
given $m$ cameras, $n$ points, and $k$ observations and we often have $k\gg
m,n$, the time complexity for building the linear system is
$\mathcal{O}(m+n+k)$ and the time complexity for each iteration of the
conjugate gradient is $\mathcal{O}(m+n+k)$. Assume we distribute the problem
to $K$ GPUs, on each GPU, the time complexity for building the linear system
is $\mathcal{O}(m+n+k/K)$ and the time complexity for each iteration of the
conjugate gradient is $\mathcal{O}(m+n+k/K)$. The ring all-reduce
communication time complexity of each conjugate gradient iteration is
$\mathcal{O}(m+n)$. In summary, the total complexity of MegBA is
$\mathcal{O}(m+n+k/K)$.

##  5 MegBA Implementation

This section describes the implementation of MegBA. There are several goals of
our implementation: (i) We want to use as many SIMD operations as possible
because both computation and memory operations on GPU are essentially SIMD-
friendly. (ii) We want to optimise the memory efficiency of MegBA, thus
avoiding memory allocation and deallocation; (iii) We want to implement the
APIs of MegBA that are fully compatible with existing popular BA libraries:
g2o and Ceres. In the following, we highlight how MegBA achieves these
implementation goals.

###  5.1 SIMD-Friendly Data Structures

_JetVector_ is a novel data structure to perform auto-differentiation over
millions of edges. Compared to conventional BA data structure: _Jet_
implemented in Ceres, _JetVector_ represents a list of Jets (i.e., Array-of-
Structure) as a single data object where Jet’s data fields: data and grad
across all items are represented as single arrays (i.e., Array-within-
Structure). When we perform mathematical operations on _JetVector_ , we will
start as many GPU threads as the elements in it, every GPU thread process one
element. Because data and grad are stored in the structure of Array-within-
Structure, the memory transactions are coalesced and make it easy to reach a
high memory throughput. The detailed structure layout of JecVector could be
found in supplementary materials.

Besides _JetVector_ , other parts of MegBA are also implemented as SIMD-
friendly data structures. The construction of linear system (Line 3 in
Algorithm 1) uses L1 cache on GPU to store Jacobian blocks in a SIMD manner.
The DPCG algorithm includes a lot of matrix/vector operations which also be
benefited from the SIMD structure. A full list of SIMD operations implemented
in MegBA can be found in supplementary materials.

###  5.2 Memory-Efficient Runtime

BA computation involves massive objects to be allocated in GPU memory. To
avoid expensive memory allocation [19], we leverage a key observation in BA
computation: The automatic differentiation works on GPU buffers that are in
the same size across BA iterations. By monitoring the sizes of GPU buffers
used in the forward pass of differentiating the BA errors, we can predict the
sizes of all memory buffers involved in future BA iterations. Based on this
observation, we can pre-allocate these memory buffers in a memory pool, thus
avoiding calling the CUDA driver to allocate memory during runtime.

###  5.3 Easy-to-use APIs

The APIs of MegBA comprises of two major components:

(i) Declaring BA problems. Following the API convention of g2o and Ceres, a BA
problem in MegBA is declared a graph that contains nodes and edges. The MegBA
nodes describe the 3D coordinates or the poses of cameras and these nodes can
be directly imported from g2o and Ceres applications. The MegBA edges are
error functions that can be written using the Eigen library [13], identical to
Ceres. A MegBA user can build a large BA problem by adding BA nodes and edges
(using the g2o-equivalent _addEdge_ and _addNode_ functions).

(ii) Choosing BA solvers. MegBA also support users to choose solvers given the
characteristics of their BA problems. The default solver is the SIMD-optimised
DPCG which can automatically use multiple GPUs. Given a small-scale BA problem
where intrinsic parallelism is not sufficient, MegBA provides users with the
CPU-optimised CHOLMOD solver [4].

##  6 Experimental Evaluation

Benchmark | Dataset | #Points | #Observations
---|---|---|---
BAL | Trafalgar-257 | 65132 | 225911
| Ladybug-1723 | 156502 | 678718
| Dubrovnik-356 | 226730 | 1255268
| Venice-1778 | 993923 | 5001946
| Final-13682 | 4456117 | 28987644
Synthesised | Synthesised-20000 | 80000 | 80000000
1DSfM | Alamo-577 | 140080 | 816891
| Ellis_Island-233 | 9210 | 20500
| Gendarmenmarkt-704 | 76964 | 268747
| Madrid_Metropolis-347 | 44479 | 195660
| Montreal_Notre_Dame-459 | 151876 | 811757
| Notre_Dame-548 | 224153 | 1172145
| NYC_Library-334 | 54757 | 211614
| Piazza_del_Popolo-336 | 29731 | 150161
| Piccadilly-2292 | 184475 | 798085
| Roman_Forum-1067 | 223844 | 1031760
| Tower_of_London-484 | 126648 | 596690
| Trafalgar-5052 | 327920 | 1266102
| Union_Square-816 | 26430 | 90668
| Vienna_Cathedral-846 | 154394 | 495940
| Yorkminster-429 | 100426 | 376980
Table 1: Dataset Statistics.

We conduct a comprehensive evaluation with MegBA. The evaluation comprises of
BAL [2], 1DSfM [38], and a large synthetic dataset modelled after a city-scale
BA application we have in production. The dataset statistic is shown in Table
1. Due to the page limit, this section only presents the results with BAL [2],
and we put the results of 1DSfM and the synthetic dataset in the supplementary
materials.

We compare MegBA with four baselines: (i) Ceres [1] (version 2.0) is the most
popular BA library that can efficiently use massive CPU cores, (ii) g2o [12]
is a lightweight CPU-based BA library, (iii) RootBA [7] is a recent CPU-based
BA library that uses Nullspace-Marginalization in place of Schur Complement,
and (iv) DeepLM [16] is the state-of-the-art GPU-based BA library (2021), and
it was shown to out-perform other popular BA libraries: STBA [42] and PBA [39]
(We provide comparison results between PBA and MegBA in the supplementary
materials).

We run experiments on a server that has 80 Intel Xeon 2.5GHz CPU cores, 8
Nvidia V100 GPUs and 320GB memory. The GPUs are inter-connected using NVLink
2.0. We use 64-bit floating points (FP64) unless otherwise specified.

(a) Venice-1778 (b) Final-13682

(c) Ladybug-1723 (d) Trafalgar-257

(e) Dubrovnik-356

Figure 2: Mean Squared Error over Time. MegBA-$X$-Y refers to $X$ GPUs while -a refers to auto-differentiation Jacobian and -m refers to analytical differentiation Jacobian. Ceres/RootBA/g2o-$X$ refers to $X$ CPU threads. | Trafalgar-257 | Ladybug-1723 | Dubrovnik-356
---|---|---|---
| MSE | Time | Mem | MSE | Time | Mem | MSE | Time | Mem
Ceres-16 | 0.434 | 8.160 | 1.659 | 0.562 | 34.50 | 2.093 | 0.393 | 116.0 | 2.550
DeepLM | 0.434 | 3.820 | 1.445 | 0.573 | 3.930 | 2.144 | 0.396 | 6.119 | 2.693
g2o-16 | 0.434 | 21.69 | 1.358 | 1.961 | 140.7 | 1.866 | 0.394 | 94.39 | 2.308
RootBA-16 | 0.433 | 3.307 | 1.468 | 0.562 | 7.050 | 2.423 | 0.393 | 78.16 | 3.942
MegBA-1-a | 0.438 | 1.364 | 1.270 | 0.560 | 0.932 | 2.450 | 0.411 | 3.640 | 3.940
MegBA-1-m | 0.438 | 1.148 | 1.010 | 0.560 | 0.774 | 1.660 | 0.411 | 3.263 | 2.480
Table 2: Small-scale experiments We only report the results of MegBA with a
GPU because the datasets in this table are small. MSE is the final Mean
Squared Error (pixels), Time is BA duration, and Mem is the memory in GB.

###  6.1 Overall Performance

We first evaluate the overall performance of MegBA, Ceres, g2o, RootBA, and
DeepLM. MegBA uses from 1 to 8 GPUs, and CPU-based algorithms use 16 threads.
We measure the Mean Squared Error (MSE) in pixels over time.

Figure 2 shows the evaluation results. In the Venice-1778 dataset (Figure
2(a)), MegBA achieves the best performance with 8 GPUs, while DeepLM can only
use a single GPU. MegBA completes with 3.34 seconds while Ceres, RootBA, g2o
uses 319.0, 73.94, and 890.6 seconds, respectively. It shows the substantial
speed-up (95.5$\times$, 22.1$\times$, and 266.6$\times$), which indicates the
benefits of fully exploiting GPUs to accelerate BA computation. For GPU-based
BA libraries, MegBA can complete with 11.96 seconds while DeepLM spent 24.44
seconds, showing the effectiveness of implementing full vectorisation for BA
on a single GPU. With more GPUs, MegBA out-performs DeepLM by 7.316$\times$,
which reflects the necessity of adopting multiple GPUs.

Thanks to the vectorisation and distributed BA designs, MegBA becomes the
state-of-the-art in the large BA dataset (i.e., Final-13682). As shown in
Figure 2(b), MegBA completes in 22.10 seconds, while DeepLM uses 149.6 seconds
(6.769$\times$ speed-up), Ceres uses 916 seconds (41.45$\times$ speed-up), g2o
uses 13161 seconds (595.5$\times$ speed-up), and RootBA uses 1427 seconds
(64.57$\times$ speed-up). In other datasets (Figure 2(c)-(e)), we observe
similar speed-up achieved by MegBA, indicating the general effectiveness of
our proposed approaches. We omit the discussion of these datasets, and their
results are reported in Table 2.

###  6.2 Scalability

| Venice-1778 | Final-13682
---|---|---
| MSE | Time | Mem | MSE | Time | Mem
Ceres-16 | 0.334 | 319.0 | 5.983 | 0.749 | 916.0 | 26.08
DeepLM | 0.333 | 24.44 | 6.256 | 0.751 | 149.6 | 14.89
g2o-16 | 0.335 | 890.6 | 5.999 | 1.061 | 13161 | 36.89
RootBA-16 | 0.337 | 73.94 | 14.14 | 0.773 | 1,427 | 263.2
MegBA-1-a | 0.333 | 11.96 | 13.68 | OOM | OOM | OOM
MegBA-2-a | 0.333 | 7.133 | 14.51 | OOM | OOM | OOM
MegBA-4-a | 0.333 | 4.767 | 16.76 | 0.748 | 28.70 | 81.03
MegBA-8-a | 0.333 | 3.340 | 22.61 | 0.748 | 22.10 | 89.74
MegBA-1-m | 0.333 | 10.92 | 7.870 | OOM | OOM | OOM
MegBA-2-m | 0.333 | 6.618 | 8.693 | 0.748 | 50.57 | 43.60
MegBA-4-m | 0.333 | 4.617 | 10.95 | 0.748 | 26.46 | 47.33
MegBA-8-m | 0.333 | 3.014 | 16.79 | 0.748 | 20.68 | 56.06
Table 3: Large-scale experiments.

Table 3 further provides the detailed experimental results to show the
scalability of MegBA, Ceres and DeepLM. In the Venice-1778 dataset, MegBA can
consistently improve its performance by adding GPUs (from 11.96 seconds to
3.34 seconds if we increase the number of GPUs from 1 to 8). In addition, the
large dataset (Final-13682) can better show the scalability of MegBA. By
increasing the number of GPUs from 4 to 8, the time can be reduced to 22.10
seconds.

###  6.3 Floating point precision

The accuracy of solving a BA problem is sensitive to the choice of floating
point precision (i.e., 32-bit vs. 64-bit floating points). We further evaluate
MegBA in all datasets with 32-bit and 64-bit floating points, and we report
the results of Venice-1778 and Final-13682 in Table 4. Other datasets show
consistent results and we omit them here. In the dataset of Final-13682, with
4 GPUs, MegBA (FP32) can complete in 4.804 seconds and MegBA (FP64) can
complete in 28.70 seconds, while both of them are reaching the same MSE. This
shows the exactness of the distributed BA algorithm in MegBA. Even with lower
precision, MegBA can reach the same MSE as double-precision; but offering
5.97$\times$ speed up, making MegBA (FP32) be the state-of-the-art in
Final-13682.

| Venice-1778 | Final-13682
---|---|---
| MSE | Time | Mem | MSE | Time | Mem
MegBA-1-a(FP32) | 0.334 | 2.620 | 8.300 | OOM | OOM | OOM
MegBA-1-a(FP64) | 0.333 | 11.96 | 13.68 | OOM | OOM | OOM
MegBA-1-m(FP32) | 0.333 | 2.065 | 4.821 | 0.750 | 11.82 | 24.51
MegBA-1-m(FP64) | 0.333 | 10.92 | 7.870 | OOM | OOM | OOM
MegBA-2-a(FP32) | 0.333 | 1.903 | 8.447 | 0.750 | 11.04 | 42.48
MegBA-2-a(FP64) | 0.333 | 7.133 | 14.51 | OOM | OOM | OOM
MegBA-2-m(FP32) | 0.333 | 1.353 | 5.541 | 0.750 | 5.133 | 25.63
MegBA-2-m(FP64) | 0.333 | 6.618 | 8.693 | 0.748 | 50.57 | 43.60
MegBA-4-a(FP32) | 0.333 | 1.680 | 10.50 | 0.749 | 4.804 | 45.28
MegBA-4-a(FP64) | 0.333 | 4.767 | 16.76 | 0.748 | 28.70 | 81.03
MegBA-4-m(FP32) | 0.334 | 1.274 | 7.598 | 0.748 | 4.279 | 28.43
MegBA-4-m(FP64) | 0.333 | 4.617 | 10.95 | 0.748 | 26.46 | 47.33
MegBA-8-a(FP32) | 0.334 | 1.622 | 16.02 | 0.748 | 8.973 | 52.15
MegBA-8-a(FP64) | 0.333 | 3.340 | 22.60 | 0.748 | 22.10 | 89.74
MegBA-8-m(FP32) | 0.333 | 1.271 | 12.99 | 0.747 | 7.582 | 35.31
MegBA-8-m(FP64) | 0.333 | 3.014 | 16.79 | 0.748 | 20.68 | 56.06
Table 4: Performance with 32-bit and 64-bit floating points.

##  7 Conclusion

We present MegBA, a novel GPU-based distributed BA library. MegBA has a set of
algorithms that enables automatically distributing BA computation to parallel
GPUs. It has a group of SIMD-optimised data structures, and a memory-efficient
runtime, making MegBA capable of fully utilising a GPU. MegBA has high-level
and compatible APIs, making it quickly become a popular open-sourced BA
library. Experimental results show that MegBA can out-perform SOTA BA
libraries by orders of magnitudes in several large-scale BA benchmarks.

## Acknowledgement

We would like to thank our colleagues at Megvii Inc., Dikai Fan, Shuxue Peng,
Yijia He, Zheng Chai, Haotian Zhang and Can Huang for their invaluable inputs.
Also, we thank Qingtian Zhu at Peking University and Xingxing Zuo at Technical
University of Munich to help with proof reading.

## References

  * [1] Agarwal, S., Mierle, K., Others: Ceres solver. http://ceres-solver.org
  * [2] Agarwal, S., Snavely, N., Seitz, S.M., Szeliski, R.: Bundle adjustment in the large. In: European conference on computer vision (2010)
  * [3] Chen, X., Ma, H., Wan, J., Li, B., Xia, T.: Multi-view 3d object detection network for autonomous driving. In: IEEE conference on Computer Vision and Pattern Recognition (2017)
  * [4] Chen, Y., Davis, T.A., Hager, W.W., Rajamanickam, S.: Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate. ACM Transactions on Mathematical Software (TOMS) (2008)
  * [5] Cieslewski, T., Choudhary, S., Scaramuzza, D.: Data-efficient decentralized visual slam. In: IEEE international conference on robotics and automation (2018)
  * [6] Demmel, N., Gao, M., Laude, E., Wu, T., Cremers, D.: Distributed photometric bundle adjustment. In: IEEE International Conference on 3D Vision (2020)
  * [7] Demmel, N., Sommer, C., Cremers, D., Usenko, V.: Square root bundle adjustment for large-scale reconstruction. In: IEEE Conference on Computer Vision and Pattern Recognition (2021)
  * [8] Demmel, N., Sommer, C., Cremers, D., Usenko, V.: Square root bundle adjustment for large-scale reconstruction. In: IEEE/CVF Conference on Computer Vision and Pattern Recognition (2021)
  * [9] Eisenstat, S.C.: Efficient implementation of a class of preconditioned conjugate gradient methods. SIAM Journal on Scientific and Statistical Computing (1981)
  * [10] Eriksson, A., Bastian, J., Chin, T.J., Isaksson, M.: A consensus-based framework for distributed bundle adjustment. In: IEEE Conference on Computer Vision and Pattern Recognition (2016)
  * [11] Fratarcangeli, M., Bradley, D., Gruber, A., Zoss, G., Beeler, T.: Fast nonlinear least squares optimization of large-scale semi-sparse problems. In: Computer Graphics Forum (2020)
  * [12] Grisetti, G., Kümmerle, R., Strasdat, H., Konolige, K.: g2o: A general framework for (hyper) graph optimization. In: IEEE International Conference on Robotics and Automation (2011)
  * [13] Guennebaud, G., Jacob, B., et al.: Eigen v3. http://eigen.tuxfamily.org (2010)
  * [14] Guo, Y., Liu, J., Li, G., Mai, L., Dong, H.: Fast and flexible human pose estimation with hyperpose. In: Proceedings of the 29th ACM International Conference on Multimedia (2021)
  * [15] Hackel, T., Savinov, N., Ladicky, L., Wegner, J.D., Schindler, K., Pollefeys, M.: Semantic3d. net: A new large-scale point cloud classification benchmark. arXiv preprint arXiv:1704.03847 (2017)
  * [16] Huang, J., Huang, S., Sun, M.: Deeplm: Large-scale nonlinear least squares on deep learning frameworks using stochastic domain decomposition. In: IEEE/CVF Conference on Computer Vision and Pattern Recognition (2021)
  * [17] Kaess, M., Johannsson, H., Roberts, R., Ila, V., Leonard, J.J., Dellaert, F.: isam2: Incremental smoothing and mapping using the bayes tree. The International Journal of Robotics Research (2012)
  * [18] Kaess, M., Ranganathan, A., Dellaert, F.: isam: Incremental smoothing and mapping. IEEE Transactions on Robotics (2008)
  * [19] Koliousis, A., Watcharapichat, P., Weidlich, M., Mai, L., Costa, P., Pietzuch, P.: Crossbow: Scaling deep learning with small batch sizes on multi-gpu servers. Proceedings of the VLDB Endowment
  * [20] Konolige, K., Agrawal, M.: Frameslam: From bundle adjustment to real-time visual mapping. IEEE Transactions on Robotics (2008)
  * [21] Levinson, J., Askeland, J., Becker, J., Dolson, J., Held, D., Kammel, S., Kolter, J.Z., Langer, D., Pink, O., Pratt, V., et al.: Towards fully autonomous driving: Systems and algorithms. In: IEEE intelligent vehicles symposium (2011)
  * [22] Li, P., Chen, X., Shen, S.: Stereo r-cnn based 3d object detection for autonomous driving. In: IEEE/CVF Conference on Computer Vision and Pattern Recognition (2019)
  * [23] Liu, H., Chen, M., Zhang, G., Bao, H., Bao, Y.: Ice-ba: Incremental, consistent and efficient bundle adjustment for visual-inertial slam. In: IEEE Conference on Computer Vision and Pattern Recognition (2018)
  * [24] Ma, X., Wang, Z., Li, H., Zhang, P., Ouyang, W., Fan, X.: Accurate monocular 3d object detection via color-embedded 3d reconstruction for autonomous driving. In: IEEE/CVF International Conference on Computer Vision (2019)
  * [25] Mai, L., Li, G., Wagenländer, M., Fertakis, K., Brabete, A.O., Pietzuch, P.: Kungfu: Making training in distributed machine learning adaptive. In: 14th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 20) (2020)
  * [26] Marquardt, D.W.: An algorithm for least-squares estimation of nonlinear parameters. Journal of the society for Industrial and Applied Mathematics (1963)
  * [27] Mayer, H.: Rpba – robust parallel bundle adjustment based on covariance information (2019)
  * [28] Ortiz, J., Pupilli, M., Leutenegger, S., Davison, A.J.: Bundle adjustment on a graph processor. In: IEEE/CVF Conference on Computer Vision and Pattern Recognition (2020)
  * [29] Park, Y., Lepetit, V., Woo, W.: Multiple 3d object tracking for augmented reality. In: IEEE/ACM International Symposium on Mixed and Augmented Reality (2008)
  * [30] Powell, M.J.: A hybrid method for nonlinear equations. Numerical methods for nonlinear algebraic equations (1970)
  * [31] Schönberger, J.L., Frahm, J.M.: Structure-from-motion revisited. In: Conference on Computer Vision and Pattern Recognition (2016)
  * [32] Schönberger, J.L., Zheng, E., Pollefeys, M., Frahm, J.M.: Pixelwise view selection for unstructured multi-view stereo. In: European Conference on Computer Vision (2016)
  * [33] Tang, C., Tan, P.: Ba-net: Dense bundle adjustment network. arXiv preprint arXiv:1806.04807 (2018)
  * [34] Triggs, B., McLauchlan, P.F., Hartley, R.I., Fitzgibbon, A.W.: Bundle adjustment—a modern synthesis. In: International workshop on vision algorithms (1999)
  * [35] Vo, M., Narasimhan, S.G., Sheikh, Y.: Spatiotemporal bundle adjustment for dynamic 3d reconstruction. In: IEEE Conference on Computer Vision and Pattern Recognition (2016)
  * [36] Wang, H.R., Lei, J., Li, A., Wu, Y.H.: A geometry-based point cloud reduction method for mobile augmented reality system. Journal of Computer Science and Technology (2018)
  * [37] Wedderburn, R.W.: Quasi-likelihood functions, generalized linear models, and the gauss—newton method. Biometrika (1974)
  * [38] Wilson, K., Snavely, N.: Robust global translations with 1dsfm. In: Proceedings of the European Conference on Computer Vision (2014)
  * [39] Wu, C., Agarwal, S., Curless, B., Seitz, S.M.: Multicore bundle adjustment. In: IEEE Conference on Computer Vision and Pattern Recognition (2011)
  * [40] Zhang, R., Zhu, S., Fang, T., Quan, L.: Distributed very large scale bundle adjustment by global camera consensus. In: 2017 IEEE International Conference on Computer Vision (2017)
  * [41] Zhao, Y., Guo, T.: Pointar: Efficient lighting estimation for mobile augmented reality. In: European Conference on Computer Vision (2020)
  * [42] Zhou, L., Luo, Z., Zhen, M., Shen, T., Li, S., Huang, Z., Fang, T., Quan, L.: Stochastic bundle adjustment for efficient and scalable 3d reconstruction. In: European Conference on Computer Vision (2020)

##  Appendix 0.A JetVector Design

To show the novel design of JetVector, we compare it with its predecessor: Jet
implemented in prior BA libraries such as Ceres [1].

A Jet object comprises value and grad (i.e., gradients). It has atomic
operators such as $+$, $-$, $*$, and $/$ which are essential for auto-
differentiation. Multiple Jet objects are contained in an Array-of-Structures
(AoS). Though effective for CPUs, the use of AoS makes it difficult to
coalesce GPU memory transactions, resulting in low efficiency in utilising
GPUs.

In contrast, JetVector is implemented as a Structure-of-Arrays which allows
GPU memory transactions to be effectively coalesced. Figure 3 compares Jet and
JetVector. As we can see, if we stack Jet in a vector, the memory addresses of
Value and Grad of different each Jet objects are not continuous. In JetVector,
however Value and Grad objects are organised in a continuous manner in the
memory, making it friendly for GPUs.

We show how does MegBA add JetVector objects in Figure 4. Supposing there are
$N$ Jets in JetVector. We launch $N$ GPU threads and each GPU thread processes
a Jet object. In this way, GPU threads will access adjacent data elements,
thus coalescing memory transactions. As a result, the memory loading/storing
operations in a warp (a warp is 32 consecutive GPU threads) can be achieved by
a 128-byte memory transaction. Notably, the memory transactions of a naive AoS
structure comprises multiple serialised 32-byte memory transactions, which can
adversely affect memory performance.

(a) Data structure of Jet

(b) Data structure of JetVector

Figure 3: Comparison between Jet and JetVector. Figure 4: Example of adding
two JetVector objects.

##  Appendix 0.B Edge-based Partitioning of BA Problems

We provide an example that shows how MegBA partitions non-zero Hessian blocks
and assigns the partitions to multiple GPUs in a load balancing manner.

Figure 5(a) shows the Jacobian structure of a BA problem. The black blocks
refer to non-zero blocks and the white blocks refer to zero blocks. It is
partitioned into two sub-matrices shown in Figure 5(c) and Figure 5(c).

Figure 6(a) is the structure of the Hessian matrix to be computed by MegBA. A
Hessian, though stored in the sparse format, is still too large to fit into a
GPU. To address this, MegBA stores a part of a Hessian in each GPU, as shown
in Figure 6(c) and Figure 6(c). Since we have $H=H_{0}+H_{1}$, it is
guaranteed that computations based on partitioned Hessian sub-matrices are
equivalent to the original computation.

(a) Jacobian $J$

(b) Partitioned Jacobian $J_{0}$ on GPU $0$

(c) Partitioned Jacobian $J_{1}$ on GPU $1$

Figure 5: A Jacobian matrix. (a) is the Jacobian of a BA problem. We partition
Jacobian into two matrix blocks (b) and (c). MegBA stores these two matrix
blocks on two GPUs.

(a) Hessian $H$

(b) Hessian $H_{0}$ on GPU $0$

(c) Hessian $H_{1}$ on GPU $1$

Figure 6: A Hessian matrix. (a) is the Hessian $H$ constructed from Jacobian
$J$ that $H=J^{T}J$. On each GPU, MegBA constructs a Hessian matrix so that
$H_{0}=J_{0}^{T}J_{0},H_{1}=J_{1}^{T}J_{1}$.

##  Appendix 0.C Experimental Results

We evaluate the performance of MegBA with the FP64, 1-GPU configuration and
compare MegBA against Ceres [1], g2o [12], RootBA [7] and DeepLM [16]. MegBA
uses the Levenberg–Marquardt algorithm and the trust-region strategy (same as
Ceres).

###  0.C.1 BAL Datasets

We summarise the results of the BAL dataset in Table 5 (Trafalgar), Table 6
(Ladybug), Table 7 (Dubrovnik), Table 8 (Venice), and Table 9 (Final). As we
can see in these tables, MegBA consistently outperforms all baselines by up to
24x. We also evaluated PBA [39] on the BAL dataset. PBA supports FP32 only, so
we evaluate MegBA(FP32) with PBA in Table 10, showing that MegBA outperforms
PBA by a large margin.

| Ceres-16 | DeepLM | g2o-16 | RootBA-16 | MegBA-1-a | MegBA-1-m
---|---|---|---|---|---|---
| MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time
Trafalgar-21 | 0.83 | 0.42 | 0.83 | 1.72 | 0.83 | 0.79 | 0.83 | 0.16 | 0.83 | 0.58 | 0.83 | 0.50
Trafalgar-39 | 0.95 | 1.02 | 0.95 | 2.63 | 0.95 | 2.52 | 0.95 | 0.38 | 0.95 | 0.88 | 0.95 | 0.81
Trafalgar-50 | 0.70 | 1.04 | 0.70 | 2.62 | 0.70 | 1.90 | 0.70 | 0.44 | 0.70 | 0.97 | 0.70 | 0.81
Trafalgar-126 | 0.62 | 3.36 | 0.62 | 3.98 | 0.62 | 8.76 | 0.62 | 8.06 | 0.62 | 1.83 | 0.62 | 1.52
Trafalgar-138 | 0.53 | 8.11 | 0.53 | 4.00 | 0.52 | 12.90 | 0.52 | 11.88 | 0.53 | 1.18 | 0.53 | 1.03
Trafalgar-161 | 0.47 | 4.76 | 0.47 | 4.11 | 0.46 | 9.04 | 0.46 | 14.64 | 0.47 | 1.12 | 0.47 | 1.11
Trafalgar-170 | 0.47 | 6.77 | 0.47 | 4.00 | 0.46 | 14.23 | 0.46 | 20.65 | 0.47 | 1.40 | 0.47 | 1.14
Trafalgar-174 | 0.47 | 5.63 | 0.46 | 4.17 | 0.45 | 23.42 | 0.46 | 19.35 | 0.46 | 1.22 | 0.46 | 1.19
Trafalgar-193 | 0.46 | 7.13 | 0.46 | 4.17 | 0.45 | 16.92 | 0.46 | 18.45 | 0.46 | 2.38 | 0.46 | 2.02
Trafalgar-201 | 0.48 | 2.83 | 0.46 | 3.72 | 0.50 | 27.17 | 0.46 | 5.18 | 0.47 | 1.30 | 0.47 | 1.22
Trafalgar-206 | 0.45 | 12.80 | 0.45 | 4.36 | 0.45 | 22.26 | 0.45 | 17.30 | 0.46 | 1.52 | 0.46 | 1.32
Trafalgar-215 | 0.45 | 9.97 | 0.45 | 4.33 | 0.45 | 24.01 | 0.45 | 12.63 | 0.45 | 1.24 | 0.45 | 1.21
Trafalgar-225 | 0.44 | 5.55 | 0.44 | 4.49 | 0.44 | 19.43 | 0.44 | 23.96 | 0.44 | 1.35 | 0.44 | 1.29
Trafalgar-257 | 0.43 | 8.16 | 0.43 | 3.82 | 0.43 | 21.69 | 0.43 | 3.307 | 0.44 | 1.36 | 0.44 | 1.15
Table 5: The results of the Trafalgar dataset. | Ceres-16 | DeepLM | g2o-16 | Rootba-16 | MegBA-1-a | MegBA-1-m
---|---|---|---|---|---|---
| MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time
Ladybug-49 | 0.42 | 1.07 | 0.42 | 2.23 | 0.42 | 2.50 | 0.42 | 0.32 | 0.42 | 0.19 | 0.42 | 0.19
Ladybug-73 | 0.37 | 1.24 | 0.37 | 2.72 | 0.37 | 2.12 | 0.37 | 0.67 | 0.37 | 0.24 | 0.37 | 0.22
Ladybug-138 | 0.70 | 2.72 | 0.70 | 3.00 | 0.70 | 31.21 | 0.82 | 0.79 | 0.70 | 0.32 | 0.70 | 0.29
Ladybug-318 | 0.48 | 3.04 | 0.48 | 2.95 | 0.48 | 65.70 | 0.48 | 3.78 | 0.48 | 0.44 | 0.48 | 0.33
Ladybug-372 | 0.55 | 4.45 | 0.55 | 2.10 | 0.55 | 83.52 | 0.55 | 2.96 | 0.55 | 0.54 | 0.55 | 0.41
Ladybug-412 | 0.50 | 4.07 | 0.50 | 3.34 | 0.49 | 86.94 | 0.50 | 2.83 | 0.49 | 0.58 | 0.49 | 0.44
Ladybug-460 | 0.53 | 6.02 | 0.53 | 3.47 | 0.53 | 243.32 | 0.53 | 5.38 | 0.53 | 0.69 | 0.53 | 0.55
Ladybug-539 | 0.55 | 8.37 | 0.55 | 3.85 | 0.55 | 65.63 | 0.55 | 5.08 | 0.55 | 0.75 | 0.55 | 0.57
Ladybug-598 | 0.55 | 8.54 | 0.55 | 3.34 | 0.54 | 293.69 | 0.55 | 4.16 | 0.55 | 0.79 | 0.55 | 0.60
Ladybug-646 | 0.55 | 9.21 | 0.55 | 3.56 | 0.56 | 200.83 | 0.55 | 5.44 | 0.55 | 0.92 | 0.55 | 0.74
Ladybug-707 | 0.57 | 15.30 | 0.56 | 4.48 | 0.56 | 283.15 | 0.56 | 6.59 | 0.56 | 0.98 | 0.56 | 0.76
Ladybug-783 | 0.53 | 7.69 | 0.53 | 5.01 | 0.53 | 539.85 | 0.52 | 15.95 | 0.53 | 0.98 | 0.53 | 0.78
Ladybug-810 | 0.52 | 9.67 | 0.52 | 4.77 | 0.52 | 391.78 | 0.52 | 13.71 | 0.52 | 1.00 | 0.52 | 0.74
Ladybug-856 | 0.51 | 13.30 | 0.51 | 4.81 | 0.51 | 114.89 | 0.51 | 8.01 | 0.51 | 1.10 | 0.51 | 0.86
Ladybug-885 | 0.51 | 15.00 | 0.51 | 5.08 | 0.51 | 308.21 | 0.51 | 8.03 | 0.51 | 1.01 | 0.51 | 0.79
Ladybug-931 | 0.55 | 15.30 | 0.55 | 4.56 | 0.55 | 220.54 | 0.55 | 11.08 | 0.55 | 1.18 | 0.55 | 0.96
Ladybug-969 | 0.55 | 17.50 | 0.54 | 5.21 | 0.55 | 197.59 | 0.54 | 7.79 | 0.54 | 1.19 | 0.54 | 0.94
Ladybug-1031 | 0.55 | 12.90 | 0.55 | 4.04 | 0.55 | 203.63 | 0.55 | 10.24 | 0.55 | 1.48 | 0.55 | 1.18
Ladybug-1064 | 0.55 | 15.40 | 0.56 | 3.17 | 0.55 | 355.78 | 0.55 | 4.47 | 0.55 | 1.23 | 0.55 | 1.03
Ladybug-1118 | 0.57 | 15.40 | 0.58 | 3.68 | 0.58 | 311.09 | 0.58 | 6.16 | 0.57 | 1.39 | 0.57 | 1.11
Ladybug-1152 | 0.56 | 13.70 | 0.56 | 3.07 | 0.56 | 488.39 | 0.56 | 7.53 | 0.56 | 1.46 | 0.56 | 1.10
Ladybug-1197 | 0.57 | 17.30 | 0.57 | 4.13 | 0.57 | 122.18 | 0.57 | 8.58 | 0.57 | 1.40 | 0.57 | 1.13
Ladybug-1235 | 0.56 | 20.30 | 0.58 | 3.34 | 0.56 | 212.62 | 0.57 | 5.89 | 0.56 | 1.30 | 0.56 | 1.06
Ladybug-1266 | 0.57 | 23.60 | 0.56 | 4.53 | 0.56 | 342.70 | 0.56 | 8.81 | 0.56 | 1.44 | 0.56 | 1.18
Ladybug-1340 | 0.57 | 26.20 | 0.57 | 4.54 | 2.02 | 453.25 | 0.57 | 9.17 | 0.57 | 1.61 | 0.57 | 1.31
Ladybug-1469 | 0.56 | 26.30 | 0.57 | 3.84 | 0.58 | 724.13 | 0.56 | 8.30 | 0.56 | 1.73 | 0.56 | 1.37
Ladybug-1514 | 0.56 | 25.20 | 0.56 | 4.69 | 1.18 | 459.05 | 0.56 | 8.61 | 0.56 | 1.59 | 0.56 | 1.34
Ladybug-1587 | 0.59 | 37.50 | 0.56 | 4.82 | 1.15 | 567.71 | 0.58 | 6.77 | 0.56 | 2.01 | 0.56 | 1.58
Ladybug-1642 | 0.58 | 36.10 | 0.56 | 3.64 | 0.76 | 127.00 | 0.56 | 16.12 | 0.56 | 1.73 | 0.56 | 1.43
Ladybug-1695 | 0.56 | 32.30 | 0.56 | 4.24 | 0.62 | 799.46 | 0.59 | 5.46 | 0.56 | 1.87 | 0.56 | 1.51
Ladybug-1723 | 0.56 | 34.50 | 0.57 | 3.93 | 1.96 | 140.69 | 0.56 | 7.05 | 0.56 | 0.93 | 0.56 | 0.77
Table 6: The results of the Ladybug dataset. | Ceres-16 | DeepLM | g2o-16 | Rootba-16 | MegBA-1-a | MegBA-1-m
---|---|---|---|---|---|---
| MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time
Dubrovnik-16 | 0.22 | 0.79 | 0.22 | 1.75 | 0.22 | 1.72 | 0.22 | 0.36 | 0.22 | 1.32 | 0.22 | 1.24
Dubrovnik-88 | 0.75 | 7.69 | 0.75 | 4.77 | 0.75 | 21.87 | 0.75 | 11.14 | 0.75 | 1.89 | 0.75 | 1.63
Dubrovnik-135 | 0.67 | 17.00 | 0.67 | 5.88 | 0.67 | 34.21 | 0.67 | 17.99 | 0.67 | 1.52 | 0.67 | 1.50
Dubrovnik-142 | 0.48 | 17.30 | 0.48 | 6.03 | 0.48 | 30.11 | 0.48 | 6.74 | 0.48 | 2.34 | 0.48 | 2.25
Dubrovnik-150 | 0.43 | 15.20 | 0.43 | 6.27 | 0.43 | 31.32 | 0.43 | 33.01 | 0.43 | 1.11 | 0.43 | 1.00
Dubrovnik-161 | 0.41 | 16.80 | 0.41 | 7.61 | 0.41 | 31.96 | 0.41 | 19.32 | 0.41 | 0.71 | 0.41 | 0.69
Dubrovnik-173 | 0.41 | 15.40 | 0.41 | 7.38 | 0.41 | 38.44 | 0.41 | 59.45 | 0.41 | 0.96 | 0.41 | 0.84
Dubrovnik-182 | 0.45 | 13.40 | 0.45 | 8.11 | 0.45 | 41.91 | 0.45 | 67.52 | 0.45 | 1.08 | 0.45 | 1.00
Dubrovnik-202 | 0.43 | 21.70 | 0.43 | 9.40 | 0.44 | 54.12 | 0.43 | 22.19 | 0.43 | 1.41 | 0.43 | 1.29
Dubrovnik-237 | 0.42 | 25.50 | 0.42 | 9.95 | 0.42 | 59.35 | 0.42 | 17.27 | 0.41 | 0.91 | 0.41 | 0.80
Dubrovnik-253 | 0.38 | 29.00 | 0.38 | 10.82 | 0.39 | 69.64 | 0.38 | 50.39 | 0.38 | 1.45 | 0.38 | 1.24
Dubrovnik-262 | 0.37 | 95.90 | 0.37 | 11.25 | 0.37 | 64.68 | 0.37 | 52.60 | 0.37 | 2.79 | 0.37 | 2.59
Dubrovnik-273 | 0.37 | 47.90 | 0.37 | 11.41 | 0.37 | 63.24 | 0.36 | 60.37 | 0.37 | 1.87 | 0.37 | 1.78
Dubrovnik-287 | 0.36 | 26.80 | 0.36 | 11.43 | 0.37 | 93.30 | 0.36 | 40.84 | 0.36 | 2.35 | 0.36 | 2.13
Dubrovnik-308 | 0.37 | 33.80 | 0.37 | 11.22 | 0.37 | 91.82 | 0.37 | 37.15 | 0.37 | 3.23 | 0.37 | 2.88
Dubrovnik-356 | 0.39 | 116.00 | 0.40 | 6.12 | 0.39 | 94.39 | 0.39 | 78.16 | 0.41 | 3.64 | 0.41 | 3.26
Table 7: The results of the Dubrovnik dataset. | Ceres-16 | DeepLM | g2o-16 | Rootba-16 | MegBA-1-a | MegBA-1-m
---|---|---|---|---|---|---
| MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time
Venice-52 | 0.75 | 7.95 | 0.75 | 4.66 | 0.68 | 17.94 | 0.68 | 2.39 | 0.73 | 1.85 | 0.73 | 1.59
Venice-89 | 0.50 | 8.47 | 0.50 | 6.41 | 0.50 | 27.53 | 0.51 | 4.51 | 0.50 | 5.05 | 0.50 | 4.85
Venice-245 | 0.84 | 31.90 | 0.87 | 5.66 | 0.87 | 123.27 | 0.85 | 10.47 | 0.84 | 3.26 | 0.84 | 3.06
Venice-427 | 0.63 | 58.20 | 0.63 | 15.46 | 0.62 | 115.24 | 0.63 | 43.12 | 0.63 | 5.17 | 0.63 | 5.17
Venice-744 | 0.51 | 109.0 | 0.51 | 27.79 | 0.51 | 659.83 | 0.51 | 43.60 | 0.50 | 6.26 | 0.50 | 6.04
Venice-951 | 0.44 | 84.30 | F | F | 0.57 | 793.15 | 0.45 | 59.45 | 0.45 | 13.80 | 0.45 | 13.64
Venice-1102 | 0.40 | 146.0 | 0.39 | 32.80 | 0.56 | 585.86 | 0.40 | 94.13 | 0.39 | 14.42 | 0.39 | 13.24
Venice-1158 | 0.42 | 33.10 | 0.43 | 12.07 | 0.48 | 581.79 | 0.36 | 51.18 | 0.36 | 6.74 | 0.36 | 6.32
Venice-1184 | 0.77 | 30.10 | 0.35 | 36.32 | 0.46 | 600.57 | 0.35 | 38.16 | 0.34 | 36.24 | 0.34 | 28.52
Venice-1238 | 0.35 | 71.40 | 0.34 | 13.20 | 0.46 | 531.25 | 0.34 | 90.36 | 0.34 | 8.36 | 0.34 | 7.75
Venice-1288 | 0.33 | 153.0 | 0.33 | 19.79 | 0.45 | 630.90 | 0.33 | 95.45 | 0.33 | 8.36 | 0.33 | 7.64
Venice-1350 | 0.34 | 52.30 | 0.34 | 10.47 | 0.44 | 551.41 | 0.34 | 53.85 | 0.33 | 6.79 | 0.33 | 6.43
Venice-1408 | 0.35 | 76.50 | 0.35 | 10.72 | 0.47 | 687.48 | 0.35 | 47.00 | 0.34 | 5.66 | 0.34 | 5.34
Venice-1425 | 0.34 | 153.0 | 0.34 | 10.80 | 0.44 | 635.88 | 0.34 | 183.28 | 0.34 | 5.67 | 0.34 | 5.28
Venice-1473 | 0.33 | 110.0 | 0.33 | 43.24 | 0.42 | 571.33 | 0.33 | 61.21 | 0.33 | 7.05 | 0.33 | 6.68
Venice-1490 | 0.33 | 66.80 | 0.33 | 14.25 | 0.33 | 571.17 | 0.33 | 123.52 | 0.32 | 5.64 | 0.32 | 5.32
Venice-1521 | 0.33 | 75.60 | 0.33 | 10.98 | 0.33 | 666.63 | 0.33 | 93.50 | 0.32 | 7.47 | 0.32 | 6.78
Venice-1544 | 0.33 | 173.0 | 0.33 | 26.59 | 0.42 | 653.50 | 0.33 | 132.69 | 0.32 | 5.71 | 0.32 | 5.33
Venice-1638 | 0.57 | 75.40 | 0.58 | 30.04 | 0.56 | 1015 | 0.57 | 54.58 | 0.58 | 17.66 | 0.58 | 16.60
Venice-1666 | 0.48 | 74.00 | 0.48 | 23.48 | 0.45 | 1033 | 0.43 | 85.30 | 0.44 | 10.75 | 0.44 | 10.08
Venice-1672 | 0.38 | 260.0 | 0.38 | 33.49 | 0.38 | 991.73 | 0.37 | 123.52 | 0.40 | 12.20 | 0.40 | 11.06
Venice-1681 | 0.37 | 48.40 | 0.34 | 46.98 | 0.34 | 1067 | 0.34 | 118.23 | 0.34 | 9.16 | 0.34 | 8.38
Venice-1682 | 0.34 | 164.0 | 0.33 | 38.92 | 0.34 | 956.35 | 0.33 | 274.34 | 0.33 | 17.39 | 0.33 | 15.78
Venice-1684 | 0.34 | 207.0 | 0.33 | 43.94 | 0.33 | 1658 | 0.33 | 202.10 | 0.33 | 11.78 | 0.33 | 11.25
Venice-1695 | 0.37 | 53.40 | 0.34 | 44.99 | 0.34 | 1648 | 0.33 | 292.24 | 0.34 | 16.16 | 0.34 | 14.48
Venice-1696 | 0.38 | 43.00 | 0.33 | 47.04 | 0.34 | 1335 | 0.33 | 243.68 | 0.34 | 7.43 | 0.34 | 7.06
Venice-1706 | 0.34 | 234.0 | 0.34 | 38.17 | 0.33 | 1845 | 0.34 | 196.48 | 0.34 | 19.72 | 0.34 | 18.59
Venice-1776 | 0.33 | 192.0 | 0.33 | 46.69 | 0.35 | 1147 | 0.34 | 62.10 | 0.33 | 10.56 | 0.33 | 10.02
Venice-1778 | 0.33 | 319.0 | 0.33 | 24.44 | 0.34 | 890.55 | 0.34 | 73.94 | 0.33 | 11.96 | 0.33 | 10.92
Table 8: The results of the Venice dataset. F stands for failed. | Ceres-16 | DeepLM | g2o-16 | Rootba-16 | MegBA-1-a | MegBA-1-m
---|---|---|---|---|---|---
| MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time
Final-93 | 0.51 | 4.10 | 0.51 | 1.56 | 0.51 | 6.22 | 0.51 | 1.84 | 0.51 | 0.78 | 0.51 | 0.67
Final-394 | 0.56 | 30.10 | 0.56 | 5.10 | 0.56 | 101.72 | 0.56 | 20.07 | 0.56 | 1.10 | 0.56 | 0.92
Final-871 | 0.62 | 115.0 | 0.62 | 26.12 | 0.62 | 481.22 | 0.62 | 76.65 | 0.62 | 8.20 | 0.62 | 7.75
Final-961 | 0.94 | 44.30 | 0.94 | 12.82 | 0.94 | 710.55 | 0.94 | 348.11 | 0.94 | 3.32 | 0.94 | 2.85
Final-1936 | 0.89 | 106.0 | 0.89 | 17.03 | 0.89 | 789.57 | 0.89 | 595.28 | 0.89 | 7.03 | 0.89 | 6.67
Final-3068 | 1.09 | 31.60 | 1.03 | 5.49 | 1.13 | 26203 | 0.98 | 72.27 | 1.01 | 2.82 | 1.01 | 2.40
Final-4585 | 0.57 | 417.0 | 0.56 | 66.13 | 0.56 | 23221 | 0.56 | 373.30 | 0.57 | 22.22 | 0.57 | 15.60
Table 9: The results of the Final dataset. | PBA(FP32) | MegBA-1-m(FP32)
---|---|---
| MSE | Time | MSE | Time
Trafalgar-257 | 0.85 | 0.57 | 0.44 | 0.11
Ladybug-1695 | 34.13 | 0.37 | 1.33 | 0.32
Dubrovnik-287 | 0.84 | 1.22 | 0.38 | 0.10
Venice-89 | 0.86 | 1.32 | 0.51 | 0.28
Final-4585 | 1.87 | 4.96 | 0.57 | 3.13
Table 10: We evaluate PBA with the BAL dataset. To make a fair comparison, we
compare PBA against the FP32 MegBA. Both PBA and FP32 MegBA suffer from
numerical instability issues in some datasets. In addition, PBA uses texture
memory, it failed to solve large-scale BA problems due to memory limitations.
We chose the largest datasets that can be solved successfully by PBA to
evaluate its performance. The result shows that MegBA outperforms PBA by a
large margin, with maximum 12.5 times speed up while achieving a lower MSE.

###  0.C.2 1DSfM Dataset

We further compare MegBA with the baselines using the 1DSfM [38] dataset. The
statistics, such as the number of observations, of the 1DSfM dataset can be
found in Table 11. Results are shown in Table 12. 1DSfM is a challenging
dataset and other algorithms can converge to local optima that have several
magnitudes of orders larger MSE than other algorithms, while MegBA always
converges to a competitive optima.

Dataset | #Points | #Observations
---|---|---
Alamo-577 | 140080 | 816891
Ellis Island-233 | 9210 | 20500
Gendarmenmarkt-704 | 76964 | 268747
Madrid Metropolis-347 | 44479 | 195660
Montreal Notre Dame-459 | 151876 | 811757
Notre Dame-548 | 224153 | 1172145
NYC Library-334 | 54757 | 211614
Piazza del Popolo-336 | 29731 | 150161
Piccadilly-2292 | 184475 | 798085
Roman Forum-1067 | 223844 | 1031760
Tower of London-484 | 126648 | 596690
Trafalgar-5052 | 327920 | 1266102
Union Square-816 | 26430 | 90668
Vienna Cathedral-846 | 154394 | 495940
Yorkminster-429 | 100426 | 376980
Table 11: The statistics of the datasets in 1DSfM. | Ceres-16 | DeepLM | g2o-16 | Rootba-16 | MegBA-1-m
---|---|---|---|---|---
| MSE | Time | MSE | Time | MSE | Time | MSE | Time | MSE | Time
Alamo | 4.8E+2 | 219.0 | 2.3E+2 | 40.4 | 7.6E+5 | 48.9 | 2.3E+2 | 31.1 | 2.2E+2 | 18.9
Ellis Island | 2.2E+3 | 6.3 | 2.3E+3 | 9.9 | 4.6E+5 | 2.7 | 1.6E+3 | 1.8 | 2.8E+3 | 0.5
Gen.markt | 2.5E+1 | 63.0 | 6.2E+1 | 31.2 | 4.4E+4 | 46.0 | 1.5E+2 | 20.8 | 1.6E+1 | 2.4
M.Metropolis | 8.9E+3 | 45.7 | 8.4E+3 | 21.2 | 2.8E+8 | 51.4 | 8.0E+3 | 15.8 | 1.1E+4 | 1.1
M.N.Dame | 8.4E+2 | 206.0 | 7.1E+2 | 73.4 | 5.0E+5 | 89.1 | 7.3E+2 | 46.4 | 7.4E+2 | 10.6
Notre Dame | 1.5E+4 | 282.0 | 1.6E+4 | 92.3 | 4.2E+8 | 68.6 | 2.3E+4 | 48.3 | 2.0E+4 | 10.8
NYC Library | 4.8E+0 | 62.7 | 1.2E+1 | 37.4 | 2.4E+0 | 23.2 | 3.0E+1 | 17.7 | 4.1E+0 | 1.4
P.del Popolo | 7.5E+1 | 26.5 | 8.4E+1 | 4.6 | 5.3E+6 | 1.4 | 3.2E+1 | 19.3 | 3.7E+1 | 1.4
Piccadilly | 4.8E+2 | 213.0 | 6.3E+2 | 68.3 | 3.9E+6 | 21.6 | 3.7E+2 | 73.5 | 3.7E+2 | 14.3
R.Forum | 6.5E+1 | 104.0 | 5.6E+1 | 88.5 | 1.1E+5 | 201.6 | 5.8E+1 | 105.8 | 5.5E+1 | 6.0
T.of London | 1.4E+3 | 81.9 | 1.8E+3 | 53.5 | 6.6E+8 | 7.1 | 1.9E+3 | 45.9 | 1.3E+3 | 14.7
Trafalgar | 1.7E+2 | 358.0 | 1.6E+2 | 112.1 | 4.8E+5 | 34.8 | 1.7E+2 | 211.3 | 1.5E+2 | 15.0
Union Square | 5.4E+2 | 24.2 | 3.9E+2 | 5.8 | 2.9E+6 | 23.1 | 2.3E+2 | 8.9 | 3.4E+2 | 0.8
V.Cathedral | 5.7E+2 | 139.0 | 5.2E+2 | 49.2 | 1.1E+5 | 51.6 | 4.9E+2 | 13.1 | 4.9E+2 | 12.2
Yorkminster | 1.2E+1 | 115.0 | 9.2E+0 | 18.9 | 4.2E+0 | 25.9 | 1.0E+2 | 29.8 | 9.0E+0 | 4.4
Table 12: Results of the 1DSfM dataset. | ceres16 | DeepLM | g2o | RootBA | MegBA-1-a | MegBA-1-m
---|---|---|---|---|---|---
Alamo | 2.16 | 1.86 | 1.96 | 4.39 | 2.88 | 1.85
Ellis Island | 1.46 | 1.53 | 1.15 | 1.06 | 0.67 | 0.65
Gen.markt | 1.72 | 1.83 | 1.42 | 1.13 | 1.36 | 1.07
M.Metropolis | 1.62 | 1.73 | 1.33 | 1.14 | 1.15 | 0.95
M.N.Dame | 2.16 | 2.40 | 1.92 | 3.85 | 2.87 | 1.85
Notre Dame | 2.50 | 2.39 | 2.29 | 5.54 | 3.88 | 2.47
NYC Library | 1.64 | 1.76 | 1.35 | 1.13 | 1.20 | 0.98
P.del Popolo | 1.57 | 1.67 | 1.29 | 1.31 | 1.02 | 0.89
Piccadilly | 2.20 | 2.47 | 2.18 | 3.30 | 2.85 | 1.85
R.Forum | 2.38 | 2.56 | 2.17 | 3.40 | 3.49 | 2.17
T.of London | 1.99 | 2.17 | 1.71 | 1.80 | 2.28 | 1.55
Trafalgar | 2.72 | 2.99 | 3.15 | 4.97 | 4.15 | 2.51
Union Square | 1.53 | 1.58 | 1.24 | 1.09 | 0.86 | 0.80
V.Cathedral | 1.96 | 2.10 | 1.67 | 1.36 | 2.00 | 1.40
Yorkminster | 1.81 | 2.02 | 1.52 | 1.19 | 1.67 | 1.23
Table 13: The memory overhead of the 1DSfM dataset. The unit is GB.

###  0.C.3 Large Synthesised Dataset

In the end, we evaluate MegBA on a larger synthetic dataset:
SynthesisedData-20000. This dataset emulates a real-world BA problem we have
in production. We use $\mathcal{U}(a,b)$ to denote a uniform distribution with
a minimum as $a$ and a maximum as $b$. We generate 20,000 cameras on a circle
of radius 8 uniformly. We add a uniform noise
$\epsilon\sim\mathcal{U}(0,0.01)$ to camera rotation and translation, as well
as a uniform noise $\epsilon\sim\mathcal{U}(0,0.5)$ to the camera intrinsic.
We generate 80, 000 points. For each point, the $x,y$ coordinates are sampled
from a uniform distribution $\mathcal{U}(-0.1,0.1)$ and the $z$ coordinate is
sampled from a uniform distribution $\mathcal{U}(-0.03,0.03)$. The $x,y$
coordinates of each point are also added noise
$\epsilon\sim\mathcal{U}(-0.1,0.1)$. Each point can be captured by 1,000
cameras, and there are 80,000,000 observations in total.

The evaluation was conducted on a GPU server with 300 GB CPU memory in total
and 32 GB GPU memory for each GPU. RootBA [7] and DeepLM [16] incur out-of-
memory (OOM) errors and fail to complete in this experiment. g2o cannot return
results in a reasonable time (1 hour) and throws an out-of-time error.

Figure 7 shows the experimental results of MegBA, Ceres, DeepLM and g2o.
DeepLM and g2o both cannot solve this problem due to memory and computation
limitations. Only MegBA and Ceres can process such a large dataset; however,
MegBA is 20$\times$ faster than Ceres, making MegBA the state-of-the-art BA
library for large BA problems (with more than 100s millions of observations).

Figure 7: Experimental results of the synthesised dataset. DeepLM throws an
OOM error after the 6th step. g2o is significantly slower than all other
baselines and we terminate its execution after 6 steps.

##  Appendix 0.D GPU Utilisation

Solving BA problems on a single GPU can suffer from computation and memory
bottlenecks. We are interested in if the computation bottleneck can be
resolved by adding more GPUs. We used nvprof to profile MegBA and calculated
its GPU utilisation. In particular, we plot the time proportion of different
computation in MegBA. We repeat the same experiment with 1, 2, 4 GPUs and
report the results in Figure 8. According to this figure, the profiling
results show that using multiple GPUs incurs only marginal communication cost
and the computation time decreases around 50%.

Figure 8: The visualisation of the profiling result returned by nvprof. We use
the dataset Venice-1778. SpMV denotes sparse matrix-vector multiplication.
Using 2 GPUs, the all-reduce communication accounts for 2.5%, while the
overall time reduces 60.6% compared with the result using 1 GPU. Using 4 GPUs,
the all-reduce communication accounts for 7.8%, while the overall time reduces
69.8% compared with the result using 2 GPUs.

##  Appendix 0.E Reconstruction Plots

Finally, we show the visualisation results of MegBA. Specifically, we show the
landmark point clouds (i.e. black dots) with cameras (i.e. red frames) for
each dataset. The point clouds are rendered by COLMAP [31][32].

Figure 9: Ladybug Figure 10: Trafalgar Square Figure 11: Dubrovnik Figure 12:
Venice Figure 13: Final