Pole Allocation for Rational Gauss Quadrature Rules for Matrix Functionals Defined by a Stieltjes Function

Abstract

This paper considers the computation of approximations of matrix functionals of form $F\left(A\right):={\mathit{v}}^{T}f\left(A\right)\mathit{v}$, where A is a large symmetric positive definite matrix, v is a vector, and f is a Stieltjes function. The functional $F\left(A\right)$ is approximated by a rational Gauss quadrature rule with poles on the negative real axis (or part thereof) in the complex plane, and we focus on the allocation of the poles. Specifically, we propose that the poles, when considered positive point charges, be allocated to make the negative real axis (or part thereof) approximate an equipotential curve. This is easily achieved with the aid of conformal mapping.

1. Introduction

We are concerned with the development of numerical methods for the approximation of expressions of the form

$$F\left(A\right):={\mathit{v}}^{T}f\left(A\right)\mathit{v},$$

(1)

where $A\in {\mathbb{R}}^{N\times N}$ is a large symmetric positive definite matrix, the vector $\mathit{v}\in {\mathbb{R}}^N\backslash \left\{\mathbf{0}\right\}$ is of unit Euclidean norm, f is a Stieltjes function, and the superscript ^T denotes transposition. Thus, f can be represented as

$$f\left(z\right)={\int }_0^{\infty }\frac{1}{t+z}d\mu \left(t\right),z\in \mathbb{C}\backslash \left(-\infty ,0\right],$$

(2)

where $\mu$ is a nonnegative distribution function defined in the interval $[0,\infty )$ such that $f\left(z\right)$ is defined; see, e.g., [1,2,3,4] for discussions and illustrations of Stieltjes functions. Examples of Stieltjes functions include

$$\begin{array}{ccc}\hfill f\left(z\right)& =& z^{-a}=\frac{sin\left(a\pi \right)}{\pi }{\int }_0^{\infty }\frac{1}{t+z}d\mu \left(t\right),\mathrm{with}d\mu \left(t\right)=t^{-a}dt,a\in (0,1),\hfill \\ \hfill f\left(z\right)& =& \frac{log(z+1)}{z}={\int }_1^{\infty }\frac{1}{t+z}d\mu \left(t\right),\mathrm{with}d\mu \left(t\right)=t^{-1}dt.\hfill \end{array}$$

When the matrix A is large, straightforward evaluation of (1) by first computing $f\left(A\right)$ may be prohibitively expensive; see Higham [5] for discussions and analyses of many methods for the evaluation of matrix functions. We are, therefore, interested in exploring the possibility of reducing A to a fairly small matrix $H_m$ by a rational Krylov subspace method and approximating $F\left(A\right)$ by evaluating $f\left(H_m\right)$. This results in a rational approximant of $F\left(A\right)$, whose poles are allocated by the user on the negative real axis. The rational approximant can be chosen to converge to zero as $z\to \infty$, similarly as $f\left(z\right)$. This paper focuses on the choice of poles. In particular, we propose to allocate the poles so that, when considered as positive point charges, they roughly make the negative real axis or part thereof an equipotential curve. This can easily be achieved with the aid of conformal mapping.

Similarly, as in [6], we express the rational approximant as a linear combination of rational orthogonal functions with specified poles and use the fact that these functions satisfy a recurrence relation with few terms only. We use a rational Gauss quadrature rule to determine an approximation of an integral of this rational approximant. The novelty of the present paper is that it addresses the allocation of the poles of the rational approximant. This issue was ignored in [6]. A nice introduction to rational Gauss quadrature is provided by Gautschi (Section 3.1.4 in [7]).

The choice of poles for rational approximants of Stieltjes functions has previously been discussed by Güttel and Knizhnerman [8], who propose a greedy algorithm, and by Massei and Robol [4], who express the approximation problem as a Zolotarev minimization problem. We will allocate the poles with the aid of conformal mapping. This approach is based on the observation that the approximation error is invariant under conformal mapping. Specifically, we allocate the poles equidistantly on the unit circle and then map the circle conformally to the negative real axis or a subset thereof. A benefit of this approach is that it is very easy to implement; in particular, it does not require the solution of a minimization problem.

This paper is organized as follows. Section 2 expresses (1) as an integral of the Stieltjes function over the positive real axis with respect to a measure that is determined by both the matrix A and vector v. This kind of transformation is described, e.g., by Golub and Meurant (Chapter 7 in [9]). Properties of orthogonal rational functions with real poles are reviewed in Section 3, and rational Gauss quadrature is discussed in Section 4. Our approach to choosing the poles is described in Section 5. A few computed examples are presented in Section 6, and concluding remarks can be found in Section 7.

2. Integration of Stieltjes Functions

We derive an integral that will be approximated by a rational Gauss quadrature rule. A similar derivation can be found in [6]. Our derivation uses the spectral factorization

$$A=U\mathsf{\Lambda }{U}^T,\mathsf{\Lambda }=\mathrm{diag}[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_N],$$

(3)

where the matrix $U\in {\mathbb{R}}^{N\times N}$ is orthogonal and the eigenvalues ${\lambda }_i$ of A are ordered according to $0<{\lambda }_1\le \cdots \le {\lambda }_N$. We remark that the spectral factorization (3) is helpful for deriving our approximation method. However, the application of the quadrature rules of this paper does not require the computation of this factorization.

