{'fraction_non_alphanumeric': 0.04264332454341244, 'fraction_numerical': 0.013404629358335775, 'mean_word_length': 4.3114382051614575, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient and Consistent Bundle Adjustment on Lidar Point Clouds

Bundle Adjustment (BA) refers to the problem of simultaneous determination of sensor poses and scene geometry, which is a fundamental problem in robot vision. This paper presents an efficient and consistent bundle adjustment method for lidar sensors. The method employs edge and plane features to represent the scene geometry, and directly minimizes the natural Euclidean distance from each raw point to the respective geometry feature. A nice property of this formulation is that the geometry features can be analytically solved, drastically reducing the dimension of the numerical optimization. To represent and solve the resultant optimization problem more efficiently, this paper then proposes a novel concept {\it point clusters}, which encodes all raw points associated to the same feature by a compact set of parameters, the {\it point cluster coordinates}. We derive the closed-form derivatives, up to the second order, of the BA optimization based on the point cluster coordinates and show their theoretical properties such as the null spaces and sparsity. Based on these theoretical results, this paper develops an efficient second-order BA solver. Besides estimating the lidar poses, the solver also exploits the second order information to estimate the pose uncertainty caused by measurement noises, leading to consistent estimates of lidar poses. Moreover, thanks to the use of point cluster, the developed solver fundamentally avoids the enumeration of each raw point (which is very time-consuming due to the large number) in all steps of the optimization: cost evaluation, derivatives evaluation and uncertainty evaluation. The implementation of our method is open sourced to benefit the robotics community and beyond.


I. INTRODUCTION
Light detection and ranging (lidar) has become an essential sensing technology for robots to achieve a high level of autonomy [1,2]. Enabled by the direct, dense, active and accurate 1  (DDAA) depth measurements, lidar sensors have the ability to build a dense and accurate 3D map of the environment in realtime and at a relatively low computation cost. These unique advantages have made lidar sensors essential to a variety of applications that require real-time, dense, and accurate 3D mapping of the environment, such as autonomous driving [3,4], unmanned aerial vehicles navigation [5]- [7], and realtime mobile mapping [8]- [10]. This trend becomes even more evident with recent developments in lidar technologies which have enabled the commercialization and mass production of lightweight and high-performance solid-state lidars at a significantly lower cost [11,12]. The central task of many lidar-based techniques, such as lidar-based odometry, simultaneous localization and mapping (SLAM), and multi-lidar calibration, is to register multiple point clouds, each measured by the lidar at different poses, into a consistent global point cloud map. However, the predominant point cloud registration methods, such as iterative closest point (ICP) [13] and its variants (e.g, generalized-ICP [14]), normal distribution transformation (NDT) [15,16], and surfel registration [17], allow registration of two point clouds only. Such a pairwise registration leads to an incremental scan registration process for an odometry system (e.g., [17]- [20]), which would rapidly accumulate drift, or a repeated pairwise registration process for 3D mapping [21] or multilidar calibration [22], which would bring dramatic computation cost. All these necessitate an efficient concurrent multiple scan registration technique.
Concurrent multiple lidar scan registration requires determining all lidar poses and the scene geometry simultaneously, a process referred to as bundle adjustment (BA) in computer vision. Compared to visual BA, which has been well-established in photogrammetry and played a fundamental role in various vital applications, including visual odometry (VO) [23]- [25], visual-inertial odometry (VIO) [26,27], 3D visual reconstruction [28,29] and multi-camera calibration [30,31], lidar BA has a similarly fundamental role but is much less mature due to two major challenges. First, lidar has a long measuring range but low resolution between scanning lines. The measured point cloud are sparsely (sometimes even not repeatedly [12]) distributed in a large 3D space, making it difficult (almost impossible) to scan the same point feature in the space across different scans. This has fundamentally prevented the use of straightforward visual bundle adjustment formulation, which is largely based on point features benefiting from the high-resolution images accurately capturing individual point features. The second challenge lies in the large number of raw points (from tens of thousands to million points) collected by arXiv:2209.08854v1 [cs.RO] 19 Sep 2022 a practical lidar sensors. Processing all these points in the lidar BA is extremely computation intensive.
In this work, we propose an efficient and consistent BA framework specifically designed for lidar point clouds. The framework follows our previous work BALM [32], which formulates the lidar BA problem based on edge and plane features that are abundant in lidar scans. The BA formulation naturally minimizes the straightforward Euclidean distance of each point in a scan to the corresponding edge or plane, while the decision variables include the lidar poses and feature (edge and plane) parameters. Furthermore, it is shown that the geometry parameters (i.e., edge and plane) can be solved analytically, leading to an optimization that depends on the lidar poses only. Since the number of geometry features is often large, elimination of these geometry features from the optimization will drastically reduce the optimization dimension (hence time).
A novel design in our proposed BA framework is the point cluster, which summarizes all points of a lidar scan associated to one feature by a compact set of parameters, point cluster coordinates. Based on the point cluster, we derive the closed-form derivatives (up to second order) of the BA optimization with respect to (w.r.t.) its decision variables (i.e., lidar poses). We prove that the formulated BA optimization and the closed-form derivatives can both be represented fully by the point cluster without enumerating the large number of individual points in a lidar scan. The removal of dependence on individual raw points drastically speeds up the evaluation of the cost function and derivatives, which further enables us to develop an efficient and consistent secondorder solver, BALM2.0, which is also released on Github to benefit the community. Our experiment video is available on https://youtu.be/MDrIAyhQ-9E.
We conduct extensive evaluations on the proposed BA method. Simulation study shows that the BA method produces consistent lidar pose estimate. Exhaustive benchmark comparison on 19 real-world open sequences shows that the BA method produces consistently higher performance (pose estimation accuracy, mapping accuracy, and computation efficiency) than other counterparts. We finally integrate the BA method in three vital lidar applications: lidar-inertial odometry, multi-lidar calibration, and global mapping, and show how their accuracy and/or computation efficiency are improved by the proposed BA.

A. Multi-view registration
The bundle adjustment problem is similar to the multi-view registration problem that has been previously researched [33]- [38]. These methods all adopt a two layer framework: the first layer estimates the relative poses of a selected set of scan pairs using the pairwise registration methods (e.g., ICP [13]); From the relative poses, the second layer constructs and solves a pose graph to obtain a maximum a posteriori estimate of all lidar poses. Such a two-layer framework decouples the raw point registration from the global pose estimation, so that each raw point registration only involves a small amount of local points contained in the two scans (instead of all scans sharing overlaps) and the pose graph optimization only involves a small amount of constraints arising from the relative poses (instead of raw points). The net effect is a significant saving of time, hence being largely used in online lidar SLAM systems [39,40]. However, the advantage in computation efficiency comes with fundamental limitation in accuracy: the pairwise scan registration only considers the overlap among two scan at a time, while the overlap is really shared by all scans and should be registered concurrently. Moreover, the pose graph optimization only considers constraints from the relative poses, while the mapping consistency indicated by the raw points are completely ignored. Consequently, it is usually difficult to produce (or even be aware of) a globally consistent map that is necessary for high-accuracy localization and mapping tasks.
Some early works in computer vision and computer graphics have proposed multi-view registration methods that directly optimize the mapping consistency from multiple range images, aiming for consistent surface modeling of 3D objects. [41] is a direct extension of the ICP method, it minimizes the Euclidean distance between a pre-known control point in one scan to all matched control points in the rest scans. Within this framework, [42] uses a quaternion representation in the optimization, and [43] extends the distance between control points to the distance between surfaces around the respective control points. More recently, Zhu et al. [44] proposes a two step registration method: the first step uses a K-means clustering to cluster points from all scans, and the second step estimates the scan poses by minimizing the Euclidean distance between each point in a cluster to the centroid. Since these methods rely on point features in the scan, they require densely populated point cloud (e.g., by depth camera) for extracting such salient point features. While this is not a problem for small object reconstruction for which these methods are designed, it is not the case for scene reconstruction where the LiDAR measurements are very sparse (sometimes even nonrepetitive) as explained above.

B. Bundle or plane adjustment
In recent years, researchers in the robotics community have shown increasing interests to address the bundle adjustment problem on (lidar) point clouds more formally. Kaess [45] exploits the plane features in the bundle adjustment and minimizes the difference between the plane measured in a scan and the plane predicted from the optimization variables: scan poses and plane parameters. This formulation was later integrated into a key-frame-based online SLAM system [46]. Since the method minimizes the plane-to-plane distance, it requires to segment each scan and estimate the contained local planes in advance. Such plane segmentation and estimation usually require dense point clouds measured by RGB-D cameras on which the work were demonstrated.
A more formal bundle adjustment method on lidar point cloud, termed as the plannar (bundle) adjustment, was later proposed in [47] which minimizes the natural Euclidean distance between each point in a scan to the plane predicted from the scan poses and plane parameters (the optimization variables). Compared with plane-to-plane distance in [45], the point-to-plane metric is faster, more accurate, and more suitable for lidar sensors, where local plane segmentation or estimation are less reliable due to sparse point clouds. Moreover, the direct use of raw points in the point-to-plane metric could also lead to a more consistent estimate of the optimization variables by considering the measurement noises in the raw points. Then, the formulated non-linear least square problem is solved by a Levenberg-Marquardt (LM) algorithm. To lower the computation load caused by the large number of points measurements associated to the same plane feature, [47] propose a reduction technique to eliminate the enumeration of individual points in the evaluation of the residual and Jacobian. Furthermore, due to the very similar structure to the visual bundle adjustment, the proposed bundle adjustment is also compatible with the Schur complement trick [48], which eliminates the plane parameters in each iteration of the LM algorithm. This plane adjustment method is largely used in many online lidar SLAM systems developed subsequently [49,50] or before [51,52].
On the other hand, Ferrer [53] exploits plane features similar to [47] and minimizes any deviation of each raw point from the plane equation. The resultant optimization cost then reduces to the minimum eigenvalue of a covariance matrix and is thus termed as the eigen-factor. The author further derived the closed-form gradient of the cost function w.r.t. to both the scan poses and plane parameters and employed a gradientbased method to solve the optimization iteratively. Due to the second order nature of the eigenvalue (as confirmed in [32], see below), the gradient method converges very slowly (requiring a few hundreds of iterations) [50,53].
Our previous work BALM [32] takes another step towards more efficient bundle adjustment. Similar to [47], BALM minimizes the natural Euclidean distance between each point in a scan to the plane (i.e., point-to-plane metric). Most importantly, BALM proved that the point-to-plane distance is essentially equivalent to the eigen-factor in [53], thus unifying the two metrics in [47] and [53]. It was further proved that all plane parameters can be analytically solved with closedform solutions in advance, hence the large number of plane parameters can be completely removed from the resultant optimization. Such feature elimination significantly lowered the optimization dimension and poses a fundamental difference from all the previous plannar adjustment methods [47,49]- [51,53] and visual bundle adjustment methods [48]. The feature elimination also removed the various issues in the plane representation from the optimization, such as normal constraint in the Hesse normal representation (n, d) [47,49,50,53], singularity issue in the closest-point (CP) representation nd [51] and over-parameterization issue in the quaternion representation [45]. A further advantage of BALM against previous method [47,51,53] is that the whole framework is naturally extendable to edge features besides plane features. To solve the resultant optimization problem, BALM [32] subsequently derived the second order derivatives of the cost function and developed a LM-like second-order solver. The developed BA method and solver are very efficient, achieving real-time sliding window optimization when integrated to LOAM [18].
A major drawback of the BALM [32] is that the evaluation of the second-order derivatives including Jacobian and Hessian requires to enumerate each individual lidar point, leading to a computational complexity of O(N 2 ) where N is the number of points [50]. Consequently, the method is hard to be used in large-scale problems where the lidar points are huge in number. This problem is partially addressed in [54], which aggregates all points associated to the same plane feature in a scan in the scan local frame. However, to ensure convergence, [54] modifies the cost function by including an extra heuristic penalty term, which is not a true representation of the map consistency. Moreover, the cost function in [54] still involves the plane feature similar to [47,49]- [51,53]. To lower the computation load caused by optimizing the large number of feature parameters, the method further fixes the feature parameters in the optimisation, which could slow down the optimisation speed.
Our BA formulation in this paper is based on BALM [32], hence inheriting the fundamental feature elimination advantage when compared to [47,49]- [51,53]. To address the computational complexity of O(N 2 ) in BALM, we propose a novel point cluster concept, which fundamentally eliminates the enumeration of each individual point in the evaluation of the cost function, Jacobian and Hessian matrix. Consequently, the computational complexity is irrelevant to both the feature dimension (similar to BALM [32]) and the point number (similar to [47,49]). The point cluster in our method is similar to the point aggregation used in [54], but the overall BA formulation is fundamentally different: (1) it minimizes the true map consistency (the point-to-plane distance) without trading off with any other heuristic penalty; and (2) it performs exact feature elimination with rigorous proof instead of empirical fixation. Based on these nice theoretical results, we develop an efficient second-order solver, termed as BALM2.0. Besides solving the nominal lidar poses, the solver also estimates the uncertainty of the estimated lidar pose by leveraging the second-order derivative information, which is another new contribution compared with existing works.

III. BUNDLE ADJUSTMENT FORMULATION AND
OPTIMIZATION In this chapter, we derive our BA formulation and optimisation. First, following [32], we formulate the BA as minimizing the the point-to-plane (or point-to-edge) distance (Sec. III-A) and show that the feature parameters can be eliminated from the formulated optimisation (Sec. III-B). Then, we introduce the point cluster in Sec. III-C, based on which the first and second order derivatives are derived in Sec. III-D. Based on these theoretical results, we present our second-order solver in Sec. III-E. Then the time complexity of the proposed BA is analyzed in the Sec. III-F . Finally, in Sec. III-G, we show how to estimate the uncertainty of the BA solution. Throughout this paper, we use notations summarized in Table I or otherwise specified in the context.

A. BA formulation
Shown in Fig. 1, assume there are M f features, each denoted by parameter π i (i = 1, ..., M f ), observed by M p The set of m × m symmetric matrices. The encapsulated "boxplus" operations on manifold.
The value of (·) expressed in lidar local frame, (·) The value of (·) expressed in global frame.
The number of features and poses, respectively. i, j, k The indexes of features, poses and points, respectively. l The index of eigenvalue and eigenvector of a matrix. p,q The indexes of (block) row and column in a matrix. e l The vector in R 4 with all elements being zeros except the l-th element being one (l ∈ {1, 2, 3, 4}). lidar poses, each denoted by T j = (R j , t j ) (j = 1, ..., M p ), the bundle adjustment refers to simultaneously determining all the lidar poses (denoted by T = (T 1 , · · · , T Mp )) and feature parameters (denoted by π = (π 1 , · · · , π M f )), such that reconstructed map agrees with the lidar measurements to the best extent. Denote c(π i , T) the map consistency due to the i-th feature, a straightforward BA formulation is In our BA formulation, we make use of plane and edge features that are often abundant in lidar point cloud and minimize the natural Euclidean distance between each measured raw lidar point and its corresponding plane or edge feature. Specifically, assume a total number of N ij lidar points are measured on the i-th feature at the j-th lidar pose, each denoted by p f ijk (k = 1, ..., N ij ). Its predicted location in the global frame is For a plane feature, it is parameterized by π i = (n i , q i ) with n i the plane normal vector and q i an arbitrary point on the plane, both in the global frame (see Fig. 2 (a)). Then, the Euclidean distance between a measured point p f ijk to the plane is n T i (p ijk − q i ) 2 . Aggregating the distance for all points observed in all poses leads to the total map consistency corresponding to this plane feature: where N i = Mp j=1 N ij is the total number of lidar points observed on the plane feature by all poses.
For an edge feature, it is parameterized by π i = (n i , q i ) with n i the edge direction vector and q i an arbitrary point on the edge, both in the global frame (see Fig. 2 (b)). Then, the Euclidean distance between a measured point p f ijk to the edge is Aggregating the distance for all points observed in all poses leads to the total map consistency corresponding to this edge feature: where N i = Mp j=1 N ij is the total number of lidar points observed on the edge feature by all poses.

B. Elimination of feature parameters
In this section, we show that in the BA optimization (1), the feature parameter π can really be solved with a closed-form solution. The key observation is that one cost item c(π i , T) depends solely on one feature parameter, so that the feature parameter can be optimized independently. Concretely, In case of a plane feature, we substitute (3) into c(π i , T): where λ l (A i ) denotes the l-th largest eigenvalue of matrix A i , u l (A i ) denotes the corresponding eigenvector, the matrix A i and vectorp i are defined as: The proof will be given in Supplementary II-A [55]. Note that the optimal solution q i in (6) is not unique, any deviation from q i along a direction perpendicular to n i will equally serve the optimal solution. However, these equivalent optimal solution will not change the plane nor the optimal cost (hence the results that follow next). Indeed, the point q i could be an arbitrary point on the plane as it is defined to be.
In case of an edge feature, we substitute (4) into c(π i , T): Again, the optimal solution q i in (8) is not unique, any deviation from q i along the direction n i will equally serve the optimal solution. However, these equivalent optimal solution will not change the edge nor the optimal cost (hence the results that follow next).
As can be seen from (6) and (8), the parameter π i for each feature, either it is a plane or edge, can be analytically solved and hence removed from the BA optimization process. Consequently, the original BA optimization in (1) reduces to where l ∈ {2, 3} and we omitted the exact number of eigenvalues in the cost function for brevity. Note that the matrix A i in (9) depends on the lidar pose T since each involved point p ijk depends on the pose (see (7) and (2)). Hence the decision variables of the resultant optimization in (9) involve the lidar pose T only, which dramatically reduces the optimization dimension (hence computation time).

C. Point cluster
With the feature parameters eliminated, another difficulty remaining in the BA optimization (9) is that the evaluation of matrix A i (and its Jacobian or Hessian necessary for developing a numerical solver) requires to enumerate every point observed at each lidar pose. Such an enumeration is extremely computationally expensive due to the large number of points in a lidar scan. In this section, we show such point enumeration can be avoided by a novel concept point cluster, which is detailed as follows.
A point cluster is a finite point set denoted by set C = {p k ∈ R 3 |k = 1, · · · , n}, the corresponding point cluster coordinate, denoted as (C), is defined as: where S 4×4 denotes the set of 4 × 4 symmetric matrix.
A point cluster can be thought as a generalized point, for which a rigid transform could be applied. Similarly, we can define rigid transformation on a point cluster as follows.
Besides rigid transformation, we also define cluster merging operation, as follows.
Definition 2. (Cluster merging) Given two point clusters with point collections C 1 = {p 1 k ∈ R 3 |k = 1, · · · , n 1 } and C 2 = {p 2 k ∈ R 3 |k = 1, · · · , n 2 }, respectively. The merged cluster, denoted by C 1 ⊕ C 2 , is defined as Next we will show that the two operations defined above can be fully represented by their point cluster coordinates.
Theorem 2. Given two point clusters C 1 and C 2 . The merged cluster satisfies Proof. See Supplementary II-C [55].  As can be seen, rigid transformation and cluster merging operations on point clusters can be represented by usual matrix multiplication and addition on the point cluster coordinates. A visual illustration of the two operations and their coordinate representations are shown in Fig. 3. These results are crucially important: Theorem 1 indicates that the point cluster can be constructed in one frame (e.g., local lidar frame) and transformed to another (e.g., the global frame) without enumerating each individual points; Theorem 2 indicates that two (and by induction more) point clusters can be further merged to form a new point cluster. A particular case of Theorem 2 is when the second point cluster contains a single point, indicating that the point cluster can be constructed incrementally as lidar points arrives sequentially.
Remark: A point set and its coordinate is not a one-to-one mapping. While it is obvious that the coordinate is uniquely determined from the point set as shown in (10), the reverse way does not hold: the point set cannot be recovered from its coordinate uniquely. Since different point sets may lead to the same coordinate, the raw points must be saved if a reclustering is needed. Now, we apply the point cluster to the BA problem concerned in this paper. To start with, we group all points on the same feature as a point cluster. For example, the point cluster for the i-th feature is C i {p ijk |j = 1, · · · , M p , k = 1, · · · , N ij }. Denote C i the coordinate of the point cluster, following (10), we obtain A key result we show now is that this point cluster coordinate C i is completely sufficient to represent the matrix A i required in the BA optimization (9). According to (7), we have where we denote A i as a function of C i since the A i is fully represented (and uniquely determined) by C i . We slightly abuse the notation here by denoting the function as A(·).
On the other hand, since the point set C i is defined on points in the global frame, the coordinate C i is dependent on the lidar pose, which remains to be optimized. To explicitly parameterize the lidar pose, we note that where p f ijk is a point represented in the lidar local frame (see (2)), C ij is the point set that is composed of all points on the i-th feature (either plane or edge) observed at the j-th lidar pose, and C fij is the same point set as C ij , but represented in the lidar local frame (see Fig. 4).
The relation between the point cluster C i and C fij shown in (20) will lead to a relation between their coordinates C i and C fij as below: As a result, the BA optimization in (9) further reduces to where the function A(·) is defined in (16). Note that the cost function in (22) only requires the knowledge of point cluster coordinate C fij without enumerating each individual points. The coordinate C fij is computed as (following (10)): which can be constructed during the feature association stage before the optimization. In particular, if the j-th pose observes no point on the i-th feature, C fij = 0 4×4 , which will naturally remove the dependence on the j-th pose for the i-th cost item as shown in (22).
Theorem 3 implies that left multiplying all poses T j , ∀j, by the same transform T 0 does not change the optimization at all. That is, the BA optimization is invariant to the change of the global reference frame, which is the well-known gauge freedom in a bundle adjustment problem.

D. First and Second Order Derivatives
As shown in the previous section, the BA optimization problem in (22) is completely equivalent to the original formulation (1), where each cost item standards for the squared Euclidean distance from a point to a plane (or edge). This squared distance is essentially a quadratic optimization, which requires the knowledge of the second-order information of the cost function for efficient solving. In this section, we derive such first and second order derivatives. Without loss of generality, we only discuss the i-th feature, which contributes a cost item in the form of with C fij ∈ R 4×4 being a pre-computed matrix (see (23)).
To derive the derivative of the cost item (25) w.r.t. the pose T j , which is an element of the Special Euclidean group SE(3), we parameterize its perturbation by a special addition, called boxplus ( -operation). For the pose vector T = (· · · , T j , · · · ), we define the the operation as below where δT (· · · , δT j , · · · ) ∈ R 6Mp with δT j (δφ j , δt j ) ∈ R 6 , ∀j, is the perturbation on the pose vector.
For a scalar function f (T) : (T 0 ) its first-order derivative and ∂f 2 (T) ∂T 2 (T 0 ) its second-order derivative, both at a chosen point of the input T 0 . The operation enables us to parameterize the input of function f (·), T, by its perturbation δT from a given point T 0 : T = T 0 δT. Since the map between T and δT is bijective if δφ j < π, ∀j, the scalar function f (T) in terms of T can be equivalently written as a function f (T 0 δT) in terms of δT. As a consequence, the derivatives of f (T) w.r.t. T at the point T 0 can be defined as the derivatives of f (T 0 δT) w.r.t. δT at zero, where the latter is a normal derivative w.r.t. Euclidean vectors: In the following discussion, we use T as the reference point to replace T 0 in the derivatives and omit it for the sake of notation simplification. Based on the derivatives defined in (28) and (29), we have the following results for the first and second-order derivatives of the cost item (25).
(2) Poses T j ∈ SE(3), j = 1, · · · , M p ; eigenvalue of A with corresponding eigenvector u l , the Jacobian and Hessian matrix of the function λ l with respect to the poses T are where Moreover, the Jacobian J and Hessian matrix H satisfy: and Proof. See Supplementary II-E [55].

Remark 1:
The results in (39) essentially implies that the cost function (25) in the BA optimization does not change along the direction where all the pose are perturbed by the same quantity w, which agrees with gauge freedom stated in Theorem 3.

Algorithm 1: LM optimization
Input: Initial poses T; Point cluster in the local frame C fij ; Compute current cost c = c(T) and the new cost c = c(T ) from (22); Output: Final optimized states T;

Remark 2:
The results in (40) implies that the Jacobian and Hessian matrices are zeros and hence their computation can be saved if any of the related poses does not observe the current feature. This sparse structure could save much computation time if a feature is observed only by a sparse set of poses.
Remark 3: The derivatives in Theorem 4 are obtained based on the pose perturbation defined in (26), which multiplies the perturbation δT on the left of the current pose (i.e., a perturbation in the global frame). If other perturbation (denoted by δT) is preferred (e.g., a perturbation in the local frame to integrate with other measurements such as IMU pre-integration), where δT = LδT with L the Jacobian between the two perturbation parameterization, its first and second order derivatives can be computed asJ = J · L andH = L T · H · L, respectively. It can be seen thatJ andH preserves a nullspace of L · W with W defined in (39).

E. Second Order Solver
The Jacobin and Hessian matrix from Theorem 4 are computed for one cost item (25) that corresponds to one feature in the space. Denote J i , H i the Jacobian and Hessian matrix for the i-th feature (or cost item), to determine the incremental update ∆T, we make use of the second order approximation of the total cost function c(T) in (22): where , ∀i, we have Jd = 0 and d T Hd = 0, ∀i, which means that any additional update along d ∈ W does not change the approximation at all. One way to resolve the gauge freedom is fixing the first pose at its initial value throughout the optimization. That is, setting ∆T 1 = 0 in (41). Then, setting the differentiation of the cost approximation in (41) w.r.t. ∆T (excluding ∆T 1 ) to zero leads to the optimal update ∆T : where we used a Levenberg-Marquardt (LM) algorithm-like method to re-weight the gradient and Newton's direction by the damping parameter µ. The complete algorithm is summarized in Algorithm 1.

F. Time Complexity Analysis
In this section, we analyze the time complexity of the proposed BA solver in Algorithm 1 and compare it with other similar BA methods, including BALM [32] and Plane Adjustment [47]. The computation time of Algorithm 1 mainly consists of the evaluation of the Jacobian J and Hessian matrix H on Line 6, and the solving of the linear equation on Line 9. For the evaluation of the Jacobian J, according to (41), it consists of M f items, each item J i requires to evaluate M p block elements, and each block g j ll requires a constant computation time according to (33). Therefore, the time complexity for the evaluation of Similarly, for the Hessian matrix, according to (41), it consists of M f items, each item H i requires to evaluate M 2 p block elements, and each block H ij requires a constant computation time according to (35). Therefore Our previous work BALM [32] also eliminated the feature parameters from the optimization, leading to an optimization similar to (9) whose dimension of Jacobian and Hessian are also 6M p . To solve the cost function, BALM [32] adopted a second order solver similar to Algorithm 1, where the linear solver on Line 9 has the same time complexity O(M 3 p ). However, when deriving the Jacobian and Hessian matrix of the cost function, the chain rule was used where the Jacobian and Hessian is the multiplication of the cost w.r.t. each point of a feature and the derivative of the point w.r.t. scan poses (see Section III. C of [32]). As a consequence, for each feature, the evaluation of Jacobian has complexity of O(N M p ) and the Hessian has O(N 2 M 2 p ), where N is the average number of points on each feature. Therefore, the time complexity including the evaluation of all features' Jacobian and Hessian matrices and the linear solver are The plane adjustment method in [47] is a direct mimic of the visual bundle adjustment, which does not eliminate the feature parameters but optimizes them along with the poses in each iteration. The resultant linear equation at each iteration is in the form of J T Jδx = b, where J is the Jacobian of the cost function w.r.t. both the pose and feature parameters. Due to a reduced residual and Jacobian technique similar to our point cluster method, the evaluation of J does not need to enumerate each point of a feature, so the time complexity of computing , which is linear to M f , the number of feature, quadratic to N , the number of point, and cubic to M p , the number of pose. The plane adjustment [47] which is similar to [47] (i.e., linear to M f , irrelevant to N , and cubic to M p ) but has less operations.

G. Covariance Estimation
Assume the solver converges to an optimal pose estimate T , it is often useful to estimate the confidence level of the estimated pose. Let T gt be the ground-true pose and δT be the difference between the optimal estimate T and the ground-true T gt , where T gt = T δT . The aim is to estimate the covariance of the error δT , denoted by Σ δT .
Ultimately, the estimation error δT is caused by the measurement noise in each raw point. Denote p gt f ijk the ground-true location of the k-th point observed on the i-th feature at the jth lidar pose, with measurement noise δp f ijk ∈ N (0, Σ p f ijk ), the measured point location, denoted by p f ijk , is Aggregating the ground-true raw points and the measured ones lead to the ground-true point cluster, denoted by C gt f ijk , and the measured point cluster, denoted by C f ijk , respectively (see (23)): where which can be constructed in advance along with the point cluster C fij during the feature associations stage.
In the following discussion, to simplify the notation, we the ground-truth, measurements, and noises of all point clusters observed on any features at any lidar poses.
Although the ground-true pose T gt and point cluster C gt f are unknown, they are genuinely the optimal solution of (22) and hence the Jacobian evaluated there should be zero, i.e., where we wrote the Jacobian as an explicit function of the pose and point clusters. Now, we approximate the left hand side of (47) by its first order approximation: Noticing that which implies Since T is the the converged solution using the measured cluster C f , they should lead to zero update, i.e., We defer the exact derivation and results of [55]. Note that the evaluation of Σ δT only requires the covariance δC fij , which has been constructed in advance according to (46), avoiding the enumeration of each raw point during the optimization.

IV. IMPLEMENTATIONS
We implemented our proposed method in C++ and tested it in Unbuntu 20.04 running on a desktop equipped with Intel i7-10750H CPU and 16Gb RAM. Since the reduced optimization problem (22) is not in a standard least square problem, which existing solvers (e.g., Google Ceres [56]) applies to, we implemented the optimization algorithm with steps and parameters described in Algorithm 1. When solving the linear equation on Line 9 at each iteration, we use the LDLT Cholesky decomposition decomposition method implemented in Eigen library 3.3.7. The termination conditions on Line 19 are iteration number below 50 (i.e., j max = 50), rotation update below 10 −6 rad, and translation update below 10 −6 m.

V. CONSISTENCY EVALUATION
This study aims to verify the consistency of the proposed BA method. That is, whether the estimated covariance Σ δT from (52) agrees with the ground-true covariance of the pose estimation error δT . As the ground-true covariance is unknown, we refer to a standard measure of consistency, the normalized estimation error squared (NEES) [57,58], which is defined below: δT δT , where δT j is the estimation error of the pose defined according to (27): where the superscript "gt" denotes the ground-true poses and (R j , t j ) denotes the estimated pose for the j-th scan. Assume the pose estimate (R j , t j ) is unbiased (i.e., E δT j = 0), if the computed covariance Σ δT is the ground-truth, we can obtain the expectation That is, if the solver is consistent, the expectation of NEES should be equal to the dimension of the optimization variable. If the expectation of NEES is far higher than the dimension, the estimator is over-confident (i.e., the computed covariance is less than the ground-truth). Conversely, it is conservative.
In practice, the expectation of NEES is evaluated by Monte Carlo method, where the NEES is computed for many runs and then averaged to produce the empirical expectation.
where η (i) is the NEES computed at the i-th Monte Carlo run.
To conduct the Monte Carlo evaluation, we simulate a 16channel lidar along a rectangular trajectory in a cuboid semiclosed space shown in Fig. 5. The size of the space is 30  To simulate realistic measurements, each point is corrupted with an independent isotropic Gaussian noise with multiple standard deviations σ p ∈ {0.05, 0.10, · · · , 1.00} m and for each value of the standard deviation σ p , we performed 100 Monte Carlo experiments, leading to a total number of 2000 experiments. In each run, we compute the optimal pose estimate from Algorithm 1 with the same parameters specified in Section IV and the covariance matrix from (52). The initial trajectory required in Algorithm 1 is obtained by perturbing the ground-true trajectory with a Gaussian noise with standard deviation δφ = 2 deg and δt = 0.1 m on each pose.To avoid unnecessary errors, we use the ground-true plane association across different scans and ignore the in-frame motion distortion in the simulation. Fig. 6 shows orientation and position errors with the corresponding 3σ bounds in one Monte Carlo experiment with σ p = 0.05 m. As can be seen, the pose estimation errors are very small and they all remain within the 3σ bounds very well, which suggests that our new method is consistent.
Furthermore, we test the consistency of our BA method under different levels of point noise, where the standard deviation of a point noise ranges from σ p = 0.05 m to 1 m. The results are shown in Fig. 7(a) for the NEES averaged over 100 runs for each noise level and in Fig. 7(b) for the average pose error. For better visualization, the average NEES is normalized by the pose dimension (i.e., 600 for 100 poses on the trajectory) in Fig. 7(a). As can be seen, the normalized average NEES is very close to one, which suggests that our method is consistent, when the point noise is up to 0.3 m. Beyond this noise level, the first order approximation in (49) no longer holds, which undermines the accuracy of the computed covariance. We should note that this noise level rarely occurs in actual lidar sensors, which are well below 0.1 m. Moreover, from Fig. 7(b), we can see that our method produces accurate pose estimation even when the point noise are unrealistically large (up to 1 m, see Fig. 7(c) for the point cloud map at this point noise level).

VI. BENCHMARK EVALUATION
In this section, we compare our method with other multiview registration methods for lidar point clouds. The experiment will be divided into two parts: Section VI-A evaluates all methods with known feature association on virtual point clouds, and Section VI-B evaluates the overall BA pipeline including both optimization solver and feature association on various real-world open datasets.
To verify the effectiveness of our method, we compare it with four state-of-the-art methods that focus on the lidar bundle adjustment (or similar) problem: Eigen-Factor (EF) [53], BALM [32], Plane Adjustment (PA) [49] and BAREG [54]. Among them, EF 3 , BALM 4 , BAREG 5 are open sourced, so we use the available implementation on Github. PA is not available anywhere, so we re-implemented it in C++. To reduce the time cost of PA, we used the reduced Jacobian and residual technique in [47] (we derived it based on the cost function in Equation (10) of [49]), which avoids the enumeration of each individual point, and the Schur complement trick, which reduces the linear equation dimension at each optimization iteration. The re-implemented PA uses the same LM algorithm pipeline as in Algorithm 1 to make a fair comparison.
For the optimization solvers, our method and the reimplemented PA use the same parameters in Algorithm 1, while EF, BALM and BAREG use their default parameters as available on their open source implementation. All methods use the same termination condition shown in Sec. IV (i.e. maximal iteration number beow 50, rotation update below 10 −6 rad, and translation update below 10 −6 m), except for EF, which we found it converges too slowly and hence set the maximal iteration steps to 1000. In addition to the open source version of BALM (denoted by BALM), which samples only three points from each plane to lower the computation load, we also evaluated another vision (denoted by BALMfull) which keeps all the points on a plane and use the same default parameters as its open sourced version. All solvers use the same initial pose trajectories detailed later.
For the feature association, we use the same adaptive voxelization proposed in our previous work BALM [32], which registers all points in the world frame (using the initial trajectory) and recursively cuts the space into smaller voxels until the voxel contains a single feature that associates the points from multiple scans. EF did not address the feature association problem and PA did not open relevant association codes, so we use the method of BALM [32] too. BAREG used a adaptive voxelization method similar to BALM but has its own implementation, which we retain. All feature associations have the same set of parameters: the root voxel size L, the maximum voxelization layer l max , the minimum number of points n min for a feature test, and the feature test thresholds γ. Exact values of these parameters are provided later.
Our proposed optimization method and the adaptive voxelization method in BALM [32] are applicable to both edge and plane features. However, we notice that in real-world point cloud data, edge features extracted based on local smoothness (e.g., [18]) are very noisy because the laser pulse emitted by lidars can barely hit an edge exactly due to the limited angular resolution. The situation is even worsened when the edge is at far or when the lidar has increased laser beam divergence, which creates many bleeding points behind an edge and degrades the edge points extraction [59]. As a consequence, we found that in practice, adding the edge features does not contribute much to the accuracy. Moreover, since the other methods, including EF, PA, and BAREG are only designed for plane features, for a fair comparison, we use only plane features in the following experiments for all methods. In fact, an edge is often created by an foreground object, which also makes a good plane feature, so dropping the edge feature does not reduce the number of constraints significantly.  1) Convergence: First we investigate the convergence performance of these methods. Fig. 8(a) shows the cost convergence for all methods in one repeat experiment of the nominal settings (i.e., M f = M p = N = 40). Since different method uses different cost function, to compare them in one figure, the cost value of each method is normalized by its initial cost and then it is re-based such that the converged cost value of all methods are aligned at the same value. As can be seen, EF converes rather slowly and requires the most number of iterations. This is because EF uses a cost function similar to ours in (9), which is essentially a quadratic optimization, but uses only the gradient information for optimization. Indeed, slow convergence of the gradient descent method on a quadratic cost function is a very typical phenomenon [60]. PA and BAREG converge fast at the beginning but slowly when approaching the final convergence value. This is because PA optimizes both the plane parameters and scan poses, leading to a very large number of optimization variables that significantly slow down the speed at convergence. For BAREG, the empirical fixation of plane parameters also causes the optimization to slow down. In contrast, BALM, BALM-full and our method eliminates the plane parameters exactly and the resultant optimization problem is only in dimension of the pose number. Further leveraging the exact Hessian information in their optimization update, BALM, BALM-full and our method converge in four iterations, which represent the fastest convergence. Fig. 8(b) shows the iteration experienced by each BA optimization method, where for each data point, the y-axis value represents how many experiments out of the total (13×3 = 130) experienced a certain iteration number indicated by the x-axis value. As can be seen, the overall trend agrees with the results in Fig. 8 (b) very well: our proposed method and BALM-full require only four or five iterations in all experiments, while BAREG and PA require up to 20 and 50 iterations, respectively. EF requires even more iterations beyond 50. The only exception is BALM, which experienced more iterations than BAREG in the worst case, the reason is that BALM samples only three points from each plane feature, which may reduce the robustness and increase the iteration steps in some cases.

