Error estimators and their analysis for CG, Bi-CG and GMRES

TL;DR

This paper analyzes the uncertainty in convergence criteria for iterative linear solvers like CG, Bi-CG, and GMRES, proposing new error estimators that improve accuracy and computational efficiency for large, ill-conditioned problems.

Contribution

It introduces novel error estimators for CG, Bi-CG, and GMRES algorithms, enhancing stability and accuracy in convergence assessment for large-scale problems.

Findings

Proposed an $ ext{O}(n)$ error estimator for Bi-CG.

Modified GMRES error estimator for increased stability.

Error estimators improve efficiency and accuracy in large, ill-conditioned systems.

Abstract

The demands of accuracy in measurements and engineering models today, renders the condition number of problems larger. While a corresponding increase in the precision of floating point numbers ensured a stable computing, the uncertainty in convergence when using residue as a stopping criterion has increased. We present an analysis of the uncertainty in convergence when using relative residue as a stopping criterion for iterative solution of linear systems, and the resulting over/under computation for a given tolerance in error. This shows that error estimation is significant for an efficient or accurate solution even when the condition number of the matrix is not large. An $\mathcal{O}(1)$ error estimator for iterations of the CG algorithm was proposed more than two decades ago. Recently, an $\mathcal{O}(k^2)$ error estimator was described for the GMRES algorithm which allows for non-symmetric linear systems as well, where $k$ is the iteration number. We suggest a minor modification in this GMRES error estimation for increased stability. In this work, we also propose an $\mathcal{O}(n)$ error estimator for A-norm and $l_{2}$ norm of the error vector in Bi-CG algorithm. The robust performance of these estimates as a stopping criterion results in increased savings and accuracy in computation, as condition number and size of problems increase.

Equations

$$\kappa \frac{\Vert r_k\Vert }{\Vert b\Vert }\ge \frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }\ge \frac{\Vert r_k\Vert }{\Vert b\Vert }\frac{1}{\kappa }$$

$$\Vert {ϵ}_k{\Vert }_A^2\approx \Vert r_0{\Vert }_2^2[(T_{k+d}^{-1}{)}_{(1,1)}-(T_k^{-1}{)}_{(1,1)}]$$

$$\Vert {ϵ}_k{\Vert }_A^2\approx \sum _{j=k}^{k+d-1}{\alpha }_j{r}_j^T{r}_j$$

$$H_k=\left(\begin{array}{cc}{H}_{k-d}& W_{k-d}\\ Y_{k-d}^T& {\tilde{H}}_{k-d}\end{array}\right)$$

$$\gamma_{k-d}=\dfrac{h_{k-d+1,k-d}\big{(}e_{k-d},H_{k-d}^{-1}e_{1}\big{)}}{1-h_{k-d+1,k-d}\big{(}e_{k-d},H_{k-d}^{-1}w_{k-d}\big{)}}$$

$$w_{k-d}=W_{k-d}{\tilde{H}}_{k-d}^{-1}{e}_1$$

$${\delta }_{k+1}=\frac{h_{k+1,k}^2}{1+h_{k+1,k}^2{t}_{kk}}$$

$$u_k={\delta }_{k+1}{t}_k$$

$$\frac{{\Vert {ϵ}_{k-d}\Vert }^2}{{\Vert r_0\Vert }^2}\approx {\gamma }_{k-d}^2{\Vert {\tilde{H}}_{k-d}^{-1}{e}_1\Vert }^2+{\Vert {\gamma }_{k-d}{H}_{k-d}^{-1}{w}_{k-d}+(e_{k-d},H_{k-d}^{-1}{e}_1)u_{k-d}\Vert }^2$$

$$\frac{{\Vert {ϵ}_{k-d}\Vert }^2}{{\Vert r_0\Vert }^2}=\underset{\text{term1}}{\underbrace{\frac{{\Vert {ϵ}_{k-d}\Vert }^2-{\Vert {ϵ}_k\Vert }^2}{{\Vert r_0\Vert }^2}}}+\underset{\text{term2}}{\underbrace{{\Vert s_k\Vert }^2}}+\underset{\text{term3}}{\underbrace{{\zeta }_k}}$$

$$\begin{array}{rl}{\zeta }_k=& 2h_{k+1,k}\left(e_k,H_k^{-1}{e}_1\right)\left((H_n^{-1}{e}_{k+1}{)}^k,H_k^{-1}{e}^1+s_k\right)\\ & -2h_{k-d+1,k-d}\left(e_{k-d},H_{k-d}^{-1}{e}_1\right)\left((H_n^{-1}{e}_{k-d+1}{)}^{k-d}-(H_k^{-1}{e}_{k-d+1}{)}^{k-d},H_{k-d}^{-1}{e}_1+s_{k-d}\right)\end{array}$$

$$\frac{{\Vert {ϵ}_{k-d}\Vert }^2}{{\Vert r_0\Vert }^2}\approx \left|{\gamma }_{k-d}^2{\Vert {\tilde{H}}_{k-d}^{-1}{e}_1\Vert }^2+{\Vert {\gamma }_{k-d}{H}_{k-d}^{-1}{w}_{k-d}+(e_{k-d},H_{k-d}^{-1}{e}_1)u_{k-d}\Vert }^2-{\Vert s_k\Vert }^2\right|$$

$$\Vert {ϵ}_k{\Vert }_A^2=\left|{ϵ}_k^{T}A{ϵ}_k\right|=\left|r_k^T(A^T{)}^{-1}{r}_k\right|=\left|r_k^T{A}^{-1}{r}_k\right|$$

