Stochastic Bundle Adjustment for Efficient and Scalable 3D Reconstruction

In constrast to the previous methods [32, 43, 44, 18, 17, 40] that partition the problem in the pre-processing stage, we present a reformulation of problem (1) to decompose the RCS into clusters inside the LM iterations.

Particularly, we consider the most general case that every single camera forms a cluster. In order to preserve all the projections $\mathcal{E}$, we apply point splitting to the physical points, as shown in Fig. 2. For a physical point $\mathbf{p}_{i}$ viewed by $v_{i}$ cameras, we split it into $v_{i}$ virtual points $\{\mathbf{p}_{i}^{j}\}_{j=1}^{v_{i}}$, each assigned to one cluster. Such a clustering will reformulate problem (1) equivalently as a new constrained quadratic programming (QP) problem as below

Here, $\Delta\mathbf{x^{\prime}}=[\Delta\mathbf{c}^{T},\Delta\mathbf{p^{\prime}}^{T}]^{T}$ is an expansion of $\Delta\mathbf{x}$ which considers the update steps for all the virtual points, and so is the Jacobian $\mathbf{J}^{\prime}=[\mathbf{J}_{\mathbf{c}},\mathbf{J}_{\mathbf{p}}^{\prime}]$. Accordingly, $\mathbf{D}^{\prime}=\textrm{diag}(\mathbf{J}^{\prime T}\mathbf{J}^{\prime})^{\frac{1}{2}}$. The noteworthy distinction between problems (7) and (1) is that the new equality constraints of Eq. 8 are imposed to enforce that the steps of the same points in different clusters are identical. For example, $\Delta{\mathbf{p}_{i}^{s}}=\Delta{\mathbf{p}_{i}^{t}}$ $(\forall s,t\in\{1,...,v_{i}\})$ for point $\mathbf{p}_{i}$ and the corresponding j-th row of $\mathbf{A}$ appears in a form as $\mathbf{a}_{j}=[0...,1,...,-1,...0]$. Since a point $\mathbf{p}_{i}$ introduces $v_{i}-1$ equations, $\mathbf{A}$ has a row number of $r=\sum_{i=1}^{n}(v_{i}-1)$. Besides, $\mathbf{A}$ is full row rank, because the rows each of which defines a unique equality constraint are linearly independent.

The constrained QP problem above can be easily solved by Lagrangian duality, which turns the problem into

where $\mathbf{H}_{\lambda}=\mathbf{J^{\prime}}^{T}\mathbf{J^{\prime}}+\lambda\mathbf{D}^{\prime T}\mathbf{D}^{\prime}$ and $\boldsymbol{\nu}=-(\mathbf{A}\mathbf{H}_{\lambda}^{-1}\mathbf{A}^{T})^{-1}\mathbf{A}\mathbf{H}_{\lambda}^{-1}\mathbf{J^{\prime}}^{T}\mathbf{f}$ are the Lagrangian multipliers. Eq. 9 is in the similar format to Eq. 3, but has an additional term $\mathbf{A}^{T}\boldsymbol{\nu}$ on the right hand side compared with Eq. 3. While $\mathbf{J^{\prime}}^{T}\mathbf{f}$ includes the gradients w.r.t. the independent virtual points, $\mathbf{A}^{T}\boldsymbol{\nu}$ acts as a correction term to ensure that the solution complies with the constraints of Eq. 8.

The ultimate benefit of the clustering-based reformulation is revealed below. By likewise applying Schur complement to Eq. 9, we have

where $\mathbf{E}^{\prime}=\mathbf{J}_{\mathbf{c}}^{T}\mathbf{J^{\prime}_{p}}$, $\mathbf{C^{\prime}}=\mathbf{J^{\prime}_{p}}^{T}\mathbf{J^{\prime}_{p}}+\lambda\textrm{diag}(\mathbf{J^{\prime}_{p}}^{T}\mathbf{J^{\prime}_{p}})$, $\mathbf{S}^{\prime}=\mathbf{B}-\mathbf{E}^{\prime}\mathbf{C}^{\prime-1}\mathbf{E}^{\prime T}$, and $[\mathbf{v}^{\prime T},\mathbf{w}^{\prime T}]^{T}=-(\mathbf{J^{\prime}}^{T}\mathbf{f}+\mathbf{A}^{T}\boldsymbol{\nu})$. Due to the fact that any two cameras do not share any common virtual points, $\mathbf{S}^{\prime}$ now becomes a block-diagonal matrix, i.e., $\mathbf{S}^{\prime}_{ij}=\mathbf{0},\forall i\neq j$. Then Eq. 10 can be equivalently decomposed into $m$ most elementary linear systems each corresponding to one camera.

