LU分解法：

最好的功能是在于它在解A*x = b，并且当有许多个A相同，b不同的式子时，效率较高。

算法思想（由VC++那本书上说，这个又叫杜利特尔Doolittle方法）：

我们将A分解成为两个子矩阵L和U，令A = L * U。其中L是下三角矩阵，U是上三角矩阵。

大致可以得到如下（Tex公式)：

[tex]\begin{pmatrix} L_{00} & 0 & 0 & 0 \\ L_{10} & L_{11} & 0 & 0 \\ L_{20} & L_{21} & L_{22} & 0 \\ L_{30} & L_{31} & L_{32} & L_{33} \end{pmatrix} \cdot \begin{pmatrix} U_{00} & U_{01} & U_{02} & U_{03} \\ 0 & U_{11} & U_{12} & U_{13} \\ 0 & 0 & U_{22} & U_{23} \\ 0 & 0 & 0 & U_{33} \end{pmatrix} = \begin{pmatrix} a_{00} & a_{01} & a_{02} & a_{03} \\ a_{10} & a_{11} & a_{12} & a_{13} \\ a_{20} & a_{21} & a_{22} & a_{23} \\ a_{30} & a_{31} & a_{32} & a_{33} \end{pmatrix}[/tex]

非Tex公式：

于是我们可以将整个线性方程化为如下形式：

A * x = (L * U) * x = L * (U * x) = b

先解出：

L * y = b

然后再解出：

U * x = y

即可得到最终解。

注：我们在储存L和U的时候，可以将两个矩阵储存在同一矩阵之中。

算法过程：

（1）A分解为L和U：

我们先将L[i,i]令为1，其中i = 0,.., N-1

对于for(j=0;j <= N-1; ++j)，

for(i=0;i <= j; ++i)

[tex]U_{ij} = a_{ij} - \sum_{k=0}^{i-1}(L_{ik}*U_{kj})[/tex]

非tex公式：U[i,j] = a[i,j] - Sum[ L[i,k]*U[k,j],{k,0,i-1} ]

for(i=j+1;i <= N-1; ++i)

[tex]L_{ij} = \frac {1}{U_{jj}}[L_{ij} - \sum_{k=0}^{j-1}(L_{ik}*U_{kj})][/tex]

