线性代数(10)——初等矩阵和矩阵可逆性(下)

矩阵的LU分解

概念

矩阵的分解实际上就是将一个矩阵分解成为几个矩阵乘积的形式。不同类型的矩阵分解有不同的目的，如矩阵的LU分解的目的是为了提高计算效率。

LU分解实际上是将一个矩阵$A$分解成为两个矩阵$L$和$U$的乘积，即 $A=L\cdot U$，

L：Lower triangle matrix 下三角矩阵
U：Upper triangle matrix 上三角矩阵

以方阵为例，演示上三角矩阵和下三角矩阵，

1. 下三角矩阵
   $$\left(\begin{array}{cccc}\ast & 0& 0& 0\\ \ast & \ast & 0& 0\\ \ast & \ast & \ast & 0\\ \ast & \ast & \ast & \ast \end{array}\right)$$
   下三角矩阵：有效信息均在主对角线及其下方，而上方的元素均为0。
2. 上三角矩阵
   $$\left(\begin{array}{cccc}\ast & \ast & \ast & \ast \\ 0& \ast & \ast & \ast \\ 0& 0& \ast & \ast \\ 0& 0& 0& \ast \end{array}\right)$$
   上三角矩阵：有效信息均在主对角线及其上方，而下方的元素均为0。

在矩阵分解中，经常是把矩阵分解成一个“单位下三角矩阵”(即主对角线元素均为1) 和一个上三角矩阵，不保证上三角矩阵的主对角线元素为1。

回顾之前的Gauss-Jordan消元法，前向过程最终是将矩阵经过初等变换生成了一个上三角矩阵。因为初等变换的过程又可以视为将一系列初等矩阵乘在矩阵左边，所以可以写为，
$$E_p\cdot ...\cdot E_3\cdot E_2\cdot E1\cdot A=U$$
等式左边只保留$A$，
$$E_1^{-1}\cdot E_2^{-1}\cdot E_3^{-1}\cdot ...\cdot E_p^{-1}\cdot E_p\cdot ...\cdot E_3\cdot E_2\cdot E1\cdot A=E_1^{-1}\cdot E_2^{-1}\cdot E_3^{-1}\cdot ...\cdot E_p^{-1}\cdot U$$
$$I\cdot A=E_1^{-1}\cdot E_2^{-1}\cdot E_3^{-1}\cdot ...\cdot E_p^{-1}\cdot U$$
因为最终想使得 $A=L\cdot U$，如果一切顺利(并不是所有矩阵都能够进行LU分解)，
$$L=E_1^{-1}\cdot E_2^{-1}\cdot E_3^{-1}\cdot ...\cdot E_p^{-1}$$

LU分解过程

具体的举一个例子进行演示，
$$\left(\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\\ 3& -3& 5\end{array}\right)⟹\left(\begin{array}{ccc}1& 2& 3\\ 0& -3& -6\\ 0& -9& -4\end{array}\right)⟹\left(\begin{array}{ccc}1& 2& 3\\ 0& -3& -6\\ 0& 0& 14\end{array}\right)$$
通过Gauss-Jordan消元法可以将原矩阵转化为上面这样的上三角矩阵。将每一步所对应的初等矩阵求逆后相乘得到单位下三角矩阵，
$$\begin{gathered} {\left(\begin{array}{ccc}1& 0& 0\\ -4& 1& 0\\ 0& 0& 1\end{array}\right)}^{-1}\cdot {\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ -3& 0& 1\end{array}\right)}^{-1}\cdot {\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& -3& 1\end{array}\right)}^{-1} \\ =\left(\begin{array}{ccc}1& 0& 0\\ 4& 1& 0\\ 0& 0& 1\end{array}\right)\cdot \left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 3& 0& 1\end{array}\right)\cdot \left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 3& 1\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 4& 1& 0\\ 3& 3& 1\end{array}\right) \end{gathered}$$
两个矩阵满足，
$$\left(\begin{array}{ccc}1& 0& 0\\ 4& 1& 0\\ 3& 3& 1\end{array}\right)\cdot \left(\begin{array}{ccc}1& 2& 3\\ 0& -3& -6\\ 0& 0& 14\end{array}\right)=\left(\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\\ 3& -3& 5\end{array}\right)$$
之所以很顺利的完成了矩阵A的LU分解的原因是因为整个消元过程中仅涉及到某一行减去另一行的某个倍数这样单一的初等操作，而不需要交换两行的位置。

