The extended global Lanczos method, Gauss–Radau quadrature, and matrix function approximation
Introduction
Let be a large symmetric matrix, a function defined on the convex hull of the spectrum of , and a vector of unit Euclidean norm. Application of steps of the (standard) Lanczos method to with initial vector yields the partial Lanczos decomposition where the matrix has orthonormal columns with initial column , is a symmetric tridiagonal matrix, and satisfies ; see, e.g., [1], [2] for details on the Lanczos algorithm. Throughout this paper denotes the th column of the identity matrix of suitable order, and the superscript stands for transposition. Generally, in applications of the Lanczos method. We assume that is small enough so that the Lanczos decomposition (1.1) with the stated properties exists. This is the generic situation. Then where the right-hand side denotes the Krylov subspace
Introduce the spectral factorization with an orthogonal matrix . Substituting the spectral factorization into gives where . The right-hand side of (1.5) can be written as a Stieltjes integral determined by a piece-wise constant non-decreasing distribution function with jumps at the . Let denote the measure associated with . Then We note for future reference that the measure defines an inner product for polynomials and of low enough degree.
Golub and Meurant [3], [4] observed that the nontrivial entries of the symmetric tridiagonal matrix in (1.1) are recursion coefficients for a sequence of orthonormal polynomials with respect to the inner product (1.7). In particular, the eigenvalues of are the zeros of the orthonormal polynomial of degree . This led Golub and Meurant [3], [4] to show that the expression is the -point Gauss quadrature rule for approximating the Stieltjes integral (1.6). Hence, where denotes the set of polynomials of degree at most . Substituting the spectral factorization of into (1.8) transforms the right-hand side into a sum of terms and shows that the eigenvalues of are the nodes and the square of the first component of normalized eigenvectors are the weights of the quadrature rule. This observation forms the basis for the Golub–Welsch algorithm [5] for the efficient computation of nodes and weights of a Gauss quadrature rule from . However, we remark that it may not be necessary to first evaluate (1.8) by determining the nodes and weights of the Gauss rule and then compute at the nodes. For some functions , it may be more efficient to compute directly by one of the methods discussed in [6]; see [7] for an illustration.
The quadrature rule (1.8) is obtained by applying steps of the (standard) Lanczos method to with initial vector . An analogous rule can be computed by applying the global Lanczos method to with initial block vector with block size . The global Lanczos method is a block Lanczos method with a particular inner product. It was first proposed and investigated by Elbouyahyaoui et al. [8] and Jbilou et al. [9]. The global block Lanczos method uses the inner product between block vectors and the induced Frobenius norm This method can be applied to approximate expressions of the form where with . The problem of estimating the trace of a matrix , without evaluating , is a classical problem in numerical linear algebra; see, e.g., [10], [11], [12], [13], [14], [15], [16], [17], [18]. The need to evaluate or approximate expressions of the type (1.11) arises in various applications, including in network analysis and when solving ill-posed problem; see, e.g., [13], [17], [19], [20], [21], [22] for discussions of these applications. An approximation of the trace of can be computed by approximating expressions (1.11) for several block vectors , , where the are block axis vectors. For notational simplicity, we here assume that is a multiple of ; see [13] for details on the computations.
The method described in this paper for approximating expressions of the form (1.11) is particularly well suited for problems with a matrix that allows efficient solution of linear systems of equations with this matrix. This includes semiseparable matrices [23] and, in particular, symmetric positive definite Toeplitz matrices [24].
This paper is organized as follows. Section 2 reviews the global Lanczos method. This method gives approximations of integrals that are exact when the integrand is a polynomial of low enough degree, analogously to (1.9). However, the error in the computed approximation may be large if the integrand cannot be well approximated by a polynomial of moderate degree. The extended global Lanczos method determines approximations of integrals that are exact when the integrand is a Laurent polynomial of low enough order. The method described in this paper may determine approximations of an integral with much higher accuracy than the (standard) global Lanczos method when the support of the measure includes points close to the origin, and the integrand has a singularity there. The extended global Lanczos method recently has been described in [25] for the special case when the numerator and denominator degrees of the Laurent polynomials that define the approximant are about the same. Section 3 extends this method to allow more general Laurent polynomials. The method is derived by applying results in [26] for the extended Lanczos method. The computations with the method exploit that, analogously to computations with the standard Lanczos method, short recursion relations can be applied in the computations with the extended global Lanczos method. Section 4 is concerned with bounding the quadrature error of Gauss–Laurent quadrature rules associated with the extended global Lanczos method. These bounds require that certain derivatives of the integrand do not change sign on the convex hull of the support of the measure. A method that can be applied to evaluate error estimates is described in Section 5. This technique has been developed by Spalević [27], [28], [29] for estimating the error in Gauss quadrature rules. We describe an extension that can be applied to estimate the quadrature error in Gauss–Laurent rules. A few computed examples are presented in Section 6 and concluding remarks can be found in Section 7.
Section snippets
The global Lanczos algorithm and Gauss quadrature
Let be a large symmetric matrix and be a block vector with . Application of steps of the global Lanczos algorithm gives the global Lanczos decomposition where the block columns are orthonormal with respect to the inner product (1.10), i.e., The global Lanczos method simplifies to the standard Lanczos method when the block size is one. Analogously to (1.2), the range of , given by
The extended global Lanczos algorithm and Gauss–Laurent rules
The extended global Lanczos method exploits that orthogonal Laurent polynomials associated with a non-negative measure on the real axis satisfy short recursion relations; see [26], [31], [32], [33], [34], [35], [36] for this and related results.
Bounds for Gauss–Laurent rules
The degree of along with its orthogonality conditions guarantees that has distinct zeros, . These zeros are eigenvalues of ; see [26, Theorem 5.1]. Assume that is differentiable. Then the Laurent–Hermite interpolation polynomial, , that interpolates and its derivative at the nodes may be constructed. An expression for the approximation error shown in [26, Theorem 5.4] can be used to derive the following error term for
Estimating the error in Gauss–Laurent rules
Spalević [27], [28], [29] describes an approach to estimate the error in Gauss quadrature rules. Related results can be found in [40]. Extensions to Gauss–Szegő quadrature rules are discussed in [41]. This section presents an extension to Gauss–Laurent rules.
The matrices that arise in Spalević’s approach to estimate the error in Gauss quadrature rules have the structure where are symmetric tridiagonal matrices with the same spectra and
Numerical examples
The computations in this section are performed using MATLAB with about significant decimal digits. The first two examples compare the performance of the standard global Lanczos method () with the global rational Lanczos methods for the cases . In all computed examples, we use Krylov subspaces of dimension , and . These dimensions are divisible by , , and , and assure that the denominator degree of the rational Lanczos methods considered increases by at least one
Conclusion
Recurrence relations for extended global Lanczos methods with the numerator degree roughly an arbitrary integer multiple of the denominator degree are described. Error bounds determined by pairs of Gauss–Laurent and Gauss–Laurent–Radau rules are discussed. These bounds apply to certain integrands. Error estimates for more general integrands are developed. These estimates generalize error estimates for Gauss quadrature rules developed by Spalević.
Acknowledgment
This research is supported in part by NSF, USA grants DMS-1720259 and DMS-1729509.
References (45)
- et al.
Some large-scale matrix computation problems
J. Comput. Appl. Math.
(1996) - et al.
Bounding matrix functionals via partial global block Lanczos decomposition
Appl. Numer. Math.
(2015) - et al.
Orthogonality and recurrence for ordered laurent polynomial sequences
J. Comput. Appl. Math.
(2010) - et al.
The extended Krylov subspace method and orthogonal Laurent polynomials
Linear Algebra Appl.
(2009) - et al.
Recursion relations for the extended Krylov subspace method
Linear Algebra Appl.
(2011) - et al.
Orthogonal laurent polynomials and strong moment theory: a survey
J. Comput. Appl. Math.
(1999) - et al.
Generalized averaged Szegő quadrature rules
J. Comput. Appl. Math.
(2017) - et al.
On the computation of Gauss quadrature rules for measures with a monomial denominator
J. Comput. Appl. Math.
(2015) - et al.
Matrix Computations
(2013) Iterative Methods for Sparse Linear Systems
(2003)