A method for using Cuda to accelerate large-scale BA parallel optimization

Abstract

This invention discloses a method for accelerating large-scale BA parallel optimization using CUDA, comprising the following steps: S1, constructing the energy residual equation of the system; S2, solving the linear equation system; S2-1, applying J... l Perform QR decomposition; S2-2 Use QR decomposition to marginalize the linear equation system; S3, Solve to obtain Δx p Δx l Substituting the original system equations, BA optimization is achieved. The above technical solution improves the numerical stability of the computation and allows for solving large-scale bundle adjustment problems using single-precision floating-point numbers, resulting in good parallelism in solving the linear equation system Ax = b. This significantly improves the parallelism of solving the linear equations, effectively increasing the efficiency of solving the linear equation system. This effectively overcomes the limitations of existing hardware resources, enabling optimization in the large-scale 3D scene reconstruction process, and greatly improving optimization efficiency, particularly for global pose and observation point optimization. The optimization results can serve as the basis for reconstructing sparse 3D point clouds of even larger-scale 3D scenes.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of inspection robot algorithm, and particularly relates to a large-scale BA parallel optimization method accelerated by Cuda. BACKGROUND

[0002] BA (Bundle Adjustment) is a core component of many visual SLAM and laser SLAM. It is one of the classic computer vision problems and has been studied by researchers for more than 60 years. Although there are different formulas, the underlying problem is always the same: given a set of approximate point (landmark) locations, observed from many different poses of cameras, what are the actual landmark locations and camera poses. One can already compute accurate 3D positions with a small number of images. However, with more images available, we will get a more complete reconstruction. With the emergence of large-scale network picture libraries, it is necessary to solve the large-scale bundle adjustment problem, that is, thousands of pictures and thousands of landmarks. The solution of large-scale linear equations has always been an important problem in scientific computing. With the continuous development of the hardware architecture of the graphics processing unit (GPU), the functions of the GPU have been derived to the field of GPU general-purpose computing. The GPU, as a co-processor of the CPU, completes large-scale intensive computing tasks, and the computing capacity is equivalent to a small-scale cluster. Compared with the cluster, the GPU has an advantage in power consumption and cost. In 2007, NVIDIA launched the CUDA parallel computing platform, which reduces the difficulty of using GPU accelerated computing, and researchers can more conveniently use GPU accelerated computing to solve problems in scientific research, which makes more and more fields begin to introduce GPU computing.

[0003] In the BA optimization problem, the core algorithm is the solution of the large-scale linear equation Ax=b. At the same time, through the analysis of the actual data, it can be known that the coefficient matrix A in the linear equation is usually a sparse matrix, and the density of the large-scale sparse matrix of the SLAM system is usually less than 1%. If the sparsity of the matrix cannot be effectively utilized, a large amount of storage and computing resources will be wasted on invalid zero elements when using a computer to process large-scale sparse matrix operations, which leads to low processing efficiency and slow speed, and part of the super large-scale sparse matrix cannot be solved by using the traditional dense matrix algorithm. SUMMARY

[0004] According to the deficiencies of the prior art, the present application provides a large-scale BA parallel optimization method accelerated by Cuda, which greatly improves the parallel degree of linear equation solving, effectively improves the solving efficiency of the linear equation system, and further improves the optimization efficiency of the global pose and observation point.

[0005] In order to solve the above technical problems, the technical scheme of the present application is as follows:

[0006] A method for using Cuda to accelerate large-scale BA parallel optimization method, comprising the following steps:

[0007] S1, constructing the energy residual equation of the system

[0008] S1-1 defining the point x in the world coordinate system p , by rotating angle R and translation transformation t, the point x p Project into pixel coordinates x p ';

[0009] S1-2 construct the pixel coordinates x p x l And the pixel coordinates x p ' of the projection constitute the system error equation r=(x l -Rx p ′-t)=(x l -x p ';

[0010] S1-3 minimize the residual equation r by Levenberg-Marquardt algorithm, have

[0011]

[0012] Wherein, r=r(x 0 ), Δx=x-x 0 ;

[0013] S2, solve the linear equation

