We explore quadratures for \({\cal F}(s)=B^T \phi (A,s) B\) where \(A\) is a symmetric, nonnegative-definite matrix in \(\mathbb {R}^{n \times n}\), \(B\) is a tall matrix in \(\mathbb {R}^{n \times p}\), and \(\phi (\cdot ,s)\) is a matrix function with parameter \(s\) [1, 2]. These formulations commonly arise in the computation of multiple-input, multiple-output transfer functions for diffusion PDEs.

We derive bounds and averaging schemes for quadrature rules for \(B^T \phi (A,s) B\) computed via the block-Lanczos algorithm, which are particularly efficient for discretizations of PDE operators with continuous spectra. Additionally, we demonstrate that these bounds and averaging schemes are applicable to parametric model reduction of dynamical systems via Galerkin projections.

1 Block-Lanczos Approximations to \(B^T \phi (A,s) B\)

We propose an approximation scheme for \({\cal F}(s)=B^T \phi (A,s) B\) leveraging the block-Lanczos algorithm [3] and its representation via Stieltjes matrix continued fractions.

The block-Lanczos recursion for the block-Lanczos vectors \(Q_{i}\in \mathbb {R}^{n\times p}\) reads

\begin{equation} \label {eq:Lanczos} A Q_i=Q_{i+1} {\boldsymbol \beta }_{i+1} +Q_{i} {\boldsymbol \alpha }_{i}+Q_{i-1} ({\boldsymbol \beta }_{i})^T, \end{equation}

with block coefficients \({\boldsymbol \alpha }_{i},{\boldsymbol \beta }_{i} \in \mathbb {R}^{p\times p}\). Using Stieltjes matrix continued fractions, we show that this block-Lanczos algorithm defines a block-Gauss quadrature approximation \({\cal F}_m(s)\) and converges monotonically for \(\phi (A,s) =(A+sI)^{-1}\) with \(s \in \mathbb {R}^+\) and \(I\) the identity matrix. We further show that a monotonically convergent block Gauss-Radau quadrature \(\tilde {\cal F}_m(s)\) can be readily defined through this Stieltjes continued fraction representation. In the literature, Gauss-Radau quadratures for symmetric matrices are often defined via rank-one updates of the Lanczos matrix in the non-block case [4, 5] or rank \(p\) updates in the block case [6].

Here we define Gauss-Radau quadrature through Stieltjes matrix continued fractions which can be written via the recursion

\begin{equation*} \mathcal {C}_j(s)=\cfrac {1}{s \hat {\boldsymbol \gamma }_j + \cfrac {1}{{\boldsymbol \gamma }_j+\mathcal {C}_{j+1}(s)}} , \end{equation*}

where \(\hat {\boldsymbol \gamma }_j, {\boldsymbol \gamma }_j \in \mathbb {R}^{p \times p}\) are symmetric positive definite matrices directly related to the block-Lanczos coefficients \({\boldsymbol \alpha }_j\) and \({\boldsymbol \beta }_j\). We show that the Gauss quadrature approximation to \(B^T (A+sI)^{-1} B\) after \(m\) iterations of block-Lanczos corresponds to \(\mathcal {C}_1(s)\), defined through the above recursion, terminated with \(\mathcal {C}_{m+1}=0\), whereas the Gauss-Radau quadrature corresponds to truncation with \(\mathcal {C}_{m+1}=\infty \).

Through Stieltjes matrix continued fractions, we demonstrate that Gauss quadrature provides a lower bound, while Gauss-Radau quadrature provides an upper bound to the matrix function. Combined with the monotonicity result, this yields an ordering of the Gauss and Gauss-Radau quadrature approximations to \(\cal F\). Given \(m\), the quadrature order, we obtain

\[ 0<{\cal F}_{m-1}(s)<{\cal F}_{m}(s)< {\cal F}(s)< \tilde {\cal F}_{m}(s)<\tilde {\cal F}_{m-1}(s) \quad \forall s\in \mathbb {R}_+, \]

where, for two symmetric matrices \(G_1, G_2\), the notation \(G_1<G_2\) indicates that \(G_2-G_1\) is positive definite.

This ordering further enables derivation of easily computable error bounds of the form

\[ \|{\cal F}- {\cal F}_m \| < \|\tilde {\cal F}_m - {\cal F}_m \|. \]

In this contribution, we discuss averaging schemes of Gauss and Gauss-Radau quadrature motivated by potential theory for Padé approximations. We show numerical examples for various \(\phi (A,s)\), where \(A\) is a graph Laplacian or discretization of an operator with a continuous spectrum (e.g., PDE operators in unbounded domains). These examples demonstrate that the derived error bound is tight in important applications, and that the averaging schemes reduce the approximation error by an order of magnitude for discretizations of operators with continuous spectra, as shown in [7]. In the next section we discuss the applicability of such Gauss-Radau bounds to parametric model reduction, a subsequent result not covered in [7].

2 Parametric Model Reduction via Galerkin Projection

The first iteration of the block-Lanczos procedure can be interpreted as a Galerkin projection of a symmetric positive definite matrix on a general basis. This insight allows us to directly apply the Gauss-Radau bound and averaging schemes to, for instance, projection-based parametric model reduction [8].

