## 2 Approximation by Gauss Quadrature  
Let \(A\) be a symmetric matrix of order \(n\) and suppose that we are interested in computing bounds for bilinear forms  
\[F_{\mathbf{u},\mathbf{v}}(A):=\mathbf{u}^{\mathrm{T}}f(A)\mathbf{v}, \tag{8}\]  
where \(\mathbf{u}\) and \(\mathbf{v}\) are given vectors and \(f\) is a smooth function defined on an interval \([a,b]\subset\mathbb{R}\) that contains the spectrum of \(A\). Since  
\[F_{\mathbf{u},\mathbf{v}}(A)=\frac{1}{4}(F_{\mathbf{u}+\mathbf{v},\mathbf{u}+ \mathbf{v}}(A)-F_{\mathbf{u}-\mathbf{v},\mathbf{u}-\mathbf{v}}(A)),\]  
we can focus on the case \(\mathbf{u}=\mathbf{v}\).  
The matrix \(A\) has the spectral decomposition \(A=Q\Lambda Q^{T}\). Then we can write  
\[F_{\mathbf{u},\mathbf{u}}(A)=\mathbf{u}^{\mathrm{T}}Qf(\Lambda)Q^{T}\mathbf{u }=\boldsymbol{\mu}^{\mathrm{T}}f(\Lambda)\boldsymbol{\mu}=\sum_{i=1}^{n}f( \lambda_{i})\mu_{i}^{2}=\int_{a}^{b}f(\lambda)\;d\mu(\lambda), \tag{9}\]  
i.e., we may regard \(F_{\mathbf{u},\mathbf{u}}(A)\) as a Stieltjes integral; see [20; 24] for further details. We approximate this integral by Gauss-type quadrature rules as follows. Let \(\mathbf{u}\) be of unit Euclidean norm. Application of \(k\) steps of the Lanczos algorithm to \(A\) with initial vector \(\mathbf{u}\) gives the \(k\times k\) symmetric tridiagonal matrix \(T_{k}\). It can be shown that \(\mathbf{e}_{1}^{\mathrm{T}}f(T_{k})\mathbf{e}_{1}\) is a \(k\)-node Gauss quadrature rule \(\mathcal{G}_{k}\) for the Stieltjes integral (9). A \((k+1)\)-node Gauss-Radau quadrature formula \(\mathcal{G}_{k+1}\) with a fixed node at \(\theta\geq\rho(A)\) for approximating the Stieltjes integral also can be defined. This discussion assumes that the Lanczos algorithm does not break down. Breakdown is very rare and allows the computations to be simplified.  
Under the assumption that the derivatives of \(f(x)\) have constant sign on the convex hull of the support of the measure, which is met by the functions (2) and (3), and the Radau node \(\theta\) is suitably chosen, pairs of Gauss and Gauss-Radau rules furnish lower and upper bounds of increasing accuracy for the quadratic form (9). For the functions (2) and (3), and the Radau node \(\theta\) chosen as described, we have  
\[\mathcal{G}_{k}\leq\mathcal{G}_{k+1}\leq F_{\mathbf{u},\mathbf{u}}(A)\leq \mathcal{G}_{k+2}\leq\mathcal{G}_{k+1}.\]  
For a user-chosen accuracy \(\tau\), we terminate the iterations with the [PERSON] algorithm when  
\[\frac{|\mathcal{G}_{k}-\mathcal{G}_{k+1}|}{|\mathcal{G}_{k}|}\leq\tau. \tag{10}\]The default value in the code is \(\tau=10^{-3}\).  
The matrix functions are applied to the tridiagonal matrixes using their spectral factorization. Thus, let \(T_{k}=W_{k}\Lambda_{k}W_{k}^{T}\) be the spectral factorization. Then \(f(T_{k})=W_{k}f(\Lambda_{k})W_{k}^{T}\). When \(f\) is the exponential function, we let \(\mu\) be the largest eigenvalue of \(T_{k}\) and evaluate  
\[\mathrm{e}^{-\mu}\mathrm{e}^{T_{k}}=\mathrm{e}^{T_{k}-\mu I}=W\mathrm{e}^{ \Lambda_{k}-\mu I}W^{T} \tag{11}\]  
instead of \(\mathrm{e}^{T_{k}}\) to avoid overflow.  
Regarding the choice of the Radau node \(\theta\), we often may let \(\theta=\|A\|_{\infty}\). Alternatively, we can use the MATLAB function eigs or the function irbleigs described in [25; 26] to determine an estimate of the largest eigenvalue \(\lambda_{1}\) of \(A\).