The major problem with the clustering based reformulation above is the excessive cost of evaluating the Lagrangian multipliers $\boldsymbol{\nu}$, because it requires the evaluation of $\mathbf{H}_{\lambda}^{-1}$. In order to make the problem practically solvable, we would like to eliminate the need to evaluate the correction term $\mathbf{A}^{T}\boldsymbol{\nu}$ of Eq. 9 by means of relaxation for problem (7).

We multiply a random binary variable $\theta_{i}$ with the i-th equality constraint of Eq. 8, which results in $\theta_{i}\mathbf{a}_{i}\Delta\mathbf{x}^{\prime}=0$. It could be interpreted that, if $\theta_{i}=1$, the constraint $\mathbf{a}_{i}\Delta\mathbf{x}^{\prime}=0$ must be satisfied; otherwise, the constraint is allowed to be violated. In this way, Eq. 8 will be relaxed into chance constraints [25], which leads to

where $\alpha\in(0,1]$ is a predefined confidence level. It means that, instead of enforcing the hard constraints, we allow them to be satisfied with a probability above $\alpha$. The advantage of the chance constrained relaxation is that we can determine the reliability level of approximation by controlling $\alpha$. The larger $\alpha$ is, the closer the chance constrained problem (11) will be to the original deterministic problem (7). One approach to problem (11) is called sampled convex program [8, 9]. It extracts $N$ independent samples $\Theta(\alpha)=\{\theta_{i}^{(n)}|i=1,...,r,n=1,...,N\}$ with a minimum sampling probability of $\alpha$ for each variable $\theta_{i}^{(n)}$ and replaces the chance constraints (12) with the sampled ones. Below we will elaborate on how problem (11) can be solved given the samples $\Theta(\alpha)$.

For a sample $\theta_{i}^{(n)}\in\Theta(\alpha)$, if $\theta_{i}^{(n)}=0$, the equality constraint $\mathbf{a}_{i}\Delta\mathbf{x}^{\prime}=0$ is dropped; if $\theta_{i}^{(n)}=1$, it enforces the equality of the steps of two virtual points, e.g., $\Delta\mathbf{p}_{j}^{s}=\Delta\mathbf{p}_{j}^{t}$. Here the virtual point $\mathbf{p}_{j}^{s}$ belongs to the single-camera cluster of camera $\mathbf{c}_{s}$ and similarly $\mathbf{p}_{j}^{t}$ to $\mathbf{c}_{t}$. Then we merge $\mathbf{p}_{j}^{s}$ and $\mathbf{p}_{j}^{t}$ into one point as shown in Fig. 2, and we call the operation point binding as opposed to point splitting introduced in Sec. 4.1. On the one hand, the point binding leads to the consequence that the Lagrangian multiplier $\nu_{i}=0$, because the equality $\Delta\mathbf{p}_{j}^{s}=\Delta\mathbf{p}_{j}^{t}$ always holds. On the other hand, the block $\mathbf{S}^{\prime}_{st}$ of $\mathbf{S}^{\prime}$ (c.f. Eq. 10) becomes nonzero, since the merged point is now shared by cameras $\mathbf{c}_{s}$ and $\mathbf{c}_{t}$. After applying point binding to all the virtual points involved in the sampled constraints, i.e., $\{\theta_{i}^{(n)}\in\Theta(\alpha)|\theta_{i}^{(n)}=1\}$, all the constraints will be eliminated and there is no need to evaluate $\boldsymbol{\nu}$ any more. Meanwhile, it will bring the cameras sharing common points into the same clusters.

Since the cameras in different clusters have no points in common after the point binding, the matrix $\mathbf{S}^{\prime}$ will appear in a block-diagonal structure which we call cluster-diagonal, as illustrated in Fig. 2. It means that each diagonal block of $\mathbf{S}^{\prime}$ corresponds to a camera cluster. In particular, we can intentionally design the sampler $\Theta(\alpha)$ in order to shape $\mathbf{S}^{\prime}$ into the desired cluster-diagonal structures, as we will present in Sec. 4.4. As a result, this structure of $\mathbf{S}^{\prime}$ still enables the decomposition of Eq. 10 into smaller independent linear systems each relating to one camera cluster. Provided that there are $l$ clusters, Eq. 10 can be equivalently re-written as

