CT重建笔记(三)——共轭梯度法
求解力学问题的有限元方法、CT迭代重建方法等都存在以下问题:
A
x
=
b
Ax=b
Ax=b
其中,
A
A
A为大型稀疏矩阵。利用系统矩阵的稀疏性可以把原问题缩减为一个小问题,从而加速求解。对于同样维度的问题,采用缩减法对于线性方程组问题不如特征值问题的效率提高明显。
我的毕业论文就是研究有限元模型的缩减,里面研究了一下这些方法。
关于缩减线性方程组问题、特征值问题的文献推荐,尤其是这名柏兆俊
老师,我觉得他太厉害了。
书籍:
[1] Saad Y. Iterative methods for sparse linear systems[M]. Philadelphia: Society for
Industrial and Applied Mathematics, 2003.
[2] Bai Z, Demmel J, Dongarra J, et al. Templates for the solution of algebraic eigenvalue
problems: a practical guide[M]. Philadelphia: Society for Industrial and Applied
Mathematics, 2000.
[3] Gutknecht M H. A Brief Introduction to Krylov Space Methods for Solving Linear
Systems[M]. Berlin: Springer, 2007.
[4] Barrett R, Berry M, Chan T F, et al. Templates for the solution of linear systems:
building blocks for iterative methods[M]. Philadelphia: Society for Industrial and Applied
Mathematics, 1994.
经典论文:
[5] Lanczos C. An Iteration Method for the solution of the Eigenvalue Problem of Linear
Differential and Integral Operators[J]. Journal of research of the National Bureau of
Standards, 1950(4).
[6] Arnoldi W E. The principle of minimized iterations in the solution of the matrix
eigenvalue problem[J]. Quarterly of Applied Mathematics, 1951(1).
求解线性方程组
线性方程组:
A
x
=
b
Ax=b
Ax=b
方阵(解的个数等于方程个数):
x
=
A
−
1
b
x=A^{-1}b
x=A−1b
方程组超定时(解的个数少于方程个数):
x
=
(
A
T
A
)
−
1
A
T
b
x=(A^TA)^{-1} A^Tb
x=(ATA)−1ATb
推导:
优化问题:
x
=
a
r
g
m
i
n
f
(
x
)
x = argmin \space f(x)
x=argmin f(x)
目标函数:
f
(
x
)
=
∣
∣
A
x
−
b
∣
∣
2
2
=
(
A
x
−
b
)
T
(
A
x
−
b
)
=
x
T
A
T
A
x
−
2
x
T
A
T
b
−
b
T
b
f(x) = ||Ax-b||_2^2 = (Ax-b)^T (Ax-b) = x^TA^TAx - 2x^TA^Tb - b^Tb
f(x)=∣∣Ax−b∣∣22=(Ax−b)T(Ax−b)=xTATAx−2xTATb−bTb
梯度:
∂
x
f
(
x
)
=
2
A
T
(
A
x
−
b
)
\partial_x f(x) = 2A^T(A x- b)
∂xf(x)=2AT(Ax−b)
由梯度为0得到:
x
=
(
A
T
A
)
−
1
A
T
b
x=(A^TA)^{-1} A^Tb
x=(ATA)−1ATb
方程组欠定时(解的个数多于方程个数):
x
=
A
T
(
A
A
T
)
−
1
b
x=A^T(AA^T)^{-1}b
x=AT(AAT)−1b
推导:
优化问题:
x
=
a
r
g
m
i
n
f
(
x
)
x = argmin \space f(x)
x=argmin f(x)
目标函数:
f
(
x
)
=
∣
∣
x
∣
∣
2
2
+
Λ
(
A
x
−
b
)
f(x) = ||x||_2^2 + \Lambda (Ax-b)
f(x)=∣∣x∣∣22+Λ(Ax−b),
Λ
=
d
i
a
g
(
λ
)
\Lambda=diag(\lambda)
Λ=diag(λ)
梯度:
∂
x
f
(
x
)
=
2
x
T
+
Λ
A
,
∂
Λ
f
(
x
)
=
A
x
−
b
\partial_x f(x) = 2x^T+\Lambda A, \partial_{\Lambda} f(x) = Ax-b
∂xf(x)=2xT+ΛA,∂Λf(x)=Ax−b
由梯度为0得到:
x
=
A
T
(
A
A
T
)
−
1
b
x=A^T(AA^T)^{-1}b
x=AT(AAT)−1b
方阵 A A A满秩,则 A A A可逆
行数大于列数的矩形阵
A
A
A列满秩,则
A
T
A
A^TA
ATA可逆
推导:
A
x
=
0
Ax=0
Ax=0,若
A
A
A的列满秩,则
N
u
l
l
(
A
)
=
0
Null(A)=0
Null(A)=0。
因为
N
u
l
l
(
A
)
=
N
u
l
l
(
A
T
A
)
Null(A)=Null(A^TA)
Null(A)=Null(ATA),所以若
A
A
A的列满秩,则
A
T
A
x
=
0
A^TAx=0
ATAx=0的唯一解为
0
0
0,得出
A
T
A
A^TA
ATA满秩。
行数小于列数的矩形阵
A
A
A行满秩,则
A
A
T
AA^T
AAT可逆
推导:
A
T
x
=
0
A^Tx = 0
ATx=0,若
A
T
A^T
AT的列满秩(
A
A
A的行满秩),则
N
u
l
l
(
A
T
)
=
0
Null(A^T)=0
Null(AT)=0。
因为
N
u
l
l
(
A
T
)
=
N
u
l
l
(
A
A
T
)
Null(A^T)=Null(AA^T)
Null(AT)=Null(AAT),所以若
A
A
A的行满秩,则
A
A
T
x
=
0
AA^Tx=0
AATx=0的唯一解为
0
0
0,得出
A
A
T
AA^T
AAT满秩。
如果
A
A
A不满秩,则
A
,
A
T
A
,
A
A
T
A,A^TA,AA^T
A,ATA,AAT都不可逆。
此时可用广义逆求解,广义逆可通过SVD分解、QR分解得到,之前采用C++ Eigen库发现SVD分解效率更高:
A
=
U
Σ
V
T
A=U\Sigma V^T
A=UΣVT,奇异值按由大到小的顺序排列
A
+
=
V
Σ
+
U
T
A^+ = V\Sigma^+ U^T
A+=VΣ+UT
重建图像:
x
=
A
+
b
x = A^+b
x=A+b
[1] 线代随笔04-ATA可逆条件. https://bourneli.github.io/linear-algebra/2016/03/03/linear-algebra-04-ATA-inverse.html
可逆、奇异、正定
只有方阵才能说是否可逆
可逆
定义在域
R
R
R上的矩阵
A
n
×
n
A_{n\times n}
An×n可逆:
A
−
1
A
=
A
A
−
1
=
I
A^{-1}A = AA^{-1} = I
A−1A=AA−1=I,
其中,
I
=
d
i
a
g
(
1
)
I=diag(1)
I=diag(1)为单位矩阵(identity matrix,unit matrix)
可逆-性质
r
a
n
k
(
A
)
=
n
rank(A)=n
rank(A)=n
A
T
i
s
i
n
v
e
r
t
i
b
l
e
A^T \space is \space invertible
AT is invertible
N
u
l
l
(
A
)
=
0
Null(A)=0
Null(A)=0
A
[
:
,
1
]
,
A
[
:
,
2
]
,
.
.
.
A
[
:
,
n
]
A[:,1],A[:,2],...A[:,n]
A[:,1],A[:,2],...A[:,n]构成
R
n
R^n
Rn的基
A
[
1
,
:
]
,
A
[
2
,
:
]
,
.
.
.
A
[
n
,
:
]
A[1,:],A[2,:],...A[n,:]
A[1,:],A[2,:],...A[n,:]构成
R
n
R^n
Rn的基
d
e
t
(
A
)
≠
0
det(A)\neq 0
det(A)=0
A
A
A的特征值不为0
奇异
若
d
e
t
(
A
)
≠
0
det(A)\neq 0
det(A)=0,则矩阵非奇异;反之,矩阵奇异。
其中,
d
e
t
(
⋅
)
det(\cdot)
det(⋅)是行列式,是一个将矩阵映射为标量的函数。
非奇异矩阵可逆。
正定
如果定义在域
R
R
R上的实对称矩阵
M
M
M,满足
x
T
M
x
>
0
,
∀
x
∈
R
n
x^TMx>0, \forall x \in R^n
xTMx>0,∀x∈Rn,则称其正定。
如果定义在域
R
R
R上的实对称矩阵
M
M
M,满足
x
T
M
x
≥
0
,
∀
x
∈
R
n
x^TMx\geq 0, \forall x \in R^n
xTMx≥0,∀x∈Rn,则称其半正定。
二次型可以写作
x
T
M
x
x^TMx
xTMx,可以通过判断
M
M
M是否正定来确定二次型是否正定。
正定-性质
正定矩阵是可逆矩阵,其逆矩阵也正定;
正定矩阵
M
M
M缩放为
r
M
,
r
>
0
rM, r>0
rM,r>0,依然正定;
两个正定矩阵
M
,
N
M,N
M,N满足:若
M
≥
N
M \geq N
M≥N,则
N
−
1
≥
M
−
1
N^{-1} \geq M^{-1}
N−1≥M−1,
M
M
M的第k个特征值大于
N
N
N的第k个特征值;
正定矩阵+正定矩阵=正定矩阵,正定矩阵+半正定矩阵=正定矩阵,半正定矩阵+半正定矩阵=半正定矩阵;
半正定矩阵的集合是凸的;
实数值函数在一点处取极值的条件:其梯度向量为0,海森矩阵为半正定。
正定矩阵的所有特征值大于0
对任意矩阵
A
A
A:
r
a
n
k
(
A
A
T
)
=
r
a
n
k
(
A
)
=
r
a
n
k
(
A
T
)
=
r
a
n
k
(
A
T
A
)
rank(AA^T)=rank(A)=rank(A^T)=rank(A^TA)
rank(AAT)=rank(A)=rank(AT)=rank(ATA)
对任意矩阵
A
A
A:
A
T
A
,
A
A
T
A^TA, AA^T
ATA,AAT至少半正定
推导
x
T
A
T
A
x
=
(
A
x
)
T
(
A
x
)
x^TA^TAx = (Ax)^T(Ax)
xTATAx=(Ax)T(Ax)
x
T
A
A
T
x
=
(
A
T
x
)
T
(
A
T
x
)
x^T AA^T x = (A^Tx)^T(A^Tx)
xTAATx=(ATx)T(ATx)
因为向量内积
u
T
u
=
Σ
i
u
i
2
≥
0
u^Tu = \Sigma_i u_i^2 \geq 0
uTu=Σiui2≥0,所以
x
T
A
T
A
x
≥
0
,
x
T
A
A
T
x
≥
0
x^TA^TAx\geq 0, x^TAA^Tx\geq 0
xTATAx≥0,xTAATx≥0。
又因为
(
A
T
A
)
T
=
A
T
A
,
(
A
A
T
)
T
=
A
A
T
(A^TA)^T=A^TA, (AA^T)^T=AA^T
(ATA)T=ATA,(AAT)T=AAT,所以两个矩阵对称。
综上,两个矩阵至少半正定。
共轭梯度法用于正定矩阵,文献[4]证明了当矩阵 A A A半正定且系统相容(至少有一个解),则可以用共轭梯度法求解。
[1] 正定矩阵. https://en.wikipedia.org/wiki/Definite_matrix
[2] ATA对称正定的证明. https://blog.csdn.net/wu_nan_nan/article/details/75097015
[3] What happens when I use a conjugate gradient solver with a symmetric positive semi-definite matrix? https://scicomp.stackexchange.com/questions/35239/what-happens-when-i-use-a-conjugate-gradient-solver-with-a-symmetric-positive-se
[4] Convergence of the Conjugate Gradient Method on Singular Systems. https://arxiv.org/abs/1809.00793.
共轭梯度法
线性方程组:
A
x
=
b
Ax=b
Ax=b,
A
A
A为正定矩阵 式1
求解式1
等价于求解优化问题:
x
∗
=
a
r
g
m
i
n
x
∈
R
n
1
2
x
T
A
x
−
x
T
b
x^* = argmin_{x\in R^n} \frac{1}{2}x^TAx - x^T b
x∗=argminx∈Rn21xTAx−xTb
共轭梯度法就是求解这个优化问题的一种方法。
共轭梯度法
更新公式:
x
⃗
k
+
1
=
x
⃗
k
+
α
k
+
1
p
⃗
k
+
1
\vec x_{k+1} = \vec x_k+\alpha_{k+1} \vec p_{k+1}
xk+1=xk+αk+1pk+1
推出:
x
⃗
k
=
x
⃗
0
+
Σ
j
=
1
k
α
j
p
j
=
x
⃗
0
+
P
k
α
⃗
\vec x_k = \vec x_0 + \Sigma_{j=1}^k \alpha_j p_j = \vec x_0 + P_k \vec \alpha
xk=x0+Σj=1kαjpj=x0+Pkα
其中,
P
k
=
(
p
⃗
1
,
.
.
.
,
p
⃗
k
)
P_k = (\vec p_1, ..., \vec p_k)
Pk=(p1,...,pk),
α
⃗
=
(
α
1
,
.
.
.
,
α
k
)
\vec \alpha = (\alpha_1, ..., \alpha_k)
α=(α1,...,αk)
∣
∣
x
k
−
x
∗
∣
∣
A
≤
2
∣
∣
x
0
−
x
∗
∣
∣
A
(
k
−
1
)
k
/
(
k
+
1
)
k
||x_k - x^*||_A \leq 2||x_0 - x^*||_A (\sqrt k - 1)^k/(\sqrt k + 1)^k
∣∣xk−x∗∣∣A≤2∣∣x0−x∗∣∣A(k−1)k/(k+1)k
其中,
k
k
k为谱条件数,
∣
∣
w
∣
∣
A
=
w
T
A
w
||w||_A=\sqrt{w^TAw}
∣∣w∣∣A=wTAw。
因此,条件数越大矩阵越病态,共轭梯度法收敛很慢。
所以一般执行共轭梯度算法之前,要先进行预处理,使新的求解问题的条件数较小。
预处理
将
A
A
A拆分:
A
=
M
−
N
A = M-N
A=M−N
其中,
M
M
M为正定矩阵,Cholesky分解
M
=
R
T
R
M=R^TR
M=RTR,
R
R
R为上三角矩阵。
预处理后得到新问题:
A
^
x
=
b
^
\hat{A} x = \hat{b}
A^x=b^
其中,
A
^
=
R
−
T
A
R
−
1
,
b
^
=
R
−
T
b
\hat{A} = R^{-T}AR^{-1}, \hat{b} = R^{-T} b
A^=R−TAR−1,b^=R−Tb,
y
=
R
x
y = Rx
y=Rx
我们希望
k
(
A
^
)
≪
k
(
A
)
k(\hat{A}) \ll k(A)
k(A^)≪k(A)
预条件矩阵举例:
M
=
d
i
a
g
(
a
11
,
.
.
.
,
a
n
n
)
M = diag(a_{11},...,a_{nn})
M=diag(a11,...,ann)
A
^
=
D
−
1
/
2
A
D
−
1
/
2
,
b
^
=
D
−
1
/
2
b
\hat{A} = D^{-1/2}AD^{-1/2}, \hat{b} = D^{-1/2}b
A^=D−1/2AD−1/2,b^=D−1/2b
[1] 共轭梯度法. https://www.cse.psu.edu/~b58/cse456/lecture20.pdf
[2] 矩阵条件数. https://leslielee.blog.csdn.net/article/details/120564282
[3] 预条件. https://math.ecnu.edu.cn/~jypan/Teaching/IMP/slides_IMP05_beamer.pdf
原文地址:https://blog.csdn.net/qq_37083038/article/details/145166579
免责声明:本站文章内容转载自网络资源,如侵犯了原著者的合法权益,可联系本站删除。更多内容请关注自学内容网(zxcms.com)!