On the backward stability of s-step GMRES
Erin Carson Yuxin Ma

SUSY Quantum Mechanics, (non)-Analyticity and … Phase Transitions
    Alexander V. Turbiner
    September 9, 2024
=================================================================



bstract.tex
ntroduction.tex
reliminary.tex
tability-sstep.tex
rthogonalization.tex
rylov-basis.tex
xperiments.tex
onclusion.tex
cknowledgments.tex

abbrvurl




§ PROPERTIES OF BCGSI+


Given X ∈ℝ^m × n with m ≥ n, in Algorithm <ref>, we present the (k+1)-th step of the BCGSI+ algorithm to compute X = QT with an orthonormal matrix Q ∈ℝ^m × n and an upper triangular matrix T ∈ℝ^n × n.
Note that MGS or any unconditionally stable QR algorithm, e.g., Householder QR or Tall-Skinny QR (TSQR), described in <cit.>, can be utilized in Line <ref>, while any backward stable QR algorithm, i.e.,

    Ŵ^(2) + Δ W^(2) = Q̂_ks+1:(k+1)sT̂^(2)  with  Δ W^(2)≤() Ŵ^(2),

can be employed in Line <ref>; see <cit.> for details.



The properties of BCGSI+ have already been studied in <cit.>.
Based on the results in <cit.>, it is easy to obtain the following lemmas.

 
    Let Q_1:js and T_1:js be computed by Algorithm <ref>.
    If assuming () κ(X_1:js) < 1, then
    
    X_1:js + Δ X_1:js = Q_1:jsT_1:js,   Δ X_i≤() X_i

    for any i ≤ js, and
    
    Q_1:jsQ_1:js -  I≤().




    Similarly to the proof of <cit.>, it is easy to verify (<ref>), since each line of Algorithm <ref> is column-wise backward stable.
    The conclusion (<ref>) is directly followed by <cit.>.


From Lemma <ref>, = (), and (<ref>) is satisfied when the Q-factor is not well-conditioned.
Then it remains to estimate  defined by (<ref>).

 
    Assume that for is+j, Q_1:is+j and T_1:is+j are computed by Algorithm <ref>.
    If
    
    Q_1:is+j-1Q_1:is+j-1 - I≤(), 
       Q_1:is+jQ_1:is+j - I > (),

    then
    
    max{Q_is+jT_is+j,is+j, T_is+j,is+j}≤() X_is+j.




    Without loss of generality, we only need to prove the case j = 1 since Algorithm <ref> is columnwise backward stable.
    From the assumption (<ref>),
    
    Q_1:is+1Q_1:is+1 - I≤() + Q_is+1Q_is+1 - I + 2Q_1:isQ_is+1.

    Note that <cit.> do not depend on <cit.> when using an unconditionally stable QR algorithm in Line <ref> of Algorithm <ref>.
    From <cit.> it holds that
    
    Q_is+1Q_is+1 - I≤().

    Then it remains to estimate Q_1:isQ_is+1, which can be bounded as
    
    Q_1:isQ_is+1   ≤Q_1:isW̃^(2)_1 (T^(2)_1,1)^-1 + Q_1:isΔW̃^(2)_1 (T^(2)_1,1)^-1
       ≤(I - Q_1:isQ_1:is) Q_1:isU_1/T^(2)_1,1
                + Q_1:isΔW̃^(2)_1/T^(2)_1,1,

    where W^(2) = W̃^(2) + ΔW̃^(2) with W̃^(2) = (I - Q_1:isQ_1:is) U.
    Together with <cit.> and the assumption (<ref>), we obtain
    
    Q_1:isQ_is+1≤()/T^(2)_1,1.

    Then together with <cit.>, it follows that
    
    W̃^(2)_1 + Δ W^(2)_1 = Q_is+1T^(2)_1,1,   Δ W^(2)_1≤(),

    and further,
    
    Q_1:isQ_is+1   ≤() Q_is+1/(I - Q_1:isQ_1:is) U_1 - Δ W^(2)_1
       ≤()/(I - Q_1:isQ_1:is) U_1 - ().

    This means that Q_1:is+1Q_1:is+1 - I≤() if 2 (I - Q_1:isQ_1:is) U_1 > () holds.
    Furthermore, the contrapositive is that 2 (I - Q_1:isQ_1:is) U_1≤() if Q_1:is+1Q_1:is+1 - I > () guaranteed by the assumption (<ref>).
    Then by <cit.>, and
    
    T_is+1,is+1 = T^(2)_1,1T^(1)_1,1 + Δ T_is+1,is+1  with  Δ T_is+1,is+1≤() T^(2)_1,1T^(1)_1,1,

    we have
    
    Q_is+1T_is+1,is+1≤Q_is+1T^(2)_1,1T^(1)_1,1 + () Q_is+1T^(2)_1,1T^(1)_1,1≤() X_is+1,

    which also implies T_is+1,is+1≤() X_is+1 by noticing Q_is+1≥ 1 - ().





