对称矩阵的三对角分解(Lanzos分解算法）-MINRES算法预热-Toy模板网

这篇具有很好参考价值的文章主要介绍了对称矩阵的三对角分解(Lanzos分解算法）-MINRES算法预热。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

这篇博客看完以后接着看下一篇博客添加链接描述专门介绍MINRES算法实现就容易了

Lanzos分解

首先介绍Lanczos分解，Lanzos把对称矩阵转换为一个三对角对称矩阵。考虑三对角对称矩阵如下，考虑正交分解
$T = Q^T A Q$
$T=\left(\begin{array}{cccccc} \alpha_1 & \beta_1 & 0 & \cdots & 0 & 0\\ \beta_1 & \alpha_2 & \beta_2 & 0 & \cdots & 0 \\ 0 & \beta_2 & \alpha_3 & \beta_3 & \cdots & 0\\ 0 & \cdots & \cdots & \cdots & \cdots & 0\\ 0 & \cdots & 0 & \beta_{n-2} & \alpha_{n-1} & \beta_{n-1}\\ 0 & 0 & \cdots & 0 & \beta_{n-1} & \alpha_n \end{array}\right)$

下面重点考虑正交矩阵 $Q$ 和三对角矩阵 $T$ 的形成，我们记 $Q=[q_1,q_2,\ldots,q_n],q_i \in R^n$ ，则根据 $Q T = A Q$ ，我们会得到下面这个等式，约定 $\beta_0 q_0 = \beta_n q_n = 0$ ：
$Aq_i = \beta_{i - 1}q_{i - 1} + \alpha_i q_i + \beta_i q_{i + 1},1 \leq i \leq n.$
我们先考虑 $\alpha_1,\beta_1$ 的确定，任意取一个向量 $q_1 \in R^n, \| q_1 \|_2 = 1$ ，则有
$\left\{\begin{array}{l} Aq_1 = \alpha_1 q_1 + \beta_1 q_2,\\ \alpha_1 = q_1^T A q_1,\\ \beta_1 = \|Aq_1 - \alpha_1 q_1\|_2,\\ q_2 = (Aq_1 - \alpha_1 q_1)/\beta_1. \end{array}\right.$
类似地，假设我们已经得到了 $[q_1,q_2,\ldots,q_k]$ 和 $[\alpha_1,\alpha_2,\ldots,\alpha_{k - 1}],[\beta_1,\beta_2,\ldots,\beta_{k-1}]$ ，下面同样可以类似得到 $\alpha_k,\beta_k,q_{k+1}$ ，
$\left\{\begin{array}{l} Aq_k = \beta_{k - 1}q_{k - 1} + \alpha_k q_k + \beta_k q_{k+1},\\ \alpha_k = q_k^T (Aq_k - \beta_{k - 1}q_{k - 1}),\\ \beta_k = \|Aq_k - \beta_{k - 1}q_{k - 1}\|_2,\\ q_{k+1} = (Aq_k - \beta_{k - 1}q_{k - 1})/\beta_k. \end{array}\right.$
引入记号 $Q_k = [q_1,q_2,\ldots,q_k]$ ，和 $T_k$ ，其中有
$T_k=\left(\begin{array}{cccccc} \alpha_1 & \beta_1 & 0 & \cdots & 0 & 0\\ \beta_1 & \alpha_2 & \beta_2 & 0 & \cdots & 0 \\ 0 & \beta_2 & \alpha_3 & \beta_3 & \cdots & 0\\ 0 & \cdots & \cdots & \cdots & \cdots & 0\\ 0 & \cdots & 0 & \beta_{k-2} & \alpha_{k-1} & \beta_{k-1}\\ 0 & 0 & \cdots & 0 & \beta_{k-1} & \alpha_k\\ \end{array}\right)$
则有 $AQ_k = Q_k T_k + \beta_{k} q_{k + 1} e_{k}^T$ ，其中 $e_k \in \boldsymbol{R}^k$ 且最后一个元素为1，其余元素为0，如果 $\beta_k \neq = 0,1 \leq k \leq n - 1$ ，则可以顺利得到对称矩阵A的三对角分解。如果求解过程中某个 $\beta_{k_0}= 0$ ，那么显然可以根据 $AQ_k = Q_k T_k$ 得到前 $k_0 - 1$ 个特征向量，如果要计算剩下的特征向量和特征值，只需要重新初始化一个 $q_1$ 即可。

MINRES算法解读

MINRES主要应用于对称不定方程求解，考虑线性方程组如下\eqref{line}所示，其中 $\in \boldsymbol{R}^{n \times n}$ 是对称矩阵。
$A x = b$
MINRES算法的出发点是寻求一个向量 $x^{(k)} \in x^{(0)} + \kappa_k(A,r_0)$ ，即
$x^{k} = \min_{x \in x^0 + \kappa_k(A,r_0)} \| Ax -b \|_2.$
其中 $\kappa_k(A,r_0) = (r_0,Ar_0,\ldots,A^{k-1}r_0)$ 形成的向量空间。

选择 $q_1 = \frac{r^0}{\|r^0\|},r^0 = b - Ax^0$ ，则 $(r^0,Ar^0,\ldots,A^{k-1}r^0)=(q_1,Aq_1,\ldots,A^{k-1}q_1)\|r^0\|_2$ ，再根据 $Aq_k = \beta_{k - 1}q_{k - 1} + \alpha_k q_k + \beta_k q_{k+1}$ 可得 $q_{k+1} \in \kappa_{k}(q_1) = \kappa_{k} (r^0)$ ，由此可以得到对于任意一个向量 $x$ 属于 $Q_k$ 张成的向量空间，则必有 $\in \kappa_{k}(r^0)$ ，反之亦然。
引入 $Q_{k+1} = [Q_k,q_{k+1}]$ ，则根据上面的推导可以得到 $AQ_k = Q_{k+1}T_{k+1,k}$ ，
$T_{k+1,k}=\left(\begin{array}{c} T_k \\ \beta_k e_k^T \end{array}\right)$

于是乎可以把原始问题转换为下面这个形式，特别注意，这里的 $e_1 \in \boldsymbol{R}^{k+1}$ ，上面的 $e_k \in \boldsymbol{R}^k$ ：
$\left\{\begin{aligned} x^{k} &= \min_{x \in x^0 + \kappa_k(A,r_0)} \| Ax -b \|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| A(x^0 + Q_k y) -b \|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| r^0 - AQ_k y \|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| r^0 - Q_{k+1}T_{k+1,k} y\|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| \|r^0\|_2 q_1 - Q_{k+1}T_{k+1,k} y\|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| \|r^0\|_2 Q_{k+1}e_1 - Q_{k+1}T_{k+1,k} y\|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| Q_{k + 1}(\|r^0\|_2 e_1 - T_{k+1,k} y)\|_2 \\ & = \min_{y \in \boldsymbol{R}^{k}} \| \|r^0\|_2 e_1 - T_{k+1,k} y\|_2. \end{aligned}\right.$

QR分解求解二范数优化问题

引入记号 $a = \|r^0\|_2$ ，原始问题进一步转化为求下面这个问题，对于这个问题，我们引入QR分解来处理。
$\min_{y \in \boldsymbol{R}^k} \| a e_1 - T_{k + 1,k} y \|_2.$
考虑矩阵 $T_{k+1,k}$ 的QR分解，为了区分上面对称矩阵 $A$ 的正交分解，这里我们假设 $T_{k+1,k} = V_{k+1}R_{k+1,k}$ ，其中 $V_{k+1} \in R^{(k+1) \times (k+1)}$ 是正交矩阵，
$R_{k+1,k}=\left(\begin{array}{cccccc} r_{0,0} & & & \cdots & \cdots & \\ 0 & r_{1,1} & & & \cdots & \\ 0 & 0 & \cdots & \cdots & & \\ 0 & 0 & \cdots & \cdots & &r_{k-1,k-1} \\ 0 & 0 & \cdots & 0 & 0 & 0\\ \end{array}\right) =\left(\begin{array}{c} R_k \\ 0^T \\ \end{array}\right)$
有了上面的介绍以后，我们直接带入问题就可以得到
$\begin{aligned} \left\|a e_1-T_{k+1, k} y\right\|_2 & \left.=\| a e_1-V_{k+1} R_{k+1, k} y\right) \|_2 \\ & =\left\|V_{k+1}\left(V_{k+1}^{\top}\left(a e_1\right)-\left[\begin{array}{c} R_k \\ 0 \end{array}\right] y\right)\right\|_2 \\ & =\left\|\left(V_{k+1}^{\top}\left(a e_1\right)-\left[\begin{array}{c} R_k y\\ 0 \end{array}\right] \right)\right\|_2 \\ & =\left\|\left[V_{k+1, k}, v_{k+1}\right]^{\top}\left(a e_1\right)-\left[\begin{array}{c} R_k y \\ 0 \end{array}\right]\right\|_2, \end{aligned}$
如果 $R_k$ 非奇异，那么显然取 $y$ 满足 $V_{k + 1,k}^{T} a e_1 - R_{k}y = 0$ 即可，此时有$y = a R_{k}^{-1} V_{k+1,k}^T e_1 $，因此根据 M I N R E S 算法得到的第$ k$次迭代向量如下：
$x^{k} = x^0 + a Q_k R_{k}^{-1} V_{k+1,k}^T e_1$

Lanzos分解

MINRES算法解读

QR分解求解二范数优化问题

T k + 1 , k T_{k+1,k} Tk+1,k的QR分解实现

我们希望可以在 $T_{k,k-1}$ 的QR分解基础上做一次Givens变化得到 $T_{k+1,k}$ 的QR分解，假设 $T_{k,k-1}$ 的QR分解如下：
$T_{k,k-1} = (G_{k-1} G_{k-2} \ldots G_1)^T R_{k,k-1} = V_{k} R_{k,k-1}$
其中
$R_{k,k - 1}= \left(\begin{array}{cccccc} r_{0,0} & & & \cdots & \cdots & \\ 0 & r_{1,1} & & & \cdots & \\ 0 & 0 & \cdots & \cdots & & \\ 0 & 0 & \cdots & \cdots & &r_{k-1,k-1} \\ 0 & 0 & \cdots & 0 & 0 & 0\\ \end{array}\right) =\left(\begin{array}{c} R_{k - 1} \\ 0^T \\ \end{array}\right)$
以及Givens变换矩阵形如下面这个形式：
$=\left(\begin{array}{ccccccc} 1 & 0 & 0 & \cdots & \cdots & 0 & 0\\ 0 & 1 & 0 & 0 & \cdots & 0 & 0\\ 0 & 0 & \cdots & \cdots & 0 & 0 & 0\\ 0 & 0 & \cdots & c & \cdots & s & 0\\ 0 & 0 & \cdots & \cdots & \cdots & 0 & 0\\ 0 & 0 & \cdots & -s & \cdots & c & 0\\ 0 & 0 & \cdots & 0 & \cdots & 0 & 1\\ \end{array}\right)$
针对这个问题使用的Givens变换 $G_i$ 仅仅改变矩阵 $i, i + 1$ 的元素，因此 $G_i$ 的具体形式如下：
$G_i=\left[\begin{array}{cccc} I_{i-1} & & & \\ & c_i & s_i & \\ & -s_i & c_i & \\ & & & I_{k-i-1} \end{array}\right]$

我们重新回忆一下 $T_{k+1,k}$ 的形状，
$\left\{\begin{aligned} T_{k+1,k} &=\left(\begin{array}{c} T_k \\ \beta_k e_k^T \end{array}\right) =\left[\begin{array}{c|c} T_{k, k-1} & 0 \\ & \vdots \\ & \beta_{k-1} \\ & \alpha_k \\ \hline 0^T & \beta_k \end{array}\right] =\left[\begin{array}{c|c} V_{k}R_{k,k-1} & 0 \\ & \vdots \\ & \beta_{k-1} \\ & \alpha_k \\ \hline 0^T & \beta_k \end{array}\right] \\ & =\left[\begin{array}{cc} V_k & 0 \\ 0 & 1 \end{array}\right]\left[\begin{array}{c|c} R_{k, k-1} & V_k^{-1}\left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ \beta_{k-1} \\ \alpha_k \end{array}\right] \\ \hline 0 & \beta_k \end{array}\right] \\ & = \left[\begin{array}{cc} V_k & 0 \\ 0 & 1 \end{array}\right] \tilde{T}_{k+1, k} \end{aligned}\right.$
而且另一方面：

$Q_k^{-1}\left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ 0 \\ \beta_{k-1} \\ \alpha_k \end{array}\right]=G_{k-1} G_{k-2} \cdots G_1\left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ \beta_{k-1} \\ \alpha_k \end{array}\right]=G_{k-1} G_{k-2}\left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ 0 \\ \beta_{k-1} \\ \alpha_k \end{array}\right]=\left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ r_{k-3,k-1} \\ r_{k-2,k-1} \\ r_{k-1,k-1} \end{array}\right]$