We introduce the vector $[{\nu }_1,{\nu }_2,\dots ,{\nu }_N]:={\mathit{v}}^{T}U$. Using the factorization (3), the expression (1) with f defined by (2) can be written as

$$\begin{array}{cc}\hfill {\mathit{v}}^{T}f\left(A\right)\mathit{v}={\mathit{v}}^{T}Uf\left(\mathsf{\Lambda }\right)U^T\mathit{v}& ={\int }_0^{\infty }{\mathit{v}}^{T}U{(tI+\mathsf{\Lambda })}^{-1}{U}^T\mathit{v}d\mu \left(t\right)\hfill \\ \\ ={\int }_0^{\infty }\sum _{i=1}^N{(t+{\lambda }_i)}^{-1}{\nu }_i^{2}d\mu \left(t\right)\hfill \\ \\ ={\int }_0^{\infty }{\int }_0^{\infty }{(t+y)}^{-1}d\nu \left(y\right)d\mu \left(t\right)\hfill \\ \\ ={\int }_0^{\infty }f\left(y\right)d\nu \left(y\right)=:\mathcal{I}\left(f\right),\hfill \end{array}$$

(4)

where the nonnegative measure $d\nu \left(y\right)$ has support at the eigenvalues ${\lambda }_i$ of A. The distribution function $\nu \left(y\right)$ can be chosen to be a nondecreasing piece-wise constant function with jumps of size ${\nu }_i^2$ at the eigenvalues ${\lambda }_i$ of A. We will approximate the integral (4) by a rational Gauss quadrature rule.

3. Recursion Relations for Orthogonal Rational Function

We first review results in [10] on recursion relations for certain orthogonal rational functions. In this reference, functions with real poles and complex conjugate poles are allowed. Here we focus on rational functions with real poles only. The number of terms in the recursion relations depends on the number of distinct poles of the rational functions as well as on the ordering of certain elementary rational basis functions to be defined below.

We introduce linear spaces of rational functions with finite distinct real poles ${\alpha }_i$ of multiplicity $k_i$,

$${\mathbb{Q}}_{i,k_i}=\mathrm{span}\left\{\frac{1}{{\left(y-{\alpha }_i\right)}^j}:j=1,2,\dots ,k_i\right\},i=1,2,\dots ,\ell .$$

(5)

In our application, the poles will be allocated on the negative real axis. Let

$$k=\sum _{i=1}^{\ell }{k}_i$$

(6)

and define the $(m+1)$-dimensional linear space

$${\mathbb{S}}_{m+1}:={\mathbb{P}}_{m+1-k}\oplus {\mathbb{Q}}_{1,k_1}\oplus \cdots \oplus {\mathbb{Q}}_{\ell ,k_{\ell }}.$$

(7)

Here $k_i$ is chosen so that $0\le k\le m-1$. The space (7) then contains linear functions. Let

$${\mathsf{\Psi }}_{m+1}=\left\{{\psi }_0,{\psi }_1,\dots ,{\psi }_m\right\}$$

(8)

denote an elementary basis for the space ${\mathbb{S}}_{m+1}$, i.e., ${\psi }_0\left(y\right)=1$ and each basis function ${\psi }_i\left(y\right)$, for $i=1,2,\dots ,m$, is one of the functions

$$y^j,\frac{1}{{\left(y-{\alpha }_i\right)}^j},$$

(9)

for some positive integers i and j. We use the notation ${\psi }_s\prec {\psi }_t$ to indicate that the basis function ${\psi }_s$ comes before ${\psi }_t$. We refer to the ordering of the basis functions (8) as natural if ${\psi }_0\left(y\right)=1$ and the remaining functions ${\psi }_j$, $j=1,2,\dots ,m$, satisfy:

• $y^p\prec y^{p+1}$ for all integers $p>0$,
• $\frac{1}{{\left(y-{\alpha }_i\right)}^p}\prec \frac{1}{{\left(y-{\alpha }_i\right)}^{p+1}}$ for all integers $p>0$ and every real pole ${\alpha }_i$.

In particular, if ${\psi }_j=y^p$ for some $p>1$, then some preceding elementary basis function ${\psi }_i$ with $i<j$ equals $y^{p-1}$. An analogous statement also holds for negative powers ${(y-{\alpha }_i)}^{-p}$ for $p>1$.

It is shown in [10] that a Stieltjes procedure for orthogonalizing elementary basis functions with respect to some inner product only requires short recursion relations when the basis functions are in natural order. A recursion relation is said to be short if the number of terms in the relation can be bounded independently of the number of orthogonal functions generated. A well-known special case is a situation when all elementary basis functions are monomials, i.e., when ${\psi }_j\left(y\right)=y^j$, $j=0,1,2,\dots$. Then the Stieltjes procedure yields a family of orthogonal polynomials ${\varphi }_0,{\varphi }_1,{\varphi }_2,\dots$ with ${\varphi }_j$ of degree j for all j. The polynomials ${\varphi }_j$ satisfy a three-term recurrence relation.

