Error Analysis of Krylov Subspace approximation Based on IDR($s$) Method for Matrix Function Bilinear Forms

The bilinear form of a matrix function, namely $\mathbf{u}^\top f(A) \mathbf{v}$, appears in many scientific computing problems, where $\mathbf{u}, \mathbf{v} \in \mathbb{R}^n$, $A \in \mathbb{R}^{n \times n}$, and $f(z)$ is a given analytic function. The Induced Dimension Reduction IDR($s$) method was originally proposed to solve a large-scale linear system, and effectively reduces the complexity and storage requirement by dimension reduction techniques while maintaining the numerical stability of the algorithm. In fact, the IDR($s$) method can generate an interesting Hessenberg decomposition, our study just applies this fact to establish the numerical algorithm and a posteriori error estimate for the bilinear form of a matrix function $\mathbf{u}^{\top} f(A) \mathbf{v}$. Through the error analysis of the IDR($s$) algorithm, the corresponding error expansion is derived, and it is verified that the leading term of the error expansion serves as a reliable posteriori error estimate. Based on this, in this paper, a corresponding stopping criterion is proposed. Numerical examples are reported to support our theoretical findings and show the utility of our proposed method and its stopping criterion over the traditional Arnoldi-based method.

💡 Research Summary

This paper presents a novel and efficient numerical method for approximating the bilinear form of a matrix function, u^T f(A) v, a computation that arises in numerous scientific and engineering applications such as option pricing, network analysis, and computational electromagnetics. The key innovation lies in the adaptation of the Induced Dimension Reduction (IDR(s)) method, originally designed for solving large-scale linear systems, to this function evaluation problem.

The authors leverage a crucial property of the IDR(s) process: its ability to generate a standard Hessenberg decomposition, A V_m = V_{m+1} \bar{H}_m, where V_m has orthonormal columns and \bar{H}_m is an upper Hessenberg matrix. This decomposition is the cornerstone of Krylov subspace methods. Using this relationship, the high-dimensional bilinear form is projected onto the Krylov subspace, yielding a low-dimensional approximation F_m = β u^T V_m f(H_m) e_1, where H_m is a section of \bar{H}_m, β = ||v||, and v_1 = v/β.

The central theoretical contribution of the work is a detailed error analysis for this IDR(s)-based approximation. Under the assumption that f is analytic on a domain containing the numerical ranges of A and H_m, Theorem 3.1 derives an infinite series expansion for the error E_m(f) = u^T f(A)v - F_m. The most significant insight from this analysis is that the leading term of this error expansion can be computed using quantities already available during the iterative IDR(s) process, specifically involving the next basis vector v_{m+1} and the Hessenberg matrix H_m. This term is proven to serve as a reliable a posteriori error estimate.

Building upon this theoretical foundation, the paper proposes a practical and efficient stopping criterion for the iterative algorithm: the iteration can be terminated once the norm of this readily available error estimate falls below a user-defined tolerance. This criterion prevents unnecessary iterations, optimizing computational resources.

Numerical experiments in Section 4 validate the theoretical findings. Tests are conducted with various functions (exponential, sine, cosine) and matrices (symmetric, nonsymmetric, sparse). The results demonstrate that the proposed a posteriori error estimate accurately tracks the true error. Furthermore, the stopping criterion performs effectively. Comparisons with the traditional Arnoldi-based method show that the IDR(s) approach can achieve comparable accuracy while offering potential advantages in terms of computational cost and storage requirements, thanks to the dimension-reduction properties of IDR(s). In conclusion, this work successfully extends the IDR(s) methodology to matrix function bilinear forms, providing a rigorous error analysis, a computationally efficient error estimator, and a practical stopping criterion, thereby offering a valuable new tool for large-scale scientific computing.