有了上述的介绍以后，我们发现只要选择一个Givens变换把 $\tilde{T}_{k+1, k}$ 右下角的 $\beta_k$ 消去即可，因此选择 $G_k$ 如下：
$G_k=\left[\begin{array}{ccc} I_{k-1} & & \\ & c_k & s_k \\ & -s_k & c_k \end{array}\right] \in \boldsymbol{R}^{(k+1) \times(k+1)}$
于是我们可以得到 $T_{k+1,k}$ 的QR分解 $T_{k+1,k} = V_{k+1}R_{k+1,k}$ ，这里假设Givens变换是矩阵 $G$ 的维度自动增长。
$V_{k+1}=\left[\begin{array}{cc} V_k & 0 \\ 0 & 1 \end{array}\right] G_k^{\top}=\left(G_k \tilde{G}_{k-1} \cdots \tilde{G}_1\right)^{\top}$
注意，由于矩阵 $T_{k+1,k}$ 是三对角矩阵，而Givens变换只改变两行两列元素，因此 $R_{k+1,k}$ 是一个上三角的带宽为3的带状矩阵，以及观察这个变换过程可以发现 $R_{k-1}$ 其实是 $R_k$ 的 $(k - 1)$ 阶顺序主子阵。

x k x^k xk的具体递推公式

