Use umfpack to solve large-scale sparse matrix equations

Recently, in the model PRT processing, we need to use the least square method to optimize the fitting of sampling points. We may know the least square method (http://en.wikipedia.org/wiki/Least_squares ), in general, in order to facilitate the calculation of the target point, a matrix is usually introduced and the following matrix equation is obtained:
Ax = B
TheXIs a parameter column vector about the target point,BIs the corresponding observed vector, andAIs the data matrix obtained based on the expected background. The linear equations in the matrix form can be obtained by the least square method.XOptimize valuation.

In most casesAIt is a large and sparse matrix. Of course, you can write a program here to solve the problem, but isn't it better to use a verified library? Umpack is a solution similarAx = BA library of such problems comes from Florida State University. You can directly. The original file needs to compile and generate Lib for makefile. If it is inconvenient, you can download the generated lib here.

The usage of the package is very simple. You can refer to the corresponding doc in the package. Here is a brief description. The simplest application is as follows:

# Include "umfpack. H "int n = 5; int AP [] = {0, 2, 5, 9, 10, 12}; int Ai [] = {0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4}; double ax [] = {2.0, 3.0, 3.0,-1.0, 4.0, 4.0,-3.0, 1.0, 2.0, 2.0, 6.0, 1.0}; Double B [] = {8.0, 45.0,-3.0, 3.0, 19.0}; Double X [5]; int main () {double * null = (double *) NULL; int I; void * symbolic, * numeric; (void) umfpack_di_symbolic (N, N, AP, AI, ax, & symbolic, null, null); (void) umfpack_di_numeric (AP, AI, ax, symbolic, & numeric, null, null); encrypt (& symbolic); (void) umfpack_di_solve (umfpack_a, AP, AI, ax, X, B, numeric, null, null); umfpack_di_free_numeric (& numeric );//... obtain the Final Solution from X and use return (0 );}

The correspondingAPrototype:

We don't need to worry about solving symboli and numeric in umfpack. Let's take a look at the Matrix Representation Methods in umfpack. Generally, for a large sparse matrix, the compression method is used to indicate better space saving. The compression method used in umfpack is also simple and straightforward (umfpack uses the matrix index of the column master ):

• APIs an integer array with the following size:AThe number of columns + 1, the first element is 0, and then each element value is the sum of the number of non-0 elements in the current column and all previous columns. For exampleAP [1] = 2The number of non-zero elements in the first column is 2;AP [2] = 5The number of non-zero elements in the second column is 5-2 = 3, and so on.
• AIInAPCorresponding to the position of non-zero elements in each column, suchAi [0], Ai [1]It indicates that the positions of two non-0 elements in the first column are 0, 1.
• AxJust likeAICorresponds to one by one, indicating the specific value of each non-zero element.
• BAndXIs the corresponding column vector, the conventional storage method.

In addition, pay attention to the FunctionUmfpack_di_solveMacros used inUmfpack_aIt is used to describe the type of matrix equation to be solved. The corresponding relationship is:
Umfpack_a: Ax = B
Umfpack_at: a' x = B
Umfpac_aat: A. 'X = B

...
For more information, see umfpack. h.

Finally, let's take a look at the efficiency: in practical applications, we can solve an equation with a matrix A scale of 4611x4611. On average, each line has about seven non-zero elements (up to 9; at least 4 ), umfpack calculates the equation for 0.025 ms, and the performance is quite good.