Similarly, when orthogonalizing elementary basis functions in the natural order with respect to a bilinear form $[f,g]=\mathcal{L}\left(fg\right)$, where $\mathcal{L}$ is a linear functional on ${\mathbb{S}}_{m+1}$, the orthogonal rational basis functions ${\varphi }_0,{\varphi }_1,{\varphi }_2,\dots$ obtained satisfy short recurrence relations, provided that $\mathcal{L}\left({\varphi }_j^2\right)\ne 0$ for all j. We will orthogonalize the elementary basis functions with respect to the bilinear form

$$\langle f,g\rangle ={\left(f\left(A\right)\mathit{v}\right)}^T\left(g\left(A\right)\mathit{v}\right).$$

(10)

This bilinear form is an inner product for all elementary basis functions of low enough order.

When the elementary basis functions are of the form (9), the number of terms required in the recurrence relations depends on how often powers of the same function appear in the sequence (8). We review two types of recursion relations from [10] that are required below. First, consider the recursion relation

$$y{\varphi }_r\left(y\right)=\sum _{j=-n_1}^{n_2}{c}_{r,r+j}{\varphi }_{r+j}\left(y\right),r=0,1,\dots ,$$

(11)

where $r-n_1$ is the largest integer smaller than r such that ${\psi }_{r-n_1}$ is a monomial if there is such a monomial (otherwise $r-n_1=0$), and $r+n_2$ is the smallest integer larger than r such that ${\psi }_{r+n_2}$ is a monomial. To introduce a real pole α we need the recursion relation

$$\frac{1}{y-{\alpha }_i}{\varphi }_r\left(y\right)=\sum _{j=-n_3}^{n_4}{c}_{r,r+j}^{\left(i\right)}{\varphi }_{r+j}\left(y\right),r=0,1,2,\dots ,$$

(12)

where $r-n_3$ is the largest integer smaller than r such that ${\psi }_{r-n_3}$ is a rational function with a pole at ${\alpha }_i$, provided that there is such a rational function (otherwise $r-n_3=0$), and $r+n_4$ is the smallest integer larger than r such that ${\psi }_{r+n_4}$ is a rational function with a pole at ${\alpha }_i$.

We introduce the vector of orthonormal rational functions,

$${\mathsf{\Phi }}_m\left(y\right):={[{\varphi }_0\left(y\right),{\varphi }_1\left(y\right),\dots ,{\varphi }_{m-1}\left(y\right)]}^T.$$

Assume that ${\psi }_m$ is a monomial and that $m-d$ is the largest integer smaller than m such that ${\psi }_{m-d}$ is a monomial. Then the recursion formulas (12) can be written in the form

$$y{\mathsf{\Phi }}_m\left(y\right)=H_m{\mathsf{\Phi }}_m\left(y\right)+\sum _{j=1}^d{h}_{m-j,m}{\varphi }_m\left(y\right){\mathit{e}}_{m+1-j},$$

(13)

where $e_j={[0,\dots ,0,1,0,\dots ,0]}^T\in {\mathbb{R}}^m$ denotes the jth axis vector. The matrix $H_m$ in (13) has the following block-diagonal structure: the matrix has $m-k-1$ square blocks along the diagonal, such that any two consecutive blocks overlap in one diagonal element. The jth block of $H_m$ is of dimension $r\times r$, where $r-2$ is the number of rational functions between consecutive monomials $y^{j-1}$ and $y^j$ for some j. In detail, the jth block of $H_m$ is the submatrix $H_m\left(r_1:r_2,r_1:r_2\right)$ (using MATLAB notation) with the entries $h_{ij}$ for $r_1\le i,j\le r_2$, where ${\psi }_{r_1}\left(y\right)=y^{j-1}$ and ${\psi }_{r_2}\left(y\right)=y^j$. Further, $r=r_2-r_1+1$. The nonvanishing entries of $H_m={\left[h_{ij}\right]}_{i,j=0}^{m-1}$ are recursion coefficients for the orthogonal rational functions ${\varphi }_i,i=0,1,\dots ,m-1$. The recursion coefficients satisfy

$$h_{i,j}=\langle y{\varphi }_i\left(y\right),{\varphi }_j\left(y\right)\rangle =\langle y{\varphi }_j\left(y\right),{\varphi }_i\left(y\right)\rangle =h_{j,i}.$$

This shows that the matrix $H_m$ depends only on the first m elementary basis functions ${\psi }_j$, $j=0,1,\dots ,m-1$. In particular, $H_m$ is independent of ${\psi }_m$. Throughout this paper, we will assume that ${\psi }_m$ is a monomial because then the zeros of ${\varphi }_m$ are eigenvalues of $H_m$.

The structure of the matrix $H_m$ is illustrated by the following example.

Example 1.

Consider the elementary basis

$$\begin{array}{cc}\hfill & \left\{1,\frac{1}{y-{\alpha }_1},y,\frac{1}{{(y-{\alpha }_1)}^2},y^2,\dots ,\frac{1}{y-{\alpha }_{\ell }},y^k\right\}.\hfill \end{array}$$

(14)

This elementary basis, together with the function ${\psi }_m\left(y\right)=y^{k+1}$, where k is defined by (6) and $m=1+2k$, satisfies the requirements of natural ordering. The basis (14) spans the space ${\mathbb{S}}_{m+1}$. The matrix $H_m$ determined by this basis has k 3 × 3 blocks along the diagonal,