$$\Vert {ϵ}_k{\Vert }_2^2={ϵ}_k^T{ϵ}_k=r_k^T(A^T{)}^{-1}{A}^{-1}{r}_k$$

$$\begin{array}{rl}{r}_k^T(A^T{)}^{-1}{r}_k-r_{k+1}^T(A^T{)}^{-1}{r}_{k+1}& =r_k^T(A^T{)}^{-1}{r}_k-(r_k-{\alpha }_{k}A{p}_k{)}^T(A^T{)}^{-1}(r_k-{\alpha }_{k}A{p}_k)\\ & ={\alpha }_k{p}_k^T{r}_k+{\alpha }_k{r}_k^T(A^T{)}^{-1}A{p}_k-{\alpha }_k^2{p}_k^{T}A{p}_k\end{array}$$

$$\begin{array}[]{ccl}r_{k}^{T}A^{-1}r_{k}-r_{k+1}^{T}A^{-1}r_{k+1}=\alpha_{k}r_{k}^{T}p_{k}+\alpha_{k}r_{k}^{T}(A^{T})^{-1}(\dfrac{r_{k}-r_{k+1}}{\alpha_{k}})-\alpha_{k}^{2}p_{k}^{T}Ap_{k}\\ \implies r_{k+1}^{T}A^{-1}r_{k+1}=-\alpha_{k}r_{k}^{T}p_{k}+r_{k+1}^{T}A^{-1}r_{k}+\alpha_{k}^{2}p_{k}^{T}Ap_{k}\end{array}$$

$$\begin{array}[]{ccl}\epsilon_{k}=\sum_{j=k}^{n}\alpha_{j}p_{j}\implies\epsilon_{k}-\epsilon_{k+d}=\sum_{j=k}^{k+d}\alpha_{j}p_{j}\\ \implies\epsilon_{k}\approx\sum_{j=k}^{k+d}\alpha_{j}p_{j}\end{array}$$

$$\Vert {ϵ}_{k+1}{\Vert }_A^2=r_{k+1}^T{A}^{-1}{r}_{k+1}\approx -{\alpha }_k{r}_k^T{p}_k+r_{k+1}^T(\sum _{j=k}^{k+d}{\alpha }_j{p}_j)+{\alpha }_k^2{p}_k^{T}A{p}_k$$

$$\begin{gathered} r_k^T(A^T{)}^{-1}{A}^{-1}{r}_k-r_{k+1}^T(A^T{)}^{-1}{A}^{-1}{r}_{k+1}=r_k^T(A^T{)}^{-1}{A}^{-1}{r}_k \\ -(r_k-{\alpha }_{k}A{p}_k{)}^T(A^T{)}^{-1}{A}^{-1}(r_k-{\alpha }_{k}A{p}_k) \end{gathered}$$

$$r_k^T(A^T{)}^{-1}{A}^{-1}{r}_k-r_{k+1}^T(A^T{)}^{-1}{A}^{-1}{r}_{k+1}=2{\alpha }_k{p}_k^T{A}^{-1}{r}_k-{\alpha }_k^2\Vert p_k{\Vert }^2$$

$$\Vert {ϵ}_{k+1}{\Vert }_2^2=r_{k+1}^T(A^T{)}^{-1}{A}^{-1}{r}_{k+1}\approx (\sum _{j=k}^{k+d}{\alpha }_j{p}_j{)}^T(\sum _{j=k}^{k+d}{\alpha }_j{p}_j)-2{\alpha }_k{p}_k^T(\sum _{j=k}^{k+d}{\alpha }_j{p}_j)+{\alpha }_k^2\Vert p_k{\Vert }^2$$

$${\kappa }_f(A,x)=\Vert A\Vert \frac{\Vert x\Vert }{\Vert Ax\Vert }$$

$${\kappa }_b(A,b)=\Vert A^{-1}\Vert \frac{\Vert b\Vert }{\Vert A^{-1}b\Vert }$$

$$\kappa =\Vert A\Vert \Vert A^{-1}\Vert$$

$$\max \left(\frac{\frac{{\chi }_k}{\Vert x\Vert }-\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }}{\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }},\frac{\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }-\frac{{\chi }_k}{\Vert x\Vert }}{\frac{{\chi }_k}{\Vert x\Vert }}\right)$$

$$\frac{\left|\frac{{\chi }_k}{\Vert x\Vert }-\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }\right|}{\min \left(\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert },\frac{{\chi }_k}{\Vert x\Vert }\right)}$$

$$U.R.(\chi )=\frac{1}{n-d}\sum _{k=0}^{n-d-1}\left(\frac{\left|\frac{{\chi }_k^2}{\Vert x{\Vert }^2}-\frac{\Vert {ϵ}_k{\Vert }^2}{\Vert x{\Vert }^2}\right|}{\min \left(\frac{\Vert {ϵ}_k{\Vert }^2}{\Vert x{\Vert }^2},\frac{{\chi }_k^2}{\Vert x{\Vert }^2}\right)}\right)$$

$$L.U.R.(\chi )=\frac{1}{n-d}\sum _{k=0}^{n-d-1}\left(\frac{\left|\frac{{\chi }_k}{\Vert x\Vert }-\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }\right|}{\min \left(\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert },\frac{{\chi }_k}{\Vert x\Vert }\right)}\right)$$

