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

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

$Q_k R_{k}^{-1}$ 的处理
假设 $D_k = [d_1,d_2, \ldots, d_k] = Q_k R_{k}^{-1}$ ，假设 $R_k$ 的第 $i$ 行第 $j$ 列元素为 $r_{i,j}$ ，
$\left\{\begin{aligned} & 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 & *& *\\ 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 & *& *\\ 0 & 0 & \cdots & \cdots & 0 & r_{k - 1,k - 1}\\ \end{array}\right) \end{aligned}\right.$
则根据 $D_k R_k = Q_k$ ，可以得到以下公式。
$\left\{\begin{aligned} & 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,\ldots k - 1,\\ \end{aligned}\right.$
从上述公式反解出 $d_i$ 关于 $q_i$ 的表达式如下：
$\left\{\begin{aligned} & 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{aligned}\right.$
由于 $Q_k = [Q_{k - 1},q_k]$ 和 $R_{k-1}$ 是 $R_k$ 的 $(k - 1)$ 阶顺序主子阵，可以得到
$\left\{\begin{aligned} D_{k} &= Q_k R_{k}^{-1} \\ & = \left[Q_{k-1}, q_k\right]\left[\begin{array}{ll} R_{k-1}^{-1} & * \\ & * \end{array}\right]\\ & =\left[D_{k-1}, d_k\right] \end{aligned}\right.$

$a V_{k+1,k}^T e_1$ 的处理
由于 $V_{k + 1} = [V_{k + 1,k},v_{k + 1}]$ ，根据下面这个公式可以有：
$V_{k + 1}^T e_1=\left(\begin{array}{c} V_{k + 1,k}^T \\ v_{k + 1}^T \\ \end{array}\right) (ae_1)=\left(\begin{array}{c} V_{k + 1,k}^T (ae_1) \\ v_{k + 1}^T (ae_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$ 和 $aV_{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,\ldots,\eta_k]$
\subsection{ $x^k$ 计算公式和残量计算}
$\bullet$ \quad 最后得到下面这个递推公式，其中 $d_k$ 可以通过公式来递推得到，
$\left\{\begin{aligned} 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{aligned}\right.$
而根据上面的信息，我们知道 $\eta_k$ 是 $a V_{k + 1}^T e_1$ 的第 $k$ 个分量，根据下面的公式发现只需要将对应的Givens变换作用在 $(a e 1)$ 即可得到 $a V_{k + 1}^T e_1$ 。
$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 残量的计算如下所示，即残量最终是等式 $aV_{k+1}e_1$ 的最后一个分量的绝对值：
$\left\{\begin{aligned} \|r_k\|_2 &= \|Ax^k - b\|_2 \\ &= \|aq_1 - Q_{k+1} T_{k+1,k} y^k\|_2 \\ &= \|Q_{k+1} (ae_1 - T_{k+1,k} y^k)\|_2 \\ &= \|(ae_1 - V_{k+1} R_{k+1,k} y^k)\|_2 \\ &= \|V_{k+1} (V_{k+1}^T ae_1 - R_{k+1,k} y^k)\|_2 \\ &= \|(V_{k+1}^T ae_1 - R_{k+1,k} y^k)\|_2 \\ &= \left\|\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]\right\|_2 \\ & = |\hat{\eta}_{k+1}| \end{aligned}\right.$

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

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

\newpage%------------------------
$\bullet \quad x^3 = x^2 + \eta_3 d_3$

$\Rightarrow \quad d_3 = (q_3 - r_{2,3} d_2 - r_{1,3} d_1)/r_{2,2},q_3 = (Aq_2 - \beta_1 q_1)/\beta_2$

$\Rightarrow \quad 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 \quad \alpha_3 = (q_3,Aq_3),\beta_3 = \|Aq_3 - \beta_2 q_2 - \alpha_3 q_3 \|_2$

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

$\Rightarrow \quad \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 \quad \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 \quad \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 \quad \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 \quad [\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} (ae_1) = G_3 \tilde{G}_2 \tilde{G}_1 (ae_1) ,$ \

$\Rightarrow \quad \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 \hat{\eta_3} \end{array}\right]$ ,

$\Rightarrow \quad \boldsymbol{\eta_3 = c_3 \hat{\eta}_3,r_{2,2} =c_3 \hat{r}_{2,2} + s_3 \beta_3,r_{0,2} = \hat{r}_{0,2},r_{1,2} = \hat{r}_{1,2},d_3 = \frac{q_3 - r_{0,2}d_1 - r_{1,2} d_2 }{r_{2,2}} }$ ,

$\Rightarrow \quad x^3 = x^2 + \eta_3 d_3$

$\left\{\begin{aligned} & \eta_1 = a c_1,\\ & d_1 = q_1/r_{0,0},\\ & q_1 = r^0/a = (Ax^0 - b)/a, a = \|r^0\|_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{aligned}\right.$

$\left\{\begin{aligned} & \eta_2 = c_2 \hat{\eta}_2,\\ & d_2 = (q_2 - r_{0,1} d_1)/r_{1,1},\\ & q_2 = (Aq_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{aligned}\right.$

$\left\{\begin{aligned} & \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 = (Aq_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{aligned}\right.$
通过前三个变量的更新过程，我们可以发现 $R_{k+1,k}$ 的最后一列其实就是 $G_k \tilde{G}_{k-1} \ldots \tilde{G}_1$ 作用在 $T_{k+1,k}$ 最后一列得到的，然后我们知道 $T_{k+1,k}$ 最后一列是 $[0,\ldots,0,\beta_{k-1},\alpha_k,\beta_k]^T$ 以及 $G_i$ 只改变第 $i, i + 1$ 行元素，所以我们有：\

$\left[\begin{array}{c} 0 \\ \vdots \\ 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 \\ \vdots \\ 0 \\ \beta_{k-1} \\ \alpha_k \\ \beta_k \end{array}\right] ,k \geq 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 \geq 3$

$\left\{\begin{aligned} & \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 = (Aq_{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{aligned}\right.$
在这里插入图片描述

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

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))))