[0014] S2-1 QR decomposition of J l , have

[0015]

[0016] S2-2 using QR decomposition edge linear equation group;

[0017] S3, the solution of Δx p , Δx l Substitute the original system equation to realize BA optimization.

[0018] As preferred, the step 2-2

[0019] Q left formula (1) in the above formula (2), have

[0020]

[0021] Q=(Q1,Q2) into the above formula (3), realize the null space edge of the system, complete the unknown Δx p , Δx lDecoupling of the solution.

[0022] As preferred, the step 2-2 further comprises

[0023] The above formula (3) is decomposed into equation groups Emin1(Δx p ), Emin2(Δx p , Δx l ) by QR decomposition, the equation group Emin1(Δx p ) is solved first to obtain Δx p , and then Δx p is brought into the equation group Emin2(Δx p , Δx l ) to obtain Δx l .

[0024] As preferred, the step S2 further comprises a step S2-3

[0025] The equation groups Emin1(Δx p ), Emin2(Δx p , Δx l ) are accelerated to be solved by using a Cuda model.

[0026] As preferred, in the step S2-1,

[0027] Let J l be a full rank matrix of m×n, m≥n, the rank of the matrix J l is n, and the matrix J l is decomposed into an m×m orthogonal matrix Q and an m×n upper triangular matrix R, since the lower (m-n) rows of R are zero, R and Q are divided into:

[0028]

[0029] where R1 is an n×n upper triangular matrix, Q1 is m×n, and Q2 is m×(m-n).

[0030] As preferred, the columns of the Q2 constitute the left null space of A, i.e.

[0031] As preferred, in the step 2-2, the above formula (3) is modified by QR decomposition of the above formula (4) and the invariance of the L2 norm under orthogonal transformation:

[0032]

[0033] As preferred, QR decomposition is performed on the above (5) to obtain

[0034] The equation group Emin1(Δx p ):

[0035]

[0036] Equation group Emin2 (Delta x p , Delta x l ) :

[0037]

[0038] The present application has the following characteristics and benefits:

[0039] By using the above technical solution, the method is algebraically equivalent to the commonly used Schur complement technique, improves the numerical stability of the calculation, and allows the use of single-precision floating-point numbers to solve large-scale bundle adjustment problems. Using zero space edge as an alternative to the traditional Schur complement algorithm makes the solution of the linear equation system Ax=b have good parallelism. At the same time, for the linear equation system after zero space edge, the linear equation solving based on GPU sparse symmetric matrix is used, which greatly improves the parallel degree of linear equation solving, effectively improves the efficiency of linear equation system solving, and does not require additional equipment to accelerate the solution, to a certain extent, reduces the cost of hardware equipment, thereby effectively getting rid of the limitation of existing hardware resources, realizing the optimization in the process of large-scale three-dimensional scene three-dimensional reconstruction, and greatly improving the optimization efficiency. In the slam industrial scene involving the present scheme, the speed of three-dimensional scene reconstruction is greatly improved, the optimization efficiency of global pose and observation points is improved, and the optimization result can be used as the basis for reconstructing a larger scale three-dimensional scene sparse three-dimensional point cloud. BRIEF DESCRIPTION OF DRAWINGS

[0040] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor.

[0041] Figure 1 The method flowchart of the present application;

[0042] Figure 2 The element distribution form of the coefficient matrix in the system residual equation;

[0043] Figure 3 The element distribution form of the coefficient matrix after edge processing in the system residual equation. DETAILED DESCRIPTION

[0044] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0045] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, the terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined with "first," "second," etc., may explicitly or implicitly include one or more of that feature. In the description of this invention, unless otherwise stated, "a plurality of" means two or more.

[0046] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art will understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0047] This invention provides a method for accelerating large-scale BA parallel optimization using CUDA, such as... Figure 1 As shown, it includes the following steps:

[0048] S1. Construct the energy residual equation of the system.

[0049] S1-1 Defines a point x in the world coordinate system. p By rotating by an angle R and translating by a translation t, point x is transformed. p Projected as pixel coordinates x p ';

