Monotonicity, bounds and extrapolation of Block-Gauss and Gauss-Radau quadrature for computing $B^T \phi (A) B$

In this paper, we explore quadratures for the evaluation of $B^T \phi(A) B$ where $A$ is a symmetric positive-definite (s.p.d.) matrix in $\mathbb{R}^{n \times n}$, $B$ is a tall matrix in $\mathbb{R}^{n \times p}$, and $\phi(\cdot)$ represents a matrix function that is regular enough in the neighborhood of $A$'s spectrum, e.g., a Stieltjes or exponential function. These formulations, for example, commonly arise in the computation of multiple-input multiple-output (MIMO) transfer functions for diffusion PDEs. We propose an approximation scheme for $B^T \phi(A) B$ leveraging the block Lanczos algorithm and its equivalent representation through Stieltjes matrix continued fractions. We extend the notion of Gauss-Radau quadrature to the block case, facilitating the derivation of easily computable error bounds. For problems stemming from the discretization of self-adjoint operators with a continuous spectrum, we obtain sharp estimates grounded in potential theory for Pad\'e approximations and justify averaging algorithms at no added computational cost. The obtained results are illustrated on large-scale examples of 2D diffusion and 3D Maxwell's equations and a graph from the SNAP repository. We also present promising experimental results on convergence acceleration via random enrichment of the initial block $B$.