超级会员免费看
快速线性求解器：Cholesky分解、QR分解与共轭梯度法
1. Cholesky分解
对于对称正定矩阵 $A$，可以进行Cholesky分解，形式为 $A = L L^T$。这里，由于矩阵的对称性，矩阵 $U$ 等于 $L$ 的转置。若要求 $L$ 的所有对角元素为正，则这种分解是唯一的。需要注意的是，此时的 $L$ 与LU分解中得到的矩阵不同，在LU分解中，所有对角元素都等于1。
我们可以通过以下方式显式地求出 $L$ 的元素 $\ell_{ij}$：
[
\begin{bmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \
\vdots & \vdots & \ddots & \vdots \
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{bmatrix}
=
\begin{bmatrix}
\ell_{11} & 0 & \cdots & 0 \
\ell_{21} & \ell_{22} & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots \
\ell_{n1} & \ell_{n2} & \cdots & \ell_{nn}
\end{bmatrix}
\begin{bmatrix}
\ell_{11} & \ell_{21} & \cdots & \ell_{n1} \
0 & \ell_{22} & \cdots & \ell_{n2} \
\vdots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & \ell_{nn}
\end{bmatrix}
]
通过使等式两边元素相等，对于第一列（$i = 1, \cdots, n$），我们可以得到 $a_{i1} = \ell_{i1}\ell_{11}$ 和 $a_{11}^2 = \ell_{11}^2$，其他列同理。以下是Cholesky算法的伪代码：
for j = 1, n
ℓ_{jj} = sqrt(a_{jj} - sum_{k = 1}^{j - 1} ℓ_{jk}^2)
for i = j + 1, n
ℓ_{ij} = (a_{ij} - sum_{k = 1}^{j - 1} ℓ_{ik}ℓ_{jk}) / ℓ_{jj}
endfor
endfor
1.1 Cholesky分解的优点
- 稳定性 ：Cholesky分解算法是稳定的，因此不需要进行选主元操作。
- 内存和运算量 ：Cholesky算法大约只需要LU分解一半的内存和一半的运算量。
- 正定性的重要性 ：矩阵的正定性对于在不进行部分选主元的情况下求出 $\ell_{ij}$ 非常重要，因为部分选主元可能会破坏矩阵 $A$ 的对称性。
- 不完全Cholesky分解 ：在某些情况下，需要进行不完全或近似的Cholesky分解，例如作为预条件子来加速迭代求解器的收敛。这可以通过将 $L$ 中与原始（假设为稀疏）矩阵 $A$ 中对应零元素的位置填充为零来实现。
2. QR分解与Householder变换
LU分解并不是对矩阵 $A$ 进行分解的唯一方法。这里介绍的Householder变换是对一般矩阵进行高效分解的基础，可将矩阵 $A$ 分解为 $A = QR$，其中 $Q$ 是正交矩阵，$R$ 是上三角矩阵。
2.1 Householder变换
Householder变换的核心是将一个全向量转换为只有一个非零元素的特殊向量。这可以通过Householder矩阵来高效实现，Householder矩阵定义为：
[
H = I - \frac{2ww^T}{w^Tw}, \quad \forall w \neq 0
]
特别地，对于向量 $\alpha e_1 = (\alpha, 0, 0, \cdots, 0)^T$，如果适当地计算向量 $w$，就可以从任意向量 $x$ 创建出该向量。具体来说，令 $w = x + \text{sign}(x_1) \cdot |x|_2 e_1$，其中 $x_1$ 是向量 $x$ 的第一个元素。这种从 $Hx$ 到 $(\alpha, 0, 0, \cdots, 0)^T$ 的变换称为Householder变换。
以下是描述Householder变换的算法伪代码：
xm = max {|x1|, |x2|, ..., |xn|}
for k = 1, n
wk = xk / xm
endfor
α = sign(w1) [w1^2 + w2^2 + ... + wn^2]^(1/2)
w1 = w1 + α
α = -α * xm
经过上述步骤，得到的期望向量为 $(\alpha, 0, 0, \cdots, 0)^T$。该操作的运算次数与 $O(n)$ 成正比，即只需要线性的工作量就能完成这个重要操作。
2.2 示例
考虑向量 $x = \begin{pmatrix} 1 \ 2 \ -3 \end{pmatrix}$，则 $x_m = 3$，$w$ 的分量的中间值为：
$w_1 = \frac{1}{3}$，$w_2 = \frac{2}{3}$，$w_3 = -1$
$\alpha$ 的中间值为：
$\alpha = + \sqrt{\frac{1}{3^2} + \frac{2^2}{3^2} + (-1)^2} = 1.2472$
更新后的值为：
$w_1 = \frac{1}{3} + 1.2472 = 1.5805$，$w_2 = \frac{2}{3}$，$w_3 = -1$
$\alpha = -1.2472 \cdot 3 = -3.7416$
所以，期望向量为 $Hx = \begin{pmatrix} -3.7416 \ 0 \ 0 \end{pmatrix}$
2.3 Householder矩阵的矩阵 - 向量乘法
与Householder矩阵 $H$ 进行矩阵 - 向量乘法是一个 $O(n)$ 操作，而一般的矩阵 - 向量乘法是 $O(n^2)$ 操作。这可以通过以下关系实现：
$Hx = x - \beta w (w^T x)$
其中 $\beta^{-1} = \frac{w^T w}{2}$ 是一个标量。在实际计算中，右侧的式子可以在一个循环内完成：
$x_i = x_i - \beta \cdot \gamma \cdot w_i, \quad i = 1, \cdots, n$
其中 $\gamma = w^T x$ 也是标量。显然，在这种情况下，我们不需要显式地构造 $H$，因为显式构造会导致 $O(n^2)$ 的操作。
2.4 QR分解的步骤
为了对一个 $n \times n$ 的方阵 $A$ 进行QR分解，我们考虑其列，并依次应用Householder变换来将每列的次对角线元素置零。这与LU分解类似，需要进行 $(n - 1)$ 个阶段。
2.4.1 第一阶段
考虑 $A$ 的第一列，确定一个Householder矩阵 $H_1$，使得：
[
H_1
\begin{bmatrix}
a_{11} \
a_{21} \
\vdots \
a_{n1}
\end{bmatrix}
=
\begin{bmatrix}
\alpha_1 \
0 \
\vdots \
0
\end{bmatrix}
]
为了确定 $H_1$，我们只需应用Householder变换算法来得到长度为 $n$ 的向量 $w_1$。经过第一阶段后，我们用 $A_1$ 覆盖 $A$，其中：
[
A_1 =
\begin{pmatrix}
\alpha_1 & a_{12}^
& \cdots & a_{1n}^
\
0 & a_{22}^
& \cdots & \cdots \
\vdots & \vdots & \ddots & \vdots \
0 & a_{n2}^
& \cdots & a_{nn}^*
\end{pmatrix}
= H_1 A
]
第一列之后的所有元素都是新的（用星号表示）。
2.4.2 第二阶段
接下来考虑更新后的矩阵 $A \equiv A_1 = H_1 A$ 的第二列，只取对角线以下的部分，得到：
[
H_2^
\begin{bmatrix}
a_{22}^
\
a_{32}^
\
\vdots \
a_{n2}^
\end{bmatrix}
=
\begin{bmatrix}
\alpha_2 \
0 \
\vdots \
0
\end{bmatrix}
]
这会得到一个长度为 $(n - 1)$ 的向量 $w_2$。该向量唯一地定义了Householder矩阵：
[
H_2^
= I - \frac{2w_2 w_2^T}{w_2^T w_2}
]
与 $H_1$ 不同，这里我们首先需要将 $H_2^
$ “扩展” 为 $H_2$，然后用 $H_2 A_1$ 覆盖 $A$，其中：
[
H_2 =
\begin{bmatrix}
1 & \cdots & 0 \
0 & H_2^*
\end{bmatrix}
]
2.4.3 第 $k$ 阶段
在QR过程的第 $k$ 阶段，通过求解：
[
H_k^
\begin{bmatrix}
a_{kk}^
\
a_{k + 1, k}^
\
\vdots \
a_{nk}^
\end{bmatrix}
=
\begin{bmatrix}
\alpha_k \
0 \
\vdots \
0
\end{bmatrix}
]
得到一个长度为 $(n - k + 1)$ 的向量 $w_k$。同样，我们用 $A_k = H_k A_{k - 1}$ 覆盖 $A$，并将 $H_k^
$ “扩展” 为：
[
H_k =
\begin{bmatrix}
I_{k - 1} & 0 \
0 & H_k^
\end{bmatrix}
]
2.5 Householder算法的效率
Householder算法的效率基于高效的乘法 $A_{k + 1} = H_{k + 1} \cdot A_k$。这个乘法不应该显式地进行，而是使用前面介绍的 $O(n)$ 矩阵 - 向量乘法算法。具体来说，我们计算：
$(I - \beta w_{n - k} w_{n - k}^T) [a_{kj}]$
其中 $\beta^{-1} = \frac{w_{n - k}^T \cdot w_{n - k}}{2}$，$a_{kj}$ 表示 $A_k$ 的列，$j = k + 1, \cdots, n - k$。确定 $H_k^*$ 需要 $(n - k)$ 次操作，因为未知向量 $w$ 的长度为 $(n - k)$。
2.6 最终结果
经过 $(n - 1)$ 个阶段后，我们得到一个上三角矩阵 $R$，其对角元素为 $\alpha_k$（$k = 1, \cdots, n$），其他元素通过以下公式计算：
$r_{ij} = a_{ij} - \gamma w_{ik}$
其中 $n \geq i, j \geq k$，且
$\gamma = \beta \sum_{i = k}^{n} w_{ik} a_{ij}$
$\beta^{-1} = (w_{n - k + 1}^T \cdot w_{n - k + 1}) / 2$
矩阵 $R$ 可以通过在每个阶段从阶数为 $(n - k + 1)$ 的 $H_k^
$ 形成一个等价的 $n$ 阶Householder矩阵 $H_k$ 来构造。具体来说，我们设置：
[
H_k =
\begin{bmatrix}
I_{k - 1} & 0 \
0 & H_k^
\end{bmatrix}
]
并计算：
$A_k = H_k A_{k - 1}$
其中 $A_1 \equiv A$。最终，上三角矩阵为：
$R = A_{n - 1} = H_{n - 1} A_{n - 2} = \cdots = H_{n - 1} H_{n - 2} \cdots H_1 A$
我们可以通过设置 $Q^T = H_{n - 1} H_{n - 2} \cdots H_1$ 来反转这个等式，由于 $H_k$ 都是正交矩阵，所以 $Q^{-1} = Q^T$。因此，$Q \cdot R = Q Q^T A \Rightarrow A = QR$。我们现在可以计算正交矩阵 $Q$，即 $Q = H_1^T H_2^T \cdots H_{n - 1}^T$。
2.7 Householder QR分解的伪代码
Begin Loop: k = 1, ..., n - 1 (number of stages)
Solve H_k^*
\begin{bmatrix}
a_{kk}^* \\
a_{k + 1, k}^* \\
\vdots \\
a_{nk}^*
\end{bmatrix}
=
\begin{bmatrix}
\alpha_k \\
0 \\
\vdots \\
0
\end{bmatrix}
(obtain w_{n - k + 1})
r_{kk} = \alpha_k
Compute \beta^{-1} = w_{n - k + 1}^T \cdot w_{n - k + 1} / 2
Zero a_{ik}: i = k + 1, ..., n
Begin Loop: j = k + 1, ..., n
\gamma_j = 0
Begin Loop: q = k, ... n
\gamma_j = \gamma_j + \beta w_{q - k + 1, k} a_{qj}
End Loop
Begin Loop: i = k, ..., n
r_{ij} = a_{ij} - \gamma_j w_{i - k + 1, k}
a_{ij} = r_{ij}
End Loop
End Loop
End Loop
2.8 示例
考虑 $3 \times 3$ 的Hilbert矩阵：
[
A =
\begin{bmatrix}
1 & \frac{1}{2} & \frac{1}{3} \
\frac{1}{2} & \frac{1}{3} & \frac{1}{4} \
\frac{1}{3} & \frac{1}{4} & \frac{1}{5}
\end{bmatrix}
]
应用Householder QR算法：
2.8.1 第一阶段（$k = 1$）
求解：
[
H_1
\begin{bmatrix}
1 \
\frac{1}{2} \
\frac{1}{3}
\end{bmatrix}
=
\begin{bmatrix}
\alpha_1 \
0 \
0
\end{bmatrix}
\Rightarrow
w_1 =
\begin{bmatrix}
2.1666 \
0.5 \
0.3333
\end{bmatrix}
]
同时，$r_{11} = \alpha_1 = -1.1666$，$\beta = \frac{2}{w_1^T w_1} = 0.3956$。
对于 $j = 2$，计算：
$\gamma_2 = \beta [w_{11} a_{12} + w_{21} a_{22} + w_{31} a_{32}] = 0.5274$
所以：
$a_{12} := r_{12} = a_{12} - \gamma_2 w_{11} = -0.6429$
$a_{22} := r_{22} = a_{22} - \gamma_2 w_{21} = 0.0696$
$a_{32} := r_{32} = a_{32} - \gamma_2 w_{31} = 0.0795$
在下一步迭代，$j = 3$，类似地计算：
$\gamma_3 = \beta [w_1 a_{13} + w_2 a_{23} + w_3 a_{33}] = 0.3615$
所以：
$a_{13} := r_{13} = a_{13} - \gamma_3 w_{11} = -0.4500$
$a_{23} := r_{23} = a_{23} - \gamma_3 w_{21} = 0.0692$
$a_{33} := r_{33} = a_{33} - \gamma_3 w_{31} = 0.0795$
同时，$a_{21} = a_{31} = 0$
2.8.2 第二阶段（$k = 2$）
求解：
[
H_2^*
\begin{bmatrix}
0.0696 \
0.0795
\end{bmatrix}
=
\begin{bmatrix}
\alpha_2 \
0
\end{bmatrix}
\Rightarrow
w_2 =
\begin{bmatrix}
2.3095 \
1.0000
\end{bmatrix}
]
同时，$r_{22} = \alpha_2 = -0.1017$，$\beta = 2 / (w_2^T \cdot w_2) = 0.315759$。
对于 $j = 3$，计算：
$\gamma_3 = \beta [w_{22} a_{23} + w_{32} a_{33}] = 0.0756$
$r_{23} = a_{23} - \gamma_3 w_{22} = -0.1053$
$r_{33} = a_{33} - \gamma_3 w_{32} = 0.0039$
2.9 计算成本
这里描述的QR分解的计算复杂度由计算矩阵 $H_k^*$ 的成本（$O(n - k)$）和从 $A_k = H_k A_{k - 1}$ 计算 $A_k$ 的成本（$O((n - k)^2)$ 操作）决定。因此，总的成本为：
[
\sum_{k = 1}^{n - 1} (n - k) + \sum_{k = 1}^{n - 1} (n - k)^2 = \frac{(n - 1)(n)}{2} + \frac{(n - 1)(n)(2n - 1)}{6} \approx \frac{n^3}{3}
]
更仔细地计算常数因子，并考虑加法和乘法，QR分解的复杂度为 $O(\frac{4}{3} n^3)$，而LU分解的复杂度为 $O(\frac{2}{3} n^3)$。
2.10 QR分解的其他特点
-
求解线性系统
：在得到矩阵 $A$ 的QR分解后，我们可以将线性系统 $Ax = b$ 求解如下：
$Ax = b \Rightarrow QRx = b \Rightarrow Q^T QRx = Q^T b \Rightarrow Rx = Q^T b$
这是一个上三角系统，可以通过回代法求解。 - 稳定性 ：基于Householder变换的QR分解总是稳定的。
- 与LU分解的比较 ：Householder QR分解的成本大约是LU分解的两倍，并且会扩展稀疏矩阵 $A$ 的带宽。相比之下，LU分解可以保持矩阵的带宽，但可能会受到数值不稳定性的影响。不过，对于大型矩阵 $A$，即使进行部分选主元，LU分解也比QR分解更高效。
-
Givens旋转
：除了Gram - Schmidt和Householder方法外，还可以使用Givens旋转矩阵来进行QR分解。Givens旋转矩阵定义为：
[
R(\theta) \equiv
\begin{bmatrix}
\cos \theta & \sin \theta \
-\sin \theta & \cos \theta
\end{bmatrix}
]
通过设置 $R(\theta) \cdot \begin{bmatrix} x \ y \end{bmatrix} = \begin{bmatrix} \sqrt{x^2 + y^2} \ 0 \end{bmatrix}$，可以得到 $\cos \theta = \frac{x}{\sqrt{x^2 + y^2}}$，$\sin \theta = \frac{y}{\sqrt{x^2 + y^2}}$。基于这个思想，可以构造一个一般的旋转矩阵来在每个迭代周期中消去一个元素，而不是像Householder变换那样消去一列。然而，Givens旋转的成本是Householder变换的两倍，是LU分解的四倍。它不用于求解方阵线性系统，而是用于最小二乘线性系统的求解和特征值求解器。
2.11 Hessenberg和三对角化约简
我们现在考虑将矩阵 $A$ 变换为一个上三角矩阵，且其第一下对角线元素非零。这种矩阵称为上Hessenberg矩阵，形式如下：
[
He \equiv
\begin{bmatrix}
* & * & * & \cdots & * \
* & * & \cdots & * \
\vdots & \vdots & \ddots & \vdots \
* & * & & * \
0 & \cdots & & *
\end{bmatrix}
]
Householder变换过程也可用于对一般矩阵进行这种变换。如果矩阵是对称的，那么得到的矩阵是三对角矩阵。
约简算法包括 $(n - 2)$ 个消元阶段。我们希望得到 $He = H \cdot A \cdot H^T$，其中矩阵 $H = H_{n - 2} \cdot H_{n - 1} \cdots H_1$。这些是Householder矩阵，通过将 $A$ 及其更新版本的子列置零来计算。例如：
[
H_1 =
\begin{bmatrix}
I_1 & 0 \
0 & H_1^
\end{bmatrix}
]
其中
[
H_1^
\begin{bmatrix}
a_{21} \
a_{31} \
\vdots \
a_{n1}
\end{bmatrix}
=
\begin{bmatrix}
\alpha_1 \
0 \
\vdots \
0
\end{bmatrix}
]
且 $A_1 = H_1 A H_1^T$ 等等。
2.12 示例
考虑 $3 \times 3$ 的Hilbert矩阵：
[
A =
\begin{bmatrix}
1 & \frac{1}{2} & \frac{1}{3} \
\frac{1}{2} & \frac{1}{3} & \frac{1}{4} \
\frac{1}{3} & \frac{1}{4} & \frac{1}{5}
\end{bmatrix}
]
将其变换为上Hessenberg矩阵。由于 $A$ 是对称的，我们期望得到一个三对角矩阵 $He$。
在第一阶段（$k = 1$），求解：
[
H_1^*
\begin{bmatrix}
\frac{1}{2} \
\frac{1}{3}
\end{bmatrix}
=
\begin{bmatrix}
\alpha \
0
\end{bmatrix}
]
得到 $w_2 = \begin{bmatrix} 2.2019 \ 0.6666 \end{bmatrix}$ 和 $\alpha = -0.6009$。
然后，使用 $a_{ij} = a_{ij} - \gamma_i w_j$ 得到 $A$ 的新元素：
[
He = H_1 A H_1^T =
\begin{bmatrix}
1 & 0 & 0 \
0 & -0.8321 & -0.5547 \
0 & -0.5547 & 0.8321
\end{bmatrix}
]
2.13 Hessenberg线性系统的求解
Hessenberg线性系统 $He \cdot x = b$ 可以通过带部分选主元的高斯消元法求解，并且保证是稳定的。包括部分选主元成本在内，其计算复杂度为 $O(n^2)$。
2.14 代码实现
以下是一个计算给定矩阵 $A$ 的上Hessenberg矩阵的函数：
void Hessenberg(SCMatrix &A){
int i,j,k,q;
double beta,gamma;
int N = A.Rows();
SCVector *x = new SCVector(N),
*w = new SCVector(N);
for(k=0;k<N-2;k++){
A.GetColumn(k,*x,k+1);
A(k+1,k) = HouseholderTrans(*x,*w);
beta = 2.0/dot(N-k-1,*w,*w);
for(i=k+2;i<N;i++)
A(i,k) = (*w)(i-k);
for(j=k+1;j<N;j++){
gamma = 0.0;
for(q=k+1;q<N;q++)
gamma += beta*(*w)(q-k-1)*A(q,j);
for(i=k+1;i<N;i++){
A(i,j) = A(i,j) - gamma*(*w)(i-k-1);
}
}
for(i=0;i<N;i++){
gamma = 0.0;
for(q=k+1;q<N;q++)
gamma += beta*(*w)(q-k-1)*A(i,q);
for(j=k+1;j<N;j++){
A(i,j) = A(i,j) - gamma*(*w)(j-k-1);
}
}
}
delete x;
delete w;
}
2.15 代码说明
- 引用传递 ：我们传递矩阵 $A$ 的引用（在参数列表中表示为 “SCMatrix &A”），因为我们希望用新的上Hessenberg矩阵替换 $A$ 中的值。如果省略 “&” 而采用值传递，那么在函数内部对矩阵 $A$ 所做的修改在函数返回调用程序时将丢失。
- 动态分配 ：在函数内部，我们动态分配了两个新的SCVector，并将它们分配给两个指针 $w$ 和 $x$。为了使用与SCVector类关联的 “( )” 运算符，我们首先使用一元运算符 “ ” 来获取指针所指向的对象。因此，表达式 $( w)$ 产生指针 $w$ 所指向的对象。表达式周围的额外括号用于确保 “*” 操作在 “( )” 运算符之前执行。
3. 预条件共轭梯度法（PCGM）
3.1 共轭梯度法（CGM）的收敛速率
对于线性系统 $Ax = b$，我们在之前已经介绍过共轭梯度法（CGM）。通过定义第 $k$ 次迭代的残差 $r_k \equiv b - Ax_k$，可以计算解和搜索方向：
$x_{k + 1} = x_k + \alpha_k p_k$
$p_{k + 1} = r_{k + 1} + \beta_k p_k$
残差也可以迭代计算：
$r_{k + 1} = r_k - \alpha_k A p_k$
通过将式（9.13）中的向量 $p_k$ 代入式（9.12），并使用残差定义，可以推导出三项递推公式：
$x_{k + 1} = (1 + \gamma_k) x_k + \alpha_k (b - Ax_k) - \gamma_k x_{k - 1}$
其中 $\gamma_k \equiv \frac{\alpha_k \beta_{k - 1}}{\alpha_{k - 1}}$，这有时被称为Rutishauser公式。对称性和正交性共同导致了我们在本书中多次看到的熟悉的三项神奇公式。
我们在之前提到过，CGM 等价于最小化一个适当定义的二次型。实际上，可以证明如果 $s$ 是精确解，那么共轭梯度迭代 $x_k$ 在由 $A$ 的幂和正交方向定义的 $k$ 维Krylov子空间上最小化范数 $|s - x|_A$。这个空间定义为：
$K_k(A, p) = \text{span} {p, Ap, A^2 p, \cdots, A^{k - 1} p}$
Rutishauser公式与Lanczos三项公式（见第10.3.6节）相似，因此CGM和Lanczos方法是相关的——它们都使用相同的Krylov子空间，并且都有三项递推公式。具体来说，可以为CGM构造一个三对角矩阵 $T_j$：
$T_j = P^T AP$
这里 $P = RD^{-1}$，其中 $D$ 是一个包含残差大小的对角矩阵，$R$ 是两个矩阵的乘积，第一个矩阵由正交搜索方向 $p_j$ 的列组成，第二个矩阵是一个主对角线为 1，主对角线上方为标量 $\beta_j$ 的双对角矩阵。
在实际应用中，即使迭代总数未达到 $n$（$n \times n$ 是对称正定矩阵 $A$ 的大小），CGM 也能产生很好的结果。我们回忆起在第4.1.7节中提到的共轭方向基本定理，它保证在 $n$ 次迭代后可以得到精确解。更具体地说，如果矩阵 $A$ 只有 $m \leq n$ 个不同的特征值，那么在 $m$ 次迭代中就可以实现收敛。舍入误差和相应的正交性损失是导致与理论结果偏差的原因，尽管舍入误差不像Lanczos方法中那样严重。
在实际应用中，我们总是使用收敛停止准则，而不是进行 $n$ 次迭代的循环。为此，收敛测试中的容差 $\epsilon$ 应该与初始残差的相对减少成比例，或者在某些情况下用右侧向量进行归一化，例如：
$|r_{k + 1}|
2 \leq \epsilon |b|_2$
这里 $\epsilon \approx 10^{-d}$，其中 $d$ 是所需的精度位数。计算 $|r
{k + 1}|
2 = (r
{k + 1}, r_{k + 1})$ 不需要额外的工作，因为这个量在 $\beta_k$ 的公式的分子中使用。需要注意的是，残差 $|r_{k + 1}|_2$ 可能不会单调递减，尽管解误差 $|s - x_k|_2$ 是单调递减的。这是因为在最小化过程中，直接目标是解的误差，而不是残差。直接最小化残差在实际中更为常见。
3.2 预条件共轭梯度法的作用
预条件共轭梯度法（PCGM）是对共轭梯度法的改进。当矩阵 $A$ 的条件数较大时，CGM 的收敛速度会变慢。预条件共轭梯度法通过引入一个预条件矩阵 $M$ 来改善这种情况。预条件矩阵 $M$ 应该满足以下条件：
- $M$ 是对称正定的。
- $M$ 容易求解，即求解 $M z = r$ 相对容易。
预条件共轭梯度法的基本思想是将原线性系统 $Ax = b$ 转化为一个等价的线性系统 $M^{-1}Ax = M^{-1}b$，使得新系统的条件数更小，从而加速收敛。
3.3 预条件共轭梯度法的步骤
以下是预条件共轭梯度法的详细步骤：
1.
初始化
：
- 选择初始猜测解 $x_0$。
- 计算初始残差 $r_0 = b - Ax_0$。
- 求解预条件系统 $M z_0 = r_0$。
- 令 $p_0 = z_0$。
2.
迭代过程
：
- 对于 $k = 0, 1, 2, \cdots$，直到满足收敛条件：
- 计算 $\alpha_k = \frac{r_k^T z_k}{p_k^T A p_k}$。
- 更新解 $x_{k + 1} = x_k + \alpha_k p_k$。
- 计算新的残差 $r_{k + 1} = r_k - \alpha_k A p_k$。
- 求解预条件系统 $M z_{k + 1} = r_{k + 1}$。
- 计算 $\beta_k = \frac{r_{k + 1}^T z_{k + 1}}{r_k^T z_k}$。
- 更新搜索方向 $p_{k + 1} = z_{k + 1} + \beta_k p_k$。
3.
收敛判断
：
- 检查是否满足收敛条件，例如 $|r_{k + 1}|
2 \leq \epsilon |b|_2$。如果满足，则停止迭代，输出 $x
{k + 1}$ 作为解；否则，继续迭代。
3.4 预条件矩阵的选择
预条件矩阵的选择是预条件共轭梯度法的关键。常见的预条件矩阵选择方法有：
-
对角预条件子
：$M$ 是一个对角矩阵，其对角元素为 $A$ 的对角元素。这种预条件子计算简单，但效果可能有限。
-
不完全Cholesky分解预条件子
：如前面所述，通过对矩阵 $A$ 进行不完全Cholesky分解得到预条件矩阵 $M$。这种预条件子在处理稀疏矩阵时效果较好。
-
ILU预条件子
：不完全LU分解预条件子，通过对矩阵 $A$ 进行不完全LU分解得到预条件矩阵 $M$。
3.5 示例代码
以下是一个简单的预条件共轭梯度法的Python示例代码：
import numpy as np
def pcg(A, b, M, x0, max_iter=1000, tol=1e-6):
x = x0
r = b - np.dot(A, x)
z = np.linalg.solve(M, r)
p = z
for k in range(max_iter):
Ap = np.dot(A, p)
alpha = np.dot(r, z) / np.dot(p, Ap)
x = x + alpha * p
r_new = r - alpha * Ap
if np.linalg.norm(r_new) < tol:
break
z_new = np.linalg.solve(M, r_new)
beta = np.dot(r_new, z_new) / np.dot(r, z)
p = z_new + beta * p
r = r_new
z = z_new
return x
# 示例矩阵和向量
A = np.array([[4, 1], [1, 3]])
b = np.array([1, 2])
M = np.diag(np.diag(A)) # 对角预条件子
x0 = np.zeros(len(b))
# 调用PCG求解
x = pcg(A, b, M, x0)
print("解:", x)
3.6 复杂度分析
预条件共轭梯度法的复杂度主要取决于矩阵 $A$ 的稀疏性和预条件矩阵 $M$ 的求解复杂度。在每次迭代中，主要的计算步骤包括矩阵 - 向量乘法 $A p$ 和预条件系统的求解 $M z = r$。如果矩阵 $A$ 是稀疏的，并且预条件矩阵 $M$ 的求解复杂度较低，那么PCGM的效率会很高。
3.7 与其他方法的比较
- 与CGM的比较 ：预条件共轭梯度法在矩阵条件数较大时比共轭梯度法收敛更快。通过引入预条件矩阵，可以改善系统的条件数，从而加速收敛。
- 与直接法的比较 ：直接法（如LU分解、Cholesky分解）可以在有限的步骤内得到精确解，但对于大型稀疏矩阵，直接法的存储和计算成本可能很高。预条件共轭梯度法是一种迭代法，对于大型稀疏矩阵，它可以在较少的迭代次数内得到近似解，并且存储成本较低。
3.8 应用场景
预条件共轭梯度法在许多领域都有广泛的应用，例如：
-
计算物理
：在求解偏微分方程时，经常会得到大型稀疏线性系统，PCGM可以有效地求解这些系统。
-
机器学习
：在某些优化问题中，需要求解线性系统，PCGM可以作为一种高效的求解方法。
-
工程计算
：在结构分析、流体力学等领域，也会遇到大型稀疏线性系统，PCGM可以提高计算效率。
4. 总结
本文介绍了几种快速线性求解器，包括Cholesky分解、QR分解和预条件共轭梯度法。这些方法在不同的场景下具有各自的优势：
| 方法 | 适用矩阵类型 | 优点 | 缺点 | 计算复杂度 |
| ---- | ---- | ---- | ---- | ---- |
| Cholesky分解 | 对称正定矩阵 | 稳定，内存和运算量约为LU分解的一半 | 要求矩阵为对称正定 | $O(\frac{2}{3}n^3)$ |
| QR分解（Householder变换） | 一般矩阵 | 稳定，可用于求解线性系统 | 成本约为LU分解的两倍，会扩展稀疏矩阵带宽 | $O(\frac{4}{3}n^3)$ |
| 预条件共轭梯度法 | 对称正定矩阵 | 对于大型稀疏矩阵收敛快，存储成本低 | 需要选择合适的预条件矩阵 | 取决于矩阵稀疏性和预条件矩阵求解复杂度 |
在实际应用中，我们应该根据矩阵的特点和问题的需求选择合适的求解方法。例如，对于对称正定矩阵，Cholesky分解是一个不错的选择；对于一般矩阵，QR分解可以保证稳定性；而对于大型稀疏对称正定矩阵，预条件共轭梯度法可能是最有效的方法。
通过合理选择和应用这些快速线性求解器，可以提高线性系统求解的效率和精度，从而为各种科学计算和工程应用提供有力的支持。
mermaid图展示PCGM的流程：
graph TD;
A[初始化: x0, r0, z0, p0] --> B[迭代开始];
B --> C[计算αk];
C --> D[更新xk+1];
D --> E[计算rk+1];
E --> F[求解Mz_k+1 = r_k+1];
F --> G[计算βk];
G --> H[更新pk+1];
H --> I{是否收敛};
I -- 是 --> J[输出解xk+1];
I -- 否 --> B;
通过以上内容，我们对快速线性求解器有了更深入的了解，希望这些知识能够帮助你在实际问题中更好地选择和应用合适的求解方法。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
