矩陣LU分解的MATLAB與C++實現

一：矩陣LU分解

矩陣的LU分解目的是將一個非奇異矩陣\(A\)分解成\(A=LU\)的形式，其中\(L\)是一個主對角線為\(1\)的下三角矩陣；\(U\)是一個上三角矩陣。

比如\(A= \begin{bmatrix} 1 & 2 & 4 \\ 3 & 7 & 2 \\ 2 & 3 & 3 \\ \end{bmatrix}\)，我們最終要分解成如下形式：

\[A=L\cdot U = \begin{bmatrix} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 2 & -1 & 1 \\ \end{bmatrix} \cdot \begin{bmatrix} 1 & 2 & 4 \\ 0 & 1 & -10 \\ 0 & 0 & -15 \\ \end{bmatrix} \]

現在主要的問題是如何由矩陣\(A\)計算得到矩陣\(L\)和\(U\)呢？我們將在下面詳細討論。

1.1 LU分解原理

首先從矩陣\(U\)入手，因為它是一個上三角矩陣，所以很容易想到高斯消元法，依次把矩陣\(A\)主對角線左下角的元素消為\(0\)就得到\(U\)了。

然后計算矩陣\(L\)，這里有個技巧，可以這樣想，正是因為有了\(L\)，所以\(U\)的左下部分才能被消為\(0\)，所以我們記錄一下把\(U\)的左下部分消為\(0\)時矩陣\(A\)每行所乘的倍數，這個減去的倍數便是\(L\)左下元素的值！

1.2 LU分解計算舉例

\[A=\begin{bmatrix} 1 & 2 & 4 \\ 3 & 7 & 2 \\ 2 & 3 & 3 \\ \end{bmatrix} \overset{(2)- \color{red}{3} \times (1)}{\underset{}{\to}} \begin{bmatrix} 1 & 2 & 4 \\ 0 & 1 & -10 \\ 2 & 3 & 3 \\ \end{bmatrix} \overset{(3)- \color{red}{2} \times (1)}{\underset{}{\to}} \begin{bmatrix} 1 & 2 & 4 \\ 0 & 1 & -10 \\ 0 & -1 & -5 \\ \end{bmatrix} \overset{(3)+ \color{red}{1} \times (2)}{\underset{}{\to}} \begin{bmatrix} 1 & 2 & 4 \\ 0 & 1 & -10 \\ 0 & 0 & -15 \\ \end{bmatrix} =U \]

在運算過程中左下相應元素減去的倍數(上面紅色的數字)便是矩陣\(L\)左下角的元素，可以得到：

\[L= \begin{bmatrix} 1 & 0 & 0 \\ \color{red}{3} & 1 & 0 \\ \color{red}{2} & \color{red}{-1} & 1 \\ \end{bmatrix}\]

1.3 計算公式總結

通用計算公式是很重要的，因為有了公式之后，編程起來就方便很多了。我們可以根據上面的推導過程整理出如下偽代碼：

\[for \text{ } i = 1 : n \hspace{6cm} \\ for \text{ } j = i : n \quad此時i為行下標，j為列下標\\ \qquad U_{ij}=A_{ij}-\sum_{k=1}^{i-1} L_{ik}U_{kj} \hspace{1cm}\\ \qquad for \text{ } x = i+1 : n \quad 此時x為行下標，i為列下標\\ \qquad L_{xi}=(A_{xi}-\sum_{k=1}^{i-1} L_{xk}U_{ki}) /U_{ii} \hspace{0cm}\\ \]

其中\(n\)為方陣的行或列長度，可以看出先計算矩陣\(U\)的第一行，再計算矩陣\(L\)的第一列，再計算矩陣\(U\)的第二行，再計算矩陣\(L\)的第二列，依此類推。

---

二：矩陣LU分解MATLAB實現

clc,clear all,close all
% 矩陣的LU分解 

%% 自己實現
A = [1 2 4;3 7 2;2 3 3] 
[n,n] = size(A);
L = eye(n,n); % L初始化為單位矩陣
U = zeros(n,n); % U初始化為零矩陣

for i = 1 : n % 根據計算公式實現
    for j = i : n
        U(i,j) = A(i,j) - sum(L(i,1 : i - 1) .* U(1 : i - 1,j)');       
    end
    for x = i + 1 : n
        L(x,i) = (A(x,i) - sum(L(x,1 : i - 1) .* U(1 : i - 1,i)')) ./ U(i,i);        
    end
end
L 
U
%% 內置函數實現

[L1,U1] = lu(A)

---

三：矩陣LU分解C++實現

#include <iostream>
#include <vector>
using namespace std;

int main()
{
	vector<vector<double>> a = { {1,2,4},{3,7,2},{2,3,3} };
	int n = a.size();
	vector<vector<double>> u(n, vector<double>(n));
	vector<vector<double>> l(n, vector<double>(n));

	for (int i = 0; i < n; i++) //初始化矩陣L和矩陣U
		for (int j = 0; j < n; j++)
		{
			u[i][j] = 0;
			if (i == j) l[i][j] = 1;
		}

	for (int i = 0; i < n; i++)
	{
		double sum = 0;
		for (int j = i; j < n; j++)
		{
			for (int k = 0; k <= i - 1; k++)
				sum += l[i][k] * u[k][j];
			u[i][j] = a[i][j] - sum; //計算矩陣U
			sum = 0;
		}

		for (int x = i + 1; x < n; x++)
		{
			for (int k = 0; k <= i - 1; k++)
				sum += l[x][k] * u[k][i];
			l[x][i] = (a[x][i] - sum) / u[i][i]; //計算矩陣L
			sum = 0;
		}
	}

	cout << "A:" << endl; //輸出矩陣A
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			printf("%.3f ", a[i][j]);
		}
		cout << endl;
	}

	cout << "L:" << endl; //輸出矩陣L
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			printf("%.3f ", l[i][j]);
		}
		cout << endl;
	}

	cout << "U:" << endl; //輸出矩陣U
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			printf("%.3f ", u[i][j]);
		}
		cout << endl;
	}


	return 0;
}