The extended global Lanczos method, Gauss–Radau quadrature, and matrix function approximation

In memory of Bill Gragg
https://doi.org/10.1016/j.cam.2020.113027Get rights and content

Abstract

The need to evaluate expressions of the form I(f) trace (WTf(A)W), where the matrix ARn×n is symmetric, WRn×k with 1kn, and f is a function defined on the convex hull of the spectrum of A, arises in many applications including network analysis and machine learning. When the matrix A is large, the evaluation of I(f) by first computing f(A) may be prohibitively expensive. In this situation it is attractive to compute an approximation of I(f) by first applying a few steps of a global Lanczos-type method to reduce A to a small matrix and then evaluating f at this reduced matrix. The computed approximation can be interpreted as a quadrature rule. The present paper generalizes the extended global Lanczos method introduced in Bentbib et al. (2018) and discusses the computation of error-bounds and error estimates. Numerical examples illustrate the performance of the techniques described.

Introduction

Let ARn×n be a large symmetric matrix, f a function defined on the convex hull of the spectrum of A, and wRn a vector of unit Euclidean norm. Application of mn steps of the (standard) Lanczos method to A with initial vector w yields the partial Lanczos decomposition AVm=VmTm+gmemT,where the matrix VmRn×m has orthonormal columns with initial column w, TmRm×m is a symmetric tridiagonal matrix, and gmRn satisfies VmTgm=0; see, e.g., [1], [2] for details on the Lanczos algorithm. Throughout this paper ei denotes the ith column of the identity matrix of suitable order, and the superscript T stands for transposition. Generally, mn in applications of the Lanczos method. We assume that m is small enough so that the Lanczos decomposition (1.1) with the stated properties exists. This is the generic situation. Then range(Vm)=Km(A,w),where the right-hand side denotes the Krylov subspace Km(A,w)span{w,Aw,,Am1w}.

Introduce the spectral factorization A=UΛUT,Λ=diag[λ1,λ2,,λn]Rn×nwith an orthogonal matrix URn×n. Substituting the spectral factorization into I(f)wTf(A)wgives I(f)=wTUf(Λ)UTw=j=1mf(λj)μj2,where μj=ejTUTw. 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 λj. Let dμ(λ) denote the measure associated with μ(λ). Then I(f)=f(λ)dμ(λ).We note for future reference that the measure dμ(λ) defines an inner product (f,g)f(λ)g(λ)dμ(λ)for polynomials f and g of low enough degree.

Golub and Meurant [3], [4] observed that the nontrivial entries of the symmetric tridiagonal matrix Tm 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 Tm are the zeros of the orthonormal polynomial of degree m. This led Golub and Meurant [3], [4] to show that the expression Gm(f)e1Tf(Tm)e1is the m-point Gauss quadrature rule for approximating the Stieltjes integral (1.6). Hence, Gm(f)=I(f)fP2m1,where P2m1 denotes the set of polynomials of degree at most 2m1. Substituting the spectral factorization of Tm into (1.8) transforms the right-hand side into a sum of m terms and shows that the eigenvalues of Tm 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 Tm. 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 f at the nodes. For some functions f, it may be more efficient to compute f(Tm) directly by one of the methods discussed in [6]; see [7] for an illustration.

The quadrature rule (1.8) is obtained by applying m steps of the (standard) Lanczos method to A with initial vector w. An analogous rule can be computed by applying the global Lanczos method to A with initial block vector WRn×k with block size 1<kn. 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 W1,W2trace(W1TW2),W1,W2Rn×k,and the induced Frobenius norm W1FW1,W112.This method can be applied to approximate expressions of the form I(f)trace(WTf(A)W),where WRn×k with 1kn. The problem of estimating the trace of a matrix f(A), without evaluating f(A), 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 f(A) can be computed by approximating expressions (1.11) for several block vectors WEm, 1mnk, where the Em=[ek(m1)+1,,ekm]Rn×kare block axis vectors. For notational simplicity, we here assume that n is a multiple of k; 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 A 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 ARn×n be a large symmetric matrix and WRn×k be a block vector with 1kn. Application of m steps of the global Lanczos algorithm gives the global Lanczos decomposition A[V1,V2,,Vm]=[V1,V2,,Vm]T̂mk+βm+1Vm+1EmT,where the block columns VjRn×k are orthonormal with respect to the inner product (1.10), i.e., Vi,Vj=1,i=j,0,ij.The global Lanczos method simplifies to the standard Lanczos method when the block size k is one. Analogously to (1.2), the range of [V1,V2,,Vm], given by j=1mγjVj,γ

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 ψτLm1,im+1Lm1,im along with its orthogonality conditions guarantees that ψτ has τ distinct zeros, {xj}j=1τ. These zeros are eigenvalues of Hτ; see [26, Theorem 5.1]. Assume that f is differentiable. Then the Laurent–Hermite interpolation polynomial, LˆL2m2,2im+1, that interpolates f and its derivative at the nodes xj may be constructed. An expression for the approximation error fLˆL2m2,2im+1 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 M̂=M1β1eτOβ1eτTαβ2e1TOβ2e1M2R(2τ+1)×(2τ+1),where M1,M2Rτ×τ are symmetric tridiagonal matrices with the same spectra and

Numerical examples

The computations in this section are performed using MATLAB with about 15 significant decimal digits. The first two examples compare the performance of the standard global Lanczos method (i=0) with the global rational Lanczos methods for the cases i=1,2,3. In all computed examples, we use Krylov subspaces of dimension τ=12,24,36,48, and 60. These dimensions are divisible by 2, 3, and 4, 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)

  • GolubG.H. et al.

    Matrices, moments and quadrature

  • GolubG.H. et al.

    Matrices, Moments and Quadrature with Applications

    (2010)
  • GolubG.H. et al.

    Calculation of Gauss quadrature rules

    Math. Comp.

    (1969)
  • HighamN.J.

    Functions of Matrices: Theory and Computation

    (2008)
  • FenuC. et al.

    Network analysis via partial spectral factorization and Gauss quadrature

    SIAM J. Sci. Comput.

    (2013)
  • ElbouyahyaouiL. et al.

    Algebraic properties of the block GMRES and block Arnoldi methods

    Electron. Trans. Numer. Anal.

    (2009)
  • JbilouK. et al.

    Oblique projection methods for multiple linear systems

    Electron. Trans. Numer. Anal.

    (2005)
  • AvronH. et al.

    Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix

    J. ACM

    (2011)
  • BaiZ. et al.

    Bounds for the trace of the inverse and the determinant of symmetric positive definite matrices

    Ann. Numer. Math.

    (1997)
  • BrezinskiC. et al.

    Estimations of the trace of powers of self-adjoint operators by extrapolation of the moments

    Electron. Trans. Numer. Anal.

    (2012)
  • BrezinskiC. et al.

    Moments of a linear operator on a Hilbert space, with applications to the trace of the inverse of matrices and the solution of equations

    Numer. Linear Algebra Appl.

    (2012)
  • GolubG.H. et al.

    Generalized cross-validation as a method for choosing a good ridge parameter

    Technometrics

    (1979)
  • Cited by (0)

    View full text