A. Virtual point cloud
2) Accuracy: Fig. 9(a)-(f) shows the statistic values of the pose estimation accuracy in terms of RMSE. In each subplot, we fix two parameters of M f , M p and N at the nominal values and change the third parameter to investigate the change of pose accuracy. Since the error of EF is much larger than the others, we used a broken y-axis to better display all the RMSE. As can be seen, the accuracy overall increases with the points per plane N (in (b) and (e)) or number of plane features M f (in (c) and (f)), both increases the number of pose constraints. No such monotonic accuracy improvement is found in (a) and (d) as the pose number increases because increasing the pose number itself does not gives more pose constraints. Relatively speaking, our proposed method and BALM-full achieves the same highest accuracy, since they essentially optimizes the same cost using the same exact Hessian information. The next best methods are BAREG and PA. While optimizing the same point to plane distance with our method (and BALM-full), PA has significantly more optimization variables, which cause a much slower convergence where the solution is still slightly premature at the preset iteration number (i.e., 50). The next accurate method is BALM, which samples only three instead of all points (as in our method, BALM-full, BAREG and PA) and hence has higher RMSE. Finally, EF has the highest RMSE due to the very slow convergence, the solution is much premature even at the preset iteration number (i.e., 1000).
3) Computation time: Finally, we show the computation time of different solvers at different feature number M f , pose number M p , and point number N . The results are shown in Fig. 9(g)-(i), respectively. As can be seen, the time consumption of all methods increases with the number of poses M p (see (g)) and plane features M f (see (i)), which is reasonable since more poses or planes leads to a higher optimization dimension or more number of cost items, respectively. In contrast, as the point number N increases, the method BALM-full increases rapidly since its time complexity involves O(N 2 ) while the rest methods (including ours) do not increase notably since they do not need to evaluate every raw point. Relatively speaking, our method achieves the lowest computation time in all cases due to the small number of iteration numbers and low time-complexity per iteration. The next efficient method is BAREG, which has very low timecomplexity per iteration due to the empirical feature parameter fixation but significantly more iteration numbers due to the same reason. Compared with our method, PA has the same time complexity per iteration as discussed in Section III-F, but requires more iterations to converge. Hence its time cost is a little higher than ours and BAREG. BALM requires more iterations than our method and each iteration, it requires to enumerate the sampled three points. Collectively, it leads to a computation time higher than our method and also BAREG and PA. The slow convergence problem is more severe in EF, leading to an even higher computation time. Finally, the most time-consuming method is BALM-full, which, although has very small iteration numbers, takes significantly high time to compute the Hessian matrix in each iteration due to enumeration of each raw point.