[0050] S1-2 constructs x p pixel coordinates x l With the pixel coordinates of the projection x p The system error equation formed is r=(x l -Rx p (x'-t)=( ... l -x p ′);

[0051] S1-3 minimizes the residual equation r using the Levenberg-Marquardt algorithm, and has

[0052]

[0053] where r = r(x 0 ), Δx = x - x 0 ;

[0054] S2, solve the linear system

[0055] S2-1, QR decomposition of J l , and

[0056]

[0057] S2-2, use the QR decomposition to edge the linear system;

[0058] S3, substitute Δx p into the original system equation to achieve BA optimization. l

[0059] In the above technical solution, in order to realize the parallel operation of solving the linear equation system Ax = b, a new formula is proposed, which depends on the zero space edge of the sign variable based on QR decomposition. This method is algebraically equivalent to the commonly used Schur complement technique, which improves the numerical stability of the calculation and allows the use of single-precision floating-point numbers to solve large-scale bundle adjustment problems. Using zero space edge as an alternative to the traditional Schur complement algorithm makes the linear equation system Ax = b solving have good parallelism. At the same time, for the linear equation system after zero space edge, the linear equation solving of sparse symmetric matrix based on GPU is used, which greatly improves the parallel degree of linear equation solving and effectively improves the efficiency of linear equation system solving, thereby effectively getting rid of the limitation of existing hardware resources, realizing the optimization in the process of large-scale three-dimensional scene three-dimensional reconstruction, and greatly improving the optimization efficiency, improving the optimization efficiency of global pose and observation points, and the optimization result can be used as the basis for reconstructing larger scale three-dimensional scene sparse three-dimensional point cloud.

[0060] Specifically, in step S1, x p is defined as follows: Assuming a very general AB optimization scenario, let x be a state vector containing all optimization variables. We can subdivide x into a pose part x p , which contains the extrinsic and possibly intrinsic camera parameters of all n images (indexed by i), and a landmark part x l , which contains the 3D coordinates of all n l landmarks (indexed by j). The sum of squared energy residuals of this system is:

[0061] ​

[0062] Then the residual equation r is minimized by Levenberg-Marquardt algorithm, that is,

[0063]

[0064] Where r = r(x 0 ), Δx = x-x 0 . At this time, the coefficient matrix of the equation group Emin(Δx p , Δx l ) is in the form as shown in the following formula (2). Figure 2

[0065] Figure 2 The light gray block in the matrix represents a non-zero element, and the dark gray block and the white block represent zero elements. As can be seen, the number of non-zero elements in the matrix J p and J l is much smaller than the total number of matrix elements, and the matrix has strong sparsity.

[0066] Further arrangement of the application, the step 2-2

[0067] Multiply Q in the above formula (2) on the left of formula (1), and have

[0068]

[0069] Again, Q=(Q1,Q2) is brought into the above formula (3), the zero space of the system is edge, and the decoupling of the unknown Δx p , Δx l is completed.

[0070] Specifically, in the step S2-1,

[0071] Let J l be a full rank matrix of m*n, m>=n, the rank of the matrix J l is n, and the matrix J l is decomposed into an m*m orthogonal matrix Q and an m*n upper triangular matrix R, since the lower (m-n) row of R is zero, R and Q are divided into:

[0072]

[0073] Where R1 is an n*n upper triangular matrix, Q1 is m*n, and Q2 is m*(m-n). It can be understood that this division of Q directly means that the columns of Q2 constitute the left null space of A, that is,

[0074] Further, in the step 2-2, the QR decomposition of the above formula (4) is used, and the invariance of L2 norm under orthogonal transformation, and the above formula (3) is modified as:​

[0075]

[0076] Further, QR decomposition is performed in (5) above to obtain

[0077] Equation group Emin1 (Δx p ):

[0078]

[0079] Equation group Emin2 (Δx p , Δx l ):

[0080]

[0081] It can be understood that, by Figure 3 It can be seen that, after QR decomposition and edge linearization, the sparsity of the matrix is greatly reduced, and the number of zero elements in the matrix block R1, is much smaller than the total number of non-zero elements, and the matrix has very weak sparsity. Usually, an N*N array is used to store an N*N matrix, but for sparse matrices, this method is not efficient. Since zero elements do not participate in the operation, they occupy a large part of the space for storage. For example: for a 30000*30000 single-precision sparse matrix with a sparsity of 5%, if it is stored in the general method, the space required is 30000*30000*sizeof(float)*8 / (1024*1024*1024) = 26.8GB. But if only the values of non-zero elements are stored, the space required is 30000*30000*5%*sizeof(float)*8 / (1024*1024*1024) = 1.34GB. From the above comparison, it can be seen that for sparse matrices, the strategy of storing non-zero elements is obviously more efficient than the general method, saving storage space and also saving calculation time.

[0082] It can be understood that, in the above technical solution, by decomposing the algebraic problem of simultaneously solving unknown Δx p and Δx l into the solution of two linear equation groups, using zero space edge linearization as an alternative to the traditional Schur complement algorithm, the linear equation group Ax = b is solved with good parallelism. At the same time, for the linear equation group after zero space edge linearization, the linear equation solving based on GPU sparse symmetric matrix is used, which greatly improves the parallel degree of linear equation solving and effectively improves the efficiency of linear equation group solving.

[0083] Further, the step 2-2 further includes

[0084] Equations Emin1 (Delta x p ), Emin2 (Delta x p , Delta x l ) are obtained by QR decomposition of the above formula (3) into the equation group p , Delta x p is obtained by solving the equation group Emin1 (Delta x p ), and then Delta x p is brought into the equation group Emin2 (Delta x l ) to obtain Delta x l .

[0085] Further, the step S2 further includes a step S2-3

[0086] The equation group Emin1 (Delta x p ), Emin2 (Delta x p , Delta x l ) is accelerated to be solved by using a Cuda model.

[0087] Specifically, the functions in the CUBLAS library are fully utilized, the SpMV method is combined, and the precondition conjugate gradient method is realized by using a GPU to solve the equation group Emin1 (Delta x p ), Emin2 (Delta x p , Delta x l ).

[0088] In the Cuda programming model, the CPU is responsible for the execution of logical transactions and serial codes, the GPU is responsible for the execution of parallel codes, and the CPU and the GPU work together to calculate. In the equation solving process, each iteration process is serial, but the SpMV operation, the vector inner product, the vector addition and subtraction and other operations are data-level parallel operations, which can be realized by using the GPU in parallel, and the remaining operations are executed by the CPU. Therefore, the specific allocation mode is that the CPU is responsible for the iteration loop and the judgment of the convergence condition, the scalar division operation, the opening and release of the GPU memory and the CPU memory, and the start of the kernel function; the GPU is responsible for the parallel calculation of the vector inner product, the vector addition and subtraction and the SpMV operation in each iteration.

[0089] The specific process of solving the equation group based on Cuda is as follows: (wherein b and x are variables in the GPU memory, b is the right end vector of the equation group, x is the initial value after the iteration, M -1 is the inverse matrix of the preconditioning matrix)

[0090] (1) In the global memory of GPU, allocate arrays r(n), p(n), Ap(n), h(n) storage space, define variables alpha, beta, error_norm, r2 in CPU memory;

[0091] (2) Parallel computing r = b - A * x on GPU

[0092] (2.1) Call cudaM emcpy DeviceToDevice, copy b to r: r = b;

[0093] (2.2) Call kernel of SpMV to compute A * x: Ap = A * x;

[0094] (2.3) Call CUBLAS library function cublasSgbmv to compute r = r - Ap;

[0095] (3) Call kernel to compute h = M -1 *r;

[0096] (4) Call cudaM emcpy DeviceToDevice, copy h to p: p = h;

[0097] (5) Call vector inner product kernel to compute the inner product of r and h and assign it to the variable error_norm:

[0098] error_norm

[0099] = (h, r)

[0100] (6) for k = 1, 2, 3, … n

[0101] (6.1) Assign value to CPU-side variable: r2 = error_norm;

[0102] (6.1) Call kernel function of SpMV to compute Ap = A * p;

[0103] (6.2) Call vector inner product kernel to compute the inner product (p, Ap);

[0104] (6.3) Compute alpha = r2 / (p, Ap) on the CPU side;

[0105] (6.4) Call CUBLAS library function to compute and update variables x, r:

[0106] x = x + alpha * p; r = r * alpha * Ap;

[0107] (6.5) Call kernel function to compute h = M-1 r;

[0108] (6.6) Call vector inner product kernel to compute inner product and assign to CPU side variable error_norm

[0109] = (h, r);

[0110] (6.7) CPU side check convergence condition: if error_norm < 1e-13 exit loop;

[0111] (6.8) CPU side compute beta = error_norm / r2;

[0112] (6.9) Call CUBLAS library function to compute and update vector p:

[0113] p = r * beta * p;

[0114] (7) Free GPU side global memory variables r, p, h, Ap

[0115] end

[0116] It can be understood that in the above technical solution, by decomposing the original algebraic problem of simultaneously solving unknown quantities Δx p and Δx l into the solution of two linear equation groups, and solving through the Cuda model, the parallel degree of linear equation solving is further improved, and the efficiency of linear equation group solving is effectively improved.

[0117] The embodiments of the present application are described in detail above with reference to the drawings, but the present application is not limited to the described embodiments. For those skilled in the art, various changes, modifications, replacements and variations of the embodiments including components can be made without departing from the principles and spirits of the present application, and still fall within the protection scope of the present application.

Claims

1. A method for accelerating large-scale business intelligence (BA) parallel optimization using CUDA, characterized in that, Includes the following steps: S1. Construct the energy residual equation of the system. S1-1 Defines a point x in the world coordinate system. p By rotating by an angle R and translating by a translation t, point x is transformed. p Projected as pixel coordinates x p '; S1-2 constructs x p pixel coordinates x l With the pixel coordinates of the projection x p The system error equation formed is r=(x l -Rx p '-t)=(x l -x p '); S1-3 minimizes the residual equation r using the Levenberg-Marquardt algorithm, and has Where r = r(x) 0 ), Δx=xx 0 ; S2. Solving the system of linear equations S2-1 to J l Perform QR decomposition, and have S2-2 uses QR decomposition to marginalize the linear equation system; S3, Solve for Δx p Δx l Substitute the original system equations to achieve BA optimization.

2. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 1, characterized in that, Step 2-2 Multiplying Q in equation (2) by formula (1) on the left, we have: Substituting Q = (Q1, Q2) into equation (3) above, we can achieve null space marginalization of the system and complete the unknown quantity Δx. p Δx l Solving for decoupling.

3. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 2, characterized in that, Step 2-2 also includes Equation (3) above is decomposed into a system of equations Emin1(Δx) using QR decomposition. p ), Emin2(Δx p Δx l First, solve the system of equations Emin1(Δx) p ), thus obtaining Δx p Then Δx p Substitute into the system of equations Emin2(Δx) p Δx l Find Δx l .

4. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 2, characterized in that, Step S2 further includes steps S2-3 Using the CUDA model, the equation system Emin1(Δx) is solved. p ), Emin2(Δx p Δx l () to accelerate the solution.

5. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 4, characterized in that, In step S2-1, Let J l Let J be an m×n full-rank matrix, where m ≥ n. l The rank is n, and the matrix J l Decompose into an m×m orthogonal matrix Q and an m×n upper triangular matrix R. Since the lower (mn) row of R is zero, partition R and Q as follows: R1 is an upper triangular matrix of size n×n, Q1 is of size m×n, and Q2 is of size m×(mn).

6. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 5, characterized in that, The columns of Q2 constitute the left null space of A, that is...

7. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 5, characterized in that, In step 2-2, using the QR decomposition of equation (4) above and the invariance of the L2 norm under orthogonal transformation, equation (3) above is transformed into:

8. The method for accelerating large-scale BA parallel optimization using CUDA according to claim 7, characterized in that, The QR decomposition obtained by equation (5) above yields: The system of equations Emin1(Δx) p ): The system of equations Emin2(Δx) p Δx l ):

Images