$$H_m=\left[\begin{array}{ccccccccccc}\ast & \ast & \ast & & & & & & & & \\ \ast & \ast & \ast & & & & & & & & \\ \ast & \ast & \ast & \ast & \ast & & & & & & \\ & & \ast & \ast & \ast & & & & & & \\ & & \ast & \ast & \ast & & & & & & \\ & & & & & \ddots & & & & & \\ & & & & & & \ast & \ast & \ast & & \\ & & & & & & \ast & \ast & \ast & & \\ & & & & & & \ast & \ast & \ast & \ast & \ast \\ & & & & & & & & \ast & \ast & \ast \\ & & & & & & & & \ast & \ast & \ast \end{array}\right].$$

Here we mark matrix entries that may be nonvanishing by ∗. If we place one more elementary rational basis function between the two consecutive monomials $y^{j-1}$ and $y^j$ in the elementary basis, then the size of the jth block of $H_m$ increases by one. Analogously, if we remove one elementary rational basis function between two consecutive monomials, then the size of the corresponding block decreases by one.

Assume that the functions ${\varphi }_j$, $j=0,1,\dots ,m$, are orthonormal with respect to the bilinear form (10). Then the vectors

$$\{{\mathit{v}}_0={\varphi }_0\left(A\right)\mathit{v},\dots ,{\mathit{v}}_j={\varphi }_j\left(A\right)\mathit{v}\}$$

form an orthonormal basis for the rational Krylov subspace

$${\mathbb{K}}_{j+1}(A,\mathit{v})=\mathrm{span}\{{\psi }_0\left(A\right)\mathit{v},{\psi }_1\left(A\right)\mathit{v},\dots ,{\psi }_j\left(A\right)\mathit{v}\},$$

(15)

for $j=0,1,\dots ,m$, with respect to the bilinear form (10). Indeed,

$${\mathit{v}}_j^T{\mathit{v}}_i={\left({\varphi }_j\left(A\right)\mathit{v}\right)}^T\left({\varphi }_i\left(A\right)\mathit{v}\right)=\langle {\varphi }_j,{\varphi }_i\rangle .$$

We note that the vectors ${\mathit{v}}_j$ and the orthogonal rational functions ${\varphi }_j$ satisfy the same recursion relations. The vectors, therefore, can be determined by the rational Lanczos process (Algorithm 1). This process is analogous to the Stieltjes-type procedure (Algorithm 3.1 in [10]) for computing an orthonormal basis for the space ${\mathbb{S}}_{m+1}$. The rational Lanczos process is based on the recursion relations (12). The norm $\parallel ·\parallel$ in Algorithm 1 is the Euclidean vector norm.

The execution of Algorithm 1 requires the solution of linear systems of equations with matrices of the form $A-{\alpha }_{i}I$, where the ${\alpha }_i$ are suitable real scalars. We assume, in this paper, that the matrix A has a structure that allows fairly efficient computation of ${(A-{\alpha }_{i}I)}^{-1}{\mathit{v}}_{r-1}$ by a direct or an iterative method. This is, for instance, the case when A is banded with a small bandwidth or when A is a symmetric positive Toeplitz matrix; see, e.g., ref. [11] for fast solvers for the latter kind of systems of equations. We note here that since A is a symmetric positive definite and ${\alpha }_i<0$, the matrix $A-{\alpha }_{i}I$ also is a symmetric positive definite. Moreover, when A is obtained by discretizing an elliptic partial differential equation, it may be possible to solve the required linear systems rapidly by a multigrid method.

Algorithm 1 The rational Lanczos process

1: Input: $\mathit{v}\in {\mathbb{R}}^N\backslash \left\{\mathbf{0}\right\}$, a sequence of matrix-valued elementary basis functions ${\psi }_0,{\psi }_1,\dots ,{\psi }_m$, and functions for evaluating matrix-vector products with A and for solving linear systems of equations with matrices of the form $A-{\alpha }_{\ell }I$. Thus, we do not explicitly form the elementary basis functions $A^j$ and ${(A-{\alpha }_{\ell }I)}^{-j}$. The sequence of elementary basis functions implicitly defines the integers j and ℓ used below.

2: Output: Orthonormal basis ${\left\{{\mathit{v}}_r\right\}}_{r=0}^m$.

3: Initialization: ${\mathit{v}}_0:=\mathit{v}/\parallel \mathit{v}\parallel ;r:=1;$

4: while  $r\le m$do

5:      if ${\psi }_r=A^j$ for some $j\in \mathbb{N}$ then

6:          $\mathit{u}:=A{\mathit{v}}_{r-1}$;

7:          for $i=\widehat{r}:r-1$ do

8:              $c_{r-1,i}:={\mathit{v}}_i^T\mathit{u}$; $\mathit{u}:=\mathit{u}-c_{r-1,i}{\mathit{v}}_i$;

9:          end for

10:         ${\delta }_r:=\parallel \mathit{u}\parallel$; ${\mathit{v}}_r:=\mathit{u}/{\delta }_r$;

11:         $r=r+1$

12:      else if ${\psi }_r={(A-{\alpha }_{\ell }I)}^{-j}$ for some $j,\ell \in \mathbb{N}$ then

13:          $\mathit{u}:={(A-{\alpha }_{\ell }I)}^{-1}{\mathit{v}}_{r-1}$;

