线性代数(1):LU分解

线性代数(1):LU分解

矩阵分解

矩阵分解(Matrix Factorizations)就是将一个矩阵用两个以上的矩阵相乘的等式来表达。而矩阵乘法涉及到数据的合成(即将两个或多个线性变换的效果组合成一个矩阵),因此可以说,矩阵分解是对数据的一种分析。

矩阵分解 LU分解
分解形式 A=LUA=LU
(LL代表下三角矩阵,UU代表上三角矩阵)
目的 提高计算效率
前提 (1)矩阵是方阵(LULU分解主要是针对方阵);
(2)矩阵是可逆的,也就是该矩阵是满秩矩阵,每一行都是独立向量;
(3)消元过程中没有00主元出现,也就是消元过程中不能出现行交换的初等变换。

LU分解

分解形式

A=LU=[1000l21100l31l3210l41l42l431][u11u12u13u140u22u23u2400u33u34000u44]A = LU = \left[\begin{matrix} 1 & 0 & 0 & 0 \\ l_{21} & 1 & 0 & 0 \\ l_{31} & l_{32} & 1 & 0 \\ l_{41} & l_{42} & l_{43} & 1\end{matrix}\right]\left[\begin{matrix} u_{11} & u_{12} & u_{13} & u_{14} \\ 0 & u_{22} & u_{23} & u_{24} \\ 0 & 0 & u_{33} & u_{34} \\ 0 & 0 & 0 & u_{44}\end{matrix}\right]

简述

LULU分解常用于求解工业和商业问题中的序列方程。它是最常见的求解线性系统 Ax=bAx=b 的方法,主要思路是:把AA 分解成一个下三角矩阵(Lower Triangular Matrix)和一个上三角矩阵(Upper Triangular Matrix),简称LULU。分解后,等式Ax=bAx=b可以写成L(Ux)=bL(Ux)=b,这样我们令y=Uxy=Ux,这样就可以通过分开求解两个等式来得到xx

Ly=bUx=yLy=b \\ Ux=y

即可以先通过Ly=bLy=b求出yy,最后再用Ux=yUx=y,求出xx

  • 本质上,LULU分解是高斯消元法的一种表达方式,为了更好地理解LULU分解,我们需要解释一下高斯消元法。

高斯消元法(Gaussian Elimination)

简述

为了求解线性系统,我们基本的策略是用相等的式子去替代要求解的式子,来让求解变得更容易。对于一个线性系统中有若干个未知数,我们需要将方程组中的一方程的未知数用含有另一未知数的代数式表示,并将其代入到另一方程中,这就消去了一未知数,得到一解;或将方程组中的一方程倍乘某个常数加到另外一方程中去,也可达到消去一未知数的目的,并最终得到线性系统的解。这一求解算法称之为高斯消元法(Gaussian Elimination)。

示例

下文将展示高斯消元法的求解过程,设需求解的线性方程组如下:

{x12x2+x3=02x28x3=85x15x3=10\begin{cases} x_1-2x_2+x_3=0 \\ 2x_2-8x_3=8\\5x_1-5x_3=10 \end{cases}

将其表达为矩阵形式为:

[1210028850510]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]

  1. 为了将多余的未知数x1x_1消除,我们需要将第三行的55消除,矩阵的第一行乘5-5倍再和第三行相加,得到:

[121002880101010]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10 \end{array}\right]

  1. 接着再进行对多余的未知数x2x_2的消除,这里将第二行乘5-5倍再与第三行相加,得到:

[12100288003030]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]

至此,我们得到了一个新的三角形矩阵,转换为线性方程组的形式为:

{x12x2+x3=0x24x3=4x3=1\begin{cases} x_1-2x_2+x_3=0 \\ x_2-4x_3=4\\x_3=-1 \end{cases}

通过简单的回代,我们最终可以得到线性方程组的解:

{x1=1x2=0x3=1[100101000011]\begin{cases} x_1=1 \\ x_2=0\\x_3=-1 \end{cases} \quad\quad\quad \left[\begin{array}{ccc|c} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & -1 \end{array}\right]

高斯消元法与LU分解之间的联系

那么为什么说LULU分解其实就是高斯消元法的一种表达方式呢?这是因为我们可以将上述示例中每一步的操作用矩阵相乘的形式(行变换的本质就是左乘矩阵)表达出来。

例如上文示例中的第一步:“为了将多余的未知数x1x_1消除,我们需要将第三行的55消除,矩阵的第一行乘5-5倍再和第三行相加” ,可以表达为:

[100010501][1210028850510]=[121002880101010]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & 0 & 1 \end{matrix}\right] \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]=\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10 \end{array}\right]

同样地,第二步:“继续进行对多余的未知数x2x_2的消除,这里将第二行乘5-5倍再与第三行相加” 可以表达为:

[100010051][121002880101010]=[12100288003030]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1 \end{matrix}\right] \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10 \end{array}\right] = \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]

