演算法：矩陣 LUP 分解

本文著筆於矩陣 LUP 分解演算法，以及利用矩陣的 LUP 分解來解線性方程組、求矩陣對應行列式的值、求逆矩陣。

對於矩陣的定義程式碼如下：

struct Matrix
{
	double dat[MAX_N][MAX_N],det,LUdat[MAX_N][MAX_N],rev[MAX_N][MAX_N];
	int dgr,pi[MAX_N];
	bool LUP_Decomposition();
	VectorColumn LUP_Solve(VectorColumn);
	void GetDeterminant();
	void GetReverse();
}
 
;

矩陣的 LUP 分解

矩陣 A 的 LUP 分解即為找到下三角矩陣 L、上三角矩陣 U、行置換矩陣 P 使得 PA = LU ，實現方法為高斯消元。

$\begin{matrix} E \end{matrix} = \begin{matrix} I \end{matrix} = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix}, e_{ij} = \begin{cases} 1 & i = j \\0 & i \neq j \end{cases}$

初始矩陣 P = E。對於矩陣 A，我們可以將它分解成這樣兩個矩陣的乘積（這裡直接借用了 Introduction to Algorithms 裡的描述）
$\begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} & \omega^T \\ \nu & A' \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ \nu / a_{11} & I_{n-1} \end{pmatrix} \begin{pmatrix} a_{11} & \omega^T \\ 0 & A' - \nu \omega^T / a_{11} \end{pmatrix}$

這樣，將這個規模為 n 的問題分解為一個規模為 n - 1 的子問題，迭代之即可求得矩陣 PA 的 LU 分解。最後得到的矩陣 U 相當於高斯消元后的結果，因為求矩陣 U 的實質是用當前規模矩陣的第一行的相應倍數減其它行，將其它行首列元素消為 0。

為了避免除數為零或過小引起精度問題，在迭代到矩陣第 i 行時，我們可以找到 j 使得 $\vert A_{j, i} \vert$ 儘可能大，其中 $i \leq j \leq n$。當 i 不等於 j 時，交換矩陣 PA 的第 i 行與第 j 行，相當於矩陣 A 不變，矩陣 P 的第 i 行與第 j 行交換。

對於矩陣的 LUP 分解程式碼如下：

#define ABS(x) \
({ \
	typeof(x) __tmp_x=(x); \
	__tmp_x<0?-__tmp_x:__tmp_x; \
})

#define EXCHANGE(a,b) \
{ \
	typeof(a) __tmp_a=(a); \
	a=b,b=__tmp_a; \
}

bool Matrix::LUP_Decomposition()
{
	int p;
	for(int i=1;i<=dgr;i++)
	{
		pi[i]=i;
		for(int j=1;j<=dgr;j++)
			LUdat[i][j]=dat[i][j];
	}
	det=1;
	for(int i=1;i<=dgr;i++)
	{
		p=i;
		for(int j=i+1;j<=dgr;j++)
			if(ABS(LUdat[j][i])>ABS(LUdat[p][i]))
				p=j;
		if(LUdat[p][i]==0)
			return false;
		if(p!=i)
		{
			det=-det;
			EXCHANGE(pi[i],pi[p]);
		}
		for(int j=1;j<=dgr;j++)
		EXCHANGE(LUdat[i][j],LUdat[p][j]);
		for(int j=i+1;j<=dgr;j++)
		{
			LUdat[j][i]/=LUdat[i][i];
			for(int k=i+1;k<=dgr;k++)
				LUdat[j][k]-=LUdat[j][i]*LUdat[i][k];
		}
	}
	return true;
}

時間複雜度 $O(n^3)$。

矩陣求解線性方程組

對於線性方程組
$\begin{cases} a_{11} x_1 + & a_{12} x_2 + \cdots + & a_{1n} x_n = & b_1 \\ a_{21} x_1 + & a_{22} x_2 + \cdots + & a_{2n} x_n = & b_2 \\ \vdots & \vdots & \vdots & \vdots \\ a_{n1} x_1 + & a_{n2} x_2 + \cdots + & a_{nn} x_n = & b_n \end{cases}$

構建矩陣 A 和列向量 B，使得
$\begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix}, \begin{matrix} B \end{matrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}$

矩陣求解線性方程組即對於矩陣 A 和列向量 B，找到列向量 X，使得