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

2023-10-28

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

背景

LU分解（LU-Factorization）

（1）Theorem on LU Factorization

（2）Cholesky Theorem on LL^t Factorization

背景

考虑线性方程组Ax=b：

$\begin{bmatrix} a_{11} &a_{12} &\dots & a_{1n}\\ a_{21} &a_{22} &\dots & a_{2n} \\ \vdots & \vdots & \ddots &\vdots \\ a_{n1} & a_{n2} &\dots & a_{nn} \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n\\ \end{bmatrix} = \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n\\ \end{bmatrix}$

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

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

$\bold{input}\ n,(a_{ij}),(b_i) \\ \bold{for}\ i=1\ \bold{to}\ n\ \bold{do} \\ \quad x_i\ \leftarrow \ (b_i-\sum_{j=1}^{i-1}a_{ij}x{j})/a_{ii}\\ \bold{end\ do}\\ \bold{output}\ (x_i)$

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

$Ax=LUx=Lz=b$

$L=\begin{bmatrix} l_{11} & 0 & \dots & 0\\ l_{21} & l_{22}& \dots &0\\ \vdots& \vdots &\ddots & \vdots\\ l_{n1}& l_{n2} & \dots &l_{nn} \end{bmatrix},\ U=\begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n}\\ 0 & u_{22}& \dots &u_{2n}\\ \vdots& \vdots &\ddots & \vdots\\ 0& 0 & \dots &u_{nn} \end{bmatrix}$

通过求解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的转置。

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

计算数学

c