因此,我们可以用一个矩阵相乘的等式来记录整个高斯消元的过程:

[100010051][100010501][1210028850510]=[12100288003030]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1 \end{matrix}\right]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & 0 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]=\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]

将三个消元矩阵合并,用EE表示则为:

EA=[100010551][1210028850510]=[12100288003030]=UEA =\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & -5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]=\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right] = U

那么我们可以得到关于AA的表达式:

A=[1210028850510]=E1U=[100010551][12100288003030]=LUA = \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right] = E^{-1}U = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right] = LU

这样一来我们就得到了AALULU分解,即将矩阵A分解为一个下三角矩阵(LL)和一个上三角矩阵(UU)的乘积。

伪代码

U = A, L = I
for j = 1 : n - 1 do
	for i = j + 1 : n do
		l_ij = u_ij / u_jj
		for k = j : n do
			u_ik = u_ik - l_ij * u_jk
		end for
	end for
end for
  • 此处是矩阵index从1开始计数。
  • Credit to http://iacs-courses.seas.harvard.edu/courses/am205/slides/am205_lec07.pdf

算法效率

接下来我们来思考一下算法效率的问题,假如我们有一个100100100*100的矩阵,并且各项均非00,那我们需要计算多少次才能得到矩阵UU呢?

对于这个问题我们先从列的角度来考虑,第一列消元运算结束以后,矩阵变为:

[100]\left[\begin{matrix} 1 & \cdots & \cdots & \cdots \\ 0 & \ddots & \cdots & \cdots \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & \ddots\end{matrix}\right]

这一步中,第一列的元素运算了100100次,而第一行共有100100个元素,于是仅第一行与第一列消元结束后,我们就计算了1002100^2次。之后我们要研究剩下的999999*99的矩阵,以此类推可知,最后的计算量为k=1nk2\sum^{n}_{k=1}{k^2}

接着,我们可以利用微积分计算得到,AA的操作次数为13n3\frac{1}{3}n^3bb的操作次数为n2n^2次。

结合上文伪代码, 因为有三重循环嵌套,因此LULU分解的复杂度为O(n3)O(n^3),其中加法操作为13n3\frac{1}{3}n^3,乘法操作为13n3\frac{1}{3}n^3,总共为23n3\frac{2}{3}n^3次,求解bbyy各需要n2n^2

前提条件

(1)矩阵是方阵(LU分解主要是针对方阵);
(2)矩阵是可逆的,也就是该矩阵是满秩矩阵,每一行都是独立向量;
(3)消元过程中没有0主元出现,也就是消元过程中不能出现行交换的初等变换。

LUD分解(LDU Decomposition

上文中我们已经得到了A=LUA=LU,我们还可以继续分解下去得到A=LDUA=LDUDD对角矩阵(Diagonal Matrix),矩阵的LDULDU分解是在LULU分解之后,把UU再次分解,目的是把UU的对角线元素都化为11,上文示例中我们得到的A=LUA=LU如下

A=[1210028850510]=LU=[100010551][12100288003030]A = \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right] = LU = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]

接下来,让UU矩阵中位于对角线上的元素都为11,就可以提出一个对角矩阵DD

LU=[100010551][12100288003030]=[100010551][1000200030][121001440011]=LDULU = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right] =\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 30 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 1 & -4 & 4 \\ 0 & 0 & 1 & -1 \end{array}\right]=LDU

LUP分解(LU Decomposition with Partial Pivoting)

我们来考虑一个矩阵$ A=\left[\begin{matrix} 0 & 1 \ 1 & 1 \end{matrix}\right],尽管矩阵A使非奇异,并且很简单,但是使用LU分解的算法会失败,因为第一个主元A_{0,0} = 0$。在高斯消元法中,我们常常通过人为的行变换来消除主元为00时的影响。例如我们需要把矩阵A的第一行和第二行进行交换得到$ A=\left[\begin{matrix} 1 & 1 \ 0 & 1 \end{matrix}\right],这样我们就能将消元进行下去。这样的操作也可以用矩阵相乘的方法来表达,那么在该示例中,即为PA=[0110][0111]=[1101]PA=\left[\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}\right]\left[\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right]=\left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right],这里的矩阵PPermutationMatrix称为置换矩阵(Permutation Matrix)。那么就有了LUPLUdecompositionwithpartialpivoting分解(LU decomposition with partial pivoting):PA = LU。这一算法的出现显著改善了LU$分解算法的稳定性。

伪代码

U = A, L = I, P = I
for j = 1 : n - 1 do
	Select i (>= j) that maximizes |u_ij|
	Interchange rows of U: u_(j,j:n) <-> u_(i,j:n)
	Interchange rows of L: l_(j,1:j-1) <-> l_(i,1:j-1)
	Interchange rows of P: p_(j,:) <-> p_(i,:)
	for i = j + 1 : n do
		l_ij = u_ij / u_jj
		for k = j : n do
			u_ik = u_ik - l_ij * u_jk
		end for
	end for
end for