Consider the parametric dynamical system

\[ (A(\rho )+s I)X(s,\rho ) = B, \quad \text {with transfer function } {\cal F}(s) = B^T X(s,\rho ), \]

where \(A(\rho )\in \mathbb {R}^{n \times n}\) is symmetric positive definite for all parameters \(\rho \) of interest, \(B\in \mathbb {R}^{n\times p}\), and the Laplace frequency \(s \in \mathbb {R}^+\). This formulation arises in inverse problems where \(A(\rho )\) represents a discretized PDE operator with PDE coefficients parametrized by \(\rho \).

Assuming for simplicity, that the Galerkin projection basis \(U\in \mathbb {R}^{n \times k}\) contains \(B\) as \(U=[B, U_0]\) and is orthogonalized \(U^TU=I_k\) (e.g. a rational Krylov subspace with a shift at \(\infty \)). Then the Galerkin approximation of \(\cal F\) is

\begin{align*} {\cal F}_{\rm Gal}(s) &= (U^TB)^T (U^TA(\rho )U+s I_k)^{-1} (U^T B) \\ {} &=\begin{bmatrix}I_p\\ 0\end {bmatrix}^T (H^{\rm ROM}(\rho )+sI_k)^{-1} \begin{bmatrix}I_p\\ 0\end {bmatrix}. \end{align*}

We can interpret \(U\) as the first Lanczos vector \(Q_1\) in equation (1) for \(i=1\) with \(Q_0=0\) and \(H^{\rm ROM}(\rho )\) as the \({\boldsymbol \alpha }_{1}\) block-Lanczos coefficents. Then we can further define the \({\boldsymbol \beta }_2\) coefficent according to the block-Lanczos recursion

\[ ({\boldsymbol \beta }_2)^T{\boldsymbol \beta }_2=U^TA^2U - (H^{\rm ROM}(\rho ))^2. \]

Then the Gauss-Radau quadrature approximation to the \(k \times k\) transfer function \({\cal F}_U=U^T(A(\rho )+sI)U\) as defined in [7] reads

\[ \tilde {\cal F}_{U}(s) = \begin {bmatrix}I_k\\ 0\end {bmatrix}^T \left ( \begin {bmatrix}H^{\rm ROM}(\rho ) & {\boldsymbol \beta }_2^T\\{\boldsymbol \beta }_2& {\boldsymbol \beta }_2 (H^{{\rm ROM}}(\rho ))^{-1} {\boldsymbol \beta }_2^T\end {bmatrix} +sI_{2k} \right )^{-1} \begin {bmatrix}I_k\\ 0\end {bmatrix}. \]

Since the Gauss Radau approximation \(\tilde {\cal F}_{U}\) is an upper bound to \({\cal F}_{U}\) their difference \(\tilde {\cal F}_{U}-{\cal F}_{U}\) is s.p.d. and any subspace projection of a s.p.d. matrix is s.p.d. Hence the leading \(p \times p\) block of \(\tilde {\cal F}_{U}\) provides an easy-to-compute error bound for the leading block of \({\cal F}_{U}\) which coincides with \({\cal F}_{\rm Gal}\) and holds for any \(U\) and \(\rho \). In this contribution, we will show that such a bound provides a meaningful tool for applications in inverse problems and greedy selection of interpolation points in parametric model reduction.

References

• [1]  Gene H. Golub and Gérard Meurant. Matrices, Moments and Quadrature with Applications. Princeton University Press, 2010.

• [2]  C. Fenu, D. Martin, L. Reichel, and G. Rodriguez. Block Gauss and anti-Gauss quadrature with application to networks. SIAM Journal on Matrix Analysis and Applications, 34(4):1655–1684, 2013.

• [3]  G.H. Golub and R. Underwood. The block Lanczos method for computing eigenvalues. In John R. Rice, editor, Mathematical Software, pages 361–377. Academic Press, 1977.

• [4]  G. López Lagomasino, L. Reichel, and L. Wunderlich. Matrices, moments, and rational quadrature. Linear Algebra and its Applications, 429(10):2540–2554, 2008. Special Issue in honor of Richard S. Varga.

• [5]  Andreas Frommer, Kathryn Lund, Marcel Schweitzer, and Daniel B. Szyld. The Radau–Lanczos method for matrix functions. SIAM Journal on Matrix Analysis and Applications, 38(3):710–732, 2017.

• [6]  Kathryn Lund. A New Block Krylov Subspace Framework with Applications to Functions of Matrices Acting on Multiple Vectors. Phd thesis, Department of Mathematics, Temple University and Fakultät Mathematik und Naturwissenschaften der Bergischen Universität Wuppertal, Philadelphia, Pennsylvania, USA, May 2018.

• [7]  Jörn Zimmerling, Vladimir Druskin and Valeria Simoncini. Monotonicity, bounds and extrapolation of Block-Gauss and Gauss-Radau quadrature for computing \(B^T \phi (A) B\), 2024; arXiv:2407.21505.

• [8]  Benner, P., Gugercin, S. & Willcox, K. A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems. SIAM Review. 57, 483-531 (2015),