14:          for $i=\widehat{r}:r-1$ do

15:              $c_{r-1,i}:={\mathit{v}}_i^T\mathit{u}$; $\mathit{u}:=\mathit{u}-c_{r-1,i}{\mathit{v}}_i$;

16:          end for

17:          ${\delta }_r:=\parallel \mathit{u}\parallel$; ${\mathit{v}}_r:=\mathit{u}/{\delta }_r$;

18:          $r=r+1$

19:      end if

20: end while

In Algorithm 1, we assume that the basis (8) satisfies conditions 1 and 2 of natural ordering. The value of $\widehat{r}$ is such that ${\psi }_{\widehat{r}}$ is a basis function with the same poles as ${\psi }_r$, with at least one of the poles of smaller order; if no such basis function ${\psi }_{\widehat{r}}$ exists, then $\widehat{r}=0$. When carrying out m steps of Algorithm 1, we obtain an orthonormal basis $\{{\mathit{v}}_0,{\mathit{v}}_1,\dots ,{\mathit{v}}_m\}$ for the rational Krylov subspace ${\mathbb{K}}_{m+1}(A,\mathit{v})$.

The matrix

$$V_m=[{\mathit{v}}_0,{\mathit{v}}_1,\dots ,{\mathit{v}}_{m-1}]\in {\mathbb{R}}^{N\times m},\mathrm{with}{\mathit{v}}_0=\mathit{v},$$

and the symmetric matrix

$$H_m={\left[h_{i,j}\right]}_{i,j=0}^{m-1}\in {\mathbb{R}}^{m\times m},h_{i,j}=\mathcal{I}\left(y{\varphi }_i{\varphi }_j\right)={\mathit{v}}_i^{T}A{\mathit{v}}_j,$$

satisfy the relation

$$A{V}_m=V_m{H}_m+\sum _{j=1}^d{h}_{m-j,m}{\mathit{v}}_m{\mathit{e}}_{m+1-j}^T,$$

(16)

where the vector ${\mathit{v}}_m\in {\mathbb{R}}^N$ is such that $V_m^T{\mathit{v}}_m=\mathbf{0}$. The subspace ${\mathbb{K}}_i$, $1\le i\le m$, is spanned by the first i columns of the matrix $V_m$. In particular, the rational Krylov subspaces are nested, i.e., ${\mathbb{K}}_1\subseteq {\mathbb{K}}_2\subseteq \dots$. The orthogonal projection of A onto the rational Krylov subspace (15) is given by

$$H_m=V_m^{T}A{V}_m.$$

We will assume that m is small enough so that the decomposition (16) with the stated properties exists. This is the generic situation. Note that the determination of the matrix $H_m$ only requires that $m-1$ steps of Algorithm 1 be carried out.

4. Rational Gauss Quadrature

This section describes properties of rational Gauss quadrature rules for the approximation of functionals (1) when f is a Stieltjes function (2) and A is a symmetric positive definite matrix. Then

$$F\left(A\right)={\mathit{v}}^{T}f\left(A\right)\mathit{v}={\int }_0^{\infty }{\mathit{v}}^T{(tI+A)}^{-1}\mathit{v}d\mu \left(t\right)=:\mathcal{I}\left(f\right).$$

(17)

It is shown in [12] that

$${\widehat{\mathcal{G}}}_m\left(f\right):={\mathit{e}}_1^{T}f\left(H_m\right){\mathit{e}}_1$$

(18)

is a rational Gauss quadrature rule for the approximation of (17). The rational Gauss quadrature rule is characterized by the property

$${\widehat{\mathcal{G}}}_m\left(f\right)=\mathcal{I}\left(f\right),\forall f\in {\mathbb{S}}_{2m},$$

where

$${\mathbb{S}}_{2m}:={\mathbb{P}}_{2m-2k}\oplus {\mathbb{Q}}_{1,2k_1}\oplus \cdots \oplus {\mathbb{Q}}_{\ell ,2k_{\ell }}$$

denotes a $2m$-dimensional linear subspace spanned by the chosen elementary rational functions.

In view of (4), the rational Gauss rule (18) can be expressed as

$${\widehat{\mathcal{G}}}_m\left(f\right)={\int }_0^{\infty }{\mathit{e}}_1^T{(tI+H_m)}^{-1}{\mathit{e}}_{1}d\mu \left(t\right),$$

(19)

where we recall that the measure $d\mu$ defines function f; the rational Gauss rule is with respect to the measure $d\nu$ in (4).

The orthonormal rational function ${\varphi }_m$ has m distinct zeros ${\left\{y_i\right\}}_{i=1}^m$ that lie in the convex hull of the support of the measure $d\nu$, and they are the eigenvalues of $H_m$ in (13); see (Theorem 2.5 in [12]). Recall that we assume ${\psi }_m$ to be a monomial. Thus the rational function ${\varphi }_m$ in (13) can be written as

$${\varphi }_m\left(y\right)=c_m\frac{{\prod }_{i=1}^m(y-y_i)}{w\left(y\right)}\in {\mathbb{S}}_{m+1},$$

where

$$w\left(y\right)=\prod _{i=1}^{\ell }{(y-{\alpha }_i)}^{k_i}$$

(20)