非tex公式：L[i,j] = 1/U[j,j]*(L[i,j] - Sum[ L[i,k]*U[k,j], {k,0,j-1} ]

（2）对于L * y = b解y得到：

非tex公式：y[0] = b[0]/L[0,0]; y[i] = 1/L[i,i]*(b[i] - Sum[ L[i,j] * y[j], {j,0,i-1} ]), i = 1, 2, ..., N-1

对于U * x = y解x得到：

非tex公式：x[N-1] = y[N-1]/U[N-1,N-1]; x[i] = 1/U[i,i]*(y[i] - Sum[ U[i,j] * x[j], {j,i+1,N-1} ], i = N-2, N-3, ..., 0

使用概述：

如果需要重复使用A矩阵来解不同b的线性方程时，我们在使用完第一次solve(b,x)之后，重新刷新b的值，再度调用solve(b,x)即可在x中获得新的解。如有不明白之处请查看本文中LUdecomposition_test1.cpp

Version 1.0

/**
*@file LUdecomposition.h
*/
/**
*LU Decomposition Introduction:
*Solves "A * x = b"
*(1) Let "L * U = A", shown as follow:
*		(	L00		0		0		0	)		(	U00	U01	U02	U03	)		(	a00	a01	a02	a03	)
*		(	L10	L11		0		0	)		(		0	U11	U12	U13	)		(	a10	a11	a12	a13	)
*		(	L20	L21	L22		0	)	*	(		0		0	U22	U23	)	=	(	a20	a21	a22	a23	)
*		(	L30	L31	L32	L33	)		(		0		0		0	U33	)		(	a30	a31	a32	a33	)
*(2)	Solve "A * x = b" using a decomposition as:
*			A * x = (L * U) * x = L * (U * x) = b
*		then
*			L * y = b
*		then solve
*			U * x = y
*(3) solve for L and U as follow:
*			L[i,0] * U[0,j] + ... = a[i,j]
*		set "a[i,i] = 1, i = 0, ..., N-1"
*		for i = 0, 1, ..., j, solve for U[i,j] as follow
*			U[i,j] = a[i,j] - Sum[L[i,k] * U[k,j], {k, 0, j-1}]
*		for i = j+1, j+2, ..., N-1 solve for L[i,j] as follow
*			L[i,j] = 1/U[j,j] * (a[i,j] - Sum[L[i,k] * U[k,j], {k, 0, j-1}])
*(4) forward substitution as follow:
*			y[0] = b[0]/L[0,0];
*			y[i] = 1/L[i,i]*(b[i] - Sum[L[i,j] * y[j], {j, 0, i-1}]), i = 1, 2, ..., N-1
*		backsubstitution as follow:
*			x[N-1] = y[N-1]/U[N-1,N-1]
*			x[i] = 1/U[i,i]*(y[i] - Sum[U[i,j] * x[j], {j, i+1, N-1}]), i = N-2, N-3, ..., 0
*(5) save the L and U in a matrix as follow:
*			(	U00	U01	U02	U03	)
*			(	L10	U11	U12	U13	)
*			(	L20	L21	U22	U23	)
*			(	L30	L31	L32	U33	)
*		 Matrix--2.3.14
*/

#ifndef LUDECOMPOSITION_H_INCLUDED
#define LUDECOMPOSITION_H_INCLUDED

/**
*@brief Object for solving linear equations A*x = b using LU decomposition, and related functions.
*/
struct LUdecomposition{
	Int n;
	MatDoub lu; ///Stores the decomposition.
	VecInt indx; ///Stores the permutation.
	Doub det_d; ///Used by det.
	LUdecomposition(MatDoub_I &a); ///Constructor. Argument is the matrix A.
	void solve(VecDoub_I &b, VecDoub_O &x); ///Solve for a single right-hand side.
	void solve(MatDoub_I &b, MatDoub_O &x); ///Solve for multiple right-hand sides.
	void inverse(MatDoub_O &a_inverse); ///Calculate matrix inverse A^(-1).
	Doub det(); ///Return determinant of A.
	void mprove(VecDoub_I &b, VecDoub_IO &x);
	MatDoub_I &aref; ///Used only by mprove.
};

/**
*@brief Constructor
*@param  a[0..n-1][0..n-1]  Input as Matrix A in A * x = b, and Output a as after LU decomposition.
*Input:
*	Given a matrix a[0..n-1][0..n-1], this routine replaces it by the LU decomposition of a
*		rowwise permutation of itself. "a" is input.
*Output:
*	"a" is arranged as in "Matrix--2.3.14" above;
*	indx[0..n-1] is an output vector that records the row permutation effected by
*		the partial pivoting
*	d is output as 1 or -1 depending on whether the number of row interchanges
*		was even or odd, respectively.
*This routine is used in combination with solve to solve linear equations or invert a matrix.
*/
LUdecomposition :: LUdecomposition(MatDoub_I &a) : n(a.nrows()), lu(a), aref(a), indx(n) {
	const Doub TINY = 1.0e-40; //A small number;
	Int i, imax, j, k;
	Doub big, temp;
	VecDoub vv(n); //vv stores the implicit scaling of each row.
	det_d = 1.0; //No row interchages yet.
	for(i = 0; i < n; ++i){
		//Loop over rows to get the implicit scaling information
		big = 0.0;
		for(j = 0; j < n; ++j){
			if((temp = abs(lu[i][j])) > big)
				big = temp;
		}//for
		if(big == 0.0) //No nonzero largest element.
			NRthrow("Singular matrix in LUdecomposition");
		vv[i] = 1.0/big;
	}//for

	for(k = 0; k < n; ++k){
		//This is the outermost k[j,j] loop.
		big = 0.0; //Initialize for the search for largest pivot element.
		for(i = k; i < n; ++i){
			temp = vv[i] * abs(lu[i][k]);
			if(temp > big){
				//Is the figure of merit for the pivot better than the best so far?
				big = temp;
				imax = i;
			}//if
		}//for
		if(k != imax){
			//Interchange rows
			for(j = 0; j < n; ++j){
				temp = lu[imax][j]; lu[imax][j] = lu[k][j]; lu[k][j] = temp;
			}//for
			det_d = -det_d; //change the parity of d
			vv[imax] = vv[k]; //interchange the scale factor.
		}//if
		indx[k] = imax;
		if(lu[k][k] == 0.0)
			lu[k][k] = TINY;
			/*
			*If the pivot element is zero, the matrix is singular (at least to the precision of
			*the algorithm). For some applications on singular matrices, it is desirable to
			*substitute TINY for zero.
			*/
			for(i = k+1; i < n; ++i){
				temp = lu[i][k] /= lu[k][k]; //Divide by the pivot element.
				for(j = k+1; j < n; ++j){
					//Innermost loop: reduce remaining submatrix.
					lu[i][j] -= temp * lu[k][j];
				}//for
			}//for
	}//for
}//LUdecomposition(MatDoub_I &a)

/**
*@brief LU Decomposition solve A * x = b for a single right-hand side.
*Solves the set of n linear equations A * x = b using the stored LU decomposition of A.
*b[0..n-1] is input as the right-hand side vector b, while x returns the solution vector x;
*b and x may reference the same vector, in which case the solution overwrites the input.
*This routine takes into account the possibility that b will begin with many zero elements,
*so it is efficient for use in matrix inversion.
*
*Input:
*@param  b[0..n-1]  Input as right-hand side vector
*Output:
*@param  x[0..n-1]  Output as solution vector x of A^(-1) * b
*/
void LUdecomposition :: solve(VecDoub_I &b, VecDoub_O &x){
	Int i, ii = 0, ip, j;
	Doub sum;
	if(b.size() != n || x.size() != n)
		NRthrow("LUdecomposition : solve bad sizes. The right-hand vector should be the same as"
				"the solution vector.");
	for(i = 0; i < n; ++i)
		x[i] = b[i];
	for(i = 0; i < n; ++i){
		/*
		*When ii is set to a positive value, it will become the index of the first
		*nonvanishing element of b. We now do the forward substitution in Information Step(4).
		*The only new small problem is to get the permutation as we go.
		*(
		*	If you read the "next_permutation" function of STL in "stl_algo.h", you will find a similar
		*	algorithm.
		*)
		*/
		ip = indx[i];
		sum = x[ip];
		x[ip] = x[i];
		if(ii != 0){
			for(j = ii-1; j < i; ++j)
				sum -= lu[i][j] * x[j];
		}//if
		else if(sum != 0.0){
			//A nonzero element was encountered, so from now we will have to do the sums in the loop above.
			ii = i+1;
		}//else if
		x[i] = sum;
	}//for
	for(i = n-1; i >= 0; --i){
		//Now we do the backsubstitution in Information Step(4).
		sum = x[i];
		for(j = i+1; j < n; ++j)
			sum -= lu[i][j] * x[j];
		x[i] = sum/lu[i][i]; //Store a component of the solution vector X.
	}//for
}//LUdecomposition :: solve(VecDoub_I &b, VecDoub_O &x)

/**
*@brief LU Decomposition solve for multiple right-hand sides.
*Solves m sets of n linear equations A * X = B using the stored LU decomposition of A.
*The matrix b[0..n-1][0..m-1] inputs the right-hand sides, while x[0..n-1][0..m-1] returns
*the solution A^(-1) * B.
*@attention b and x may reference the same matrix, in which case the solution overwrites the input.
*
*Input:
*@param  b[0..n-1][0..m-1]  Input as right-hand side matrix
*Output:
*@param  x[0..n-1][0..m-1]  Output as solution matrix x of A^(-1) * B
*/
void LUdecomposition :: solve(MatDoub_I &b, MatDoub_O &x){
	Int i, j, m = b.ncols();
	if(b.nrows() != n || x.nrows() != n || b.ncols() != x.ncols())
		NRthrow("LUdecomposition :: solve bad sizes. The right-hand matrix should be the same as"
				"the solution matrix.");
	VecDoub xx(n);
	for(j = 0; j < m; ++j){
		//Copy and solve each column in turn.
		for(i = 0; i < n; ++i)
			xx[i] = b[i][j];
		solve(xx,xx);
		for(i = 0; i < n; ++i)
			x[i][j] = xx[i];
	}//for
}//LUdecomposition :: solve(MatDoub_I &b, MatDoub_O &x)

/**
*@brief Using the stored LU decomposition, return in a_inverse the matrix inverse, A^(-1).
*@attention the matrix inverse, A^(-1) is in a_inverse.
*/
void LUdecomposition :: inverse(MatDoub_O &a_inverse){
	Int i,j;
	a_inverse.resize(n,n);
	for(i = 0;i < n; ++i){
		for(j = 0; j < n; ++j)
			a_inverse[i][j] = 0.;
		a_inverse[i][i] = 1.;
	}
	solve(a_inverse, a_inverse);
}

/**
*@brief Using the stored LU decomposition, return the determinant of the matrix A.
*@return the determinant of the matrix A
*/
Doub LUdecomposition :: det(){
	Doub dd = det_d;
	for(Int i = 0; i < n; ++i)
		dd *= lu[i][i];
	return dd;
}


#endif // LUDECOMPOSITION_H_INCLUDED

以下为对LUdecomposition.h头文件的测试文件，仅供参照：

/**
*@file LUdecomposition_test.cpp
*/
#include <iostream>
#include "nr3.h"
#include "LUdecomposition.h"
using namespace std;

int main(){
	MatDoub_IO a,b,x;
	int rows, cols;
	//INPUT DATA:
	cout << "Please input rows, cols and the matrix data for A: " << endl;
	cin >> rows >> cols;
	a.resize(rows, cols);
	readMatrix_console(a);

	//Do the Linear Solve; a and b will not be altered.
	{//Start a temporary scope.
		LUdecomposition lu(a);
		//Read matrix B1:
		cout << "Please input rows, cols and the matrix data for B: " << endl;
		cin >> rows >> cols;
		b.resize(rows, cols);
		x.resize(rows, cols);
		readMatrix_console(b);
		lu.solve(b,x);
		//PRINT:
		cout << "Gives the solution matrix: " << endl;
		printMatrix_console(x);
		//Read matrix B2:
		cout << "Please input rows, cols and the matrix data for B: " << endl;
		cin >> rows >> cols;
		b.resize(rows, cols);
		x.resize(rows, cols);
		readMatrix_console(b);
		//PRINT:
		cout << "Gives the solution matrix: " << endl;
		printMatrix_console(x);
	}//destroy alu.

	return 0;
}

Version 2.0

mprove方法的实现：

功能：

改进已知近似解的精度，既可用于求一般方程组的高精度解，也可以用于病态方程组的求解，复杂度O(n^2)。

算法原理：

对于A * x = b（记为①式），其中x表示精确解。你使用上述LU-decomposition的solve方法得到一个近似解x+△x，于是如果你将x+△x与矩阵A相乘，得到的将是一个近似的右端项b+△b，如：“A * (x+△x) = (b+△b)”，记为②式。

我们将②式减去①式，得到“A * △x = △b”记为③式。

我们将②式带入③式，替换△b，则得到：A * △x = A * (x+△x) - b（记为④式），代码中我们将A * (x+△x) - b记为sdp。

我们对④式再进行LUdecomposition分解，即代码中的solve(r, r)，得到△x。

最后将△x从x中减去即可。

测试方法：

对于由上lu.solve(b,x)解出的值x解向量，加上一个扰动值（可以由随机数生成）。

然后执行lu.mprove(b,x);进行解的优化，最后再输出，看是否得到正确的解。

提醒：

由于mprove的复杂度只有O(n^2)，相比于solve的复杂度O(n^3)要小很多，建议在每次执行solve之后都执行1~2次的mprove，可根据精度要求增加次数。

/**
*@brief	Improve a solution vector x[0..n-1] of the linear set of equations A * x = b.
*@param	b	the right-hand side vector as input
*@param	x Input as the original solution, Output as the improved solution.
*/
void LUdecomposition::mprove(VecDoub_I &b, VecDoub_IO &x){
	Int i, j;
	VecDoub r(n);
	for(i=0; i < n; ++i){
		//Calculate the right-hand side, accumulating the residual in higher precision.
		Ldoub sdp = -b[i];
		for(j=0; j < n; ++j)
			sdp += (Ldoub)aref[i][j] * (Ldoub)x[j];
		r[i] = sdp;
	}//for
	solve(r,r); //Solve for the error term.
	for(i = 0; i < n; ++i) //Subtract it from the old solution.
		x[i] -= r[i];
} //mprove(VecDoub_I &b, VecDoub_IO &x)

• Relative Posts
• Global Relative Posts