如果一个矩阵在消元的过程中需要交换两行的位置，则这个矩阵是无法进行LU分解的。

LU分解的意义

LU分解的提出目的是加速线性系统求解的效率，假设在解一个线性系统 $Ax=b$，

1. LU分解的复杂度为 $O(0.5n^3)-O(n^3)$
   主对角线下的元素大约为 $0.5n^2$个。对于这些元素在高斯消元的过程中都要消为0，整个过程是用矩阵的一行加减$k$倍的另一行。在加减的过程中，每次加减都是对一行的$n$个元素进行操作，所以整个事件的消耗为$O(0.5n^3)$。
2. 将矩阵$A$进行LU分解后，整个线性系统就变成了 $L\cdot U\cdot x=b$。令$U\cdot x=y$，则整个线性系统转化为$L\cdot y=b$。
3. 已知$L$和$b$，可以求出$y$，求$y$的过程的复杂度约为 $O(0.5n^2)-O(n^2)$。因为对于一个标准下三角矩阵而言，通过它和$b$求取一个矩阵是一个十分简单的过程。
4. $U\cdot x=y$，目前$U$和$y$都已知，求取$x$就变得十分容易，大概也需要复杂度为$O(0.5n^2)-O(n^2)$的操作。因为上三角矩阵具有与下三角矩阵相同的性质。

整体上使用了复杂度为 $O(0.5n^3)+2\cdot O(0.5n^2)$ 的操作。

对比传统的求解线性系统的形式，

1. 传统形式所需操作的复杂度为$O(2\cdot n^3)$，因为增广矩阵中有 $2n^2$个元素，某一行加减k倍的另一行将矩阵化为行最简形式时，每一步都是在对一行中的n个元素进行操作，故复杂度为$O(2n^3)$。
2. $x=A^{-1}\cdot b$，乘法过程所需操作的复杂度为 $O(n^2)$

整体上传统形式的复杂度为 $O(2n^3)+O(n^2)$。

代码实现

# -*- coding:utf-8 -*-
from .matrix import Matrix
from .vector import Vector
from .global import EPSILON


def lu(matrix):
	assert matrix.row_num == matrix.col_num()
	n = matrix.row_num()
	# 将matrix以行向量的形式存储
    A = [matrix.row_vector(i) for i in range(n)]
    # 定义下三角矩阵L
    L = [[1.0 if i==j else 0.0 for i in range(n)] for j in range(n)]

	# 矩阵A的高斯消元过程
    for i in range(n):
        # 判断A[i][i]位置上的元素是否可以是主元
        # 如果是0，则无法进行LU分解
        if A[i][i] < EPSILON:
            return None, None
        else:
            for j in range(i+1, n):
                k = self.A[j][i] / self.A[i][i]
                A[j] -= k * A[i]
                L[j][i] = p

    return Matrix(L), Matrix([A[i].underlying_list() for i in range(n)])

非方阵LU分解

之前讲解的LU分解都是针对方阵而言，具体形式可以表示为，

一个$n$阶方阵，其分解得到的 $L$矩阵和 $U$矩阵都是 $n$阶的方阵。但是，实际上LU分解可以用于非方阵，如，

若$A$是一个行数小于列数的矩阵，则分解的结果如上；反之，如果$A$是一个行数大于列数的矩阵，则分解结果可以表示为如下形式，