整理上文，我们有

$\left\{\begin{aligned} & A = Q^T T Q ,\\ & r^0 = b - A x^0 = \|r^0\|_2 q_1 = a q_1,\\ & Q_k = [q_1,q_2,\ldots, q_k], \\ & x^k = x^0 + Q_k y,\quad \min_{y \in R^k} \|b - Ax^0 - AQ_k y\| = \min_{y \in R^k} \|r^0 - AQ_k y\|, \\ & \min_{y \in R^k} \|r^0 - AQ_k y\| = \min_{y \in R^k} \|aQ_k e_1 - AQ_k y\| = \min_{y \in R^k} \|aQ_k e_1 - Q_k T_{k+1,k} y\|, \\ & \min_{y \in R^k} \|ae_1 - T_{k+1,k}y\|_2 = \min_{y \in R^k} \|ae_1 - V_{k+1}R_{k+1,k}y\|_2, y = aR_{k}^{-1}V_{k+1,k}^T e_1, \\ & T_{k+1,k}=\left(\begin{array}{c} T_k \\ \beta_k e_k^T \end{array}\right),V_{k+1}= [V_{k + 1,k}, v_{k + 1}], R_{k+1,k}=\left(\begin{array}{c} R_k \\ 0 \end{array}\right)\\ & x^{k} = x^0 + Q_k (a R_{k}^{-1} V_{k+1,k}^T e_1),\\ \end{aligned}\right.$ 文章来源地址https://www.toymoban.com/news/detail-437592.html