Implementation and Validation of Distributed Bundle Adjustment for Super Large Scale Datasets

Abstract

We propose a distributed bundle adjustment (DBA) method using the Levenberg-Marquardt (LM) algorithm for super large-scale datasets. Most of the existing methods partition the global map to small ones and conduct bundle adjustment in the submaps. In order to fit the parallel framework, they use approximate solutions instead of the LM algorithm. However, those methods often give sub-optimal results. Different from them, we utilize the LM algorithm to conduct global bundle adjustment where the formation of the reduced camera system (RCS) is actually parallelized and executed in a distributed way. To store the large RCS, we compress it with a block-based sparse matrix compression format (BSMC), which fully exploits its block feature. The BSMC format also enables the distributed storage and updating of the global RCS. The proposed method is extensively evaluated and compared with the state-of-the-art pipelines using both synthetic and real datasets. Preliminary results demonstrate the efficient memory usage and vast scalability of the proposed method compared with the baselines. For the first time, we conducted parallel bundle adjustment using LM algorithm on a real datasets with 1.18 million images and a synthetic dataset with 10 million images (about 500 times that of the state-of-the-art LM-based BA) on a distributed computing system.

Appendices

Appendix A Accessing Efficiency of the BSMC and CSR Format

To evaluate the accessing efficiency of a compression format, we consider two operations. (1) Given an element or a block in the compressed structure, find its column and row identities in the full matrix; (2) Given the column and row identities of an element or a block in the full matrix, find the corresponding block in the compressed structure.

1.1 Appendix A.1 The first Operation

Given the ID of an element in the first array of the CSR structure, the column ID in the full matrix can be identified immediately by the second array as shown in Figure 15 (left). If the elements are continuously accessed, the row ID in the full matrix can be obtained by row counting. Otherwise, the row ID must be interpolated from the third array. Use a dichotomy algorithm, the maximum searching times is $\log _2nc$ since the size of the third array is nc where n is the block number and c is the block width.

Given the identity of a block in the first array of the BSMC structure, the column ID in the full matrix is firstly identified by the second array as shown in Figure 15 (right). Similar to CSR format, the row identity can be obtained by row counting if the blocks are continuously accessed, or the row identity must be searched in the third array with the maximum searching times of $\log _2n$ since the size of the third array is the block number n. Note that one access in CSR structure only gets one element but it gets a block which includes $c\times c$ elements in BSMC structure.

1.2 Appendix A.2 The Second Operation

Given the row and column IDs of a element in full matrix, to find its corresponding position in the CSR structure, we first locate the start and end column IDs in the second array from the third array by the row ID in the full matrix as shown in Figure 15 (left), then the column ID of the element in the first array can be searched between the start and end column IDs. The maximum searching times for each element is $\log _2cs_i$. To find a block, we only have to find the beginning element for each row of this block. Therefore the total searching times for finding a block is equal to $\sum ^c_{i=0}\log _2cs_i$ where $s_i$ is the number of non-zero blocks in row i.

Given the row and column identities of a block in the full matrix, to find its corresponding position in the BSMC structure, we first locate the start and end column IDs of the block in the second array from the third array by the row identity as shown in Figure 15 (right), and then the column ID in the first array is searched between the start and end column IDs with maximum searching times equal to $\log _2s_i$. Once the first element of the block is determined, the element of the whole block is identified. Apparently, for this operation, the BSMC is more efficient than the CSR.

Fig. 15

Accessing the data in CSR (left) and BSMC (right) format

Table 7 Accessing efficiency of the BSMC and CSR formats

Table 7 shows the searching times of these two operations for the BSMC and CSR formats. As shown in Table 7, if the elements or blocks in the compressed structure are accessed continuously, the first operation will be immediately performed for the two formats. In the real case, for matrix-vector product, the elements or blocks in the compressed structure are continuously accessed. However, for matrix updating or aggregation, the second operation is frequently applied. The accessing efficiency of BSMC for both the first and second operations is better than the CSR; thus, it’s strongly recommended for the compression of block-based sparse matrices, such as the Hessian and RCS.

Fig. 16

The jacobian matrix where $\textbf{A}_(i,j)$ and $\textbf{B}_(j,i)$ are the camera and landmark part of the jacobian generated by the edge of camera i and landmark j

Appendix B Formation of the RCS

