线性方程组求解

目录

• 容易求解的线性方程组
  • 对角矩阵
  • 下三角矩阵
  • 上三角矩阵
  • 正交矩阵
  • 排列矩阵
• 因式分解求解方法
  • 求解多个右边项的方程组
• $\mathrm{LU, Cholesky}$ 和 $\mathrm{LDL^T}$因式分解
  • Gauss消元法
  • Cholesky 因式分解
    • 稀疏矩阵的Cholesky因式分解
  • $LDL^T$ 因式分解
• 分块消元和Schur补
  • 逆矩阵引理
• 代码

《Convex Optimization》

数值解这么走下去，却不好好弄弄关于线性方程组的求解，总感觉很别扭，既然《凸优化》也很详细地介绍了这一块东西，我就先跳过别的把这一块整一整吧。

容易求解的线性方程组

先讨论$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$

$$\left [ \begin{array}{cccc} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{array} \right] \left [ \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right] = \left [ \begin{array}{c} b_1 \\ b_2 \\ \vdots \\ b_n \end{array} \right]$$

所以方程的解即为：

$$x_1 := b_1 / a_{11} \\ x_2 := (b2 - a_{21}x_1 / a_{22} \\ \vdots \\ x_n := (b_n - a_{n1}x_1 - a_{n2}x_2 - \cdots - a_{n,n-1}x_{n-1} / a_{nn}$$

上三角矩阵

下三角矩阵采用的是前向代入算法，而上三角矩阵采用的是后向代入或者称为回代算法。情况，或者说推导是类似的，这里不多赘述。

正交矩阵

正交矩阵$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_{ij}= \left \{ \begin{array}{ll} 1 & j = \pi_i \\ 0 & 其他情况 \end{array} \right .$$

于是可以得到：

$$Ax = (x_{\pi_1}, \ldots, x_{\pi_n})$$

排列矩阵的逆矩阵就是$A^T$，由此可知排列矩阵是正交矩阵。

因式分解求解方法

求解$Ax = b$的基本途径是将$A$表示为一系列非奇异矩阵的乘积：

$$A = A_1 A_2 \cdots A_k$$

因此：

$$x = A^{-1}b = A_k^{-1} A_{k-1}^{-1}\cdots A_1^{-1}b$$

我们可以从右往左一步一步地来求解。

求解多个右边项的方程组

假设我们需要求解方程组：

$$Ax_1 = b_1, Ax_2 = b_2, \cdots, Ax_m = b_m$$

求解这m个问题等价于：

$$X = A^{-1}B$$

其中：

$$X = [x_1, x_2, \cdots, x_m] \in \mathbb{R}^{n \times m} \quad B=[b_1, b_2,\cdots, b_m] \in \mathbb{R}^{n \times m}$$

$\mathrm{LU, Cholesky}$ 和 $\mathrm{LDL^T}$因式分解

每一个非奇异分解$A \in \mathbb{R}^{n \times n}$都可以因式分解为：

$$A = PLU$$

其中$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$列为：

$$[0, 0, \ldots, 1, -a_{r+1,r}^{r-1}/a_{r,r}^{r-1}, -a_{r+2, r}^{r-1}/a_{r,r}^{r-1}, \ldots, -a_{n,r}^{r-1}/a_{r,r}^{r-1}]^T$$

于是，$A_r = N_r A_{r-1}$的$a_{jr}=0,j>r$，显然，如果顺利的话（因为可能出现$a_{rr}^{r-1}=0$的情况），进行$n-1$步消元后，矩阵就化为上三角矩阵了。

$$N_{n-1} \cdots N_2 N_1 A_0 = A_{n-1}, N_{n-1} \cdots N_2 N_1 b_0 = b_{n+1}$$

于是：

$$A_0 = N_1^{-1}N_2^{-1}\cdots N_{n-1}^{-1} A_{n-1} = N A_{n-1}$$

其中$N$为单位下三角矩阵（下三角矩阵的逆为下三角矩阵，下三角矩阵的乘积为下三角矩阵）。值得一提的是，如果这种分解存在，那么它是唯一的。另外，在《代数特征值问题》一书中，给出了$L和U$各个元素的显示表达式。
接下来，我们再讨论一下如何应对$a_{rr}^{r-1}=0$的情况。我们有一个最初的假设，即$A$是满秩的，虽然这个条件并非必要的（如果没有这个条件，那么就需要在最后判断是否有解）。

$$A_r = \left [ \begin{array}{ll} A_{r,r} & A_{r, n-r} \\ A_{n-r, r} & A_{n-r, n-r} \end{array} \right]$$

经过$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$为例：

$$\begin{array}{ll} A_3 &= N_3I_{3,3'}N_2I_{2,2'}N_1I_{1,1'}\\ & =N_3I_{3,3'}N_2(I_{3,3'}I_{3,3'})I_{2,2'}N_1(I_{2,2'}I_{3,3'}I_{3,3'}I_{2,2'})I_{1,1'}A_0 \\ & =N_3(I_{3,3'}N_2I_{3,3'})(I_{3,3'}I_{2,2'}N_1I_{2,2'}I_{3,3'})(I_{3,3'}I_{2,2'}I_{1,1'}A_0 )\\ &= \widetilde{N}_{3}\widetilde{N}_{2}\widetilde{N}_{1}P^TA_0 \end{array}$$

于是

$$A_0 = P\widetilde{N}A_3$$

这也是最开始的$A = PLU$的由来。$\widetilde{N}$是下三角矩阵的证明比较简单，这里便不给出证明了。另外值得说明的一点是，我们对于$t$的选择，这么选择的原因是出于数值的稳定性（保证$N_r$的元素的绝对值都小于$1$）

Cholesky 因式分解

如果$A \in \mathbb{R}^{n \times n}$是对称正定矩阵，那么它可以因式分解为：

$$A = LL^T$$

其中$L$是下三角非奇异矩阵，对角元素均为正数。这种分解可以看成$LU$分解的一种特例，不多赘述。

稀疏矩阵的Cholesky因式分解

当$A$是对称正定稀疏矩阵时，通常可以因式分解为：

$$A = PLL^TP^T$$

举个例子便于理解：

$$A = \left [ \begin{array}{ll} 1 &u^T \\ u & D \end{array} \right] = \left [ \begin{array}{ll} 1 &0 \\ u & L \end{array} \right] \left [ \begin{array}{ll} 1 & u^T \\ 0 & L^T \end{array} \right]$$

其中$D = LL^T$，如果$D$为正对角矩阵，那么$L$的对角线元素便直接可以获得了。

$LDL^T$ 因式分解

每个非奇异对称矩阵$A$都能因式分解为：

$$A = PLDL^TP^T$$

其中$P$是排列矩阵，$L$是对角均为正数的下三角矩阵，$D$是块对角矩阵，对角块为$1 \times 1$和$2 \times 2$的非奇异矩阵。
这地方就不做分析了，因为自己没怎么细看这部分过。

分块消元和Schur补

将$x \in \mathbb{R}^n$分成俩块：

$$x = \left [ \begin{array}{c} x_1\\ x_2\\ \end{array} \right]$$

其中$x_1 \in \mathbb{R}^{n_1},x_2 \in \mathbb{R}^{n_2}$。
那么线性方程组可以这样表示：

$$\left [ \begin{array}{ll} A_{11} & A_{12} \\ A_{21} & A_{22} \end{array} \right] \left [ \begin{array}{l} x_1 \\ x_2 \end{array} \right] = \left [ \begin{array}{l} b_1 \\ b_2 \end{array} \right]$$

其中$A_{11} \in \mathbb{R}^{n_1 \times n_1},A_{22} \in \mathbb{R}^{n_2 \times n_2}$且假设$A_{11}$可逆。
由第一个方程可以获得：

$$x_1 = A_{11}^{-1} (b_1 - A_{12} x_2)$$

代入第二个方程可以得到：

$$(A_{22} - A_{21} A_{11}^{-1} A_{12}) x_2 = b_2 - A_{21} A_{11}^{-1} b_1$$

注意到，$S = A_{22}-A_{21}A_{11}^{-1}A_{12}$即$A_{11}$的Schur补。
由上面的启发，我们可以先计算$x_2$再计算$x_1$，虽然这种方法对于稠密的无结构矩阵而言没有什么优点，但如果要消去的变量对于的子矩阵容易因式分解，这种方法会很有效。

逆矩阵引理

分块消元法的想法是先消去部分变量，然后求解包含这些变量的Schur补的小方程组。同样的想法可以反向应用：如果讲某个矩阵视为Schur补，就可以引入新变量，然后形成并求解一个大方程组。很多情况下这样做没有好处，因为我们最终要求解一个更大的方程组。但是，如果所形成的大方程组具有可以利用的特殊结构，引入新变量就可能导致更加有效的求解方法。最经常利用的是可以从大方程组中消去另一部分变量的情况。
假设有下面的线性方程组：

$$(A + BC)x = b$$

其中$A \in \mathbb{R}^{n \times n}$非奇异，$B \in \mathbb{R}^{n \times p}$,$C \in \mathbb{R}^{p \times n}$。我们引入新变量$y=Cx$，并将方程组重新写成

$$Ax + By = b, \quad y = Cx$$

即：

$$\left [ \begin{array}{ll} A & B \\ C & -I \end{array} \right ] \left [ \begin{array}{l} x \\ y \end{array} \right ] = \left [ \begin{array}{l} b \\ 0 \end{array} \right ]$$

注意到$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