is a polynomial of degree k with k defined by (6). The scalars ${\alpha }_i$ are the poles of the rational approximant; see the beginning of Section 3. Moreover, $c_m$ is a scalar. A remainder term for the rule (19) is derived in [6].

5. Allocation of Poles

The formulas in the previous section are valid independently of the location of the distinct negative real poles ${\alpha }_i$ and their multiplicities. However, their location and multiplicities affect the size of the quadrature error

$${\mathcal{E}}_m\left(f\right):=|F\left(A\right)-{\widehat{\mathcal{G}}}_m\left(f\right)|.$$

(21)

Example 5.3 reported in [6] indicates that when the rational approximant of the Stieltjes function f has more distinct poles, the rational Gauss rule (19) may give a smaller quadrature error (21) than when the rational approximant is of the same order with fewer distinct poles. However, no discussion on the allocation of the poles is provided in [6]. Here we propose to allocate the poles so that, when considered positive point charges, they approximately make the negative real axis an equipotential curve. The allocation of poles that satisfy this property is very easy to determine with the aid of conformal mapping. Computed examples in Section 6 illustrate the performance of this pole allocation technique.

We remark that using many distinct poles may not be advantageous because each distinct pole ${\alpha }_i$ requires the solution of a linear system of equations of the form

$$(A-{\alpha }_{i}I)\mathit{z}=\mathit{w}$$

for some vector w; see the discussion in Section 3. When the systems of equations are solved by factorization of the system matrix, each distinct pole requires the computation of a factorization. It follows that the computing time required to achieve a quadrature error (21) that is smaller than a specified tolerance may be smaller when using a rational approximant with fewer distinct poles than when using a rational approximant of the same order with more distinct poles. The optimal choice of the number of distinct poles depends on the structure of A, i.e., on how demanding the computation of the factorization of the matrix $A-{\alpha }_{i}I$ is. This section considers the allocation of ℓ distinct poles; cf. (5).

We set the poles to be images under a conformal mapping of equidistant points on the unit circle. These poles are easy to compute. The application of conformal mappings is justified by the fact that the approximation errors of analytic functions by rational functions are invariant under conformal mapping because conformal mappings map equipotential curves onto equipotential curves.

Consider the problem of approximating an analytic function outside the unit circle by a rational function with poles on the unit circle. If no particular properties of the analytic function are known, then it is natural to allocate the poles equidistantly on the unit circle; see, e.g., Walsh (Chapter 9 in [13]) for a discussion on rational approximation and equipotential curves. Let the unit circle be in the z-plane. We can map the unit circle to the interval $[-1,1]$ by the Joukowski map

$$w=\frac{1}{2}(z+1/z);$$

see, e.g., Henrici (Chapter 5 in [14]). As z traverses the unit circle, w traverses the interval $[-1,1]$ twice. We now map w to the semi-infinite interval $[-\infty ,\alpha ]$ with $\alpha \le 0$ by the Moebius transformation

$$\zeta =\frac{w-1}{w+1}+\alpha .$$

(22)

Then $\zeta =-\infty$ when $w=-1$, and $\zeta =\alpha$ for $w=1$. When z traverses the top half of the unit circle counterclockwise, starting at $z=1$ and ending at $z=-1$, $\zeta$ traverses the interval $[-\infty ,\alpha ]$ starting at $\alpha$. Similarly, when z traverses the bottom half of the unit circle clockwise, $\zeta$ also traverses the interval $[-\infty ,\alpha ]$ starting at $\alpha$.

To allocate ℓ distinct poles on the interval $[-\infty ,\alpha ]$, we first define the ℓ equidistant points

$$z_j=exp(\pi i(j-1/2)/\ell ),i=\sqrt{-1},j=1,2,\dots ,\ell ,$$

on the upper unit semicircle. The Joukowski map yields the points

$$w_j=\frac{1}{2}(z_j+1/z_j),j=1,2,\dots ,\ell ,$$

on the interval $[-1,1]$. Finally, using (22) gives the distinct poles

$${\alpha }_j=\frac{w_j-1}{w_j+1}+\alpha ,j=1,2,\dots ,\ell ,$$

(23)

on $[-\infty ,\alpha ]$. We will use these poles in the computed examples of the following section.

6. Computed Examples

This section illustrates the application of rational Gauss rules to Stieltjes functions of a symmetric positive definite matrix. We illustrate the pole allocation described in Section 5 and compare it with ad hoc pole allocation. The computations were carried out using MATLAB R2017b on a 64-bit MacBook Pro personal computer; arithmetic was performed with about 15 significant decimal digits. We determine the quadrature error by explicitly evaluating the functionals (1). This limits the size of the matrices A that we can consider.

Example 2.

Consider the approximation of the functional

$$F\left(A\right):={\mathit{v}}^T{A}^{-1/2}\mathit{v},$$

(24)

where $A\in {\mathbb{R}}^{1000\times 1000}$ is a symmetric positive definite Toeplitz matrix with first row $[1,1/2,\dots ,1/1000]$, and $\mathit{v}={[1/\sqrt{1000},\dots ,1/\sqrt{1000}]}^T\in {\mathbb{R}}^{1000}$. The smallest eigenvalue of A is ${\lambda }_1=0.3863$, and the largest one is ${\lambda }_{1000}=12.1259$. The value of $F\left(A\right)$ is about $0.2897$. We compute approximations of (24) by rational Gauss rules. The computations require the solution of linear system of equations with symmetric positive definite Toeplitz matrices $A-{\alpha }_{i}I$, where ${\alpha }_i<0$ are poles. Ammar and Gragg [11] describe fast algorithms for the solution of systems of equations with this kind of matrix.

