- 前言
- LU直接分解法 — Doolittle_Decomposition
- 代码实现
- test example
- 三对角矩阵的特殊解法
- 代码实现
- test example
- 迭代法解线性方程组
- Jacobi迭代法
- Gauss-Seidel迭代法
首先很多内容都跟矩阵论的部分重合了，相关的关于LU分解、LDV分解、LU解线性方程组、求行列式、求逆等都可以在矩阵论专栏中对应查看，这里仅补充一些其他方法。
LU直接分解法 — Doolittle_Decomposition
- 道立特分解属于一种先设定好U矩阵第一行和L矩阵第一列的值，便可以开始迭代求解L和U的其他值的算法，分解结果如下（以三阶矩阵为例）：
[ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ] = [ 1 0 0 l 21 1 0 l 31 l 32 1 ] [ u 11 u 12 u 13 0 u 22 u 23 0 0 u 33 ] \begin{bmatrix} a_{11} & a_{12} & a_{13} \ a_{21} & a_{22} & a_{23} \ a_{31} & a_{32} & a_{33} \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0 \ l_{21} & 1& 0 \ l_{31} & l_{32} & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \ 0 & u_{22} & u_{23} \ 0 & 0 & u_{33} \end{bmatrix} ⎣⎡a11a21a31a12a22a32a13a23a33⎦⎤=⎣⎡1l21l3101l32001⎦⎤⎣⎡u1100u12u220u13u23u33⎦⎤ - 算法流程（先求解U的对应行再求解L的对应列）：
- 求解U的第i行第j个元素： u i j = a i j − ∑ k = 1 i − 1 l i k ∗ u k j u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}*u_{kj} uij=aij−∑k=1i−1lik∗ukj
- 求解L的第j列第i个元素： l j i = ( a j i − ∑ k = 1 i − 1 l j k ∗ u k i ) / u i i l_{ji}=(a_{ji}-\sum_{k=1}^{i-1}l_{jk}*u_{ki})/u_{ii} lji=(aji−∑k=1i−1ljk∗uki)/uii
@staticmethod
def domain_row_transform(arr, eps=1e-6, Test=False):
"""domain_row_transform 选取主元并进行行变换
Args:
arr ([np.darray]): [input matrix]
eps ([float]): numerical threshold.
test (bool, optional): [show checking information]. Defaults to False.
Returns:
[A]: [resort one of arr]
Author: Junno
Date: 2022-04-27
"""
assert len(arr.shape) == 2
A = np.copy(arr)
M, N = A.shape
for i in range(M):
for j in range(i, N):
if Test:
print('During process in row:{}, col:{}'.format(i, j))
if sum(abs(A[i:, j])) > eps:
zero_index = []
non_zero_index = []
for k in range(i, M):
if abs(A[k, j]) > eps:
non_zero_index.append(k)
else:
zero_index.append(k)
non_zero_index = sorted(
non_zero_index, key=lambda x: abs(A[x, j]))
sort_idnex = non_zero_index+zero_index
if Test:
print('Sorting index for {} cols:\n'.format(i), sort_idnex)
if sort_idnex[0] != i:
A[[sort_idnex[0], i], :] = A[[i, sort_idnex[0]], :]
if Test:
print('After resorting\n', A)
return A
@classmethod
def Doolittle_Decomposition(self, target, domain_selection=True, test=False):
"""Doolittle_Decomposition equal to LU_Decomposition
Args:
target ([np.darray]): [input matrix]
domain_selection (bool, optional): [operate domain selection on target]. Defaults to True.
test (bool, optional): [show testing information]. Defaults to False.
Returns:
[resort input, L, U] if domain_selection=True else [L,U]
Author: Junno
Date: 2022-04-27
"""
assert len(target.shape) == 2
arr = np.copy(target)
M, N = arr.shape
L = np.zeros_like(arr, dtype=arr.dtype)
U = np.zeros_like(arr, dtype=arr.dtype)
if domain_selection:
arr = self.domain_row_transform(arr)
for i in range(M):
if test:
print(i)
print(L, U, sep='\n')
for j in range(i, N):
if test:
print('U', i, j, arr[i, j], L[i, :i], U[:i, j])
U[i, j] = arr[i, j]-L[i, :i]@U[:i, j]
if i == j:
L[i, j] = 1
else:
if test:
print('L', j, i, arr[j, i],
L[j, :i], U[:i, i], U[i, i])
L[j, i] = (arr[j, i]-(L[j, :i]@U[:i, i]))/U[i, i]
if domain_selection:
print('return resort input matrix, L, U')
return arr, L, U
else:
return L, U
- Tips：为避免大数吃小数的问题和舍入误差问题，加入了选主元的选项，可以对原输入矩阵进行行变换，不改变结果的正确性，可以自行选择
A=np.random.randn(5,5)
A
>>>
array([[-1.50948856, -1.61033161, -0.20671825, 0.4286645 , -0.2572608 ],
[-1.28615783, 0.87991375, 1.71080037, -1.61832726, 1.52242031],
[-0.98092031, -0.32273876, 2.20229521, 0.853808 , -2.15168542],
[-2.02003405, 0.05052931, -1.01116616, 0.26437957, -1.65870163],
[-0.83327297, 0.32495716, -1.02878403, -0.96607898, -0.59284746]])
L,U=Doolittle_Decomposition(A,domain_selection=False,test=False)
L
>>>
[[ 1. 0. 0. 0. 0. ]
[ 4.688 1. 0. 0. 0. ]
[ 5.484 -3863.106 1. 0. 0. ]
[ 8.449 -8279.836 2.143 1. 0. ]
[ -8.763 10132.862 -2.623 -11.842 1. ]]
U
>>>
[[ -0.165 0.618 -0.885 -0.515 -1.011]
[ 0. 0.001 3.434 2.532 6.176]
[ 0. 0. 13269.089 9782.795 23865.142]
[ 0. 0. 0. 0.216 -4.26 ]
[ 0. 0. 0. 0. -43.798]]
print(L@U)
>>>
array([[-1.50948856, -1.61033161, -0.20671825, 0.4286645 , -0.2572608 ],
[-1.28615783, 0.87991375, 1.71080037, -1.61832726, 1.52242031],
[-0.98092031, -0.32273876, 2.20229521, 0.853808 , -2.15168542],
[-2.02003405, 0.05052931, -1.01116616, 0.26437957, -1.65870163],
[-0.83327297, 0.32495716, -1.02878403, -0.96607898, -0.59284746]])
-
对于一般的线性方程组，做高斯消去法的乘法复杂度大概为 O ( N 3 ) O(N^3) O(N3)， 当N超过100时就比较吃力了。在很多物理仿真的求解计算中，最常见的问题便是求解线性方程组，但是它们的维度一般都成千上万。但同时也存在一个比较有用的特性，即稀疏特性，现在有很多先进的算法来求解这些问题。这里我们介绍一种特殊情况下，即对角占优下三对角线性方程组的快速解法——Thomas分解或者追赶法，它的复杂度可以与N同阶，为 O ( 3 N − 2 ) O(3N-2) O(3N−2)。
-
对角占优的三对角方程组形如下：
[ b 1 c 1 0 ⋯ 0 a 2 b 2 c 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 a n − 1 b n − 1 c n − 1 0 0 0 a n b n ] [ x 1 x 2 ⋮ x n ] = [ f 1 f 2 ⋮ f n ] \begin{bmatrix} b_{1} & c_{1} & 0 & \cdots & 0 \ a_{2} & b_{2} & c_{2} & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & a_{n-1} & b_{n-1} & c_{n-1}\ 0& 0& 0& a_{n} & b_{n} \end{bmatrix} \begin{bmatrix} x_1\ x_2\ \vdots\ x_{n} \end{bmatrix} = \begin{bmatrix} f_1\ f_2\ \vdots\ f_{n} \end{bmatrix} ⎣⎢⎢⎢⎢⎡b1a2000c1b2⋱000c2⋱an−10⋯⋯⋱bn−1an000cn−1bn⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡f1f2⋮fn⎦⎥⎥⎥⎤
-
对角占优条件：
- ∣ b 1 ∣ > ∣ c 1 ∣ > 0 |b_1| > |c_1| >0 ∣b1∣>∣c1∣>0
- ∣ b i ∣ ≥ ∣ a i ∣ + ∣ c i ∣ , a i , c i ≠ 0 , ( i = 2 , 3 , . . . , n − 1 ) |b_i| \geq |a_i| +|c_i|, a_i,c_i \neq0,(i=2,3,...,n-1) ∣bi∣≥∣ai∣+∣ci∣,ai,ci=0,(i=2,3,...,n−1)
- ∣ b n ∣ > ∣ a n ∣ > 0 |b_n| > |a_n|>0 ∣bn∣>∣an∣>0
-
类似于LU分解解线性方程组的做法，尝试将三对角矩阵分解为上下三角矩阵，形如：
A = [ b 1 c 1 0 ⋯ 0 a 2 b 2 c 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 a n − 1 b n − 1 c n − 1 0 0 0 a n b n ] = L U = [ α 1 0 0 ⋯ 0 γ 2 α 2 0 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 γ n − 1 α n − 1 0 0 0 0 γ n α n ] [ 1 β 1 0 ⋯ 0 0 1 β 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 0 1 β n − 1 0 0 0 0 1 ] A=\begin{bmatrix} b_{1} & c_{1} & 0 & \cdots & 0 \ a_{2} & b_{2} & c_{2} & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & a_{n-1} & b_{n-1} & c_{n-1}\ 0& 0& 0& a_{n} & b_{n} \end{bmatrix}\ =LU=\begin{bmatrix} \alpha_{1} & 0 & 0 & \cdots & 0 \ \gamma_{2} & \alpha_{2} & 0 & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & \gamma_{n-1} & \alpha_{n-1} & 0\ 0& 0& 0& \gamma_{n} & \alpha_{n} \end{bmatrix} \begin{bmatrix} 1 & \beta_{1} & 0 & \cdots & 0 \ 0 & 1 & \beta_{2} & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & 0 &1 & \beta_{n-1}\ 0& 0& 0& 0 & 1 \end{bmatrix} A=⎣⎢⎢⎢⎢⎡b1a2000c1b2⋱000c2⋱an−10⋯⋯⋱bn−1an000cn−1bn⎦⎥⎥⎥⎥⎤=LU=⎣⎢⎢⎢⎢⎡α1γ20000α2⋱0000⋱γn−10⋯⋯⋱αn−1γn0000αn⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡10000β11⋱000β2⋱00⋯⋯⋱10000βn−11⎦⎥⎥⎥⎥⎤ -
比较待定系数，可以逐步求解出各系数，（推导过程详见教材第七章P184）此时 A x = f Ax=f Ax=f等价于两个三角方程组 L y = f Ly=f Ly=f和 U x = y Ux=y Ux=y，得到追赶法求解过程：
- 计算 β i \beta_i βi： β 1 = c 1 / b 1 , β i = c i / ( b i − a i β i − 1 ) , i = ( 2 , 3 , . . . , n − 1 ) \beta_1=c_1/b_1,\beta_i=c_i/(b_i-a_i \beta_{i-1}), i=(2,3,...,n-1) β1=c1/b1,βi=ci/(bi−aiβi−1),i=(2,3,...,n−1)
- 解 L y = f Ly=f Ly=f： y 1 = f 1 / b 1 , y i = ( f i − a i y i − 1 ) / ( b i − a i β i − 1 ) , i = ( 2 , 3 , . . . , n ) y_1=f_1/b_1,y_i=(f_i-a_i y_{i-1})/(b_i-a_i \beta_{i-1}), i=(2,3,...,n) y1=f1/b1,yi=(fi−aiyi−1)/(bi−aiβi−1),i=(2,3,...,n)
- 解 U x = y Ux=y Ux=y： x n = y n , x i = y i − β i x i + 1 , ( i = n − 1 , n − 2 , . . . , 2 , 1 ) x_n=y_n,x_i=y_i-\beta_i x_{i+1}, (i=n-1,n-2,...,2,1) xn=yn,xi=yi−βixi+1,(i=n−1,n−2,...,2,1)
- 实际计算时只需要传入两个次对角线和主对角线的三个向量值便可以求解，内存消耗更小！
def check_SDD_vec(a, b, c, eps=1e-4, test=False):
"""check_SDD_vec 向量形式的严格对角占优条件判断
Args:
a ([np.darray]): [lower-secondary diagonal]
b ([np.darray]): [diagonal]
c ([np.darray]): [upper-secondary diagonal]
eps ([type], optional): [threshold]. Defaults to 1e-4.
test (bool, optional): [show testing information]. Defaults to False.
Returns:
[bool]: [whether or not]
Author: Junno
Date: 2022-04-28
"""
assert len(a) == len(c) and len(a)+1 == len(b)
N = len(b)
a_ = np.concatenate(([eps], a), axis=0)
c_ = np.concatenate((c, [eps]), axis=0)
if test:
print(a_, c_)
sum_a_c = np.abs(a_)+np.abs(c_)
b_ = np.abs(b)
for i in range(N):
if b_[i] < sum_a_c[i]:
return False
return True
def check_SDD_M(arr,eps=1e-4):
"""check_SDDM 判断输入矩阵是否满足严格对角占优条件
Args:
arr ([np.darray]): [input matrix]
eps ([type], optional): [threshold]. Defaults to 1e-4.
Returns:
[bool]: [whether or not]
Author: Junno
Date: 2022-04-28
"""
assert arr.shape[0]==arr.shape[1]
return check_SDD_vec(np.diagonal(arr,offset=-1),np.diagonal(arr,offset=0),np.diagonal(arr,offset=1),eps=eps)
def Thomas_Linear_Equation_Solve(a, b, c, d, eps=1e-4, test=False):
"""Thomas_Linear_Equation_Solve
Args:
a ([list or np.darray]): [lower-secondary diagonal]
b ([list or np.darray]): [diagonal]
c ([list or np.darray]): [upper-secondary diagonal]
d ([list or np.darray]): [b of Ax=b]
eps ([type], optional): [threshold]. Defaults to 1e-4.
test (bool, optional): [show testing information]. Defaults to False.
Returns:
x [list]: [solution of Ax=b]
Author: Junno
Date: 2022-04-28
"""
assert len(a) == len(c) and len(a)+1 == len(b)
if not check_SDD_vec(a, b, c, eps=eps):
print('The answer may be wrong because the input are not meet SDD condition completely')
N = len(b)
if test:
print(a, b, c, d)
beta = np.zeros_like(b)
x = np.zeros_like(b)
y = np.zeros_like(b)
beta[0] = c[0]/b[0]
y[0] = d[0]/b[0]
for i in range(1, N):
if i < N-1:
beta[i] = c[i]/(b[i]-a[i-1]*beta[i-1])
y[i] = (d[i]-a[i-1]*y[i-1])/(b[i]-a[i-1]*beta[i-1])
for i in range(N-1, -1, -1):
if i == N-1:
x[i] = y[i]
else:
x[i] = y[i]-beta[i]*x[i+1]
return x
A=np.array([[3,1,0,0],[1,2,1,0],[0,1,3,2],[0,0,3,4]],dtype=np.float32).reshape((4,4))
A
>>>
array([[3., 1., 0., 0.],
[1., 2., 1., 0.],
[0., 1., 3., 2.],
[0., 0., 3., 4.]], dtype=float32)
# create testing example
x=np.random.randn(4,1)
print(x)
>>>
[[-0.35812746]
[ 0.44621921]
[-0.45669147]
[ 1.85002579]]
f=A@x
print(f)
>>>
[[-0.62816316]
[ 0.07761949]
[ 2.77619638]
[ 6.03002876]]
Thomas_Linear_Equation_Solve(np.diagonal(A,offset=-1),np.diagonal(A,offset=0),np.diagonal(A,offset=1),f,test=False)
>>>
array([-0.35812744, 0.44621915, -0.45669138, 1.8500258 ], dtype=float32)
迭代法解线性方程组
-
对于大型的稀疏矩阵，有时候迭代求解方法在效率和成本上较直接求解更优。对于方程组 A x = b Ax=b Ax=b，可以通过变形得到 x = B x + f x=B x+f x=Bx+f，再通过迭代方法，逐步逼近真实解， x k + 1 = B x k + f x^{k+1}=B x^k +f xk+1=Bxk+f。其中，B称为迭代传递矩阵。
-
如果 lim k → ∞ x ( k ) \lim_{k \rightarrow \infty}x^{(k)} limk→∞x(k)存在，称迭代过程收敛。引进误差向量， ϵ k + 1 = x k + 1 − x ∗ \epsilon^{k+1}=x^{k+1}-x^* ϵk+1=xk+1−x∗，递推可以得到， ϵ k = B x k − 1 = B k ϵ 0 \epsilon^{k}=B x^{k-1}=B^k \epsilon^0 ϵk=Bxk−1=Bkϵ0。即若想要迭代收敛，需要知道B满足什么情况时有 B k → 0 ( k → ∞ ) B^k \rightarrow 0 \,\, (k \rightarrow \infty) Bk→0(k→∞)。
-
迭代法基本定理：当 ρ ( B ) < 1 \rho(B) <1 ρ(B)<1时，上述迭代法在任意初始向量 x 0 x_0 x0及任意f开始迭代后均收敛。其中， ρ ( B ) \rho(B) ρ(B)表示B的谱半径。可以由Jordan分解简单证明得到，可以参考教材P206-207。
-
由Jordan的证明过程还可以得到，当 ρ ( B ) \rho(B) ρ(B)越小时，迭代收敛的速度越快，若给定求解精度为 1 0 − s 10^{-s} 10−s，则至少需要的迭代次数可以由以下式子确定：
{ [ ρ ( B ) ] k ≤ 1 0 − s k ≥ s ln 10 − ln ρ ( B ) \begin{cases} [\rho(B)]^k \leq 10^{-s}\ k \geq \frac{s \ln{10}}{-\ln{\rho(B)}}\ \end{cases} {[ρ(B)]k≤10−sk≥−lnρ(B)sln10 -
由于任意阶矩阵范数都是谱半径的上届，所以可以有一些容易计算的矩阵范数得到收敛性的大概估计。设某一矩阵范数表示为 ∣ ∣ B ∣ ∣ v = q < 1 ||B||_v=q <1 ∣∣B∣∣v=q<1，第k步的误差估计可以表示为：
∣ ∣ x ∗ − x k ∣ ∣ v ≤ q 1 − q ∣ ∣ x k − x k − 1 ∣ ∣ v ||x^*-x^k||_v \leq \frac{q}{1-q} ||x^k-x^{k-1}||_v ∣∣x∗−xk∣∣v≤1−qq∣∣xk−xk−1∣∣v
- 设有方程组
∑
j
=
1
n
a
i
j
x
j
=
b
i
\sum_{j=1}^n a_{ij}x_j=b_i
∑j=1naijxj=bi，记为
A
x
=
b
Ax=b
Ax=b，将A分解为
A
=
D
+
L
+
U
A=D+L+U
A=D+L+U，格式如下：
D = [ a 11 a 22 ⋱ a n n ] D=\begin{bmatrix} a_{11} & & & \ & a_{22} & & \ & &\ddots& \ & & &a_{nn} \ \end{bmatrix} D=⎣⎢⎢⎡a11a22⋱ann⎦⎥⎥⎤
L = [ 0 a 21 0 ⋮ ⋱ ⋱ a n 1 ⋯ a n , n − 1 0 ] L=\begin{bmatrix} 0 & & & \ a_{21} & 0 & & \ \vdots & \ddots &\ddots& \ a_{n1} & \cdots & a_{n,n-1} &0 \ \end{bmatrix} L=⎣⎢⎢⎢⎡0a21⋮an10⋱⋯⋱an,n−10⎦⎥⎥⎥⎤
U = [ 0 a 12 ⋯ a 1 n 0 ⋱ ⋮ ⋱ a n − 1 , n 0 ] U=\begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \ & 0& \ddots & \vdots \ & &\ddots&a_{n-1,n} \ & & &0 \ \end{bmatrix} U=⎣⎢⎢⎢⎡0a120⋯⋱⋱a1n⋮an−1,n0⎦⎥⎥⎥⎤
- 对于方程组
∑
j
=
1
n
a
i
j
x
j
=
b
i
\sum_{j=1}^n a_{ij}x_j=b_i
∑j=1naijxj=bi，可以改写为
x
i
=
1
a
i
i
(
b
i
−
∑
j
=
1
,
j
≠
i
n
a
i
j
x
j
)
,
(
i
=
1
,
2
,
.
.
.
,
n
)
x_i=\frac{1}{a_{ii}}(b_i-\sum_{j=1,j \neq i}^n a_{ij}x_j),(i=1,2,...,n)
xi=aii1(bi−j=1,j=i∑naijxj),(i=1,2,...,n)若是以
x
=
B
x
+
f
x=Bx+f
x=Bx+f，则可以对照写出矩阵形式：
{ B = − D − 1 ( L + U ) f = D − 1 b \begin{cases} B=-D^{-1}(L+U)\ f=D^{-1}b\ \end{cases} {B=−D−1(L+U)f=D−1b
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）