MINRES（极小残差算法）求解线性系统详细解读

本博客参考了添加链接描述这篇知乎
先看我这篇博客介绍添加链接描述

$Q_k{R}_k^{-1}$的处理
假设$D_k=[d_1,d_2,\dots ,d_k]=Q_k{R}_k^{-1}$，假设$R_k$的第$i$行第$j$列元素为$r_{i,j}$，
$\{\begin{array}{rl}& R_{k+1,k}=\left(\begin{array}{cccccc}{r}_{0,0}& r_{0,1}& r_{0,2}& \cdots & \cdots & 0\\ 0& r_{1,1}& r_{1,2}& r_{1,3}& \cdots & 0\\ 0& 0& \cdots & \cdots & \ast & \ast \\ 0& 0& \cdots & \cdots & 0& 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)\\ & R_k=\left(\begin{array}{cccccc}{r}_{0,0}& r_{0,1}& r_{0,2}& \cdots & \cdots & 0\\ 0& r_{1,1}& r_{1,2}& r_{1,3}& \cdots & 0\\ 0& 0& \cdots & \cdots & \ast & \ast \\ 0& 0& \cdots & \cdots & 0& r_{k-1,k-1}\end{array}\right)\end{array}$
则根据$D_k{R}_k=Q_k$，可以得到以下公式。
$\{\begin{array}{rl}& r_{0,0}{d}_1=q_1,\\ & r_{0,1}{d}_1+r_{1,1}{d}_2=q_2,\\ & \sum _{s=0}^{i-1}{r}_{s,i-1}{d}_{s+1}=q_i,i=3,4,5,\dots k-1,\end{array}$
从上述公式反解出$d_i$关于$q_i$的表达式如下：
$\{\begin{array}{rl}& d_1=\frac{q_1}{r_{0,0}},\\ & d_2=(q_2-r_{0,1}{d}_1)/r_{1,1},\\ & d_i=(q_i-r_{i-3,i-1}{d}_{i-2}-r_{i-2,i-1}{d}_{i-1})/r_{i-1,i-1},\end{array}$
由于$Q_k=[Q_{k-1},q_k]$和$R_{k-1}$是$R_k$的$(k-1)$阶顺序主子阵，可以得到
$\{\begin{array}{rl}{D}_k& =Q_k{R}_k^{-1}\\ & =\left[Q_{k-1},q_k\right]\left[\begin{array}{ll}{R}_{k-1}^{-1}& \ast \\ & \ast \end{array}\right]\\ & =\left[D_{k-1},d_k\right]\end{array}$

$a{V}_{k+1,k}^T{e}_1$的处理
由于$V_{k+1}=[V_{k+1,k},v_{k+1}]$，根据下面这个公式可以有：
$a{V}_{k+1}^T{e}_1=\left(\begin{array}{c}{V}_{k+1,k}^T\\ v_{k+1}^T\end{array}\right)(a{e}_1)=\left(\begin{array}{c}{V}_{k+1,k}^T(a{e}_1)\\ v_{k+1}^T(a{e}_1)\end{array}\right)=\left(\begin{array}{c}{\eta }^{(k)}\\ {\hat{\eta }}_{k+1}\end{array}\right)$
也就是说$a{V}_{k+1,k}^T{e}_1$和$a{V}_{k+1}^T{e}_1$的前$k$个元素相同，我们记${\eta }^{(k)}=a{V}_{k+1,k}^T{e}_1$，修改以后如上所示，因此转换完以后的二范数问题最优解\eqref{yk}的解是
$y^k=a{R}_k^{-1}{V}_{k+1,k}^T{e}_1=R_k^{-1}{\eta }^{(k)}.$
由于对$T_{k+1,k}$做QR分解的时候，利用Givens变换处理，根据上面推导过程得到的$V_k$其实也是$V_{k+1}$的$k$阶顺序主子阵，因此$V_{k+1}^T{e}_1$的前$k-1$个元素相同。
$V_{k+1}=\left[\begin{array}{cc}{V}_k& 0\\ 0& 1\end{array}\right]$
也就是说${\eta }^{(k)}=[{\eta }^{(k-1)},{\eta }_k]$，最后一直推导就有${\eta }^{(k)}=[{\eta }_1,{\eta }_2,\dots ,{\eta }_k]$
\subsection{$x^k$计算公式和残量计算}
$\bullet$ \quad 最后得到下面这个递推公式，其中$d_k$可以通过公式来递推得到，
$\{\begin{array}{rl}{x}^k& =x^0+Q_k(a{R}_k^{-1}{V}_{k+1,k}^T{e}_1)\\ & =x^0+D_k{\eta }^{(k)}\\ & =x^0+\left[D_{k-1},d_k\right]\left[\begin{array}{c}{\eta }^{(k-1)}\\ {\eta }_k\end{array}\right]\\ & =x^0+D_{k-1}{\eta }^{(k-1)}+{\eta }_k{d}_k\\ & =x^{k-1}+{\eta }_k{d}_k.\end{array}$
而根据上面的信息，我们知道${\eta }_k$是$a{V}_{k+1}^T{e}_1$的第$k$个分量，根据下面的公式发现只需要将对应的Givens变换作用在$(ae1)$即可得到$a{V}_{k+1}^T{e}_1$。
$a{V}_{k+1}^T{e}_1=V_{k+1}^T(ae1)={\left(G_k{\tilde{G}}_{k-1}\cdots {\tilde{G}}_1\right)}^{\top }(ae1).$
$\bullet$ \quad 残量的计算如下所示，即残量最终是等式$a{V}_{k+1}{e}_1$的最后一个分量的绝对值：
$\{\begin{array}{rl}\Vert r_k{\Vert }_2& =\Vert A{x}^k-b{\Vert }_2\\ & =\Vert a{q}_1-Q_{k+1}{T}_{k+1,k}{y}^k{\Vert }_2\\ & =\Vert Q_{k+1}(a{e}_1-T_{k+1,k}{y}^k){\Vert }_2\\ & =\Vert (a{e}_1-V_{k+1}{R}_{k+1,k}{y}^k){\Vert }_2\\ & =\Vert V_{k+1}(V_{k+1}^{T}a{e}_1-R_{k+1,k}{y}^k){\Vert }_2\\ & =\Vert (V_{k+1}^{T}a{e}_1-R_{k+1,k}{y}^k){\Vert }_2\\ & ={\Vert \left[\begin{array}{c}{\eta }^{(k)}\\ {\hat{\eta }}_{k+1}\end{array}\right]-\left[\begin{array}{c}{R}_k{y}_k\\ 0\end{array}\right]\Vert }_2\\ & =|{\hat{\eta }}_{k+1}|\end{array}$

算法流程:前三次迭代过程解析

为了更方便理解算法流程，我们先考虑$x^k,k=1，2，3$的更新过程。

\newpage%------------------------
$\bullet x^3=x^2+{\eta }_3{d}_3$

$\Rightarrow d_3=(q_3-r_{2,3}{d}_2-r_{1,3}{d}_1)/r_{2,2},q_3=(A{q}_2-{\beta }_1{q}_1)/{\beta }_2$

$\Rightarrow T_{4,3}=\left(\begin{array}{c}{T}_3\\ {\beta }_3{e}_3^T\end{array}\right)=\left(\begin{array}{ccc}{\alpha }_1& {\beta }_1& 0\\ {\beta }_1& {\alpha }_2& {\beta }_2\\ 0& {\beta }_2& {\alpha }_3\\ 0& 0& {\beta }_3\end{array}\right)$

$\Rightarrow {\alpha }_3=(q_3,A{q}_3),{\beta }_3=\Vert A{q}_3-{\beta }_2{q}_2-{\alpha }_3{q}_3{\Vert }_2$

先用Givens变换${\tilde{G}}_1{T}_{4,3}$

$\Rightarrow \left[\begin{array}{cccc}{c}_1& s_1& 0& 0\\ -s_1& c_1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left(\begin{array}{ccc}{\alpha }_1& {\beta }_1& 0\\ {\beta }_1& {\alpha }_2& {\beta }_2\\ 0& {\beta }_2& {\alpha }_3\\ 0& 0& {\beta }_3\end{array}\right)=\left[\begin{array}{ccc}{r}_{0,0}& {\hat{r}}_{0,1}& {\hat{r}}_{0,2}^{(1)}\\ 0& {\hat{r}}_{1,1}& {\hat{r}}_{1,2}^{(1)}\\ 0& {\beta }_2& {\alpha }_3\\ 0& 0& {\beta }_3\end{array}\right]$,

$\Rightarrow {\hat{r}}_{0,2}^{(1)}=s_1{\beta }_2,{\hat{r}}_{1,2}^{(1)}=c_1{\beta }_2$

再用Givens变换${\tilde{G}}_2{\tilde{G}}_1{T}_{4,3}$，从上面更新$x^2$的过程我们知道${\tilde{G}}_2$会把$[{\hat{r}}_{1,1},{\beta }_2{]}^T$这部分变成$[r_{1,1},0{]}^T$

$\Rightarrow \left[\begin{array}{cccc}1& 0& 0& 0\\ 0& c_2& s_2& 0\\ 0& -s_2& c_2& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{r}_{0,0}& {\hat{r}}_{0,1}& {\hat{r}}_{0,2}^{(1)}\\ 0& {\hat{r}}_{1,1}& {\hat{r}}_{1,2}^{(1)}\\ 0& {\beta }_2& {\alpha }_3\\ 0& 0& {\beta }_3\end{array}\right]=\left[\begin{array}{ccc}{r}_{0,0}& {\hat{r}}_{0,1}& {\hat{r}}_{0,2}\\ 0& {\hat{r}}_{1,1}& {\hat{r}}_{1,2}\\ 0& 0& {\hat{r}}_{2,2}\\ 0& 0& {\beta }_3\end{array}\right]$

Givens变换把$[{\hat{r}}_{2,2},{\beta }_3{]}^T$这部分变成$[r_{2,2},0{]}^T$，选择$\gamma =\frac{{\hat{r}}_{2,2}}{{\beta }_3},s_3=\frac{1}{\sqrt{1+{\gamma }^2}},c_3=\gamma s_3$

$\Rightarrow \left[\begin{array}{ccc}{r}_{0,0}& r_{0,1}& r_{0,2}\\ 0& r_{1,1}& r_{1,2}\\ 0& 0& r_{2,2}\\ 0& 0& 0\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& c_3& s_3\\ 0& 0& -s_3& c_3\end{array}\right]\left[\begin{array}{ccc}{r}_{0,0}& {\hat{r}}_{0,1}& {\hat{r}}_{0,2}\\ 0& {\hat{r}}_{1,1}& {\hat{r}}_{1,2}\\ 0& 0& {\hat{r}}_{2,2}\\ 0& 0& {\beta }_3\end{array}\right]$,

$\Rightarrow [{\eta }_1,{\eta }_2,{\eta }_3]={\eta }^{(3)},\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\\ {\hat{\eta }}_4\end{array}\right]=V_4^T(a{e}_1)=G_3{\tilde{G}}_2{\tilde{G}}_1(a{e}_1),$ \

$\Rightarrow \left[\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\eta }_3\\ {\hat{\eta }}_4\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& c_3& s_3\\ 0& 0& -s_3& c_3\end{array}\right]\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\\ {\hat{\eta }}_3\\ 0\end{array}\right]=\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\\ c_3{\hat{\eta }}_3\\ -s_3\widehat{{\eta }_3}\end{array}\right]$,

$\Rightarrow {\mathbf{\eta }}_3={\mathbf{c}}_3{\hat{\mathbf{\eta }}}_3,{\mathbf{r}}_{2,2}={\mathbf{c}}_3{\hat{\mathbf{r}}}_{2,2}+{\mathbf{s}}_3{\mathbf{\beta }}_3,{\mathbf{r}}_{0,2}={\hat{\mathbf{r}}}_{0,2},{\mathbf{r}}_{1,2}={\hat{\mathbf{r}}}_{1,2},{\mathbf{d}}_3=\frac{{\mathbf{q}}_3-{\mathbf{r}}_{0,2}{\mathbf{d}}_1-{\mathbf{r}}_{1,2}{\mathbf{d}}_2}{{\mathbf{r}}_{2,2}}$,

$\Rightarrow x^3=x^2+{\eta }_3{d}_3$

$\{\begin{array}{rl}& {\eta }_1=a{c}_1,\\ & d_1=q_1/r_{0,0},\\ & q_1=r^0/a=(A{x}^0-b)/a,a=\Vert r^0{\Vert }_2,\\ & r_{0,0}={\alpha }_1{c}_1+{\beta }_1{s}_1,\\ & x^1=x^0+{\eta }_1{d}_1,\\ & {\hat{\eta }}_2=-a{s}_1.\end{array}$

$\{\begin{array}{rl}& {\eta }_2=c_2{\hat{\eta }}_2,\\ & d_2=(q_2-r_{0,1}{d}_1)/r_{1,1},\\ & q_2=(A{q}_1-{\alpha }_1{q}_1)/{\beta }_1,\\ & r_{0,1}={\hat{r}}_{0,1}=c_1{\beta }_1+s_1{\alpha }_2,\\ & r_{1,1}=c_2{\hat{r}}_{1,1}+s_2{\beta }_2,{\hat{r}}_{1,1}=-s_1{\beta }_1+c_1{\alpha }_2,\\ & x^2=x^1+{\eta }_2{d}_2,\\ & {\hat{\eta }}_3=-s_2{\hat{\eta }}_2.\end{array}$

$\{\begin{array}{rl}& {\eta }_3=c_3{\hat{\eta }}_3,\\ & d_3=(q_3-r_{0,2}{d}_1-r_{1,2}{d}_2)/r_{2,2},\\ & q_3=(A{q}_2-{\beta }_1{q}_1-{\alpha }_2{q}_2)/{\beta }_2,\\ & r_{0,2}=s_1{\beta }_2,r_{1,2}=c_1{\beta }_2,\\ & r_{2,2}=c_3{\hat{r}}_{2,2}+s_3{\beta }_3,{\hat{r}}_{2,2}=-s_2{\hat{r}}_{1,2}^{(1)}+c_2{\alpha }_3,{\hat{r}}_{1,2}^{(1)}=c_1{\beta }_2,\\ & x^2=x^1+{\eta }_2{d}_2,\\ & {\hat{\eta }}_3=-s_3{\hat{\eta }}_3.\end{array}$
通过前三个变量的更新过程，我们可以发现$R_{k+1,k}$的最后一列其实就是$G_k{\tilde{G}}_{k-1}\dots {\tilde{G}}_1$作用在$T_{k+1,k}$最后一列得到的，然后我们知道$T_{k+1,k}$最后一列是$[0,\dots ,0,{\beta }_{k-1},{\alpha }_k,{\beta }_k{]}^T$以及$G_i$只改变第$i,i+1$行元素，所以我们有：\

$\left[\begin{array}{c}0\\ ⋮\\ 0\\ r_{k-3,k-1}\\ r_{k-2,k-1}\\ r_{k-1,k-1}\\ 0\end{array}\right]=G_k{\tilde{G}}_{k-1}{\tilde{G}}_{k-2}\left[\begin{array}{c}0\\ ⋮\\ 0\\ {\beta }_{k-1}\\ {\alpha }_k\\ {\beta }_k\end{array}\right],k\ge 3$ \

$\left[\begin{array}{c}{r}_{k-3,k-1}\\ r_{k-2,k-1}\\ r_{k-1,k-1}\\ 0\end{array}\right]=G_k{\tilde{G}}_{k-1}{\tilde{G}}_{k-2}\left[\begin{array}{c}0\\ {\beta }_{k-1}\\ {\alpha }_k\\ {\beta }_k\end{array}\right]=G_k{\tilde{G}}_{k-1}\left[\begin{array}{c}{r}_{k-3,k-1}\\ {\hat{\beta }}_{k-1}\\ {\alpha }_k\\ {\beta }_k\end{array}\right]=G_k\left[\begin{array}{c}{r}_{k-3,k-1}\\ r_{k-2,k-1}\\ {\hat{\alpha }}_k\\ {\beta }_k\end{array}\right]=\left[\begin{array}{c}{r}_{k-3,k-1}\\ r_{k-2,k-1}\\ r_{k-1,k-1}\\ 0\end{array}\right],k\ge 3$

$\{\begin{array}{rl}& {\eta }_k=c_k{\hat{\eta }}_k,\\ & d_k=(q_k-r_{k-3,k-1}{d}_{k-2}-r_{k-2,k-1}{d}_{k-1})/r_{k-1,k-1},\\ & q_k=(A{q}_{k-1}-{\beta }_{k-2}{q}_{k-2}-{\alpha }_{k-1}{q}_{k-1})/{\beta }_{k-1},\\ & r_{k-3,k-1}=s_{k-2}{\beta }_{k-1},\\ & r_{k-2,k-1}=c_{k-1}{\hat{\beta }}_{k-1}+s_{k-1}{\alpha }_k,{\hat{\beta }}_{k-1}=c_{k-2}{\beta }_{k-1},\\ & r_{k-1,k-1}=c_k{\hat{\alpha }}_k+s_k{\beta }_k,{\hat{\alpha }}_k=-s_{k-1}{\hat{\beta }}_k+c_{k-1}{\alpha }_k,\\ & x^k=x^{k-1}+{\eta }_k{d}_k,\\ & {\hat{\eta }}_{k+1}=-s_k{\hat{\eta }}_k.\end{array}$

详细代码以及介绍参考本人知乎添加链接描述

import numpy as np
import time
N = 1000
A = np.zeros([N,N])
penalty = 10
for i in range(N):
    for j in range(i + 1,N):
        A[i,j] = np.random.rand(1)
        A[j,i] = A[i,j]
    A[i,i] = N + penalty*abs(np.random.rand(1))#主对角占优，确保非奇异
#print(A)
x = np.random.rand(N,1)
b = A@x

def MINRES(A,b,x0,N,eps):

    r = b - A@x0
    a = np.linalg.norm(r)
    # q_old = q_{k-1},q_new = q_{k},c_old = c_{k-2},c_mid = c_{k-1},c_new = c_k
    q_old = 0*r.copy();q_new = r/a
    d_old = 0*r.copy();d_mid = 0*r.copy()
    xi_old = a;beta_old = 0
    r_old = 0;r_mid = 0
    c_old = 0;c_mid = 0
    s_old = 0;s_mid = 0
    ls = np.zeros([N,N])
    for k in range(N):
        w = A@q_new - beta_old*q_old
        alpha = np.dot(w.T,q_new)[0,0]
        w = w - alpha*q_new
        beta_new = np.linalg.norm(w)
        
        if k == 0:
            alpha_hat = alpha
        elif k == 1:
            r_mid = c_mid*beta_old + s_mid*alpha
            alpha_hat = -s_mid*beta_old + c_mid*alpha
        else:
            r_old = s_old*beta_old
            beta_hat = c_old*beta_old
            
            r_mid = c_mid*beta_hat + s_mid*alpha
            alpha_hat = -s_mid*beta_hat + c_mid*alpha
        if abs(alpha_hat) > abs(beta_new):
            gamma = beta_new/alpha_hat
            c_new = 1.0/np.sqrt(1 + gamma**2);s_new = c_new*gamma
        else:
            gamma = alpha_hat/beta_new
            s_new = 1.0/np.sqrt(1 + gamma**2);c_new = s_new*gamma
        
        r_new = c_new*alpha_hat + s_new*beta_new
        
        xi_new = -s_new*xi_old
        xi_old = c_new*xi_old
        
        d_new = (q_new - r_old*d_old - r_mid*d_mid)/r_new
        x0 += xi_old*d_new
        print('res:%.2e,real res:%.2e'%(abs(xi_new),np.linalg.norm(b - A@x0)))
        ls[:,k:k + 1] = q_new
        if abs(xi_new) < eps:
            break
        else:
            xi_old = xi_new
            q_old = q_new.copy()
            q_new = w.copy()/beta_new
            
            beta_old = beta_new
            c_old = c_mid
            c_mid = c_new
            s_old = s_mid
            s_mid = s_new
            
            r_old = r_mid
            r_mid = r_new
            
            d_old = d_mid.copy()
            d_mid = d_new.copy()
            
            #print(w,beta_new)
            
            
            
    if abs(xi_new) < eps:
        print('success')
    else:
        print('fail')
    return x0,ls
min_t0 = time.time()
x0 = np.random.rand(N,1)

eps = 1e-7 
x_p,ls = MINRES(A,b,x0,N,eps) 
min_ela = time.time() - min_t0
print('time:%.2f,err:%.2e'%(min_ela,max(abs(x_p - x))))

min_t0 = time.time()
xp = np.linalg.solve(A,b)
min_ela = time.time() - min_t0
print('time:%.2f,err:%.2e'%(min_ela,max(abs(xp - x))))