B. Real-world datasets
In this experiment, we conduct benchmark comparison on three real-world datasets. The first dataset is "Hilti" [61] which is a handheld SLAM dataset including indoor and outdoor environments. We use the lidar data collected by Ouster OS0-64 in the dataset. The ground-true lidar pose trajectory is captured by a total station or motion capture ssytem. The second dataset "VIRAL" [62] is collected on an unmanned aerial vehicle (UAV) equipped with two 16-channel OS1 lidars. One lidar is horizontal and the other is vertical. We will use the horizontal one in this experiment. It used a Leica Nova MS60 MultiStation to track a crystal prism on the UAV to provide ground-true positions. The last dataset "UrbanLoco" [63] is collected by a car driving on an urban streets. The lidar is Velodyne HDL 32E and the ground truth is given by the Novatel SPAN-CPT, a navigation system incorporating Real Time Kinematic (RTK) and precisional IMU measurements.
Two preprocessing are performed for all sequences: motion compensation and scan downsample. To compensate the points distortion caused by continuous lidar movements within a scan, we run a tightly-coupled lidar-inertial odometry, FAST-LIO2 [19], which estimates the IMU bias (and other state variables) and compensates the point motion distortion in realtime. We kept all points in a scan whose distortion has been compensated by FAST-LIO2 and discard the odometry output. The processed data are then downsampled from the original 10 Hz in all sequences to 2 Hz. This is because the BA methods need to process all scans at once, a 10Hz scan rate causes prohibitively high computation load for all BA methods. The downsampling is also similar to the keyframe selection in common SLAM frameworks.
We compare our method with EF, BALM, PA, and BAREG. Noticing that the computation time of BALM-full is prohibitively high due to the extremely large number of lidar points, we hence remove it from the benchmark comparison. For the rest methods, their solver parameters are kept the same for all sequences with values detailed before. For feature association, the parameters are l max = 3, n min = 20, γ = 1 25 . For the parameter L, it is set to 1 m for "Hilti" and "VIRAL" and 2 m for "UrbanLoco".
In addition to the multi-view registration methods, we also compare with classic pairwise registration methods, including ICP, GICP, and NDT offered in PCL library. We run the pairwise registration methods in an incremental manner, where each new scan is registered and merged to the map incrementally. To constrain the computation time, in each new scan registration, only the last 20 scans are kept in the map. The pose estimation from the ICP is then used as the initial trajectory for the feature association and optimization of the BA methods, including EF, BALM, PA, BAREG, and ours. 1) Accuracy: Table II shows the ATE results. As can be seen, our method consistently achieves the best results in all 19 sequences. The next accurate method is PA and BAREG, followed by BALM and EF. This trend is in great agreement with the results on virtual point-cloud shown in Section VI-A-2 with explanations detailed therein. In particular, our method achieves an accuracy within a few centimeters in all sequences of "Hilti" and "VIRAL", with only one exception (i.e., UzhArea2), which will be analyzed later. The centimeter level accuracy achieved by our method is at the same level of lidar point noises. Moreover, using only lidar measurements, our method achieved an average accuracy of 4.2 cm on all VIRAL dataset sequences, which outperforms the accuracy 4.7 cm reported in VIRAL-SLAM [62] that fuses all data from stereo camera, IMU, lidar, and UWB. The accuracy on "UrbanLoco" is lower (analyzed later) than other datasets, but still outperforms the other BA methods. Finally, we can notice that the BA methods (i.e., EF, BALM, PA, BAREG and ours) generally outperforms the pairwise registration methods (i.e., ICP, GICP and NDT) due to the full consideration of multi-view constraints.
Now we investigate the performance degradation on "Ur-banLoco" and the sequence UzhArea2 in "Hilti" more closely. For the "UrbanLoco" dataset, we found that the RTK groundtruth had some false sudden jumps, which contributes the large ATEs. This sudden jump may be caused by tall buildings in the crowded urban area which lowers the quality of the groundtruth. For the sequence UzhArea2, we register the point cloud with the ground-true pose trajectory and compare it with the point cloud registered with our BA method in Fig. 10. As can be seen, with the ground-true pose, points on the side wall are very blurry and points on the wall form a plane with standard deviation up to 15.3 cm (Fig. 10(d)); with the ground-true translation but with rotations optimized by our BA method, the points on the side wall are much more ordered and form an apparent plane of standard deviation 6.8 cm (Fig. 10(e)); with poses fully optimized by our method, the points are even more consistent and the standard deviation is 1.7 cm (Fig.  10(f)). From these results, we suspect that the ground-truth may be affected by some unknown errors (e.g., marker position change during the data collection). Indeed, we found similar problem on this sequence also occurred in other works [64]. Moreover, the standard deviation of 1.7cm achieved by our method is exactly the ranging accuracy of the lidar sensor, which confirms that our method achieves a mapping accuracy at the noise level of raw points as if the sensor motion had not occurred.
2) Mapping quality: A significant advantage of the BA method is the direct optimization of the map accuracy. To evaluate the map quality without a ground-true map, we adopt a method proposed by Anton et al. [65]. The method cuts the space into small cells and then counts the number of cells that lidar points occupy. The less the occupied cells, the higher the map quality. This indicator is intuitive: if points from different scans are registered accurately, they should agree with each other to the best extent, hence occupying the minimum possible number of cells. Based on this method, Table III presents the number of occupied cell with size 0.1 m. To better show the difference among different methods, the number of occupied cells are subtracted by our method for each sequence. We show the number of occupied cells by our method and the difference value of other methods. As can be seen, our method consistently achieved the best performance in all sequences and the next best is PA and BAREG. This trend also agrees with the ATE results very well.
3) Computation time: Finally, we compare the computation time. Since the pairwise registration methods, including ICP, GICP, and NDT, perform repetitive incremental registration at each scan reception, its computation time is very different from the BA methods that perform batch optimization on all scans at once. Therefore, we only compare the computation time of BA methods. Table IV shows the total computation time of optimization. As can be seen, our method consumes the least computation time, about one sixth of the BAREG, one eighth of PA, one tenth of BALM, and one seventieth of EF. The overall trend agree with the results on vitual point cloud in Section VI-A-3 with explanations detailed therein.