$$U.R.(r)=\frac{1}{n-d}\sum _{k=0}^{n-d-1}\left(\frac{\left|\frac{\Vert r_k{\Vert }^2}{\Vert b{\Vert }^2}-\frac{\Vert {ϵ}_k{\Vert }^2}{\Vert x{\Vert }^2}\right|}{\min \left(\frac{\Vert {ϵ}_k{\Vert }^2}{\Vert x{\Vert }^2},\frac{\Vert r_k{\Vert }^2}{\Vert b{\Vert }^2}\right)}\right)$$

$$L.U.R.(r)=\frac{1}{n-d}\sum _{k=0}^{n-d-1}\left(\frac{\left|\frac{\Vert r_k\Vert }{\Vert b\Vert }-\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert }\right|}{\min \left(\frac{\Vert {ϵ}_k\Vert }{\Vert x\Vert },\frac{\Vert r_k\Vert }{\Vert b\Vert }\right)}\right)$$

Full text

AMS Codes: 65F10; 65G99.

1 Introduction

Solving a system of linear equations in the form $Ax=b$ is a ubiquitous requirement in science and engineering (where $A$ is a given matrix, $x$ and $b$ are the unknown and known vectors respectively; $x$ $\in\mathbb{R}^{n}$ and $b$ $\in\mathbb{R}^{n}$ if $A$ $\in\mathbb{R}^{n\times n}$). Iterative methods like Conjugate Gradient (CG) [7], Bi-Conjugate Gradient (Bi-CG) [4], and Generalized Minimal Residual (GMRES) [14], mostly known as Krylov subspace methods, are commonly used to solve large linear systems. While using such methods, iterations are to be stopped when the norm of the error ($\epsilon_{k}=x-x_{k}$) is less than a desired tolerance, where $x$ is the unknown solution to the linear system and $x_{k}$ is its iterative approximation at the $k^{th}$ iteration. Since the actual error is unknown, relative residue ($\frac{\lVert r_{k}\rVert}{\lVert b\rVert}$) is considered as a stopping criteria where $r_{k}=b-Ax_{k}$ is the residual vector at $k^{th}$ iteration. Such stopping criteria are useful when the system is well-conditioned but they can become unreliable for even a moderate increase in the relative condition number of $A$ ($\kappa\gtrsim 10$, where $\kappa=\mathinner{\!\left\lVert A\right\rVert}_{2}\mathinner{\!\left\lVert A^{-1}\right\rVert}_{2}$), depending on the choice of initial approximation. It is well known that relying on relative residue may lead one to stop the iterations too early when the norm of error is still much larger than tolerance, or not stop early enough with too many floating point operations done for the required accuracy ([12]; [7]). Moreover, when condition number of the matrix is large, the residual of algorithms like Bi-CG need not show monotonic behaviour and oscillate [1] while the actual error might still be (however slowly) converging, and similarly, oscillations of errors in the GMRES algorithm are possible when the residue converges slowly. Note that the $l_{2}$ norm of the relative errors and residues can differ significantly as given below:

[TABLE]

Even when the iterative algorithms are used with preconditioners that reduce the condition number of the matrix, a reduction of the condition number of the backward problem to be solved is not guaranteed, as required for the faster convergence (see section 3.3.1). This is observed with matrices of larger dimensions or higher condition numbers, making the preconditioners ineffective in many cases. Moreover, most of the preconditioning techniques often fail when matrices are non-symmetric and indefinite [13]. Thus, for even $\kappa(A)\gtrsim 10^{2}$, either the accuracy or the efficiency of computation is degraded by the above conundrum. The requirements of precision in measurements and engineering today renders both the size and condition number of most problems large; making accurate stopping and restarting criteria indispensable in ensuring computational efficiency of solvers. This motivated methods to estimate norms of the error in iterative solvers. Approximations of the errors in $x_{k}$ for an iterative method were suggested in general more than two decades ago [2] [3]. Error estimators for the CG algorithm were also presented [11], and for solving non-symmetric linear systems using Full-Orthogonalization (FOM) and GMRES methods, formulas for estimation of errors have been suggested more recently [10]. We suggest a minor modification to this GMRES estimator to increase its stability and precision. We also derive an efficient estimator of error in solving non-symmetric linear systems using Bi-CG.

In the second part of this work, we present an analysis highlighting the significance of these error estimators for CG, Bi-CG and GMRES methods. We confirm that there is indeed less certainty in relying on relative residual as an estimator of the true relative error, in comparison to using an error estimator. We quantify this ‘uncertainty’ and show that its expected value is proportional to condition numbers i.e. forward or backward condition numbers or their weighted averages depending on the algorithm used. Numerical results averaged over a large number of linear systems were used to demonstrate the theorem on this uncertainty in convergence. The under and over computation due to the uncertainty is demonstrated, and the gains due to the estimators are relatively straight-forward to infer. Our analysis shows that these error estimators are robust and increase the efficiency/accuracy of computing notably. Section 2 briefly introduces related work and the previous research in the estimation of error in CG and GMRES. Section 3 presents the modification to the GMRES estimator, the proposed estimates of error in Bi-CG, and the analysis of error estimates. Section 4 discusses the experimental results on the modified GMRES error estimates, proposed BiCG error estimates, and the performance of error estimates as an efficient stopping criteria for these iterative methods of solving linear systems.

2 Related work

2.1 Error estimators for CG algorithm