The formation of RCS is a major step in BA. It dominates the most part of the computation complexity. In this section, we look into the details of the formation of a RCS and try to parallelize it. To be specific, given a dataset with 5 cameras and 3 landmarks, the Jacobian is shown in Figure 16, and the Hessian is obtained by normalizing the Jacobian as shown in Figure 17. The damping is not shown in the figure, but it does not affect the structure of the Hessian.

The computation of non-zero blocks of the Hessian are shown in equation (B1), (B2) and (B3).

$$\begin{aligned} \textbf{A}_i^T \textbf{A}_i=\sum _{j=0}^{p_i}\textbf{A}_{i,j}^T\textbf{A}_{i,j} \end{aligned}$$

(B1)

Where $p_i$ is the number of landmarks seen by camera i.

$$\begin{aligned} \textbf{B}_i^T \textbf{B}_i=\sum _{i=0}^{q_j}\textbf{B}_{j,i}^T\textbf{B}_{j,i} \end{aligned}$$

(B2)

Where $q_j$ is the number of cameras seeing landmark j.

$$\begin{aligned} \textbf{P}_{i,j}=\textbf{A}_{i,j}^T \textbf{B}_{j,i} \end{aligned}$$

(B3)

According to equation in the main text, conduct Shur complement to the Hessian, we get a matrix $\textbf{Q}=\textbf{WV}^{-1} \textbf{W}^T$ as shown in Figure 19 (left). Then the RCS is computed by equation $\textbf{R}=\textbf{U}-\textbf{Q}$ as shown in Figure 19 (right). So, to compute Shur complement is actually to compute the matrix $\textbf{Q}$. Each block $\textbf{Q}_{i,j}$ in matrix $\textbf{Q}$ is computed according to equation (B4). Note that the block $\textbf{Q}_{i,j}$ could be zero blocks if there is none common point between camera i and camera j.

Fig. 17

The Hessian matrix

Fig. 18

Structure of matrix $\textbf{Q}_k$. The non-zero blocks are not identified since k is unknown

$$\begin{aligned} \textbf{Q}_{i,j}= & \sum _{k=0}^{t_{i,j}}\textbf{Q}_{i,j,k} \end{aligned}$$

(B4)

$$\begin{aligned} \textbf{Q}_{i,j,k}= & \textbf{P}_{i,k}(\textbf{B}_k^T\textbf{B}_k)^{-1}\textbf{P}_{j,k}^T \end{aligned}$$

(B5)

Where k is the landmark ID, $t_{i,j}$ is the number of common points between camera i and camera j.

$$\begin{aligned} \textbf{Q}=\sum _{k=0}^L\textbf{Q}_k \end{aligned}$$

(B6)

Where the matrix $\textbf{Q}_k$ is a component of the RCS, L is the number of 3D points, and k is the ID of 3D points.

Conduct Shur complement to the 3D points, each of which generates a matrix $\textbf{Q}_k$ where k is the ID of the landmark, as shown in Figure 18. The matrix $\textbf{Q}$ is the sum of $\textbf{Q}_k$ as shown in equation (B6). This indicates that the Shur complement trick can be done point by point.

According to the Jacobian shown in figure 16, the 3D point 1 is seen in the camera 1, 2 and 3, the 3D point 2 is seen in the camera 2, 3 and 4, and the 3D point 3 is see in camera 3, 4 and 5, thus the actual structure of matrix $\textbf{Q}_1$, $\textbf{Q}_2$ and $\textbf{Q}_3$ are shown in Figure 20. The number and positions of non-zero blocks are determined by the number and identities of cameras which see this 3D point. For an instance, if a landmark k is seen in m cameras, then the number of non-zero blocks in $\textbf{Q}_k$ is $m^2$, and positions of non-zero blocks are determined by the identities of these cameras.

The formation of an RCS can be summarized as follows:

(1) For each 3D point k, compute its Jacobian and Hessian, and then conduct the Shur complement to obtain $\textbf{Q}_k$;

(2) Update the RCS with $\textbf{Q}_k$;

(3) Stop until all the 3D points are completed.

Fig. 19

The structure of the matrix $\textbf{Q}$ (left) and $\textbf{R}$ (right). The filled blocks are non-zero blocks and the unfilled ones are zero blocks

Fig. 20

The structure of matrix $\textbf{Q}_1$, $\textbf{Q}_2$ and $\textbf{Q}_3$ generated by landmarks 1, 2 and 3 respectively

According to the above procedures, the formation of RCS can be divided into parallel tasks, with each consisting of a group of 3D points. Those tasks are independent to each other; therefore, they can be executed in a distributed way.