VII. APPLICATIONS
Bundle adjustment is the central technique of many lidarbased applications. In this section, we show how our bundle adjustment method can effectively improve the accuracy or computation efficiency of three vital applications: lidar-inertial odometry, multi-lidar calibration, and global mapping.
A. Lidar-inertial odometry with sliding window optimization A local bundle adjustment in a sliding window of keyframes has been widely used in visual odometry and proved to be very effective in lowering the odometry drift [26,27,66,67]. Similar idea could apply to the lidar-inertial odometry based  11. Overview of the LIO system with BA and IMU preintegration on our BA method. To demonstrate the effectiveness, we design a lidar-inertial odometry system as shown in Fig. 11. The system is divided into two parts: the EKF-based frontend, which provides initial yet timely pose estimation (EKF odometry), and the BA-based back-end, which refines the pose estimation in a sliding window (Local mapping). The EKF design is similar to FAST-LIO2 [19] which compensates the motion distortion in the incoming lidar scans and performs EKF propagation and update on the state manifold. The sliding window optimization performs a local BA among the most recent 20 scans considering constraints of IMU preintegration [58] and constraints from plane features co-visible among all scans in the window (Section III). After the convergence of the local BA, points in the local window are built into a k-d tree to register the next incoming scan. Furthermore, the oldest scan is removed from window and its contained points are merged into the global point cloud map. We compare our proposed system with one state-of-art lidarinertial odometry, FAST-LIO2 [19], which performs incremen- tal pairwise scan registration via GICP. The comparison is conducted on "utbm" , "uclk", and "nclt" dataset that evaluated by FAST-LIO2, so the results of FAST-LIO2 are directly read from the original paper [68]. As can be seen in Table V, benefiting from the abundant multi-view constraints in the local sliding window, our system consistently outperforms FAST-LIO2 in terms of accuracy with considerable margins except the sequence "nclt4". The improvements in odometry accuracy confirm the effectiveness of the local BA optimization. The improvement in accuracy comes with increased computation costs as shown in the last two columns of Table V, where for FAST-LIO2, we record the time of each scan-to-map registration and for our LIO system, the time consumption is divided into two parts: one is the front-end scan-to-local map registration and the back-end local BA optimization. As can be seen, our system has a little more time consumption in frontend than FAST-LIO2 because of the building of a k-d tree, which takes a constant time overhead, while FAST-LIO2 uses a more efficient incremental k-d tree structure. In addition, our system requires an average of 73 ms in the back-end local BA optimization, which ensures a 10 Hz running frequency for both the front-end and back-end. Since the back-end runs in parallel in a separate thread, the overall odometry latency of our system is only 7 ms more than that of FAST-LIO2.