§ EXAMPLE FOR CLASSICAL S-STEP GMRES WITH BCGSPIPI+


 
    We construct the linear system Ax = b, where A is a 20-by-20 random matrix with κ(A) = 10^10 generated using the MATLAB command  and .
    The vector b is selected as the right singular vector corresponding to the fourth largest singular value, and the initial guess x_0 is the zero vector.
    
    For this specific linear system, the relative backward error Ax - b/Ax + b of the solution computed by using standard GMRES with CGSPIPI+, namely BCGSPIPI+ with s = 1, is approximately 10^-8.
    In contrast, using standard GMRES with CGSI+ results in an error of around 10^-16, as illustrated in Figure <ref>.
    This difference occurs because CGSPIPI+ cannot generate a nearly orthonormal basis V̂_1:20 when the condition number of r̂   Ŵ_1:19 exceeds approximately 10^8, implying () κ^2(r̂   Ŵ_1:19) > 1.
    Under this situation, the return value of V̂_20 is a  vector.
    Thus, it is not possible to obtain a more accurate solution than x^(19), whose backward error is approximately 10^-8.
    
    A similar result occurs for s-step GMRES with s = 2.
    The relative backward error using BCGSPIPI+ is approximately 10^-5, while for s-step GMRES with BCGSI+, the error is around 10^-11, as illustrated in Figure <ref>.











§ PROOF OF LEMMA <REF>



    First, we aim to prove (<ref>) by induction.
    For the base case, from (<ref>) with i=1, there exists Y_1:s such that B_1:s = V_1:s Y_1:s.
    Assume that B_(i-1)s+1:is = V_(i-1)s+1:is Y_(i-1)s+1:is holds for all i ≤ j - 1.
    Then our aim is to prove that it holds for j.
    Since V_1:(j-1)s is orthonormal and B_(j-1)s+1:js is the Q-factor of (I - V_1:(j-1)s V_1:(j-1)s)^2 K_(j-1)s+1:js, we obtain V_1:(j-1)s B_(j-1)s+1:js = 0.
    Together with the above assumptions on i ≤ j - 1 and (<ref>), there exists Y_(j-1)s+1:js such that
    
    B_1:js = V_1:jsY_1:(j-1)s    V_1:(j-1)s B_(j-1)s+1:js
     0     Y_(j-1)s+1:js
            = V_1:jsY_1:(j-1)s    0 
     0     Y_(j-1)s+1:js,

    which gives (<ref>) by induction on j.

    Then we will bound κ (B̂_1:ks).
    By the definition of ω_B_i, we only need to consider the off-diagonal blocks B̂_(i-1)s+1:isB̂_(j-1)s+1:js.
    From (<ref>) and dropping the quadratic terms, it holds that
    
    B̂_(i-1)s+1:isB̂_(j-1)s+1:js
         ≤Δ B_(i-1)s+1:isV̂_(j-1)s+1:jsỸ_(j-1)s+1:js
            + Ỹ_(i-1)s+1:isV̂_(i-1)s+1:isV̂_(j-1)s+1:jsỸ_(j-1)s+1:js
            + Δ B_(j-1)s+1:jsV̂_(i-1)s+1:isỸ_(i-1)s+1:is
         ≤Δ B_(i-1)s+1:isỸ_(j-1)s+1:js
        + Ỹ_(i-1)s+1:isV̂_(i-1)s+1:isV̂_(j-1)s+1:jsỸ_(j-1)s+1:js
            + Δ B_(j-1)s+1:jsỸ_(i-1)s+1:is,

    which implies that
    
    ∑_i, j = 1; i ≠ j^k    B̂_(i-1)s+1:isB̂_(j-1)s+1:js^2 
    ≤    3 ∑_i, j = 1; i ≠ j^k (Δ B_(i-1)s+1:is^2 Ỹ_(j-1)s+1:js^2 
        + Ỹ_(i-1)s+1:is^2 V̂_(i-1)s+1:isV̂_(j-1)s+1:js^2 Ỹ_(j-1)s+1:js^2 
        + Δ B_(j-1)s+1:js^2 Ỹ_(i-1)s+1:is^2) 
    ≤    6k ∑_i = 1^k (Δ B_(i-1)s+1:is^2 Ỹ_(j-1)s+1:js^2) 
                + 3 max_i(Ỹ_(j-1)s+1:js^4) ω_k^2.

    Together with B̂_(i-1)s+1:is≤√(s) (1 + ω_B_i) and
    
    Ỹ_(i-1)s+1:is≤B̂_(i-1)s+1:is + Δ B_(i-1)s+1:is/1 - V̂_(i-1)s+1:isV̂_(i-1)s+1:is - I,

    we obtain
    
    B̂_1:ksB̂_1:ks - I^2 
          = ∑_i=1^k B̂_(i-1)s+1:isB̂_(i-1)s+1:is - I^2 
                + ∑_i, j = 1; i ≠ j^k B̂_(i-1)s+1:isB̂_(j-1)s+1:js^2 
         ≤∑_i=1^k ω_B_i^2
                + 6k ∑_i = 1^k (Δ B_(i-1)s+1:is^2 Ỹ_(j-1)s+1:js^2) + 3 max_i(Ỹ_(j-1)s+1:js^4) ω_k 
         ≤∑_i=1^k ω_B_i^2
                + 24 ks ∑_i = 1^k (Δ B_(i-1)s+1:is^2) 
                + 48 s^2 ω_k^2

    by dropping the quadratic terms.
    This implies that, from the assumption (<ref>),
    
    B̂_1:ksB̂_1:ks - I≤∑_i=1^k ω_B_i
                + 5 √(ks)∑_i = 1^k Δ B_(i-1)s+1:is + 7 s ω_k≤1/2.

    Thus, (B̂_1:ks) can be bounded by
    
    (B̂_1:ks) ≥ 1 - B̂_1:ksB̂_1:ks - I≥1/2

    and further
    
    κ (B̂_1:ks) ≤∑_i B̂_(i-1)s+1:is/(B̂_1:ks)≤ 2 √(∑_i = 1^k s  (1 + ω_B_i))≤ 2 √(n) + √(s).