We will use the rational Krylov subspace

$$\begin{array}{cc}\hfill {\mathbb{K}}_m(A,\mathit{v})=& \mathrm{span}\{\mathit{v},A\mathit{v},{(A-{\alpha }_{1}I)}^{-1}\mathit{v},A^2\mathit{v},\dots ,{(A-{\alpha }_{1}I)}^{-k_1}\mathit{v},A^{k_1+1}\mathit{v},\hfill \\ \hfill & {(A-{\alpha }_{2}I)}^{-1}\mathit{v},A^{k_1+2}\mathit{v},\dots ,{(A-{\alpha }_{2}I)}^{-k_2}\mathit{v},A^{k_1+k_2+1}\mathit{v},\dots ,\hfill \\ \hfill & {(A-{\alpha }_{\ell }I)}^{-1}\mathit{v},\dots ,{(A-{\alpha }_{\ell }I)}^{-k_{\ell }}\mathit{v},A^{k+1}\mathit{v}\},\hfill \end{array}$$

where k is determined by (6).

To calculate the exact value of (24), we evaluate ${\mathit{v}}^T{A}^{-1/2}\mathit{v}$ by first computing the matrix square root $A^{1/2}$ with the MATLAB functionsqrtmand then solving a linear system of equations with the matrix A. The rational Gauss rule is evaluated as

$${\mathit{e}}_1^T{H}_m^{-1/2}{\mathit{e}}_1,$$

where $H_m^{-1/2}{\mathit{e}}_1$ is determined by first computing the matrix $H_m^{1/2}$ using the MATLAB functionsqrtmand then solving a linear system of equations with the matrix $H_m$.

The first row of

Table 1. Example 2: Errors in rational Gauss rules for approximating (24) with A is a symmetric positive definite Toeplitz matrix. For the first row, the distinct poles are ${\alpha }_i\in \{-1,-2\}$, each of multiplicity two. The second row shows the error when the poles are allocated with the aid of conformal mapping using (23) with $\alpha =0$ and $\ell =2$. Then the distinct poles are ${\alpha }_i\in \{-0.1716,-5.8284\}$.

m	${\mathcal{E}}_{\mathit{m}}\left(\mathit{f}\right)$
10	$\phantom{-}3.42\times {10}^{-11}$
10	$\phantom{-}8.19\times {10}^{-13}$

displays the error in the approximation of (24) obtained by rational Gauss rules with the somewhat arbitrarily chosen distinct poles

$${\alpha }_i\in \{-1,-2\}$$

(25)

of multiplicity two. On the second row, Table 1 shows the error in the approximation obtained with rational Gauss rules with poles defined by (23) with $\alpha =0$ and $\ell =2$. Thus, the poles used are of two multiplicity two. The distinct poles are ${\alpha }_1=-0.1716$ and ${\alpha }_2=-5.8284$. Table 1 shows the latter poles to yield a smaller error than the poles (25).

We now turn to the use of four distinct poles. The first row of

Table 2. Example 2: Errors in rational Gauss rules for approximating (24) with A a symmetric positive definite Toeplitz matrix. In the first row, the poles are ${\alpha }_i\in \{-1/2,-1,-3/2,-2\}$ of multiplicity one. The second row shows the error when the poles are allocated by using (23) with $\alpha =0$ and $\ell =4$. The poles are ${\alpha }_i\in \{-0.0396,-0.4465,-2.2398,-25.2741\}$ of multiplicity one.

m	${\mathcal{E}}_{\mathit{m}}\left(\mathit{f}\right)$
10	$\phantom{-}1.15\times {10}^{-11}$
10	$\phantom{-}2.70\times {10}^{-13}$

displays the errors in approximations determined by rational Gauss rules with the arbitrarily chosen distinct poles

$${\alpha }_i\in \{-1/2,-1,-3/2,-2\}$$

(26)

of multiplicity one. Row two of Table 2 displays the approximation of the error obtained with the four poles ${\alpha }_i\in \{-0.0396,-0.4465,-2.2398,-25.2741\}$ of multiplicity one, which is defined by (23) with $\alpha =0$ and $\ell =4$. The latter poles can be seen to give a smaller approximation error than the equidistant poles (26).

Example 3.

We would like to calculate an approximation of the functional

$$F\left(A\right):={\mathit{v}}^{T}log(A+I)A^{-1}\mathit{v},$$

(27)

where $A\in {\mathbb{R}}^{1000\times 1000}$ is a symmetric positive definite Toeplitz matrix with a first row $[3,3/2,\dots ,3/1000]$. The vector $\mathit{v}\in {\mathbb{R}}^{1000}$ is the same as in Example 1. The smallest eigenvalue of A is ${\lambda }_1=1.1589$, and the largest one is ${\lambda }_{1000}=36.3776$. The value of $F\left(A\right)$ is approximately $0.1009$.

The first row of