One of the most commonly used methods for solving linear systems with a real symmetric positive definite matrix is the Conjugate Gradient (CG) algorithm. CGQL (CG-Qaudrature Lanczos) algorithm [6, 11, 12] provides bounds of A-norm of error i.e $\mathinner{\!\left\lVert x-x_{k}\right\rVert}_{A}$ at each iteration of CG by using Lanczos and quadrature rules with a delay of $d$ iterations in the estimation. Bounds to the square of the A-norm of error at iteration $k$ of CG are given by [6]:

[TABLE]

Here, $T_{k}$ is the tridiagonal matrix produced by the Lanczos algorithm at the $k^{th}$ iteration whose entries can be computed from the coefficients of CG algorithm equally well. $(T_{k}^{-1})_{(1,1)}$ represents the $(1,1)$ element of the matrix $T_{k}^{-1}$. Also $\epsilon$ and $r$ are the error and residual vectors of the CG algorithm respectively. The essence of CGQL lies in computing difference between $(1,1)$ elements of the inverse of two tridiagonal matrices generated from a Lanczos algorithm with the same starting vectors as the CG algorithm. However, from a computational point of view, other equivalent formulas for approximation of error in CG were proposed [8], including a $LDL^{T}$ factorization of the corresponding Jacobi matrix available in CG [12]. Estimates of the $l_{2}$ norm of the error in the CG algorithm were later proposed [9], including extensions to preconditioned versions [18]. It should also be noted that with a delay in estimation $d$, where $d>1$, considering $\mathinner{\!\left\lVert\epsilon_{k+d}\right\rVert}_{A}\ll\mathinner{\!\left\lVert\epsilon_{k}\right\rVert}_{A}$, the approximation of the A-norm of error at iteration $k$ is also given using CG coefficients $\alpha_{j}$, as:

[TABLE]

2.2 Error estimator for GMRES algorithm

Generalized minimal residual (GMRES) algorithm is an iterative method for solving linear systems even for a general non-symmetric matrix [14]. Let $V_{k}$ be a matrix whose columns are orthonormal basis vectors $v_{j},\ j=1\dots k$, of Krylov subspace $\mathcal{K}_{k}(A,r_{0})$, where $r_{0}$ is initial residual. The iterates of GMRES are defined as $x_{k}=x_{0}+V_{k}z_{k}$ where vector $z_{k}$ has weights for the orthonormal basis vectors $v_{j}$ in the matrix $V_{k}$. Maintaining the notations [10] for estimates of norm of error in FOM and GMRES at the $k^{th}$ iteration of GMRES, the Hessenberg matrix $H_{k}$ resulting from the iterative orthogonal transformation of the system matrix can be decomposed block-wise as:

[TABLE]

Here, the matrices $H_{k-d}$ and $\tilde{H}_{k-d}$ are square and upper Hessenberg and $W_{k-d}$ is generally a full rectangular matrix. Let $h_{i,j}$ be the $(i,j)$ element of the Hessenberg matrix, $e_{k}$ is the $k^{th}$ column of an identity matrix, the vector $t_{k}$ be the last column of \big{(}H_{k}^{T}H_{k}\big{)}^{-1} with $t_{k,k}$ as its last element, then variables $\gamma_{k-d},w_{k-d},$ $\delta_{k+1}$, and $u_{k}$ can be defined as the following,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The error estimate $\epsilon_{k-d}$ at $(k-d)^{th}$ iteration of GMRES can be aproximated by [10]:

[TABLE]

Note that the above estimation involves only one dense matrix-vector multiplication and requires $\mathcal{O}(k^{2})$ arithmetic operations every iteration, and this is negligible for all $k\ll n$. For non-converging problems, as $k\longrightarrow n$, the total additional cost due to the estimator may approach a third of the cost of the GMRES iterations, but this is likely to happen only rarely in applications with a large $n$, and thus the estimator can provide significant benefits on the average.

3 Methods

3.1 Proposed modification to the GMRES error estimator

The above estimate of error in GMRES can also be rewritten as term-1 in the following:

[TABLE]

where $s_{k}=(e_{k},H_{k}^{-1}e_{1})u_{k}$ and

[TABLE]

where superscripts $k$ and $k-d$ have been introduced for the vector components, indicating the appropriate truncated dimensions.

The R.H.S of (4) was divided into three terms where the third term given by $\zeta_{k}$ is assumed to be negligible. Also, $H_{n}^{-1}$ required in the evaluation of $\zeta_{k}$ cannot be computed at the $k^{th}$ iteration. Note that for a large $d$ which is a delay in estimation, the second term $\mathinner{\!\left\lVert s_{k}\right\rVert}^{2}$ also becomes negligible. However $d$ is kept small in practice ($\approx$ 10), hence in such cases, this second term can significantly contribute to the value of error estimate. We propose to include this second term in the estimates of the error, and here, an absolute operation in (5) is necessary to preserve the positive estimate of the $l_{2}$ norm of error using the truncated expression. Including $s_{k}$ in (3) and considering the absolute values,

[TABLE]

.

3.2 Proposed estimator for A-norm and $l_{2}$ norm of errors in Bi-Conjugate gradient Algorithm (BiCG)

Whereas A-norm of error in CG is minimized monotonically with the number of iterations, A-norm of error in BiCG can oscillate and might not even be defined when A is an indefinite matrix. Hence, we consider the absolute values of the bi-linear products of the residues of BiCG as shown below.

[TABLE]

Note than when $A$ $\in\mathbb{R}^{N\times N}$ and $r$ $\in\mathbb{R}^{N}$ ; $r_{k}^{T}(A^{T})^{-1}r_{k}$ is a scalar quantity that is always positive when $A$ is positive definite. Moreover the square of $l_{2}$ norm of error is given by:

[TABLE]

If $A=A^{T}$, A-norm and $l_{2}$ norm of error are given by $r_{k}^{T}A^{-1}r_{k}$ and $r_{k}^{T}A^{-2}r_{k}$ respectively. In the following sections, estimates of A-norm and $l_{2}$ norm of error at every iteration of BiCG are derived. BiCG is included as Algorithm 1 in the appendix for reference.

3.2.1 Estimation of A-norm of error

By writing the consecutive difference between A-norm of error at iteration $k$ and $k+1$, the following relation is obtained when A is a non-symmetric matrix (See Appendix for BiCG Algorithm):

[TABLE]

In the above equation, the first and third terms of R.H.S can be trivially computed using iterates of the Bi-CG Algorithm. Second term involves computation of $A^{-1}$ hence we further reduce it. As $r_{k+1}=r_{k}-\alpha_{k}Ap_{k}$ and $Ap_{k}=\dfrac{r_{k}-r_{k+1}}{\alpha_{k}}$, the following relation is obtained:

[TABLE]

Error ($\epsilon_{k}=x-x_{k}$) is given by $A^{-1}r_{k}$. Also, $\epsilon_{k}$ can be written as a weighted sum of search directions from $k$ to $n$ given below in (8).

[TABLE]

If $\epsilon_{k+d}$ denotes the error at $k+d$ iteration and $\epsilon_{k+d}\ll\epsilon_{k}$ for some $d>0$, then further terms of the series can be neglected. The error vector in (8) can provide an estimate of A-norm of error after inducing a delay of $d$ iterations:

[TABLE]

When A is positive definite, the above expression is always positive and thus provides a lower bound for the square of A-norm of error. Equation (9) is derived directly from Bi-CG and when $A=A^{T}$, this estimator for $A$-norm is identical to the earlier proposed CGQL $A$-norm estimator (B). Since the matrix-vector product $Ap_{k}$ is known in terms of the residual vectors, the above evaluation costs only $2n$ arithmetic operations, compared to the $\mathcal{O}(n^{2})$ arithmetic operations in each Bi-CG iteration. The summation of $d$ terms in the above implies a latency of $d$ iterations to begin with, when an estimate of the error is not available.

3.2.2 Estimation of $l_{2}$ norm of error

By writing the consecutive difference between $l_{2}$ norm of error at iteration $k$ and $k+1$ of Bi-CG, the following relation is obtained when A is a non-symmetric matrix:

[TABLE]

Rearranging (10),

[TABLE]

Using (8) and $\epsilon_{k}=A^{-1}r_{k}$, the above equation can be rewritten as:

[TABLE]

The evaluation involves only $2n$ arithmetic operations every iteration, and for the first $d$ iterations an estimate of the error is not available. This estimation adds a negligible computational overhead compared to the $\mathcal{O}(n^{2})$ cost of every Bi-CG iteration, especially as $n$ increases.

3.3 Analysis of estimators

3.3.1 Condition number of the problem

Condition numbers of problems enable us to estimate the expected accuracy of computed result. In this case, the relative condition number for computing $b$ given $A$ and $x$ (forward problem), and condition number for computing $x$ from $A$ and $b$ (backward problem) are given by (13) and (14) respectively [19]. They help us infer the relative change in the output for a relative change in inputs of the problem representing round-offs, thus constraining the accuracy of the output. The norm ($\mathinner{\!\left\lVert.\right\rVert}$), unless specified, is considered to be the $l_{2}$ norm.

[TABLE]

[TABLE]

$\kappa_{f}(A,x)$ and $\kappa_{b}(A,b)$ are fixed for a given problem, $Ax=b$. Condition number of the matrix is given by the product of $\kappa_{f}(A,x)$ and $\kappa_{b}(A,b)$,

[TABLE]

3.3.2 Evaluation metrics

The following is the relative error in our estimates of the relative error at the $k^{th}$ iteration:

[TABLE]

where $\chi_{k}$ is the estimate of norm of error at $k^{th}$ iteration and $x$ is the solution to the problem $Ax=b$. The maximum of the two terms inside the brackets encapsulate the fact that an error estimate can either be greater than or less than the true relative error at any iteration, and it can effect a corresponding degree of under-computation or an over-computation. Equation (15) can also be concisely rewritten as,

[TABLE]

While (16) quantifies the uncertainty using the relative residue or an error estimator at a particular iteration, Uncertainty Ratios and Linear Uncertainty Ratios i.e. $U.R.(\chi)$ and $L.U.R.(\chi)$, are defined to represent average uncertainty in estimating true relative error over $n-d$ iterations, as follows:

[TABLE]

[TABLE]

where $d$ is the delay in estimation in terms of iterations. When the relative residual is used for stopping the iterative solution of a linear system, given a tolerance in relative error, the uncertainty is correspondingly:

[TABLE]

[TABLE]

In our subsequent analysis, we present a theorem on an upperbound of the expectation of $U.R.(r)$ and $L.U.R.(r)$ over converging problems having fixed singular values and constant forward condition number, we show that $\operatorname{\mathbb{E}}(U.R.(\chi))\leq\operatorname{\mathbb{E}}(U.R.(r))$, and conclude that use of an error estimator in such cases as a stopping criterion gives increased savings and accuracy in computation on average, owing to their robustness.

3.3.3 Theorem on Expectation of $U.R.$ and $L.U.R.$

Theorem 1**.**