where $[\Delta\mathbf{c}_{1}^{T},...,\Delta\mathbf{c}_{l}^{T}]^{T}=\Delta\mathbf{c}$ and $[\mathbf{b}_{1}^{T},...,\mathbf{b}_{l}^{T}]^{T}=\mathbf{v}^{\prime}-\mathbf{E^{\prime}}\mathbf{C^{\prime}}^{-1}\mathbf{w}^{\prime}$. After Eqs. 13 are evaluated, we substitute $\Delta\mathbf{c}$ into Eq. 9 to give

where $\mathbf{J^{\prime}_{p}}=[\mathbf{J}_{\mathbf{p}^{1}},...,\mathbf{J}_{\mathbf{p}^{l}}]$ and $\Delta\mathbf{p}^{\prime}=[\Delta\mathbf{p}^{1},...,\Delta\mathbf{p}^{l}]$ include the Jacobians and virtual point steps w.r.t. the $l$ clusters respectively and we omit $\mathbf{D}^{\prime}$ for ease of notation. To give a uniform step for a physical point, we equalize the steps of its virtual points in different clusters by solving the linear system below in place of Eq. 14:

Since $\sum_{i=1}^{l}\mathbf{J}_{\mathbf{p}^{i}}=\mathbf{J_{p}}$, Eq. 15 gives the same solution as the point steps in Eq. 6: $\Delta\mathbf{p}=\mathbf{C}^{-1}(\mathbf{w}-\mathbf{E}^{T}\Delta\mathbf{c})$.

So far we have presented how an update step of the camera and point parameters is determined approximately by STBA. Besides this, we keep the other components of the LM algorithm unchanged [30]. For reference, we detail the full algorithm in the supplementary material.

The chance constrained relaxation in the last section effectively decomposes the RCS, but leads to approximate solutions with decreased feasibility due to the random constraint sampling. Below, we provide an empirical analysis on the effect of the approximation and present a conditional correction step to remedy the approximation error.

The LM algorithm is known to be the interpolation of the Gauss-Newton and gradient descent methods, depending on the trust region radius controlled by the damping parameter $\lambda$. When $\lambda$ is small and the LM algorithm behaves more like the Gauss-Newton method, the approximation induced by STBA in Sec. 4.2 is admissible, in that problem (1) itself is derived from the first order approximation of the error function $\mathbf{f}(\mathbf{x})$. And the LM algorithm can automatically contract the trust region when the approximation leads to the increase of the objective.

When $\lambda$ is large, i.e., the trust region is small, the LM algorithm is closer to the gradient descent method, which gives a step towards the steepest descent direction defined by the right hand side of Eq. 9, i.e., $-(\mathbf{J^{\prime}}^{T}\mathbf{f}+\mathbf{A}^{T}\boldsymbol{\nu})$. However, a problem with STBA is that the correction term $\mathbf{A}^{T}\boldsymbol{\nu}$ is eliminated approximately by the chance constrained relaxation. As a consequence, the derived step would deviate from the steepest descent direction and thus hamper the convergence. Therefore, we propose to recover $\mathbf{A}^{T}\boldsymbol{\nu}$ to remedy the deviation in such a case. Especially, when $\lambda$ is large enough, the matrix $\mathbf{H}_{\lambda}$ in Eq. 9 will be dominated by the diagonal terms, so that we can approximate $\mathbf{H}_{\lambda}$ by $\textrm{diag}(\mathbf{H}_{\lambda})$. After that, $\mathbf{A}^{T}\boldsymbol{\nu}=-\mathbf{A}^{T}(\mathbf{A}\mathbf{H}_{\lambda}^{-1}\mathbf{A}^{T})^{-1}\mathbf{A}\mathbf{H}_{\lambda}^{-1}\mathbf{J^{\prime}}^{T}\mathbf{f}$ can be evaluated efficiently because of the sparsity of $\mathbf{A}$. Since the approximation $\mathbf{H}_{\lambda}\approx\textrm{diag}(\mathbf{H}_{\lambda})$ is not accurate unless $\lambda$ is large, in practice, we enable the steepest descent correction particularly when $\lambda\geq 0.1$. Fig. 3 visualizes the effect of the correction.

The chance constrained relaxation in Sec. 4.2 necessitates an effective random constraint sampler $\Theta(\alpha)$. Among many of the possible designs, we propose a viable implementation named stochastic graph clustering in this section.

The design of the clustering method considers the following requirements. First, the sampler should be randomized with respect to the chance constraints (12). Since (12) indicates that the expectation $\mathrm{E}(\theta_{i})$ should have $\mathrm{E}(\theta_{i})\geq\alpha$ ($i=1,...,r$), the upper bound of the confidence level $\alpha$ is defined as $\min_{i=1}^{r}\mathrm{E}(\theta_{i})$. Therefore, the sampler should sample as many constraints as possible on average to increase the upper bound of $\alpha$. Second, the random sampler is intended to partition the cameras into small independent clusters so that Eqs. 13 can be solved efficiently.