B. Multiple-lidar Calibration
With the ability of concurrent optimization of multiple lidar poses, our BA method can be applied to multi-lidar extrinsic calibration. We consider the problem in [69], which aims to calibrate the extrinsic of multiple solid-state lidars shown in Fig. 12. Due to the very small FoV, these lidars have very small or even no FoV overlap. To create co-visible features, the vehicle is rotated for one cycle, during which a set of point cloud scans are collected by all lidars. Due to the rotation, it further introduces unknown lidar poses, in addition to the extrinsic, to estimate.
To formulate the simultaneous localization and extrinsic calibration problem, we choose a reference lidar as the base and calibrate the extrinsic of the rest lidars relative to it. Denote G R T j , j = 1, · · · , M p , the pose of the reference lidar at j-th scan in a rotation, R L k T the extrinsic of the k-th lidar relative to the base lidar. Then, following (25), the cost item corresponding to the i-th plane feature is where C r fij , C k fij are the point cluster of the reference lidar and the k-th lidar, respectively, observed on the i-th feature at the j-th scan. Following Theorem 4, the first and second order derivatives of (54) can be obtained by chain rules, whose details are omitted here due to limited space.
We demonstrate the advantage of the proposed BA approach in multi-lidar calibration with the latest state-of-theart calibration methods based on ICP [22] and BALM [69]. The ICP-based method [22] optimizes the extrinsic parameter R L k T and base lidar pose G R T j by repeatedly registering the point cloud by each lidar at each scan to the rest lidar points until convergence, whereas in both BALM-based method [69] and our method, the extrinsic parameter R L k T and base lidar pose G R T j are optimized concurrently. In [22] the feature correspondence searching is conducted by a k-d tree data structure, while in [69] and our method, this is resolved by adaptive voxelization proposed in BALM [32]. To restrain the computation time, the point number in each voxel of [69] have been down-sampled to four, whereas in our method, all feature points are used.
We test these methods on two lidar setups in the system shown in Fig. 12: lidars with small field-of-view (FoV) overlap (MID-100 self calibration) and lidars without FoV overlap (AVIA and MID-100 calibration). Since the extrinsic between the internal lidars of MID-100 (i.e., L1 L0 T, L1 L2 T) is known fromt he manufacturer's in-factory calibration, it is thus served as the ground truth in both setups. Moreover, in the MID-100 self calibration setup, the middle L 1 is chosen as the reference lidar to calibrate the extrinsic of adjacent lidars, i.e., L1 L0 T, L1 L2 T. This lidar setup is tested with data collected in both [22] and [69] under two scenes which contributes to  Fig. 13. Extrinsic calibration results of the ICP-based [22], BALM-based [69] and our proposed method in two experiment setups (with small or without FoV overlap). eight calibration results. In the MID-100 and AVIA calibration setup, the L 3 is selected as the reference lidar to calibrate the extrinsic of each internal lidar of MID-100, i.e., L3 L0 T, L3 L1 T and L3 L2 T. We then calculate the relative pose L1 L0 T, L1 L2 T from the calibrated results and compare them with the ground truth. This lidar setup is tested with data collected under two scenes in [69] since no AVIA lidar was used in [22], which contributes another four calibration results.
The total twelve independent calibration results are illustrated in Fig. 13 and Fig. 14. It is seen the our proposed method outperforms the ICP-based method [22] especially in translation. This is due to the concurrent optimization of all poses and extrinsic at the same time, which leads to full convergences within a few iterations. In contrast, the repetitive pairwise ICP registration leads to very slow convergence. The optimization did not fully converge after the maximum iteration number (max iter=40), a phenomenon detailed in [69]. Furthermore, since all raw points on a feature are utilized in our method, the accuracy of our proposed work is also slightly improved when compared with the BALM-based method [69], which sample a few points on a feature to restrain the computation time. The averaged time consumption in each step of these methods have been summarized in Table VI. It is seen our proposed work and [69] has significantly shortened the calibration time in each step compared with the ICP-based method. This is due to the use of adaptive voxelization which saves a great amount of time in k-d tree build and nearest neighbor search used in [22]. Compared with [69], our proposed BA has further increased the computation speed due to the use of point cluster technique avoiding the enumerating of individual points in the BALM-based method [69].