For linear problems of non-singular real matrices A with a set of fixed singular values, fixed forward condition number $\kappa_{f}(A,x)>>1$, fixed backward condition number $\kappa_{b}(A,b)>>1$,

[TABLE]

&

[TABLE]

where $n$ is the dimension of the matrix $A$ and $\alpha$ is the average probability of relative residual over predicting the relative error of solution i.e. $\alpha=\dfrac{1}{n-d}\sum_{k=0}^{n-d-1}P_{cond}\left(\dfrac{\lVert r_{k}\rVert}{\lVert b\rVert}\geq\dfrac{\lVert\epsilon_{k}\rVert}{\lVert x\rVert}\right)$, where $P_{cond}$ refers to conditional probability with conditions stated at the start of the theorem. The sample space for conditional probability and conditional expectation is generated by sampling left-singular matrix $U$ and right-singular matrix $V$ from uniformly distributed orthogonal matrices and unit vector $\hat{b}$ from Gaussian distribution.

Proof.

Let the uncertainty ratio at $k^{th}$ iteration be

[TABLE]

and

[TABLE]

Using the definition in (15) for $U.R.^{(k)}(r)$, we can write,

[TABLE]

where $\operatorname{\mathbb{E}}_{1}$ refers to \operatorname{\mathbb{E}}\left(\bullet\Biggr{|}\dfrac{\lVert r_{k}\rVert}{\lVert b\rVert}\geq\dfrac{\lVert\epsilon_{k}\rVert}{\lVert x\rVert}\right) and $\operatorname{\mathbb{E}}_{2}$ refers to \operatorname{\mathbb{E}}\left(\bullet\Biggr{|}\dfrac{\lVert r_{k}\rVert}{\lVert b\rVert}\leq\dfrac{\lVert\epsilon_{k}\rVert}{\lVert x\rVert}\right).

[TABLE]

[TABLE]

Consider $S_{k}=\dfrac{\lVert r_{k}\rVert}{\lVert\epsilon_{k}\rVert}\dfrac{\lVert x\rVert}{\lVert b\rVert}$, thus, $S_{k}\in\left\{\dfrac{1}{\kappa_{b}(A,b)},\kappa_{f}(A,x)\right\}$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Based on fixed conditions of expectation, $\dfrac{\lVert x\rVert}{\lVert b\rVert}$ is constant,

[TABLE]

By Using definition of $l_{2}$ norm of matrix,

[TABLE]

Similarly,

[TABLE]

Hence, using above bounds and $\alpha$ as defined in theorem 1,

[TABLE]

[TABLE]

Similarly, we can derive the following bounds for $L.U.R(r)$:

[TABLE]

∎

3.3.4 Note on Delay-based Error Estimators

For delay-based estimators, $\chi_{k}^{2}\approx\mathinner{\!\left\lVert\epsilon_{k}\right\rVert}^{2}-\mathinner{\!\left\lVert\epsilon_{k+d}\right\rVert}^{2}$, and

[TABLE]

On the other hand, the relative residual has a strict lower bound of $\dfrac{1}{\kappa(A)}\dfrac{\mathinner{\!\left\lVert\epsilon\right\rVert}}{\mathinner{\!\left\lVert x\right\rVert}}$ and a strict upper bound of $\kappa(A)\dfrac{\mathinner{\!\left\lVert\epsilon\right\rVert}}{\mathinner{\!\left\lVert x\right\rVert}}$, for a given matrix $A$. this allowed us to provide conditional upper and lower bounds when $\mathinner{\!\left\lVert\epsilon_{k+d}\right\rVert}^{2}<\mathinner{\!\left\lVert\epsilon_{k}\right\rVert}^{2}$. Hence, use of delay based error estimators is reliable when an iterative method like GMRES or BiCG has a decreasing $l_{2}$ norm of error for most iterations. It can be argued that value of $d$ can be increased to an extent such that $\mathinner{\!\left\lVert\epsilon_{k+d}\right\rVert}^{2}<\mathinner{\!\left\lVert\epsilon_{k}\right\rVert}^{2}$ for all the iterations, but in practice, a constraint ($d\ll n$) is required due to low latency requirements. The above conclusions are valid without any restriction on the dimension of the matrices.

4 Numerical results

4.1 Modified $l_{2}$ norm error estimator for GMRES

In this experiment, we empirically highlight the advantages of including term-2 in the approximation of the GMRES error, as proposed in equation (5). We conducted experiments on non-symmetric positive definite high condition number matrices (generated via eigen decomposition with eigen values spread exponentially between $1$ and $\kappa(A)$) in which we compared the values of different error estimators and relative residue on each iteration.

In Figure 1, we compare relative $l_{2}$ norm of error estimator proposed by Meurant [10], the modified $l_{2}$ norm error estimator proposed in equation (5), the relative residual and the true relative error at each iteration. We observe that including term-2 in equation (5) helped in avoiding unstable overshoots in estimation for some of the iterations as seen in Figure 1.

It is also observed that GMRES error estimators may behave erratically when numerical precision of computing system is exhausted as seen in Figure 2. The exhaustion of numerical precision can lead to near singularity of Hessenberg matrix $H_{k}$ formed during the Arnoldi iteration which can cause large errors in $H_{k}^{-1}$ and its functions. Note that $\mathinner{\!\left\lVert s_{k}\right\rVert}$ is also a function of $H_{k}^{-1}$ and can be used as trigger to predict this exhaustion of numerical precision as seen in Figure 2 making this error estimate a robust stopping criterion.

4.2 A-norm and $l_{2}$ norm error estimators for Bi-CG

