线性方程组求解
数值解这么走下去，却不好好弄弄关于线性方程组的求解，总感觉很别扭，既然《凸优化》也很详细地介绍了这一块东西，我就先跳过别的把这一块整一整吧。
容易求解的线性方程组
先讨论\(Ax = b\)很容易求解的情况，即\(A\)为满秩的方阵，方程有唯一的解。
对角矩阵
\(a_{ii}x_i = b_i \Rightarrow x_i = b_i / a_{ii}， a_{ii} \neq 0\)
其中\(a_{ij}\)为矩阵\(A\)的第\(i\)行，第\(j\)列元素，下同。
下三角矩阵
下三角矩阵，即\(a_{ij}=0, j > i\)
所以方程的解即为：
上三角矩阵
下三角矩阵采用的是前向代入算法，而上三角矩阵采用的是后向代入或者称为回代算法。情况，或者说推导是类似的，这里不多赘述。
正交矩阵
正交矩阵\(A \in \mathbb{R}^{n \times n}\)的条件是\(A^{T}A = I\)，即\(A^{-1}=A^T\)，所以方程的解是\(x = A^Tb\)。如果\(A\)具有特殊的结构，可以进一步简化运算。
排列矩阵
令\(\pi = (\pi_1, \ldots,\pi_n)\)为\((1, 2, \ldots, n)\)的一种排列。相应的排列矩阵\(A \in \mathbb{R}^{n \times n}\)定义为：
于是可以得到：
排列矩阵的逆矩阵就是\(A^T\)，由此可知排列矩阵是正交矩阵。
因式分解求解方法
求解\(Ax = b\)的基本途径是将\(A\)表示为一系列非奇异矩阵的乘积：
因此：
我们可以从右往左一步一步地来求解。
求解多个右边项的方程组
假设我们需要求解方程组：
求解这m个问题等价于：
其中：
\(\mathrm{LU, Cholesky}\) 和 \(\mathrm{LDL^T}\)因式分解
每一个非奇异分解\(A \in \mathbb{R}^{n \times n}\)都可以因式分解为：
其中\(P \in \mathbb{R}^{n \times n}\)是排列矩阵，\(L \in \mathbb{R}^{n \times n}\)是单位下三角矩阵，而\(U \in \mathbb{R}^{n \times n}\)是非奇异上三角矩阵。
Gauss消元法
我们定义\(A_0 = A, A_1, A_2, \ldots, A_{n-1}\)表示第\(r\)步消元后的系数矩阵。相应的，我们设计一个第\(r\)步消元的初等矩阵\(N_r\)，这个矩阵的除了第\(r\)列外，与单位矩阵无异，第\(r\)列为：
于是，\(A_r = N_r A_{r-1}\)的\(a_{jr}=0,j>r\)，显然，如果顺利的话（因为可能出现\(a_{rr}^{r-1}=0\)的情况），进行\(n-1\)步消元后，矩阵就化为上三角矩阵了。
于是：
其中\(N\)为单位下三角矩阵（下三角矩阵的逆为下三角矩阵，下三角矩阵的乘积为下三角矩阵）。值得一提的是，如果这种分解存在，那么它是唯一的。另外，在《代数特征值问题》一书中，给出了\(L和U\)各个元素的显示表达式。
接下来，我们再讨论一下如何应对\(a_{rr}^{r-1}=0\)的情况。我们有一个最初的假设，即\(A\)是满秩的，虽然这个条件并非必要的（如果没有这个条件，那么就需要在最后判断是否有解）。
经过\(r\)步消元后（假设顺利进行了），那么\(A_{n-r, r}\)为\(0\)矩阵，\(A_{r,r}\)为上三角矩阵。现在，如果\(A_{n-r, n-r}\)的首元素\(a_{r+1, r+1}\)为0，而且\(t = \arg \max \{|a_{i,r+1}||i>r+1\}\)。注意\(a_{t, t+1}\neq 0\)，否则就与我们的满秩条件相矛盾了。当然，如果撇去假设，真的出现了这种情况，我们只需让\(N_{r+1}=I\)即可，即跳过这一次。最后，我们这一次选择的变换是\(N_{r+1}I_{r+1, t}\)。其中\(I_{t+1, t}\)是指第\(r+1\)行与\(t\)行交换的初等矩阵。
为了便于说明，我们以\(n=4\)为例：
于是
这也是最开始的\(A = PLU\)的由来。\(\widetilde{N}\)是下三角矩阵的证明比较简单，这里便不给出证明了。另外值得说明的一点是，我们对于\(t\)的选择，这么选择的原因是出于数值的稳定性（保证\(N_r\)的元素的绝对值都小于\(1\)）
Cholesky 因式分解
如果\(A \in \mathbb{R}^{n \times n}\)是对称正定矩阵，那么它可以因式分解为：
其中\(L\)是下三角非奇异矩阵，对角元素均为正数。这种分解可以看成\(LU\)分解的一种特例，不多赘述。
稀疏矩阵的Cholesky因式分解
当\(A\)是对称正定稀疏矩阵时，通常可以因式分解为：
举个例子便于理解：
其中\(D = LL^T\)，如果\(D\)为正对角矩阵，那么\(L\)的对角线元素便直接可以获得了。
\(LDL^T\) 因式分解
每个非奇异对称矩阵\(A\)都能因式分解为：
其中\(P\)是排列矩阵，\(L\)是对角均为正数的下三角矩阵，\(D\)是块对角矩阵，对角块为\(1 \times 1\)和\(2 \times 2\)的非奇异矩阵。
这地方就不做分析了，因为自己没怎么细看这部分过。
分块消元和Schur补
将\(x \in \mathbb{R}^n\)分成俩块：
其中\(x_1 \in \mathbb{R}^{n_1},x_2 \in \mathbb{R}^{n_2}\)。
那么线性方程组可以这样表示：
其中\(A_{11} \in \mathbb{R}^{n_1 \times n_1},A_{22} \in \mathbb{R}^{n_2 \times n_2}\)且假设\(A_{11}\)可逆。
由第一个方程可以获得：
代入第二个方程可以得到：
注意到，\(S = A_{22}-A_{21}A_{11}^{-1}A_{12}\)即\(A_{11}\)的Schur补。
由上面的启发，我们可以先计算\(x_2\)再计算\(x_1\)，虽然这种方法对于稠密的无结构矩阵而言没有什么优点，但如果要消去的变量对于的子矩阵容易因式分解，这种方法会很有效。
逆矩阵引理
分块消元法的想法是先消去部分变量，然后求解包含这些变量的Schur补的小方程组。同样的想法可以反向应用：如果讲某个矩阵视为Schur补，就可以引入新变量，然后形成并求解一个大方程组。很多情况下这样做没有好处，因为我们最终要求解一个更大的方程组。但是，如果所形成的大方程组具有可以利用的特殊结构，引入新变量就可能导致更加有效的求解方法。最经常利用的是可以从大方程组中消去另一部分变量的情况。
假设有下面的线性方程组：
其中\(A \in \mathbb{R}^{n \times n}\)非奇异，\(B \in \mathbb{R}^{n \times p}\),\(C \in \mathbb{R}^{p \times n}\)。我们引入新变量\(y=Cx\)，并将方程组重新写成
即：
注意到\(A+BC\)是大矩阵中\(-I\)的Schur补。容易看出，当\(A,B,C\)相当稀疏，而\(A+BC\)稀疏性很差的时候，解大方程组或许比原来的更加有效。
代码
import numpy as np
class LinearEqu: # 要求矩阵A为满秩方阵
def __init__(self, A, b):
self.m, self.n = A.shape
assert self.n == len(b), "the dimensions don't match"
assert self.m == self.n, "full-rank and row-column equal matrix required"
self.A = np.array(A, dtype=float)
self.b = np.array(b, dtype=float)
@property
def rank(self):
"""返回矩阵的秩"""
return np.linalg.matrix_rank(self.A)
@property
def extendrank(self):
"""返回[A, b]的秩"""
b = self.b.reshape(-1, 1)
return np.linalg.matrix_rank(np.hstack((self.A, b)))
@property
def diagonal(self):
assert self.rank == self.extendrank, "No solution"
assert self.rank == self.n, "A is not a full-rank matrix"
"""
下面这部分是对矩阵A对角性质的考察，但是想到，万一我只是希望利用一下
对角元素呢，所以这部分引掉。
index = np.fromfunction(lambda i, j: i!=j, (self.n ,self.n))
remain = self.A[index] == 0.
if not np.all(remain):
raise TypeError("matrix A is not diagonal...")
"""
diag_A = np.diag(1 / np.diag(self.A))
return diag_A @ b
@property
def low_triangle(self):
"""对下三角矩阵求解"""
assert self.rank == self.extendrank, "No solution"
assert self.rank == self.n, "A is not a full-rank matrix"
index = np.fromfunction(lambda i,j: i < j, (self.n, self.n))
remain = self.A[index] == 0.
if not np.all(remain): #这部分我们直接给出了检查
raise TypeError("matrix A is not low-triangle...")
x = np.zeros(self.n, dtype=float)
for i in range(self.n):
if not i:
x[i] = self.b[i] / self.A[i, i]
else:
residual = self.A[i, :i] @ x[:i]
x[i] = (self.b[i] - residual) / self.A[i, i]
return x
@property
def up_triangle(self):
"""对上三角形矩阵求解"""
assert self.rank == self.extendrank, "No solution"
assert self.rank == self.n, "A is not a full-rank matrix"
index = np.fromfunction(lambda i,j: i > j, (self.n, self.n))
remain = self.A[index] == 0.
if not np.all(remain): #这部分我们直接给出了检查
raise TypeError("matrix A is not up-triangle...")
x = np.zeros(self.n, dtype=float)
for i in range(self.n):
if not i:
x[self.n-1] = self.b[-1] / self.A[-1, -1]
else:
k = self.n - 1 - i
residual = self.A[k, k+1:] @ x[k+1:]
x[k] = (self.b[k] - residual) / self.A[k, k]
return x
@property
def orthogonal(self):
"""正交矩阵"""
assert self.rank == self.extendrank, "No solution"
assert self.rank == self.n, "A is not a full-rank matrix"
"""
我们的确可以给出检查，只需：
if np.sum(np.abs(self.A @ self.A.T - np.diag(np.ones(self.n)))) > 1e-5:
raise TypeError("A is not orthogonal matrix...")
因为会存在浪费计算的问题，这里就引掉吧。
"""
return self.A.T @ self.b
@property
def gauss(self):
"""利用高斯消元法求解"""
assert self.rank == self.extendrank, "No solution"
assert self.rank == self.n, "A is not a full-rank matrix"
def find_max(A, r):
vector = A[r:, r]
max_pos = np.argmax(np.abs(vector)) + r
return max_pos
A = np.array(self.A, dtype=float)
b = np.array(self.b, dtype=float)
for r in range(self.n - 1):
max_pos = find_max(A, r) #寻找最大的点
vector = np.array(A[r]) #替换 这么做的原因是多维ndarray似乎不支持a,b=b,a
A[r] = A[max_pos]
A[max_pos] = vector
b[r], b[max_pos] = b[max_pos], b[r]
N_r = np.diag(np.ones(self.n, dtype=float))
N_r[r:, r] = -np.array(A[r:, r]) / A[r, r]
N_r[r, r] = 1.
A = N_r @ A #更新A
b = N_r @ b #更新b
temp = LinearEqu(A, b)
return temp.up_triangle