C. Global BA on Large-Scale Dataset
In this section, we show that our proposed BA method could also be used to globally refine the quality of a large-scale lidar point cloud. It is becoming popular to use pose graph  optimization (PGO) as the back-end to enhance the overall lidar odometry precision in SLAM system [70]. The optimal poses are obtained by minimizing the summed error between the relative transformation of these poses and that estimated in the front-end odometry. One drawback of the PGO is that it cannot reinforce the map quality directly. A divergence in point cloud mapping would occur if the estimation of the relative poses near the loop is ill. We show that the mapping quality (and odometry accuracy) could be further improved with our BA method even after PGO is performed. We choose to validate on three sequences from the KITTI dataset [71]. The initial lidar odometry is provided by the stateof-the-art SLAM algorithm MULLS [70] with loop closure function enabled. We directly feed this odometry to our BA algorithm and concurrently optimize the entire poses. The optimization results (RMSE of the odometry error) are summarized in Table VII. It is seen that even with loop closure function, some divergence still exist in the area near the loop of the point cloud due to ill estimations. With our proposed global BA refinement, the accuracy in odometry is further improved and the divergence in point cloud is eliminated (see closeup in Fig. 15).

VIII. DISCUSSION
Here we discuss the efficiency, accuracy, and extendability of the proposed bundle adjustment method.