In this experiment, we want to empirically evaluate the behaviour of delay-based BiCG error estimators derived in section 3.2 when compared to relative residue, in problems with a high condition number. We considered non-symmetric matrices from the Matrix Market111https://math.nist.gov/MatrixMarket/ dataset (’sherman5’) with a random right hand side vector $b$.

In Figure 4, we compare estimated relative A-norm of the error, the relative residual, and the true relative A-norm of the error at each iteration. Figure 3 shows the comparison between estimated relative $l_{2}$ norm error, the true relative $l_{2}$ norm of the error, and the relative residual. We observe that the proposed BiCG error estimators performs consistently across iterations in both cases, where the error estimator is reliable as an indicator of the actual error for stopping.

4.3 Experimental results on performance of error estimators on randomized set of matrices

The objective of the experiments222Code freely available at https://github.com/CSPL-IISc/Error-Linsolve is to asses the predictive performance of error estimators on large randomized set of matrices. This study is useful in understanding the variability in average behaviour of estimators with respect to different aspects of problem. In previous sections, we had defined the Uncertainty Ratio as the average error in any of the predictors of true relative error that may be used for stopping.

In this experiment, we measure the Uncertainty Ratio of estimator and residual for a given forward condition number and backward condition number for each problem in randomized set. The randomized dataset was chosen to consist of 10,000 problems of high condition numbers of dimension $100$ on which the evaluation metrics $L.U.R(r)$ and $L.U.R(\chi)$ are measured with $d=10$ for delay in estimation. We consider two kinds of non-symmetric matrices i.e. positive definite and general matrices in the dataset. The right hand side vectors $b$ are random vectors from a unit normal distribution of same dimension.

Figure 5 and 6 show the behaviour of Linear Uncertainty Ratio ($L.U.R.$), defined by (18) and (20), with the condition number of the problem (both forward and backward) for BiCG and GMRES algorithms respectively. Here each circle (blue for BiCG and red for GMRES) represents the uncertainty ratio corresponding to the relative residual for a particular problem i.e for a particular matrix $A$ and a right hand side vector $b$ with $x_{0}$ sampled from standard Gaussian distribution. Similarly, each (black) cross represents the uncertainty ratio corresponding to the relative error estimator.

We observe that there is a large increase in L.U.R.(r) of relative residual in high condition number problems, whereas the L.U.R.($\chi$) for the relative error estimator remains bounded, indicating the robustness of BiCG and GMRES error estimators. It is also interesting to note that Uncertainty Ratio is dependent on both forward and backward condition number of the problem. For BiCG, the $L.U.R.(r)$ increases with the forward condition number of the problem, while $L.U.R.(r)$ increases with the backward condition number in case of the GMRES algorithm. Note that the $\alpha$ values in the result of the theorem, that determine the increasing (or decreasing) relationship with the forward (or backward) condition number, can vary across algorithms. Since the same dataset was used for testing the metrics on both Bi-CG and GMRES algorithms, the scales on x-axis differ in the two figures (and note $\kappa_{f}(A,x)*\kappa_{b}(A,b)=\kappa$). While Bi-CG has much larger uncertainty ratios because of its oscillatory covergence behavior, the expected under and over computation for a required tolerance in error are comparable for both Bi-CG and GMRES, as shown in the next section.

For further demonstration of the above conclusions even in matrices of practical applications, non-symmetric matrices of varied condition numbers and dimensions were collected from Matrix Market and Matrix Sparse Collection data set and were tested with corresponding random normalized right hand side vectors. Figure 7 shows the relation between Uncertainty Ratio and condition numbers (forward and backward) of the corresponding problems. We observe the significant increase in the uncertainty ratio of the relative residual with the increase in the condition number of the problem while the uncertainty ratio of the error estimator (for both BiCG and GMRES) stays almost stagnant as the condition number of the problem increases, depicting error estimators as a robust stopping criteria. Though we mostly focus on the evaluation metrics for BiCG and GMRES algorithms, one should keep into account that the performance of the CG estimator is implicitly described by the Bi-CG estimator in the case of symmetric positive definite matrices.

4.4 Under and over computation with BiCG and GMRES estimators

While using iterative algorithms, one would like to stop the computation at error tolerances depending on the requirements of the applications. Since there is a larger uncertainty when stopping using the relative residue, as shown in Figures 5, 6 and 7, an under or over computation is expected for any desired error tolerance in solution. One might stop the iterations too early (under computation) or deploy too many iterations (over computation) due to this uncertainty in error of the solution, resulting in loss of efficiency or accuracy. In this experiment, the objective is to evaluate the degree of over or under computation due to the uncertainty in error when using the relative residual as an indicator for stopping.

We define the relative iterations as $I$ which are either saved or lost for a desired error tolerance as:

[TABLE]

For any tolerance such as $10^{-4}$, $I$ signifies the fraction of iterations over or under computed while using relative residual as a stopping criteria as compared to while using error estimator as stopping criteria. $Iter_{\frac{\mathinner{\!\left\lVert r\right\rVert}}{\mathinner{\!\left\lVert b\right\rVert}}}$ is the number of iterations required by the relative residual ($\frac{\mathinner{\!\left\lVert r\right\rVert}}{\mathinner{\!\left\lVert b\right\rVert}}$) to reach the desired tolerance, and we have a similar variable for the iterations required by the error estimate. Figure 8 shows the plot of $I$ at error tolerances of $10^{-1},10^{-2},10^{-3}$ and $10^{-4}$ with increasing $L.U.R.(r)$ averaged over the matrix problems used in Section 4.3 for both the algorithms. It should be noted that the value of $I=1$ means there is a 100 percent difference in the iterations required by $\frac{\mathinner{\!\left\lVert r\right\rVert}}{\mathinner{\!\left\lVert b\right\rVert}}$ and $\frac{\mathinner{\!\left\lVert\chi\right\rVert}}{\mathinner{\!\left\lVert x\right\rVert}}$ to reach a desired tolerance. We observe the increase in $I$ with the uncertainty ratio of relative residual which signifies the importance of error estimators while using iterative algorithms.

