文章目录
数值解这么走下去,却不好好弄弄关于线性方程组的求解,总感觉很别扭,既然《凸优化》也很详细地介绍了这一块东西,我就先跳过别的把这一块整一整吧。
容易求解的线性方程组
先讨论 A x = b Ax = b Ax=b很容易求解的情况,即 A A A为满秩的方阵,方程有唯一的解。
对角矩阵
a
i
i
x
i
=
b
i
⇒
x
i
=
b
i
/
a
i
i
,
a
i
i
≠
0
a_{ii}x_i = b_i \Rightarrow x_i = b_i / a_{ii}, a_{ii} \neq 0
aiixi=bi⇒xi=bi/aii,aii̸=0
其中
a
i
j
a_{ij}
aij为矩阵
A
A
A的第
i
i
i行,第
j
j
j列元素,下同。
下三角矩阵
下三角矩阵,即
a
i
j
=
0
,
j
>
i
a_{ij}=0, j > i
aij=0,j>i
[
a
11
0
⋯
0
a
21
a
22
⋯
0
⋮
⋮
⋱
⋮
a
n
1
a
n
2
⋯
a
n
n
]
[
x
1
x
2
⋮
x
n
]
=
[
b
1
b
2
⋮
b
n
]
\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]
⎣⎢⎢⎢⎡a11a21⋮an10a22⋮an2⋯⋯⋱⋯00⋮ann⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡b1b2⋮bn⎦⎥⎥⎥⎤
所以方程的解即为:
x
1
:
=
b
1
/
a
11
x
2
:
=
(
b
2
−
a
21
x
1
/
a
22
⋮
x
n
:
=
(
b
n
−
a
n
1
x
1
−
a
n
2
x
2
−
⋯
−
a
n
,
n
−
1
x
n
−
1
/
a
n
n
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}
x1:=b1/a11x2:=(b2−a21x1/a22⋮xn:=(bn−an1x1−an2x2−⋯−an,n−1xn−1/ann
上三角矩阵
下三角矩阵采用的是前向代入算法,而上三角矩阵采用的是后向代入或者称为回代算法。情况,或者说推导是类似的,这里不多赘述。
正交矩阵
正交矩阵 A ∈ R n × n A \in \mathbb{R}^{n \times n} A∈Rn×n的条件是 A T A = I A^{T}A = I ATA=I,即 A − 1 = A T A^{-1}=A^T A−1=AT,所以方程的解是 x = A T b x = A^Tb x=ATb。如果 A A A具有特殊的结构,可以进一步简化运算。
排列矩阵
令
π
=
(
π
1
,
…
,
π
n
)
\pi = (\pi_1, \ldots,\pi_n)
π=(π1,…,πn)为
(
1
,
2
,
…
,
n
)
(1, 2, \ldots, n)
(1,2,…,n)的一种排列。相应的排列矩阵
A
∈
R
n
×
n
A \in \mathbb{R}^{n \times n}
A∈Rn×n定义为:
A
i
j
=
{
1
j
=
π
i
0
其
他
情
况
A_{ij}= \left \{ \begin{array}{ll} 1 & j = \pi_i \\ 0 & 其他情况 \end{array} \right .
Aij={10j=πi其他情况
于是可以得到:
A
x
=
(
x
π
1
,
…
,
x
π
n
)
Ax = (x_{\pi_1}, \ldots, x_{\pi_n})
Ax=(xπ1,…,xπn)
排列矩阵的逆矩阵就是
A
T
A^T
AT,由此可知排列矩阵是正交矩阵。
因式分解求解方法
求解
A
x
=
b
Ax = b
Ax=b的基本途径是将
A
A
A表示为一系列非奇异矩阵的乘积:
A
=
A
1
A
2
⋯
A
k
A = A_1 A_2 \cdots A_k
A=A1A2⋯Ak
因此:
x
=
A
−
1
b
=
A
k
−
1
A
k
−
1
−
1
⋯
A
1
−
1
b
x = A^{-1}b = A_k^{-1} A_{k-1}^{-1}\cdots A_1^{-1}b
x=A−1b=Ak−1Ak−1−1⋯A1−1b
我们可以从右往左一步一步地来求解。
求解多个右边项的方程组
假设我们需要求解方程组:
A
x
1
=
b
1
,
A
x
2
=
b
2
,
⋯
 
,
A
x
m
=
b
m
Ax_1 = b_1, Ax_2 = b_2, \cdots, Ax_m = b_m
Ax1=b1,Ax2=b2,⋯,Axm=bm
求解这m个问题等价于:
X
=
A
−
1
B
X = A^{-1}B
X=A−1B
其中:
X
=
[
x
1
,
x
2
,
⋯
 
,
x
m
]
∈
R
n
×
m
B
=
[
b
1
,
b
2
,
⋯
 
,
b
m
]
∈
R
n
×
m
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}
X=[x1,x2,⋯,xm]∈Rn×mB=[b1,b2,⋯,bm]∈Rn×m
L U , C h o l e s k y \mathrm{LU, Cholesky} LU,Cholesky 和 L D L T \mathrm{LDL^T} LDLT因式分解
每一个非奇异分解
A
∈
R
n
×
n
A \in \mathbb{R}^{n \times n}
A∈Rn×n都可以因式分解为:
A
=
P
L
U
A = PLU
A=PLU
其中
P
∈
R
n
×
n
P \in \mathbb{R}^{n \times n}
P∈Rn×n是排列矩阵,
L
∈
R
n
×
n
L \in \mathbb{R}^{n \times n}
L∈Rn×n是单位下三角矩阵,而
U
∈
R
n
×
n
U \in \mathbb{R}^{n \times n}
U∈Rn×n是非奇异上三角矩阵。
Gauss消元法
我们定义
A
0
=
A
,
A
1
,
A
2
,
…
,
A
n
−
1
A_0 = A, A_1, A_2, \ldots, A_{n-1}
A0=A,A1,A2,…,An−1表示第
r
r
r步消元后的系数矩阵。相应的,我们设计一个第
r
r
r步消元的初等矩阵
N
r
N_r
Nr,这个矩阵的除了第
r
r
r列外,与单位矩阵无异,第
r
r
r列为:
[
0
,
0
,
…
,
1
,
−
a
r
+
1
,
r
r
−
1
/
a
r
,
r
r
−
1
,
−
a
r
+
2
,
r
r
−
1
/
a
r
,
r
r
−
1
,
…
,
−
a
n
,
r
r
−
1
/
a
r
,
r
r
−
1
]
T
[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
[0,0,…,1,−ar+1,rr−1/ar,rr−1,−ar+2,rr−1/ar,rr−1,…,−an,rr−1/ar,rr−1]T
于是,
A
r
=
N
r
A
r
−
1
A_r = N_r A_{r-1}
Ar=NrAr−1的
a
j
r
=
0
,
j
>
r
a_{jr}=0,j>r
ajr=0,j>r,显然,如果顺利的话(因为可能出现
a
r
r
r
−
1
=
0
a_{rr}^{r-1}=0
arrr−1=0的情况),进行
n
−
1
n-1
n−1步消元后,矩阵就化为上三角矩阵了。
N
n
−
1
⋯
N
2
N
1
A
0
=
A
n
−
1
,
N
n
−
1
⋯
N
2
N
1
b
0
=
b
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}
Nn−1⋯N2N1A0=An−1,Nn−1⋯N2N1b0=bn+1
于是:
A
0
=
N
1
−
1
N
2
−
1
⋯
N
n
−
1
−
1
A
n
−
1
=
N
A
n
−
1
A_0 = N_1^{-1}N_2^{-1}\cdots N_{n-1}^{-1} A_{n-1} = N A_{n-1}
A0=N1−1N2−1⋯Nn−1−1An−1=NAn−1
其中
N
N
N为单位下三角矩阵(下三角矩阵的逆为下三角矩阵,下三角矩阵的乘积为下三角矩阵)。值得一提的是,如果这种分解存在,那么它是唯一的。另外,在《代数特征值问题》一书中,给出了
L
和
U
L和U
L和U各个元素的显示表达式。
接下来,我们再讨论一下如何应对
a
r
r
r
−
1
=
0
a_{rr}^{r-1}=0
arrr−1=0的情况。我们有一个最初的假设,即
A
A
A是满秩的,虽然这个条件并非必要的(如果没有这个条件,那么就需要在最后判断是否有解)。
A
r
=
[
A
r
,
r
A
r
,
n
−
r
A
n
−
r
,
r
A
n
−
r
,
n
−
r
]
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]
Ar=[Ar,rAn−r,rAr,n−rAn−r,n−r]
经过
r
r
r步消元后(假设顺利进行了),那么
A
n
−
r
,
r
A_{n-r, r}
An−r,r为
0
0
0矩阵,
A
r
,
r
A_{r,r}
Ar,r为上三角矩阵。现在,如果
A
n
−
r
,
n
−
r
A_{n-r, n-r}
An−r,n−r的首元素
a
r
+
1
,
r
+
1
a_{r+1, r+1}
ar+1,r+1为0,而且
t
=
arg
max
{
∣
a
i
,
r
+
1
∣
∣
i
>
r
+
1
}
t = \arg \max \{|a_{i,r+1}||i>r+1\}
t=argmax{∣ai,r+1∣∣i>r+1}。注意
a
t
,
t
+
1
≠
0
a_{t, t+1}\neq 0
at,t+1̸=0,否则就与我们的满秩条件相矛盾了。当然,如果撇去假设,真的出现了这种情况,我们只需让
N
r
+
1
=
I
N_{r+1}=I
Nr+1=I即可,即跳过这一次。最后,我们这一次选择的变换是
N
r
+
1
I
r
+
1
,
t
N_{r+1}I_{r+1, t}
Nr+1Ir+1,t。其中
I
t
+
1
,
t
I_{t+1, t}
It+1,t是指第
r
+
1
r+1
r+1行与
t
t
t行交换的初等矩阵。
为了便于说明,我们以
n
=
4
n=4
n=4为例:
A
3
=
N
3
I
3
,
3
′
N
2
I
2
,
2
′
N
1
I
1
,
1
′
=
N
3
I
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
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
N
~
2
N
~
1
P
T
A
0
\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}
A3=N3I3,3′N2I2,2′N1I1,1′=N3I3,3′N2(I3,3′I3,3′)I2,2′N1(I2,2′I3,3′I3,3′I2,2′)I1,1′A0=N3(I3,3′N2I3,3′)(I3,3′I2,2′N1I2,2′I3,3′)(I3,3′I2,2′I1,1′A0)=N
3N
2N
1PTA0
于是
A
0
=
P
N
~
A
3
A_0 = P\widetilde{N}A_3
A0=PN
A3
这也是最开始的
A
=
P
L
U
A = PLU
A=PLU的由来。
N
~
\widetilde{N}
N
是下三角矩阵的证明比较简单,这里便不给出证明了。另外值得说明的一点是,我们对于
t
t
t的选择,这么选择的原因是出于数值的稳定性(保证
N
r
N_r
Nr的元素的绝对值都小于
1
1
1)
Cholesky 因式分解
如果
A
∈
R
n
×
n
A \in \mathbb{R}^{n \times n}
A∈Rn×n是对称正定矩阵,那么它可以因式分解为:
A
=
L
L
T
A = LL^T
A=LLT
其中
L
L
L是下三角非奇异矩阵,对角元素均为正数。这种分解可以看成
L
U
LU
LU分解的一种特例,不多赘述。
稀疏矩阵的Cholesky因式分解
当
A
A
A是对称正定稀疏矩阵时,通常可以因式分解为:
A
=
P
L
L
T
P
T
A = PLL^TP^T
A=PLLTPT
举个例子便于理解:
A
=
[
1
u
T
u
D
]
=
[
1
0
u
L
]
[
1
u
T
0
L
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]
A=[1uuTD]=[1u0L][10uTLT]
其中
D
=
L
L
T
D = LL^T
D=LLT,如果
D
D
D为正对角矩阵,那么
L
L
L的对角线元素便直接可以获得了。
L D L T LDL^T LDLT 因式分解
每个非奇异对称矩阵
A
A
A都能因式分解为:
A
=
P
L
D
L
T
P
T
A = PLDL^TP^T
A=PLDLTPT
其中
P
P
P是排列矩阵,
L
L
L是对角均为正数的下三角矩阵,
D
D
D是块对角矩阵,对角块为
1
×
1
1 \times 1
1×1和
2
×
2
2 \times 2
2×2的非奇异矩阵。
这地方就不做分析了,因为自己没怎么细看这部分过。
分块消元和Schur补
将
x
∈
R
n
x \in \mathbb{R}^n
x∈Rn分成俩块:
x
=
[
x
1
x
2
]
x = \left [ \begin{array}{c} x_1\\ x_2\\ \end{array} \right]
x=[x1x2]
其中
x
1
∈
R
n
1
,
x
2
∈
R
n
2
x_1 \in \mathbb{R}^{n_1},x_2 \in \mathbb{R}^{n_2}
x1∈Rn1,x2∈Rn2。
那么线性方程组可以这样表示:
[
A
11
A
12
A
21
A
22
]
[
x
1
x
2
]
=
[
b
1
b
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]
[A11A21A12A22][x1x2]=[b1b2]
其中
A
11
∈
R
n
1
×
n
1
,
A
22
∈
R
n
2
×
n
2
A_{11} \in \mathbb{R}^{n_1 \times n_1},A_{22} \in \mathbb{R}^{n_2 \times n_2}
A11∈Rn1×n1,A22∈Rn2×n2且假设
A
11
A_{11}
A11可逆。
由第一个方程可以获得:
x
1
=
A
11
−
1
(
b
1
−
A
12
x
2
)
x_1 = A_{11}^{-1} (b_1 - A_{12} x_2)
x1=A11−1(b1−A12x2)
代入第二个方程可以得到:
(
A
22
−
A
21
A
11
−
1
A
12
)
x
2
=
b
2
−
A
21
A
11
−
1
b
1
(A_{22} - A_{21} A_{11}^{-1} A_{12}) x_2 = b_2 - A_{21} A_{11}^{-1} b_1
(A22−A21A11−1A12)x2=b2−A21A11−1b1
注意到,
S
=
A
22
−
A
21
A
11
−
1
A
12
S = A_{22}-A_{21}A_{11}^{-1}A_{12}
S=A22−A21A11−1A12即
A
11
A_{11}
A11的Schur补。
由上面的启发,我们可以先计算
x
2
x_2
x2再计算
x
1
x_1
x1,虽然这种方法对于稠密的无结构矩阵而言没有什么优点,但如果要消去的变量对于的子矩阵容易因式分解,这种方法会很有效。
逆矩阵引理
分块消元法的想法是先消去部分变量,然后求解包含这些变量的Schur补的小方程组。同样的想法可以反向应用:如果讲某个矩阵视为Schur补,就可以引入新变量,然后形成并求解一个大方程组。很多情况下这样做没有好处,因为我们最终要求解一个更大的方程组。但是,如果所形成的大方程组具有可以利用的特殊结构,引入新变量就可能导致更加有效的求解方法。最经常利用的是可以从大方程组中消去另一部分变量的情况。
假设有下面的线性方程组:
(
A
+
B
C
)
x
=
b
(A + BC)x = b
(A+BC)x=b
其中
A
∈
R
n
×
n
A \in \mathbb{R}^{n \times n}
A∈Rn×n非奇异,
B
∈
R
n
×
p
B \in \mathbb{R}^{n \times p}
B∈Rn×p,
C
∈
R
p
×
n
C \in \mathbb{R}^{p \times n}
C∈Rp×n。我们引入新变量
y
=
C
x
y=Cx
y=Cx,并将方程组重新写成
A
x
+
B
y
=
b
,
y
=
C
x
Ax + By = b, \quad y = Cx
Ax+By=b,y=Cx
即:
[
A
B
C
−
I
]
[
x
y
]
=
[
b
0
]
\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 ]
[ACB−I][xy]=[b0]
注意到
A
+
B
C
A+BC
A+BC是大矩阵中
−
I
-I
−I的Schur补。容易看出,当
A
,
B
,
C
A,B,C
A,B,C相当稀疏,而
A
+
B
C
A+BC
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