Concretely, the stochastic graph clustering operates over a camera graph $\mathcal{G}_{c}=(\mathcal{C},\mathcal{E}_{c})$, where the weight $w_{ij}$ of an edge $e_{ij}\in\mathcal{E}_{c}$ between cameras $\mathbf{c}_{i}$ and $\mathbf{c}_{j}$ is equal to the number of points covisible by the two cameras. At the beginning, each camera forms an individual cluster as formulated in Sec. 4.1. Next, if $\mathbf{c}_{i}$ and $\mathbf{c}_{j}$ are joined, a number of $w_{ij}$ pairs of virtual points viewed by $\mathbf{c}_{i}$ and $\mathbf{c}_{j}$ will be merged. Therefore, $w_{ij}$ equality constraints will be satisfied. In order to join as many virtual points as possible while yielding a cluster structure of $\mathcal{G}_{c}$, we aim at finding a clustering that maximizes the modularity below inspired by [7]: $Q=\frac{1}{2s}\sum_{e_{ij}\in\mathcal{E}_{c}}\delta(\nu_{i},\nu_{j})\left(w_{ij}-\frac{k_{i}k_{j}}{2s}\right),$ where $s=\sum_{e_{ij}\in\mathcal{E}_{c}}w_{ij}$ is the total sum of edge weights, $k_{i}=\sum_{j}w_{ij}$ is the sum of weights of edges incident to camera $\mathbf{c}_{i}$, and $\nu_{i}$ denotes the cluster of $\mathbf{c}_{i}$. $\delta(\nu_{i},\nu_{j})=1$ if $\nu_{i}=\nu_{j}$ and $0$ otherwise. The modularity $Q\in[-1,1]$ measures the density of connections inside clusters as opposed to those across clusters [7]. Therefore, a larger modularity generally indicates that more virtual points are merged inside clusters. Maximizing the modularity is NP-hard [34], but Louvain’s algorithm [7] provides a greedy strategy which greedily joins the two clusters giving the maximum increase in modularity in a bottom-up manner. It can be efficiently applied to large graphs since its complexity is shown to be linear in the node number on sparse data [7].

However, Louvain’s algorithm [7] is deterministic due to its greedy nature. To ensure that every pair of virtual points is likely to be merged, we instead join clusters randomly according to a probability distribution defined based on the modularity increments [12], which is

where $\beta>0$ is a scaling parameter. Two neighboring clusters $N_{x}$ and $N_{y}$ are more likely to join together if it leads to a larger modularity increment $\Delta Q(N_{x},N_{y})$. In order to limit the sizes of the sub-problems of STBA, we stop joining clusters if their sizes exceed $\Gamma$. In Fig. 4, we visualize the stochastic clustering results.

Following previous works [15, 24], we evaluate the solvers with Performance Profiles over the total of 25 problems of 1DSfM [37] and KITTI [19]. We obtain the SfM results for 1DSfM by COLMAP [35]¹¹1Since one of the image sets Union Square has only 10 reconstructed images, we replace it with another public image set ArtsQuad. and the SLAM results for KITTI by stereo ORB-SLAM2 [31], while disabling the final bundle adjustment (BA). Since the SfM/SLAM results are generally accurate because the pipeline uses repeated BA for robust reconstruction, we make the problems more challenging by adding Gaussian noise to the points and camera centers following [17, 21, 32]. We report the number of images, tracks and projections of the datasets as well as the typical cluster number of STBA in Table 1 & 1.

First of all, we give a brief introduction of performance profiles [15]. Given a problem $p\in\mathcal{P}$ and a solver $s\in\mathcal{S}$, let $F(p,s)$ denote the final objective the solver $s$ attained when solving problem $p$. Then, for a number of solvers in $\mathcal{S}$, let $F^{*}(p)=\min_{s\in\mathcal{S}}F(p,s)$ denote the minimum objective the solvers $\mathcal{S}$ attained when solving problem $p$. Next, we define an objective threshold for problem $p$ which is $F_{\tau}(p)=F^{*}(p)+\tau(F_{0}(p)-F^{*}(p)),$ where $F_{0}(p)$ is the initial objective and $\tau\in(0,1)$ is the pre-defined tolerance determining how close the threshold is to the minimum objective. After this, we measure the efficiency of a solver $s$ by computing the time it takes to reduce the objective to $F_{\tau}(p)$, which is denoted by $T_{\tau}(p,s)$. And the most efficient solver is the one who takes the minimum time, i.e., $\min_{s\in\mathcal{S}}T_{\tau}(p,s)$.