Table 3. Example 3: Errors in rational Gauss rules for the approximating of (27) with A a symmetric positive definite Toeplitz matrix. In the first row, the poles are ${\alpha }_1=0$ of multiplicity two and ${\alpha }_2=-1/4$ of multiplicity one. In the second row, the poles are allocated using (23) with $\alpha =-1$ and $\ell =2$. The pole ${\alpha }_1=-1.1716$ is of multiplicity two, and the pole ${\alpha }_2=-6.8284$ is of multiplicity one.

m	${\mathcal{E}}_{\mathit{m}}\left(\mathit{f}\right)$
8	$\phantom{-}6.66\times {10}^{-11}$
8	$\phantom{-}4.09\times {10}^{-13}$

shows the error for the rational Gauss quadrature rule with somewhat arbitrarily chosen poles ${\alpha }_1=0$ of multiplicity two and ${\alpha }_2=-1/4$ of multiplicity one. The second row displays the error for the rational Gauss rule with poles ${\alpha }_1=-1.1716$ of multiplicity two and ${\alpha }_2=-6.8284$ of multiplicity one. The latter poles are determined by (23) with $\alpha =-1$ and $\ell =2$. We observe that the approach of choosing poles with the aid of conformal mapping provides a more accurate result.

The first and second rows of

Table 4. Example 3: Errors in rational Gauss rules for approximating (27) with A a symmetric positive definite Toeplitz matrix. In the first row, the poles are ${\alpha }_1=0$ and ${\alpha }_2=-1$ of multiplicity two. In the second row, the poles are allocated with (23) with $\alpha =-1$ and $\ell =2$. This gives the distinct poles ${\alpha }_1=-1.1716$ and ${\alpha }_2=-6.8284$ of multiplicity two.

m	${\mathcal{E}}_{\mathit{m}}\left(\mathit{f}\right)$
10	$\phantom{-}1.60\times {10}^{-13}$
10	$\phantom{-}1.29\times {10}^{-15}$

display the errors in rational Gauss rules with somewhat arbitrarily chosen poles ${\alpha }_i\in \{0,-1\}$ of multiplicity two and with poles ${\alpha }_i\in \{-1.1716,-6.8284\}$ of multiplicity two, respectively. The latter poles are defined by (23) with $\alpha =-1$ and $\ell =2$. It can be seen that the error in the second row is smaller than the error in the first row.

Example 4.

In this example, we determine an approximation of

$$F\left(A\right):={\mathit{v}}^T\left(\pi {(I+\sqrt{A})}^{-1}\right)\mathit{v},$$

(28)

where the matrix $A\in {\mathbb{R}}^{1000\times 1000}$ and vector $\mathit{v}\in {\mathbb{R}}^{1000}$ are the same as in Example 1. The value of $F\left(A\right)$ is approximately $0.7053$. We note that

$$f\left(z\right)=\frac{\pi }{1+\sqrt{z}}={\int }_0^{\infty }\frac{1}{t+z}d\mu \left(t\right),d\mu \left(t\right)=\frac{\sqrt{t}}{1+t}dt,$$

is a Stieltjes function.

Table 5. Example 4: Errors in rational Gauss rules for computing approximations of (28) when A is a symmetric positive definite Toeplitz matrix. In the first row, we choose four equidistant poles ${\alpha }_i\in \{0,-2,-4,-6\}$ of multiplicity one. The second row displays the approximation error when the poles are allocated by (23) with $\alpha =0$ and $\ell =4$; the poles are ${\alpha }_i\in \{-0.0396,-0.4465,-2.2398,-25.2741\}$ of multiplicity one. The poles on the third row are obtained by using (23) with $\alpha =0$ and $\ell =2$. Then the distinct poles are ${\alpha }_i\in \{-0.1716,-5.8284\}$ of multiplicity two.

m	${\mathcal{E}}_{\mathit{m}}\left(\mathit{f}\right)$
10	$\phantom{-}2.49\times {10}^{-12}$
10	$\phantom{-}1.01\times {10}^{-13}$
10	$\phantom{-}2.68\times {10}^{-13}$

displays the errors in approximations determined by rational Gauss rules. In the first row, we choose four somewhat arbitrary equidistant poles, ${\alpha }_i\in \{0,-2,-4,-6\}$, each of multiplicity one. In the second row, the poles are ${\alpha }_i\in \{-0.0396,-0.4465,-2.2398,-25.2741\}$ of multiplicity one, and the poles in the third row are ${\alpha }_i\in \{-0.1716,-5.8284\}$ of multiplicity two. The poles in the second and third rows are allocated using (23). It can be seen that the errors in the second and third rows are of a smaller magnitude than the error in the first row.

7. Conclusions

This paper discusses the evaluation of expressions of the form (1), where A is a symmetric positive definite matrix and f is a Stieltjes function. We approximate f by a rational function with user-specified poles and integrate the rational function by a rational Gauss quadrature rule. The rational function is expressed as a linear combination of rational orthogonal functions, which can be generated with recursion formulas with few terms. The paper focuses on the choice of poles of the rational function. They are allocated on the negative real axis and chosen to be images of equidistant points on the upper unit semicircle under a suitable conformal mapping. The computed examples show this approach of allocating the poles to give a smaller quadrature error than poles chosen in an ad hoc manner.