数值分析——LU分解（LU Factorization）

本系列整理自博主21年秋季学期本科课程 数值分析I 的编程作业，内容相对基础，参考书: David Kincaid, Ward Cheney - Numerical Analysis Mathematics of Scientific Computing (2002, Americal Mathematical Society)

目录

背景

LU分解（LU-Factorization）

辅助部分

Doolittle分解

Cholesky分解

定理

（1）Theorem on LU Factorization

（2）Cholesky Theorem on LL^t Factorization

---

背景

考虑线性方程组Ax=b：

易知，若矩阵A为下三角或上三角矩阵时，容易解出x

如A为下三角矩阵，则通过向前迭代算法（forward substitution）求解：

因此，若A可以分解为下三角矩阵L和上三角矩阵U的乘积:

 

通过求解Lz=b，再求解Ux=z，得到方程组的解x，从而线性方程组求解问题化为如何对矩阵A做LU分解的问题。

LU分解（LU-Factorization）

对于一个矩阵A可以有无数种LU分解方式，常见的有三种：

（1）Doolittle分解：L为单位下三角阵，即L对角元皆为1

（2）Crout分解：U为单位上三角阵，即U对角元皆为1

（3）Cholesky分解：U=L的转置，有L,U对角元相等

以下给出Doolittle分解和Cholesky分解算法。

辅助部分

#include<iostream>
#include<stdio.h>
#include<math.h>
#include<iomanip>
using namespace std;
#define N 7//matrix dimension
void Print(double M[N][N]) {//print matrix M
	cout << " ";
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			cout << setw(13) << M[i][j];
			if (j == N-1) {
				cout << " " << endl;
				cout << " ";
			}
		}
	}
}
void Multi(double L[N][N], double U[N][N]) {
	double M[N][N] = { 0 };//M=LU
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <N; j++) {
			for (int k = 0; k < N; k++) {
				M[i][j] += L[i][k] * U[k][j];
			}
		}
	}
	cout << " LU" << endl;
	Print(M);
}

Doolittle分解

void DoolittleFactorization(double A[N][N]) {//apply Doolittle factorization on matrix A
	double L[N][N] = { 0 }, U[N][N] = { 0 };
	for (int k = 0; k < N; k++) {
		if (k == 0) {
			for (int j = 0; j < N; j++) {
				U[0][j] = A[0][j];
				L[j][0] = A[j][0] / U[0][0];
			}
		}
		else {
			L[k][k] = 1;
			for (int i = k; i < N; i++) {
				for (int j = 0; j < k; j++) {
					U[k][i] += L[k][j] * U[j][i];
					if (i + 1 < N) {
						L[i + 1][k] += L[i + 1][j] * U[j][k];
					}
				}
				U[k][i] = (A[k][i] - U[k][i]) / L[k][k];
				if (i + 1 < N) {
					L[i + 1][k] = (A[i + 1][k] - L[i + 1][k]) / U[k][k];
				}
				
			}
		}
		
	}
	cout << "  Doolitte Factorization：" << endl;
	cout << " L" << endl;
	Print(L);
	cout << endl;
	cout << " U" << endl;
	Print(U);
	cout << endl;
	Multi(L, U);
	cout << endl;
}

Cholesky分解

void CholeskyFactorization(double A[N][N]) {
	double L[N][N] = { 0 }, U[N][N] = { 0 };
	for (int k = 0; k < N; k++) {
		if (k == 0) {
			L[0][0] = U[0][0] = sqrt(A[0][0]);
			for (int m = 1; m < N; m++) {
				L[m][0] = U[0][m] = A[0][m] / L[0][0];
			}
		}
		else {
			for (int i = k; i < N; i++) {
				for (int j = 0; j < k; j++) {
					U[k][i] += L[k][j] * U[j][i];
				}
				if (i == k) {
					L[i][k] = U[k][i] = sqrt(A[k][i] - U[k][i]);
				}
				else {
					L[i][k] = U[k][i] = (A[k][i] - U[k][i]) / L[k][k];
				}
			}
		}
	}
	cout << "  Cholesky Factorization：" << endl;
	cout << " L" << endl;
	Print(L);
	cout << endl;
	cout << " U" << endl;
	Print(U);
	cout << endl;
	Multi(L, U);
    cout << endl;
}

定理

（1）Theorem on LU Factorization

若n*n阶矩阵A的n个顺序主子式非奇异，则A有LU分解。

（2）Cholesky Theorem on LL^t Factorization

若A 实值，对称，正定，则A有唯一的分解A=LL^t，其中L为对角为正数的下三角矩阵，L^t表示矩阵L的转置。