A. Efficiency
Our method achieved lower computation time than other state-of-the-art counterparts. The efficiency of our method are attributed to three inter-related and rigorously-proved techniques that make fully use of the problem nature and lidar point cloud property. The first technique is the solving of feature parameters in a closed-form before the BA optimization. It allows the feature parameters to be removed from the optimization, which fundamentally reduces the optimization dimension to the dimension of the pose only, a phenomenon that did not exist before in visual bundle adjustment problem. The second technique is a second-order solver which fits the quadratic cost function naturally and leads to fast convergence in the iterative optimization. This is enabled by the analytical derivation of the closed-form Jacobian and Hessian matrices of the cost function. The third technique is the point cluster, which enables the aggregation of all raw points without enumerating each individual point in neither of the cost evaluation, derivatives evaluation, or uncertainty evaluation. Collectively, these three techniques lead to an BA optimization with much lower dimension and a time complexity in each iteration, which achieves the lowest time consumption when compared with other second order methods (e.g., BALM [32] and plane adjustment [47]).

B. Accuracy
Benefiting from the point cluster technique, our proposed method is able to exploit the information of all raw point measurements, achieving high pose estimation accuracy (a few centimeters) at the level of lidar measurement noise. Optimization from the raw lidar points also enables the developed method to estimate the uncertainty level of the estimated pose, which may be useful when this information is further fused with measurements from other sensors (e.g., IMU sensors). Moreover, by minimizing the Euclidean distance from each raw point to the corresponding feature, our method can reinforce the map consistency in a more direct manner than conventional pose graph optimization. While at a higher computation cost (due to the more complete consideration of features co-visible in multiple scans), it considerably improves the mapping accuracy which is important for mapping applications. Due to this reason, our method is particularly useful for accuracy refinements from a baseline pose trajectory that can be obtained by an odometry or a pose graph optimization module. The second order optimization provides very fast convergence when the solution is near to the optimal value, preventing premature solutions.

C. Extendability
As a basic technique for multiple scan registration, our proposed method can be easily be integrated with other formality of data, such as images and IMU measurements, by incorporating visual bundle adjustment factors and IMU preintegration factors [58] in the optimization. Moreover, besides the frame-based pose trajectory, which attaches each frame an independent pose to estimate, our method can also work with other forms of pose trajectories, such as continuoustime trajectories based on Splines [9,72] or Gaussian Process models [73,74], which have the capability to compensate the in-frame motion distortion. According to the chain rules, the derivatives of the BA cost with respect to the trajectory parameters will consist of two parts: the first is the derivative of the BA cost with respect to the pose of each point cluster as derived in this paper, and the second part is the derivatives of the pose with respect to the trajectory parameters, which depends on the specific trajectories being used.

IX. CONCLUSION
This paper proposed a novel bundle adjustment method for lidar point cloud. The central of the proposed method is a point cluster technique, which aggregates all raw points into a compact set of parameters without enumerating each individual point. The paper showed how the bundle adjustment problem can be represented by the point cluster and also derived the analytical form of the Jacobian and Hessian matrices based on the point cluster. Based on these derivations, the paper developed a second-order solver, which estimates both the pose and the pose uncertainty. The developed BA method is open sourced to benefit the community.
Besides the technical developments, this paper also made some theoretical contributions, including the formalization of the point cluster and its operations, revealing of the invariance property of the formulated BA optimization, the proof of null space and sparsity of the derived Jacobian and Hessan matrices, and the time complexity analysis of the proposed BA method and its comparison with others. These theoretical results serve the foundation of our developed BA techniques.
The proposed methods and implementations were extensively verified in both simulation and real-world experiments, in terms of consistency, efficiency, accuracy, and robustness. In all evaluations, the proposed method achieved consistently higher accuracy while consuming significantly lower computation time. This paper further demonstrated three applications of the BA techniques, including lidar-inertial odometry, multilidar calibration, and high-accuracy mapping. In all applications, the adoption of BA method could effectively improve the accuracy or the efficiency.
In the future, we would like to incorporate the BA method more tightly to the above applications and beyond. This would require more thorough considerations of many practical issues, such as point cloud motion compensation, removal of dynamic objects, tightly-fusion with other formality of sensor data (e.g., IMU, camera) and module (e.g., loop closure). I. LEMMAS Lemma 1. For a scalar x ∈ R and a matrix A ∈ S 3×3 which depends on x, we have the two following conclusions.
where λ l (l = 1, 2, 3) denotes the l-th largest eigenvalue and u l is the corresponding eigenvector.
Proof. Since the matrix A(x) is symmetric, its singular value decomposition is, where Λ(x) = diag (λ 1 (x), λ 2 (x), λ 3 (x)) consists of all the eigenvalues and U(x) = u 1 (x) u 2 (x) u 3 (x) is an orthonormal matrix consisting of the eigenvectors. Therefore, Both sides take the derivative of x, Since Denote D(x) ≜ U(x) T ∂U(x) ∂x . Since U(x)U(x) T = I, differentiating both sides with respect to x leads to, It is seen that D(x) is a skew symmetric matrix whose diagonal elements are zeros. Moreover, since Λ(x) is diagonal, the last two items of the right side of (6) sum to zero on diagonal positions. Only considering the diagonal elements in (6) leads to which yields the first conclusion. Now we aims to prove the second one. In (6), ∂Λ(x) ∂x is diagonal matrix and thus for the off-diagonal, k-th row, l-th column, element (k ̸ = l), where D k,l x is the k-th row, l-th column element in the skew symmetric D(x) and satisfy D k,l Since D(x) ≜ U(x) T ∂U(x) ∂x , we have ∂U(x) ∂x = U(x)D(x). Taking the l-th column on both sides leads to where D :,l x ∈ R 3 represents the l-th column of D(x). Finally, substituting U(x) = [u 1 (x) u 2 (x) u 3 (x)] and (9) into (10), we obtain which yields the second conclusion.