The method Performance Profiles regards that the solver $s$ solves the problem $p$ if $T_{\tau}(p,s)\leq\alpha\min_{s\in\mathcal{S}}T_{\tau}(p,s)$, where $\alpha\in[1,\infty)$. Therefore, if $\alpha=1$, only the most efficient solver is thought to solve the problem, while if $\alpha\rightarrow\infty$, all the solvers can be seen to solve the problem. Finally, the performance profile of the solver $s$ is defined w.r.t. $\alpha$ over the whole problem set $\mathcal{P}$ as $\rho(s,\alpha)=100*\frac{|\{p\in\mathcal{P}|T_{\tau}(p,s)\leq\alpha\min_{s\in\mathcal{S}}T_{\tau}(p,s)\}|}{|\mathcal{P}|}.$ It is basically the percentage of problems solved by $s$ and is non-decreasing w.r.t. $\alpha$.

We plot the performance profiles of the solvers in Fig. 5. To verify the benefits of using stochastic clustering, herein we also compare STBA with its variant STBA-fixed which uses a fixed clustering as previous methods [24, 17, 40]. When $\tau$ is equal to 0.1 and 0.01, our STBA is able to solve nearly 100% of the problems for any $\alpha$, because it always reaches the objective threshold $F_{\tau}(p)$ with less time than LM and DL methods by a factor of more than 5. This is mainly attributed to the reduced per-iteration cost as we can see from the convergence curves in Fig. 6. When $\tau$ becomes 0.001 and the threshold $F_{\tau}(p)$ is harder to achieve, STBA is less efficient but still performs on par with LM-iterative and LM-sparse when $\alpha<3$ and better than DL for any $\alpha$. On the contrary, the performance of STBA-fixed drops drastically when $\tau$ decreases to 0.001. The performance change for STBA and STBA-fixed when $\tau$ decreases is mainly caused by the fact that the clustering methods come with the price of slower convergence near the stationary points [17, 40]. Compared with the full second-order solvers such as LM and DL, clustering methods only utilize the second order information within clusters. However, as opposed to STBA-fixed which uses fixed clustering, STBA has mitigated the negative effect of clustering by introducing stochasticity so that different second-order information can be utilized to boost convergence as the clustering changes. Despite the slower convergence rate, the benefit of STBA that it can reduce most of the loss with the lowest time cost (e.g., 99% loss reduction with only 1/5 time of the counterparts when $\tau=0.01$ in Fig. 5(b)) is still supposed to be highlighted, especially for the real-time SLAM applications where bundle adjustment is called repeatedly to correct the drift.

To evaluate the scalability, we conduct distributed experiments on the Large-Scale (LS) dataset, which includes the urban scenes of four cities named LS-1, LS-2, LS-3 and LS-4. We run a distributed SfM program of our own to produce initial sparse reconstructions of the four scenes and add Gaussian noise with a standard deviation of 3 meters to the camera centers and points. Then we compare our distributed STBA against the state-of-the-art distributed bundle adjustment framework DBACC [40].

We visualize the sparse reconstructions and the convergence curves of the four scenes in Fig. 7 and report the statistics in Table 2. As we can see from Fig. 7, STBA achieves faster convergence rates than DBACC by an order of magnitude. The main cause of the gap is that DBACC, which is based on the ADMM formulation [5], has to take inner iterations to solve a new minimization problem in every ADMM iteration. Although we have set the maximum inner iteration number to merely 10, DBACC still takes many more Jacobian and RCS evaluations and thus has much longer iterations than STBA by an order of magnitude, as shown in Table 2. Besides, as opposed to DBACC, our STBA is free of too many hyper-parameters.

Here we perform an ablation study on steepest descent correction proposed in Sec. 4.3 to validate its efficacy. We compare our STBA with its variant called STBA* which does not use the correction. We run STBA, STBA* and LM-sparse on 1DSfM and KITTI. Since steepest descent correction is designed particularly for a small trust region, we set the lower bound of the damping parameter $\lambda$ to 0.1 in the experiments. We observe that by using the correction, STBA consistently achieves a faster convergence than STBA* and performs on par with LM-sparse on all the scenes. Visualizations of the sample convergence curves w.r.t. the iterations are shown in Fig. 8, where STBA and LM-sparse have very close convergence curves. It manifests that steep descent correction indeed facilitates the correction of the approximation errors of the STBA steps and hence boosts the convergence.