5 Conclusions

The significance of error estimators for efficient stopping (or restarting) are clearly evident for problems with condition numbers $\kappa\gtrsim 100$, and is emphasized by numerical examples and an analysis of the expected uncertainty in convergence when using the norm of the residual vector for stopping, for any given tolerance in error. Results showing the reduction of uncertainty ratio while using error estimators were presented for both Bi-CG and GMRES algorithms. It was also highlighted that the proposed estimator for Bi-CG converges to the earlier proposed estimator for the CG algorithm in the case of symmetric positive-definite matrices. Further, a simple analysis of the change in computation due to the estimators was presented along with numerical results of large number of problems of manageable dimensions. While Bi-CG has much larger uncertainty ratios because of its oscillatory behavior (and this can be improved by using extended versions of the algorithm [4] [17] [16] [20] [5]), the expected under and over computation for a required tolerance in error are comparable for both Bi-CG and GMRES, as expected. Based on the results discussed in the previous sections, we conclude that the estimate for the A-norm or the $l_{2}$ norm of the error should be implemented into software realization of iterative solvers, instead of using only the relative norm of the residual vector as a criterion for stopping.

Author declarations:

No funding was received for conducting this study. The authors have no relevant financial or non-financial interests to disclose.

Appendix A Bi-CG and its relation to non-symmetric Lanczos algorithm

Bi-Conjugate Gradient (Bi-CG) algorithm is an extension to the Conjugate Gradient (CG) algorithm which is used to solve a system of linear equations and it is suited for a non-symmetric matrix.

Bi-CG can also be derived from the non-symmetric Lanczos algorithm. Let $v_{1}$ and $\tilde{v_{1}}$ be the given starting vectors to the non-symmetric Lancozs algorithm (such that $\lVert v_{1}\rVert=1$ and $(v_{1},\tilde{v_{1}})=1)$, the two three term recurrences which help in forming two bi-orthogonal subspaces can be as follows:

For k=1,2 $\dots$

[TABLE]

The coefficient $w_{k}$ being computed as $w_{k}=(\tilde{v_{k}},Av_{k})$. The other coefficients $\eta_{k}$ and $\tilde{\eta_{k}}$ are chosen (provided ($\tilde{z_{k}},v_{k})=0$) such that $\eta_{k}\tilde{\eta_{k}}=(\tilde{z_{k}},z_{k})$ and the new vectors at step $k+1$ are given by:

[TABLE]

The relationship between the bi-orthogonal subspaces and matrix $A$ can be written in the form of a non-symmetric tri-diagonal matrix form (under the condition $\tilde{V_{k}}^{T}AV_{k}=T_{k}$) as:

[TABLE]

Appendix B Relating BiCG error estimator to CGQL

The difference between two consecutive A-norm of error (and A-measure for matrices that are not SPD) in case of a Conjugate Gradient algorithm at iteration $k$ and $k+1$ can be given by:

[TABLE]

Hence by inducing a delay of $d$ iterations we can easily compute A-norm of error.

[TABLE]

For A-norm estimation of error in BiCGQL algorithm we derived the following results:

[TABLE]

By using $r_{k+1}=r_{k}-\alpha_{k}Ap_{k}$ the above result can also be written as:

[TABLE]

For a symmetric matrix $A=A^{T}$ the above equation can be further written as:

[TABLE]

Further, using $p_{i}^{T}r_{j}=0$ for $i\neq j$ in CG, and also $r_{k+1}=r_{k}-\alpha_{k}Ap_{k}$ and thus substituting $Ap_{k}=\dfrac{r_{k}-r_{k+1}}{\alpha_{k}}$, we reduce the above as:

[TABLE]

Equation 21 is equivalent to the $A$-norm estimation in CGQL Algorithm [6], and thus the proposed $A$-norm estimator for BiCG is equivalent to the CGQL estimator when $A=A^{T}$.

Appendix C Relations between non-symmetric Lanczos tridiagonal Matrix ($T_{k}$) and residual vectors ($r$ and $\tilde{r}$) of BiCG algorithm

A direct relationship between $T_{k}$ , $r_{k}$ and $\tilde{r_{k}}$ can be given by:

[TABLE]

[TABLE]

Subtracting 22 and 23 we get:

[TABLE]

As $\tilde{r_{k}}^{T}A^{-1}r_{k}-\tilde{r}_{k+1}^{T}A^{-1}r_{k+1}=\alpha_{k}\tilde{r}_{k}^{T}r_{k}$ for a non-symmetric matrix in BiCG:

[TABLE]

Thus knowing the relation between two consecutive $T_{k}$ inverse first elements we can relate it with the estimator we developed for A-norm of error and the approach of CGQL Algorithm. However relating Non-Symmetric Lanczos and BiCG through Quadrature based methods will involve Complex Gaussian Quadratures [15]. Hence we have followed an equivalent but direct approach of estimation using the relations of Bi-CG explicitly.