Lemma 2. Given
(2) Poses T j ∈ SE(3), j = 1, · · · , M p ; and a matrix function A(C) ≜ 1 N P − 1 N 2 vv T ; (4) Two constant eigenvectors u k , u l ∈ R 3 , then the first and second order derivatives of u T k A(T)u l w.r.t. T are: where g q kl ∈ R 1×6 , H p,q kl ∈ R 6×6 , ∀p,q ∈ {1, · · · , M p }, are and Proof. Partition the matrix C as then Therefore, Since where we omitted the input argument of A for the sake of notation simplicity. Since N = Mp j=1 N j is a constant number that is irrelevant to the pose T, perturbing the input of A, which is the pose vector T, by δT yields Based on the definition of ⊞ on SE(3) in (27), we define, When δϕ q is small, which is indeed the case for the purpose of derivative computation, we have Substituting (32) into w ql (δT q ), we obtain wherew ql , δΦ ql and J w ql can be simplified Substituting (33) into (30) and keeping terms up to the second order lead to To compute Q q kl in 1 2 δT T q Q q kl δT q , note that a T ⌊δϕ⌋ 2 b = δϕ T ⌊a⌋⌊b⌋δϕ, ∀a, b, δϕ ∈ R 3 , we havē Similarly, Summing (37) and (38) and extending the δϕ into δT: For 1 2 δT T q N p,q kl δT q and 1 2 δT T p M p,q kl δT p , where D q,p = D T p,q . Thus, extending the δϕ into δT, we obtain It can be seen that (36) is quadratic w.r.t. δT, so we cast it into the following standard form u T k A(δT)u l = 1 2 δT T · H kl · δT + g kl · δT + r kl , (44) where g kl is g kl = · · · g q kl · · · ∈ R 1×6Mp (45) where D p,q = C p FC q and F = 0 3×1 0 3×1 0 1×3 1 . Since The common part in (47) is Therefore, Additionally, it is seen that if C q = 0 4×4 , g q kl = 0 1×6 , ∀k, l. For H kl in (44), partition it as where By using similar method in (48), we have The value of Q q kl is in (39). Thus, the matrix H p,q kl is which yields the solution. Additionally, it can be seen that if C q = 0 4×4 or C p = 0 4×4 , we have D p,q = C p FC q = 0 4×4 , K q kl = 0 3×3 and then H p,q kl = 0 6×6 , ∀k, l.

II. PROOF OF THEOREMS
A. Proof of formula (6) and (8) Proof. The variable to be optimized is π i = (n i , q i ) and the cost function is where h i = n i for plane feature and h i = (I − n i n T i ) for edge feature. The dimensions of h i may be different for these two features but it has no influence on following derivation.
As can be seen, the inner optimization on q i is a standard quadratic optimization problem. So, the optimum q ⋆ i can be solved by setting the derivative to zero: where N i = To solve for the optimal parameter n i in the above optimization problem, we discuss the case of plane and edge features separately, as follows.
1) Plane feature: h i = n i .
where A i is defined in (7) and is a symmetric matrix. Performing Singular Value Decomposition (SVD) of A i where with where the minimum value λ 3 is reached when m 3 = 1, i.e., m ⋆ = 0 0 1 T and n ⋆ i = U i m ⋆ = u 3 . Therefore, the optimal cost is λ 3 (A i ) and the optimum solution is n ⋆ = u 3 (A i ) and q ⋆ =p i .
2) Edge feature: h i = I − n i n T i , then h i h T i = h i and Therefore, the optimal cost is λ 2 (A i ) + λ 3 (A i ) and the optimum solution is n ⋆ = u 1 (A i ) and q ⋆ =p i .

B. Proof of Theorem 1
For the point collections C = {p k ∈ R 3 |k = 1, · · · , n}, its point cluster is The rigid transformation of C by pose T = R t 0 3×1 1 is whose point cluster is which yields the solution. □

C. Proof of Theorem 2
For two point collections C 1 = {p 1 k ∈ R 3 |k = 1, · · · , n 1 } and C 2 = {p 2 k ∈ R 3 |k = 1, · · · , n 2 }, their point clusters are The merge of C 1 and C 2 is whose point cluster is which yields the solution. □ D. Proof of Theorem 3 Since A T 0 CT T 0 and A(C) are similar by transformation R 0 , they have the same eigenvalue. □ Thus, Since Au l = λ l u l , Therefore, For the proof of δT T HδT = 0, from (97) where We will divide (109) into two parts to discuss. For the first part, and Then, from (103) and Au l = λ l u l , Mp q=1 K q ll = ⌊Au l ⌋⌊u l ⌋ + ⌊u l ⌋⌊Au l ⌋ = 2λ l ⌊u l ⌋ 2 .
Now, substituting the results in (114) and (116) Therefore, where due to ⌊u k ⌋u l = −⌊u l ⌋u k and u T k ⌊u l ⌋ = −u T l ⌊u k ⌋ r kl = λ l u T l ⌊u k ⌋ + λ k u T k ⌊u l ⌋) T λ l u T l ⌊u k ⌋ + λ k u T k ⌊u l ⌋) = − λ 2 l ⌊u k ⌋u l u T l ⌊u k ⌋ − λ l λ k ⌊u k ⌋u l u T k ⌊u l ⌋ − λ 2 k ⌊u l ⌋u k u T k ⌊u l ⌋ − λ l λ k ⌊u l ⌋u k u T l ⌊u k ⌋ = − (λ l − λ k ) 2 ⌊u l ⌋u k u T k ⌊u l ⌋.

F. Derivation of pose covariance
The quantity δC fij can be obtained by substituting (43) into the definition of C gt fij in (44) and retaining only the first order items: To derive δC f with involving any tensor, we parameterize the matrix C fij by a column vector c fij , which consists of the independent elements in C fij : c fij = vec(C fij ) ≜ e T 1 C fij e 1 e T 1 C fij e 2 e T 1 C fij e 3 e T 2 C fij e 2 e T 2 C fij e 3 e T 3 C fij e 3 e T 1 C fij e 4 e T 2 C fij e 4 e T 3 C fij e 4 T ∈ R 9 (132) where vec(·) : S 4×4 → R 9 maps a symmetric matrix to its column vector representation, e l ∈ R 4 (l ∈ {1, 2, 3, 4}) is a vector with all zero elements except for the l-th element being one. Note that the constant N in the 4-th row, 4-th column of C fij is not contained in c fij since it is a constant number independent of the noise. Correspondingly, noises in C fij becomes the noise of c fij as below: δc fij = e T 1 δC fij e 1 e T 1 δC fij e 2 e T 1 δC fij e 3 e T 2 δC fij e 2 e T 2 δC fij e 3 e T 3 δC fij e 3 e T 1 δC fij e 4 e T 2 δC fij e 4 e T 3 δC fij e 4 T (133) where S P = [I 3×3 , 0 3×1 ], E kl = e k e T l + e l e T k ∈ S 4×4 , k, l ∈ {1, 2, 3, 4}, and With the column representation of each C f ijk contained in δC f can now be computed as To derive the quantity ∂J T (T ⋆ ,C f ij ) ∂c f ij , we give two lemmas which are useful for subsequent deduction.
E kl = e k e T l + e l e T k (k ̸ = l) e l e T l (k = l) Proof. For the k-th row, l-th column element of C, denoted by C k,l , we have ∂Cw ∂C k,l = E kl w ∈ R 4 .
U l w = −⌊u l ⌋ 0 0 u l w 1:3 w 4 = ⌊w 1:3 ⌋u l w 4 u l = ⌊w 1:3 ⌋ w 4 I 3×3 u l Thus, ∂U l w ∂u l = ⌊w 1:3 ⌋ w 4 I 3×3 = g 2 (w) With these two lemmas, next we will continue to derive ∂J T (T ⋆ ,C f ij ) ∂c f ij . Theorem 4 gives the Jacobian for one cost item, to distinguish the different cost item, we add a subscript ν to (30) to denotes the ν-th item and replace C j with the actual point cluster notation C fνj corresponding to the ν-th item, leading to: J ν = · · · J p ν · · · ∈ R 1×6Mp (138) The total Jacobian is hence Next, we calculate the partial derivative ∂(J p ) T ∂c f ij . Note that in the summation of (140), only the i-th summation term (i.e., ν = i) is related to c fij , hence, Since C i = Mp ν=1 T ν C fiν T T ν , u il is the eigenvector associated to the l-th largest eigenvalue of matrix A(C i ), and , the derivative with respect to c fij on the right hand side of (141) consists of four terms from U il , C fip (only when p = j), C i , and u il , respectively. Combining with Lemma 3 and Lemma 4, we have Only the term ∂u il ∂c f ij is unknown in this formula. To compute it, we apply the second conclusion in Lemma 1 to all components of c fij and use the notation trick similar to (81): From Lemma 2, Substitute it into L p ij : and where which can be computed beforehand without enumerating each raw point in the run time.