求解力学问题的有限元方法、CT迭代重建方法等都存在以下问题：
$Ax=b$
其中，$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).

求解线性方程组

线性方程组：$Ax=b$
方阵（解的个数等于方程个数）：$x=A^{-1}b$

方程组超定时（解的个数少于方程个数）：$x=(A^{T}A{)}^{-1}{A}^{T}b$
推导:
优化问题：$x=argminf(x)$
目标函数：$f(x)=||Ax-b|{|}_2^2=(Ax-b{)}^T(Ax-b)=x^T{A}^{T}Ax-2x^T{A}^{T}b-b^{T}b$
梯度：${\partial }_{x}f(x)=2A^T(Ax-b)$
由梯度为0得到：$x=(A^{T}A{)}^{-1}{A}^{T}b$

方程组欠定时（解的个数多于方程个数）：$x=A^T(A{A}^T{)}^{-1}b$
推导:
优化问题：$x=argminf(x)$
目标函数：$f(x)=||x|{|}_2^2+\Lambda (Ax-b)$，$\Lambda =diag(\lambda )$
梯度：${\partial }_{x}f(x)=2x^T+\Lambda A,{\partial }_{\Lambda }f(x)=Ax-b$
由梯度为0得到：$x=A^T(A{A}^T{)}^{-1}b$

方阵$A$满秩，则$A$可逆

行数大于列数的矩形阵$A$列满秩，则$A^{T}A$可逆
推导:
$Ax=0$，若$A$的列满秩，则$Null(A)=0$。
因为$Null(A)=Null(A^{T}A)$，所以若$A$的列满秩，则$A^{T}Ax=0$的唯一解为$0$，得出$A^{T}A$满秩。

行数小于列数的矩形阵$A$行满秩，则$A{A}^T$可逆
推导:
$A^{T}x=0$，若$A^T$的列满秩（$A$的行满秩），则$Null(A^T)=0$。
因为$Null(A^T)=Null(A{A}^T)$，所以若$A$的行满秩，则$A{A}^{T}x=0$的唯一解为$0$，得出$A{A}^T$满秩。

如果$A$不满秩，则$A,A^{T}A,A{A}^T$都不可逆。
此时可用广义逆求解，广义逆可通过SVD分解、QR分解得到，之前采用C++ Eigen库发现SVD分解效率更高：
$A=U\Sigma V^T$，奇异值按由大到小的顺序排列
$A^{+}=V{\Sigma }^{+}{U}^T$
重建图像：$x=A^{+}b$

[1] 线代随笔04-ATA可逆条件. https://bourneli.github.io/linear-algebra/2016/03/03/linear-algebra-04-ATA-inverse.html

可逆、奇异、正定

只有方阵才能说是否可逆

可逆
定义在域$R$上的矩阵$A_{n\times n}$可逆：
$A^{-1}A=A{A}^{-1}=I$，
其中，$I=diag(1)$为单位矩阵（identity matrix，unit matrix）
可逆-性质
$rank(A)=n$
$A^{T}isinvertible$
$Null(A)=0$
$A[:,1],A[:,2],...A[:,n]$构成$R^n$的基
$A[1,:],A[2,:],...A[n,:]$构成$R^n$的基
$det(A)\ne 0$
$A$的特征值不为0

奇异
若$det(A)\ne 0$，则矩阵非奇异；反之，矩阵奇异。
其中，$det(\cdot )$是行列式，是一个将矩阵映射为标量的函数。
非奇异矩阵可逆。

正定
如果定义在域$R$上的实对称矩阵$M$，满足$x^{T}Mx>0,\forall x\in R^n$，则称其正定。
如果定义在域$R$上的实对称矩阵$M$，满足$x^{T}Mx\ge 0,\forall x\in R^n$，则称其半正定。
二次型可以写作$x^{T}Mx$，可以通过判断$M$是否正定来确定二次型是否正定。
正定-性质
正定矩阵是可逆矩阵，其逆矩阵也正定；
正定矩阵$M$缩放为$rM,r>0$，依然正定；
两个正定矩阵$M,N$满足：若$M\ge N$，则$N^{-1}\ge M^{-1}$，$M$的第k个特征值大于$N$的第k个特征值；
正定矩阵+正定矩阵=正定矩阵，正定矩阵+半正定矩阵=正定矩阵，半正定矩阵+半正定矩阵=半正定矩阵；
半正定矩阵的集合是凸的；
实数值函数在一点处取极值的条件：其梯度向量为0，海森矩阵为半正定。
正定矩阵的所有特征值大于0

对任意矩阵$A$：
$rank(A{A}^T)=rank(A)=rank(A^T)=rank(A^{T}A)$

对任意矩阵$A$：
$A^{T}A,A{A}^T$至少半正定
推导
$x^T{A}^{T}Ax=(Ax{)}^T(Ax)$
$x^{T}A{A}^{T}x=(A^{T}x{)}^T(A^{T}x)$
因为向量内积$u^{T}u={\Sigma }_i{u}_i^2\ge 0$，所以$x^T{A}^{T}Ax\ge 0,x^{T}A{A}^{T}x\ge 0$。
又因为$(A^{T}A{)}^T=A^{T}A,(A{A}^T{)}^T=A{A}^T$，所以两个矩阵对称。
综上，两个矩阵至少半正定。

共轭梯度法用于正定矩阵，文献[4]证明了当矩阵$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.

共轭梯度法

线性方程组：$Ax=b$，$A$为正定矩阵 式1

求解式1等价于求解优化问题：
$x^{\ast }=argmi{n}_{x\in R^n}\frac{1}{2}{x}^{T}Ax-x^{T}b$
共轭梯度法就是求解这个优化问题的一种方法。

共轭梯度法
更新公式：${\vec{x}}_{k+1}={\vec{x}}_k+{\alpha }_{k+1}{\vec{p}}_{k+1}$
推出：
${\vec{x}}_k={\vec{x}}_0+{\Sigma }_{j=1}^k{\alpha }_j{p}_j={\vec{x}}_0+P_k\vec{\alpha }$
其中，$P_k=({\vec{p}}_1,...,{\vec{p}}_k)$，$\vec{\alpha }=({\alpha }_1,...,{\alpha }_k)$

$||x_k-x^{\ast }|{|}_A\le 2||x_0-x^{\ast }|{|}_A(\sqrt{k}-1{)}^k/(\sqrt{k}+1{)}^k$
其中，$k$为谱条件数，$||w|{|}_A=\sqrt{w^{T}Aw}$。
因此，条件数越大矩阵越病态，共轭梯度法收敛很慢。
所以一般执行共轭梯度算法之前，要先进行预处理，使新的求解问题的条件数较小。

预处理
将$A$拆分：
$A=M-N$
其中，$M$为正定矩阵，Cholesky分解$M=R^{T}R$，$R$为上三角矩阵。

预处理后得到新问题：
$\hat{A}x=\hat{b}$
其中，$\hat{A}=R^{-T}A{R}^{-1},\hat{b}=R^{-T}b$，$y=Rx$
我们希望$k(\hat{A})\ll k(A)$

预条件矩阵举例：
$M=diag(a_{11},...,a_{nn})$
$\hat{A}=D^{-1/2}A{D}^{-1/2},\hat{b}=D^{-1/2}b$

[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