On symmetrizing the ultraspherical spectral method for self-adjoint problems

TL;DR

This paper introduces a symmetrization technique for the ultraspherical spectral method applied to self-adjoint problems, resulting in symmetric, banded discretizations and an adaptive spectral decomposition algorithm.

Contribution

It presents a novel symmetrization mechanism and an adaptive spectral decomposition algorithm for self-adjoint operators within the ultraspherical spectral framework.

Findings

Discretizations become symmetric and banded after symmetrization

The adaptive spectral decomposition effectively analyzes self-adjoint operators

Applications demonstrate the advantages of the symmetrizer and decomposition methods

Abstract

A mechanism is described to symmetrize the ultraspherical spectral method for self-adjoint problems. The resulting discretizations are symmetric and banded. An algorithm is presented for an adaptive spectral decomposition of self-adjoint operators. Several applications are explored to demonstrate the properties of the symmetrizer and the adaptive spectral decomposition.

Figures

Tables

Table 1. Table 1: Driscoll’s  [ 3 ] first four computed generalized eigenvalues of Eq. ( 47 ) with constant coefficients.

$n$	$\sqrt{{\lambda }_n^{\mathrm{Driscoll}}}$		${\lambda }_n^{\mathrm{Driscoll}}$
$0$	2.	266 742 077	5.	138 119 644
$1$	6.	923 297 244	47.	932 044 73
$2$	13.	977 668 83	195.	375 225 9
$3$	22.	819 562 61	520.	732 437 7

Table 2. Table 2: The first four computed eigenvalues of Eq. ( 47 ) in 256 256 256 -bit extended precision.

$n$	${\lambda }_n$
$0$	5.	138 119 642 550 555 649 415 412 286 723 460 005 700 130 780 549 503 937 354 845 212 533 811 …
$1$	47.	932 044 726 525 070 683 741 246 754 121 353 113 209 603 021 890 092 332 934 543 893 845 719 …
$2$	195.	375 225 987 311 298 618 562 101 588 390 534 216 626 466 473 664 548 029 401 788 843 545 607 …
$3$	520.	732 437 610 822 993 012 167 361 990 447 367 561 584 023 509 960 784 102 391 822 269 243 419 …

Equations

$$\mathcal{L}u:=\left[(-\mathcal{D}{)}^N(p_N{\mathcal{D}}^N)+(-\mathcal{D}{)}^{N-1}(p_{N-1}{\mathcal{D}}^{N-1})+\cdots +p_0\right]u=\lambda wu,\mathcal{B}u=0,$$

$$\mathcal{L}u=f,\mathrm{with}\mathcal{B}u=c,$$

$$\left(\begin{array}{c}\mathcal{B}\\ \mathcal{L}\end{array}\right)u=\lambda \left(\begin{array}{c}0\\ \mathcal{C}\end{array}\right)u,$$

$$-{\mathcal{D}}^{2}u=\lambda u,u(\pm 1)=0.$$

$$u(x)=\sum _{n=0}^{\infty }{u}_n(1-x^2){\tilde{P}}_n^{(1,1)}(x).$$

$${\tilde{P}}_n^{(\alpha ,\beta )}(x):=\sqrt{\frac{(2n+\alpha +\beta +1)\Gamma (n+\alpha +\beta +1)n!}{2^{\alpha +\beta +1}\Gamma (n+\alpha +1)\Gamma (n+\beta +1)}}{P}_n^{(\alpha ,\beta )}(x).$$

$$\mathcal{D}\left(\begin{array}{ccc}(1-x^2){\tilde{P}}_0^{(1,1)}(x)& (1-x^2){\tilde{P}}_1^{(1,1)}(x)& \cdots \end{array}\right)=\left(\begin{array}{ccc}{\tilde{P}}_0^{(0,0)}(x)& {\tilde{P}}_1^{(0,0)}(x)& \cdots \end{array}\right)\left(\begin{array}{ccccc}0& -\sqrt{2}& & & \\ & & -\sqrt{6}& & \\ & & & -\sqrt{12}& \\ & & & & \ddots \end{array}\right),$$

$$-\mathcal{D}\left(\begin{array}{ccc}{\tilde{P}}_0^{(0,0)}(x)& {\tilde{P}}_1^{(0,0)}(x)& \cdots \end{array}\right)=\left(\begin{array}{ccc}{\tilde{P}}_0^{(1,1)}(x)& {\tilde{P}}_1^{(1,1)}(x)& \cdots \end{array}\right)\left(\begin{array}{cccc}0& & & \\ -\sqrt{2}& & & \\ & -\sqrt{6}& & \\ & & -\sqrt{12}& \\ & & & \ddots \end{array}\right).$$

$$-{\mathcal{D}}^2=\left(\begin{array}{cccc}{d}_0& & & \\ & d_1& & \\ & & d_2& \\ & & & \ddots \end{array}\right),\mathrm{where}{d}_n=(n+1)(n+2).$$

$$\mathcal{M}[1-x^2]=\left(\begin{array}{cccc}{a}_0& 0& b_0& \\ 0& a_1& & \ddots \\ b_0& & a_2& \\ & \ddots & & \ddots \end{array}\right),\mathrm{where}{a}_n=\frac{2(n+1)(n+2)}{(2n+1)(2n+5)},\mathrm{and}{b}_n=-\sqrt{\frac{(n+1)(n+2)(n+3)(n+4)}{(2n+3)(2n+5{)}^2(2n+7)}}.$$

$$\left(\begin{array}{ccc}{\rho }_0& {\rho }_1& \cdots \end{array}\right)=\left(\begin{array}{ccc}{\varphi }_0& {\varphi }_1& \cdots \end{array}\right)A.$$

$$\left(\begin{array}{ccc}{\psi }_0& {\psi }_1& \cdots \end{array}\right)=\left(\begin{array}{ccc}{\varphi }_0& {\varphi }_1& \cdots \end{array}\right)Q,$$

$$\left(\begin{array}{ccc}{\rho }_0& {\rho }_1& \cdots \end{array}\right)=\left(\begin{array}{ccc}{\psi }_0& {\psi }_1& \cdots \end{array}\right)R.$$

$$\left(\begin{array}{ccc}{\sigma }_0& {\sigma }_1& \cdots \end{array}\right)R^{\ast }=\left(\begin{array}{ccc}{\psi }_0& {\psi }_1& \cdots \end{array}\right).$$

$$\left(\begin{array}{ccc}{\psi }_0& {\psi }_1& \cdots \end{array}\right)Q^{\ast }=\left(\begin{array}{ccc}{\varphi }_0& {\varphi }_1& \cdots \end{array}\right),$$

$$\left(\begin{array}{ccc}{\sigma }_0& {\sigma }_1& \cdots \end{array}\right)A^{\ast }=\left(\begin{array}{ccc}{\varphi }_0& {\varphi }_1& \cdots \end{array}\right).$$

$${\rho }_n(x)=\sqrt{\frac{(n+1)(n+2)}{(2n+1)(2n+3)}}{\tilde{P}}_n(x)-\sqrt{\frac{(n+1)(n+2)}{(2n+3)(2n+5)}}{\tilde{P}}_{n+2}(x)=(1-x^2){\tilde{P}}_n^{(1,1)}(x),$$

$$u(x)=\sum _{n=0}^{\infty }{\psi }_n(x)u_n=\left(\begin{array}{ccc}{\psi }_0& {\psi }_1& \cdots \end{array}\right)\left(\begin{array}{c}{u}_0\\ u_1\\ ⋮\end{array}\right)={\mathbf{\psi }}^{\top }\mathbf{u}.$$

$$\mathcal{L}{\mathbf{\psi }}^{\top }\mathbf{u}={\mathbf{\psi }}^{\top }L\mathbf{u}=\lambda w{\mathbf{\psi }}^{\top }\mathbf{u}=\lambda {\mathbf{\psi }}^{\top }M\mathbf{u},$$

$${\mathbf{\sigma }}^{\top }{R}^{\ast }LR\mathbf{v}=\lambda {\mathbf{\sigma }}^{\top }{R}^{\ast }MR\mathbf{v}.$$

$$\mathcal{L}{\mathbf{\rho }}^{\top }\mathbf{v}=\lambda w{\mathbf{\rho }}^{\top }\mathbf{v}.$$

$${\mathbf{\sigma }}^{\top }{A}^{\ast }{C}^{-1}{L}_{B}A\mathbf{v}=\lambda {\mathbf{\sigma }}^{\top }{A}^{\ast }{M}_{B}A\mathbf{v}.$$

$${\Vert \mathbf{u}\Vert }_{{\ell }_{\lambda }^2}^2=\sum _{n=0}^{\infty }{\left|u_n\right|}^2(n+1{)}^{2\lambda }.$$

$${\Vert A\Vert }_{{\ell }^2\to {\ell }^2}:=\underset{{\Vert \mathbf{x}\Vert }_{{\ell }^2}\le 1}{\sup }{\Vert A\mathbf{x}\Vert }_{{\ell }^2},$$

$$\mathcal{B}=\left(\begin{array}{cccc}{b}_{0,0}& b_{0,1}& b_{0,2}& \cdots \\ b_{1,0}& b_{1,1}& b_{1,2}& \cdots \\ ⋮& ⋮& ⋮& \cdots \\ b_{2N-1,0}& b_{2N-1,1}& b_{2N-1,2}& \cdots \end{array}\right)\in {\mathbb{C}}^{2N\times \infty },$$

$$\left(\begin{array}{cccc}{b}_{0,k+1}& b_{0,k+2}& \cdots & b_{0,k+2N}\\ b_{1,k+1}& b_{1,k+2}& \cdots & b_{1,k+2N}\\ ⋮& ⋮& \ddots & ⋮\\ b_{2N-1,k+1}& b_{2N-1,k+2}& \cdots & b_{2N-1,k+2N}\end{array}\right)\left(\begin{array}{c}{a}_{k+1,k}\\ a_{k+2,k}\\ ⋮\\ a_{k+2N,k}\end{array}\right)=-\left(\begin{array}{c}{b}_{0,k}\\ b_{1,k}\\ ⋮\\ b_{2N-1,k}\end{array}\right).$$

$$c_{i,0}(n+1{)}^{{\lambda }_i}\le \left|b_{i,n}\right|\le c_{i,1}(n+1{)}^{{\lambda }_i},i=0,1,$$

$$\left(\begin{array}{c}{a}_{k+1,k}\\ \text{leavevmode}\\ a_{k+2,k}\end{array}\right)=\left(\begin{array}{c}\frac{b_{0,k+2}{b}_{1,k}-b_{1,k+2}{b}_{0,k}}{b_{0,k+1}{b}_{1,k+2}-b_{0,k+2}{b}_{1,k+1}}\\ \text{leavevmode}\\ \frac{b_{1,k+1}{b}_{0,k}-b_{0,k+1}{b}_{1,k}}{b_{0,k+1}{b}_{1,k+2}-b_{0,k+2}{b}_{1,k+1}}\end{array}\right).$$

$$\left|b_{0,k+1}{b}_{1,k+2}-b_{0,k+2}{b}_{1,k+1}\right|$$

$$\ge c_{0,0}{c}_{1,0}\left[(k+2{)}^{{\lambda }_0}(k+3{)}^{{\lambda }_1}+(k+3{)}^{{\lambda }_0}(k+2{)}^{{\lambda }_1}\right],$$

$${\Vert A\Vert }_{{\ell }^2\to {\ell }^2}\le \underset{k\in {\mathbb{N}}_0}{\sup }\max \left(1+\left|a_{k+1,k}\right|+\left|a_{k+2,k}\right|,\left|a_{k+2,k}\right|+\left|a_{k+2,k+1}\right|+1\right)<\infty .$$

Full text

Keywords: spectral methods, symmetric-definite and banded eigenvalue problems, infinite-dimensional linear algebra.

1 Introduction

Let $H$ be a separable Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and norm $\left\|{\cdot}\right\|$. We are interested in the solution of the self-adjoint linear differential eigenvalue problem with self-adjoint boundary conditions posed on a subset of the real line, $D\subset\mathbb{R}$. The $2N^{\rm th}$ order problem in self-adjoint form reads:

[TABLE]

where the usual conditions on the coefficients $p_{i}$ and $w$ are required for sensible solution; at the very least, $p_{N}\neq 0$, and $w\geq 0$ and $w=0$ only on a set of Lebesgue measure zero. In this form, $\mathcal{L}$ is a self-adjoint linear differential operator in $H=L^{2}(D,w(x){\rm\,d}x)$ provided it is accompanied by $\mathcal{B}:H\to\mathbb{C}^{2N}$, a set of $2N$ self-adjoint boundary conditions. We will also assume that $p_{i}\in\mathbb{P}_{m}$, that is, the variable coefficients are at most degree-$m$ polynomials. The self-adjoint linear differential eigenvalue problem is an important problem in applied mathematics, with many scientific contributions furthering our collective understanding, too many to all be cited.

Given the quotient Hilbert space $H_{\mathcal{B}}:=\{u\in H:\mathcal{B}u=0\}$, a classical approach calls for a Ritz–Galerkin discretization for the trial and test functions taken from an $n$-dimensional subspace $V_{n}\subset H_{\mathcal{B}}$. Indeed, when $V_{n}=\mathbb{P}_{n-1}$ and an orthonormal polynomial basis satisfying the boundary conditions is utilized, the discrete realization of $\langle u,\mathcal{L}v\rangle$ is automatically symmetric (hermitian). Unfortunately, the resulting discretizations are dense. Thus data-sparsity in problem formulations, such as Eq. (1) with finite-degree polynomial coefficients, is ignored when transformed into a computational problem, resulting in $\mathcal{O}(n^{3})$ complexity algorithms [1] for approximate eigenvalues and eigenfunctions.

Another classical approach is a pseudospectral method implemented by collocation [2, 3]. Apart from the boundary conditions, differentiation matrices resulting from spectral collocation have many desirable properties such as (skew) centrosymmetry. The traditional wisdom is that, symmetry notwithstanding, even-ordered spectral differentiation matrices have real and distinct eigenvalues [4]. One of the main draws to collocation is the ease of implementation; variable coefficients are implemented as diagonal scalings. Some formulations are even symmetric for Sturm–Liouville problems [5]. Unfortunately, the discretizations are also dense and join the Ritz–Galerkin discretization behind the $\mathcal{O}(n^{3})$ computational barrier.

A contemporary approach to the solution of linear differential equations with polynomial coefficients is the Olver–Townsend ultraspherical spectral method [6]. Compared with the finite-dimensional classical schemes, it is an infinite-dimensional Petrov–Galerkin scheme that utilizes multiple ultraspherical polynomial bases to represent a linear differential operator such as the left-hand side of Eq. (1) as sparsely as possible. This spectral method is based on the facts that: differentiation of Chebyshev polynomials of the first kind is diagonal when represented in terms of ultraspherical polynomials $C_{n}^{(\lambda)}$ of parameter $\lambda=1$; subsequent differentiations are also diagonal when represented in higher order ultraspherical polynomial bases; conversion from a representation in low order ultraspherical polynomials to higher order ones is banded; and, polynomial coefficient multiplication is banded. Crucially, the bandedness implies that for linear differential equations of the form:

[TABLE]

representations of the solution are achieved in $\mathcal{O}(m^{2}n)$ complexity, and the degree $n$ may be determined adaptively through a running estimate of the forward error without restarting the partial $QR$ factorization of $\mathcal{L}$. When posing Eq. (1) as an eigenvalue problem, one approach in adaptations of the ultraspherical spectral method to polar [7] and spherical [8] geometries is to impose:

[TABLE]

where $\mathcal{C}$ is the operator that would convert the representation of the solution in the ultraspherical polynomial basis of the domain of $\mathcal{L}$ to the basis for its range. Unfortunately, even for a self-adjoint linear differential eigenvalue problem with self-adjoint boundary conditions, symmetry appears to be sacrificed for the banded sparsity111In exceptional circumstances, such as when $\mathcal{B}$ is zero apart from the first few columns and the nonzero part prepended to $\mathcal{L}$ renders the block operator discretization self-adjoint, symmetry may still be within reach. Many common boundary conditions, however, are genuinely infinite-dimensional. This destruction of symmetry has real consequences for the spectra of finite truncations. Whereas Chebyshev and Legendre second-order differentiation matrices by collocation have real spectra [9], it is easy to show by counter-example that finite truncations of ultraspherical discretizations of self-adjoint operators may produce complex spectra..

In this work, we describe a procedure to formulate self-adjoint problems such as Eq. (1) in order that symmetry may transcend the discretization and that data-sparsity in the problem is preserved in the form of a banded discretization. Thus, we describe a procedure to sparsify the Ritz–Galerkin scheme or to symmetrize the ultraspherical spectral method. With symmetric-definite banded spectral discretizations, $\mathcal{O}(mn^{2})$ complexity partial diagonalization is possible through finite-dimensional truncation and numerical linear algebra. We also present an adaptive infinite-dimensional spectral decomposition of self-adjoint linear differential operators that deflates resolved eigenpairs on-the-fly and for certain problems, achieves linear complexity partial spectral decompositions. There are many advantages of preserving symmetry. Symmetry ensures that the general Rayleigh quotient iteration [10, 11] converges cubically rather than quadratically. And for time-stepping schemes, symmetry precludes rogue eigenvalues from being scattered spuriously in the complex plane. There are also advantages of preserving the data-sparsity via a small bandwidth. This enables spectral analysis of infinite-dimensional discretizations via, e.g. Gerschgorin disks [1] and its generalizations [12]. One disadvantage of the ultraspherical spectral method is that it applies to the strong formulation of the problem, preventing the resolution of weak eigenfunctions unless they are as smooth as the strong ones.

There are a handful of spectral schemes that lead to symmetric-definite and banded discretizations in the literature. These include the Fourier spectral method with periodic boundary conditions [13] and the Hermite function and Tanh–Chebyshev spectral methods [14] as two examples of many orthogonal systems on $L^{2}(\mathbb{R})$. There are also examples with nontrivial boundary conditions as well. In [15], Shen finds symmetric and banded spectral Legendre-based discretizations of constant coefficient second- and fourth-order problems. He repeats the procedure with Chebyshev-based discretizations [16] and while the efficiency is comparable, the symmetry is lost. Banded symmetry is also common in the finite element method [17]. Through examples, we will show precisely when we should anticipate classical and other orthogonal polynomial-based discretizations for self-adjoint problems lead to self-adjoint and banded discretizations, even on non-compact domains and with piecewise-defined bases.

Polynomial approximation theory suggests that we should seek to numerically represent the variable coefficients in Eq. (1) as Chebyshev polynomial expansions [18]. An adaptation of the fast multipole method (FMM) [19] accelerates the Chebyshev–Legendre [20], ultraspherical–ultraspherical [21], and Jacobi–Jacobi [22, §3.3] connection problems to linear complexity in the degree. Provided the Jacobi parameters are not too large, the acceleration of the connection problem enables Jacobi polynomial expansions with nearly the same rapidity as Chebyshev polynomial expansions. Thus in the problems in this paper, the variable coefficients in Eq. (1) are in fact polynomials, but in practice they could very well be good numerical approximations to functions in a particular space. How such approximations, in particular of the leading coefficient $p_{N}$, affect the spectrum is beyond the scope of this report.

2 A well-known model problem

It is instructive to begin with a well-known model problem for illustrative purposes.

Consider the self-adjoint Sturm–Liouville problem [23]:

[TABLE]

Suppose we represent eigenfunctions in terms of the weighted Jacobi polynomial expansion:

[TABLE]

Here, we use the tilde to denote orthonormalization of the orthogonal polynomials with respect to their orthogonality measure. In particular, the normalized Jacobi polynomials are related to the standard normalization through [24, §18.3]:

[TABLE]

It is clear that our representation satisfies the Dirichlet boundary conditions. Furthermore, differentiation follows from:

[TABLE]

and (negative) second differentiation from:

[TABLE]

Ultimately, our representation of the negative second derivative is diagonal, when interpreted in the appropriate bases:

[TABLE]

A representation of the proposed eigenfunction in the basis of the expansion of $-\mathcal{D}^{2}u$ can be realized through the symmetric positive-definite pentadiagonal multiplication operator:

[TABLE]

In this problem, a symmetric-definite banded Petrov–Galerkin discretization is nearly immediate. The key element is choosing an appropriate polynomial basis for $H_{\mathcal{B}}$, and the discretization unfolds naturally.

This problem is well-known and difficult to generalize with the transcendent symmetry because the Dirichlet conditions are so easily enforced by weighting a reasonable polynomial basis by $1-x^{2}$ that the underlying mechanisms are hidden. To understand the mechanisms at work, we will deconstruct the discretization in abstract terms.

3 Polynomial bases satisfying the boundary conditions

Let $\{\phi_{n}\}_{n=0}^{\infty}$ be an orthonormal polynomial basis for $H$ with $\deg(\phi_{n})=n$. Basis recombination [15, 16, 25, 26, 27, 28, 29] is the idea of recombining a usually contiguous subset of a basis to satisfy additional constraints.

For linearly independent boundary conditions in $\mathcal{B}$, let $\{\rho_{n}\}_{n=0}^{\infty}$ be the recombinations of the orthonormal basis $\{\phi_{n}\}_{n=0}^{\infty}$ that annihilate $\mathcal{B}$, guaranteeing $\mathcal{B}\rho_{n}=0$. Note that $\deg(\rho_{n})$ is not necessarily $n$ as in the previous problem, but we will assume that $\rho_{n}$ are nevertheless degree-graded in the generalized sense that there exists a positive integer $\nu$ such that $\deg(\rho_{n+\nu})>\deg(\rho_{n})$ is always true. This assumption suggests the existence of an infinite-dimensional lower-triangular and banded conversion such that:

[TABLE]

With the conversion in hand, a banded $QR$ factorization reveals the isometry mapping the orthonormal polynomial basis for $H$ to the orthonormal polynomial basis $\{\psi_{n}\}_{n=0}^{\infty}$ for $H_{\mathcal{B}}$:

[TABLE]

and the auxiliary resulting conversion:

[TABLE]

We note that a banded $QR$ factorization of $A\in\mathbb{C}^{(n+2N)\times n}$ requires $\mathcal{O}(n)$ operations and either factor may be applied to a vector with the same cost. This shows that all three of the aforementioned conversions run in linear complexity without the effort of a dense Gram–Schmidt procedure.

Finally, define the polynomials $\{\sigma_{n}\}_{n=0}^{\infty}$ by:

[TABLE]

And since:

[TABLE]

we gather all three adjoint conversions:

[TABLE]

In order to keep the discussion concrete, we relate the so-defined polynomials to the simple canonical problem of § 2; we identify the orthonormal polynomial basis for $H=L^{2}([-1,1])$ as the normalized Legendre polynomials, $\phi_{n}(x)=\tilde{P}_{n}(x)$. The recombination to satisfy the Dirichlet boundary conditions are in fact weighted Jacobi polynomials:

[TABLE]

by [24, §18.7.2, 18.9.8], and furthermore an orthonormal basis for $H_{\mathcal{B}}$ is $\psi_{n}(x)=(1-x^{2})\tilde{P}_{n}^{(2,2)}(x)$. We may also identify the auxiliary polynomials $\sigma_{n}(x)=\tilde{P}_{n}^{(1,1)}(x)$, though this identification is serendipitous rather than necessary, and in the general scenario these polynomials are almost never identifiable as classical orthogonal polynomials.

To simplify the notation, we will use a shorthand for an expansion in a basis and a vector of coefficients:

[TABLE]

We will now prove that the Petrov–Galerkin scheme for Eq. (1) is symmetric-definite and banded, where the solution is represented by $u(x)=\boldsymbol{\rho}^{\top}\boldsymbol{v}$ and the residuals are estimated using $\boldsymbol{\sigma}$.

3.1 The Petrov–Galerkin scheme is self-adjoint

Having identified that $\boldsymbol{\psi}$ form an orthonormal polynomial basis for $H_{\mathcal{B}}$, it follows that the Ritz–Galerkin scheme is self-adjoint:

[TABLE]

where $L=L^{*}\in\mathbb{C}^{\infty\times\infty}$ and $M=M^{*}\in\mathbb{C}^{\infty\times\infty}$ is positive-definite. We let $\boldsymbol{u}=R\boldsymbol{v}$ and we use Eq. (14) to represent $\boldsymbol{\psi}$ in terms of $\boldsymbol{\sigma}$. Then we arrive at the Petrov–Galerkin scheme:

[TABLE]

Clearly, both $R^{*}LR$ and $R^{*}MR$ are self-adjoint and the latter is also positive-definite.

3.2 The Petrov–Galerkin scheme is banded

Let us consider the Petrov–Galerkin scheme:

[TABLE]

Combining Eqs. (11) and (16), we know that $w\boldsymbol{\rho}^{\top}=w\boldsymbol{\phi}^{\top}A=\boldsymbol{\phi}^{\top}M_{B}A=\boldsymbol{\sigma}^{\top}A^{*}M_{B}A$, where $M_{B}=M_{B}^{*}$ is also positive-definite and banded by executing the Clenshaw–Smith algorithm with the orthonormal polynomials’ Jacobi operator [30, §4.1]. Similarly, $\mathcal{L}\boldsymbol{\rho}^{\top}=\mathcal{L}\boldsymbol{\phi}^{\top}A$, and by the ultraspherical spectral method [6], there exists a banded discretization of $\mathcal{L}$, namely $L_{B}$, mapping (normalized) Legendre polynomials to (normalized) Jacobi polynomials with integer parameters in proportion to the highest derivatives. Applying the inverse (banded) conversion operator, and further utilizing $\boldsymbol{\phi}^{\top}=\boldsymbol{\sigma}^{\top}A^{*}$, we arrive at:

[TABLE]

At first glance, this formula does not appear to guarantee a banded discretization, for only in special circumstances are the inverses of banded matrices banded themselves [31]. To show that it is indeed banded, we note that the conversion matrix $C$ is upper-triangular and banded, which means that its inverse will not extend the bandwidth of $L_{B}A$ any lower. Since $A^{*}$ is also upper-triangular and since the Petrov–Galerkin scheme is self-adjoint, we conclude that it must be banded above as well.

Remark 3.1**.**

This proof does not imply that $C^{-1}L_{B}$ is self-adjoint; this operator may be strictly upper-triangular as in the model problem! But the composition with $A^{*}$ and $A$ on either side ensures the final scheme is self-adjoint. 2. 2.

The choice of normalization of the polynomial basis that satisfies the boundary conditions amounts to a congruence transformation, therefore does not affect the symmetry. On the other hand, the operator $C^{-1}L_{B}$ must be posed on the orthonormal basis rather than any other normalization.

3.3 Analysis of the conversions $A$ and $R$

The orthogonal polynomial basis for $H_{\mathcal{B}}$ has already been described by Livermore [28], who calls them Galerkin orthogonal polynomials. Livermore presents analytic and asymptotic results for the polynomials to suggest that they are reasonable bases for the constrained Hilbert space, as they mimic many of the properties of the original orthogonal polynomials for $H$.

Here, our analysis of the basis recombination is more abstract.

Definition 3.2**.**

Let $\ell_{\lambda}^{2}\subset\mathbb{C}^{\infty}$ be the discrete Sobolev space equipped with the norm:

[TABLE]

When $\lambda=0$, it is convenient to drop it from the notation.

Colloquially, we will refer to the conversion matrices $A$ as annihilators in relation to $\mathcal{B}A=0$. For infinite-dimensional annihilators $A\in\mathbb{C}^{\infty\times\infty}$ from Eq. (11), we wish to show that the induced $\ell^{2}$-norm:

[TABLE]

is bounded. For if $A:\ell^{2}\to\ell^{2}$ is a bounded operator, then the new polynomials $\{\rho_{n}\}$ form a basis for the constrained Hilbert space $H_{\mathcal{B}}$, and furthermore the $QR$ factorization exists in the strong operator topology [32, Theorem 31]. This argument does not prove that the annihilator is numerically well-conditioned, as this also depends on a bound for the norm of the inverse of the principal finite-dimensional section of $A$.

In general:

[TABLE]

and the matrix $A$ is uniquely determined if it is unit lower-triangular and banded and if it annihilates $\mathcal{B}$. Thus the nontrivial part of the $k^{\rm th}$ column of $A$ must satisfy:

[TABLE]

A proof of the existence and uniqueness of a solution to Eq. (26) as well as an upper bound on $\left\|{a_{:,k}}\right\|_{\ell^{2}}$ would require a bound on the norm of the inverse of every matrix defined by Eq. (26). Below, we prove a simple result for practical applications.

Theorem 3.3**.**

Let $\mathcal{B}\in\mathbb{R}^{2\times\infty}$ be a pair of real linear functionals such that:

[TABLE]

where $0\leq\lambda_{0},\lambda_{1}<\infty$ and $0<c_{0,0},c_{0,1},c_{1,0},c_{1,1}<\infty$, and furthermore assume that the coefficients of one functional are positive and the other alternate in sign. Then $\left\|{A}\right\|_{\ell^{2}\to\ell^{2}}<\infty$ and the infinite-dimensional $QR$ factorization exists.

Proof.

The solution of the linear system in Eq. (26) is:

[TABLE]

As one of the functionals is alternating in sign:

[TABLE]

and it follows that for every $\lambda_{0}$ and $\lambda_{1}$, both sub-diagonal entries are absolutely bounded, independent of $k$. By an estimate on the singular values of a matrix [33]:

[TABLE]

Since $A$ is bounded, the $QR$ factorization exists and $\left\|{A}\right\|_{\ell^{2}\to\ell^{2}}=\left\|{R}\right\|_{\ell^{2}\to\ell^{2}}<\infty$. ∎

For concreteness, consider the case of separated boundary conditions on the interval $[-1,1]$. Legendre polynomials are positive at the right endpoint and alternate in sign at the other [24, §18.6]. Since they and their derivatives at the endpoints are bounded above and below by algebraic powers of the degree, a pair of separated boundary conditions of the form $a_{\pm}u(\pm 1)+b_{\pm}u^{\prime}(\pm 1)=0$, with the conditions $a_{+}b_{+}>0$ and $a_{-}b_{-}<0$, produce bounded conversion operators $A$ and $R$.

An alternative to contiguous basis recombination of the form $\phi_{n}+a_{n+1,n}\phi_{n+1}+a_{n+2,n}\phi_{n+2}$ is the recombination with only the first two elements, $\phi_{n}+\hat{a}_{0,n}\phi_{0}+\hat{a}_{1,n}\phi_{1}$, which has also been discussed elsewhere [15, 28]. Not only are such recombinations not banded, but they are also not bounded in the above scenario, which demonstrates the superiority of the contiguous approach.

3.4 Implementation details

The solution of the linear system in Eq. (26) is performed by a pivoted $LU$ factorization. The number of boundary conditions is usually quite small, say $N\leq 10$, so that instability is rarely an issue in practice. Although we anticipate algebraic growth in the functionals that realize the boundary conditions, if $b_{2N-1,n}\leq(n+1)^{2N}$, then overflow in IEEE $64$-bit floating-point arithmetic only begins after $n\geq 2^{\frac{512}{N}}$, which is rather large for spectral methods on modern computers. This method is referred to as the standard method.

At the very least, self-adjoint boundary conditions $\mathcal{B}$ must be linearly independent. But linear independence in the infinite-dimensional sense does not imply that each and every linear system in Eq. (26) is soluble. There are a few pathological combinations of boundary conditions and bases that require careful attention. To not overburden the reader, we describe the implementation details for boundary conditions on orthogonal polynomials with the standard normalization.

3.4.1 The basis satisfies some of the boundary conditions

One polynomial basis for $L^{2}([-1,1])$ consists of $\{1,x\}\oplus\{P_{n}(x)-P_{n-2}(x)\}_{n=2}^{\infty}$, and it is not unreasonable to use this basis to satisfy Dirichlet and Neumann boundary conditions for a fourth-order eigenvalue problem, $u(\pm 1)=u^{\prime}(\pm 1)=0$. In fact, every element of the basis for $n\geq 2$ already satisfies the Dirichlet boundary conditions. The boundary conditions are:

[TABLE]

Most of this $\mathcal{B}$’s finite sections in use in Eq. (26) are row-rank deficient.

3.4.2 A column of zeros appears

For a second-order problem, consider the self-adjoint and separated boundary conditions:

[TABLE]

If we pose these boundary conditions on the Legendre polynomials, then a column of zeros appears precisely when $\nu=n$. This is because the degree-$n$ Legendre polynomial already satisfies the boundary conditions. In fact, analytical formulas for the Legendre recombination of the Robin boundary conditions, $a_{\pm}u(\pm 1)+b_{\pm}u^{\prime}(\pm 1)=0$, are presented by Doha and Bhrawy [29, §2.1]:

[TABLE]

where:

[TABLE]

We stress that this analytical method fails in the neighbourhood of $n=\nu$, for the reasons described above. Surely, $\mathcal{B}P_{\nu}\equiv 0,$ but this creates a column of zeros in the linear functionals, rendering the adjacent linear systems column-rank deficient. Rank-deficient boundary conditions were also overlooked for ultraspherical and Jacobi polynomial recombinations [25, 26, §4.1].

Both of these pathological scenarios are soluble by doubling the lower bandwidth of $A$. This leads to our so-called pathological method, which returns a minimum $2$-norm solution for the $k^{\rm th}$ column of $A$ using the matrix $\mathcal{B}_{0:2N-1,k+1:k+4N}$ and the same right-hand side as in the standard method. As the systems may be rank deficient, we solve them using a Moore–Penrose pseudoinverse [1] rather than a complete orthogonal factorization. We note that both of these scenarios are pathological but almost never appear. This can be argued because there are a countably infinite subset of boundary conditions, with Lebesgue measure zero, such as those described in Eq. (33) of an uncountably infinite set in general, since $a_{\pm}$ and $b_{\pm}$ range over the reals.

Although we do not have a formal proof that doubling the lower bandwidth ensures the existence of $A$ for $2N^{\rm th}$ order problems, our rationale is that the $2N$ boundary conditions are linear combinations of the $4N$ endpoint conditions $u^{(\nu)}(\pm 1)$, where $\nu=0,\ldots,2N-1$, that are linearly independent in the infinite-dimensional sense. And although most boundary conditions require only $2N$ lower bands, some may require $4N$. It is not inconceivable, however, that there may be other pathological boundary conditions not soluble by this method; after all, $\mathcal{B}$ is infinite-dimensional.

The difficulty in the implementation of $C^{-1}L_{B}$ is the appearance of the inverse of an upper-triangular banded matrix. Although inverses of banded matrices are efficiently represented by semi-separable matrices [34, Chap. 14], only the lower-triangular part of $C^{-1}L_{B}$ is required, as the upper-triangular part is populated by reflection. One computationally inexpensive scheme to obtain $\operatorname{tril}(C^{-1}L_{B})$ creates a partial inverse of $C$ that is also upper-triangular and banded, whose upper bandwidth is necessarily and sufficiently the sum of the upper bandwidth of $C$, the lower bandwidth of $L_{B}$, and the lower bandwidth of $A$. Another is to solve the banded linear systems for the lower triangular entries of each and every column of $C^{-1}L_{B}$. This approach may be more stable for larger bandwidth problems.

3.5 Skew-symmetry

It is common for a second-order self-adjoint problem to be factored in terms of the product of skew-symmetric differentiation matrices on either side of a positive-definite multiplication operator representing the non-negative variable coefficient. Skew-symmetry of the differentiation matrices then results in semi-discretizations of equations of evolution that are “stable by design” [35]. For problems on $L^{2}([-1,1])$ such as [36]:

[TABLE]

the boundary condition is represented by the linear functional on normalized Legendre polynomials:

[TABLE]

and our pathological method suggests we expand in the recombined polynomial basis:

[TABLE]

It is left as an exercise for the reader to confirm the banded skew-symmetry of the resulting differentiation matrix and the banded symmetric positive-definiteness of the accompanying conversion between trial and test functions.

4 Applications

The ultraspherical spectral method has been implemented in the software package ApproxFun.jl [37]. We have contributed to this package by adding the mechanism which symmetrizes the formally self-adjoint linear differential operators. Below, we show a variety of different applications that may help the reader in determining the applicability of the method. As the symmetrizing mechanism has been implemented in software, our discussion regarding the literal discretizations is rather kurt with the exception of one situation where it is described completely.

Numerically, it is impossible to determine whether or not a linear differential operator and its boundary conditions are formally self-adjoint. In the context of our symmetrizer and due to Eq. (22), the entries of the symmetric-definite discretization are populated in the subdiagonal part of the operator. If the boundary conditions were not self-adjoint, then the entries in superdiagonals above the banded part of the operator would not be relatively near the machine precision. If the operator were not self-adjoint, then the entries above the main diagonal would not be conjugates of the associated entries below the main diagonal. Since linear differential operators may be unbounded, even determining a heuristic tolerance for either one of the conditions of bandedness or symmetry would be challenging. Thus, in contrast with the solution of linear differential equations, there is an additional burden on the user to ensure that her/his eigenproblem is formally self-adjoint.

4.1 Asymptotics of the spectra of confined quantum anharmonic oscillators

The Schrödinger equation [38] is a celebrated Sturm–Liouville eigenvalue problem:

[TABLE]

In [39], Gaudreau et al. derive a reversed WKB series for the eigenvalues of the Schrödinger equation on the real line with the quartic anharmonic potential $V(x)=\omega x^{2}+x^{4}$ of the form:

[TABLE]

where:

[TABLE]

and more terms are available.

When the domain is truncated to $[-10,10]$ and Eq. (40) is accompanied by Dirichlet boundary conditions, Weyl’s law [40] states that we should anticipate the eigenvalues with absolute value greater than the potential to approach those of a square well potential:

[TABLE]

where $\mu([-10,10])=20$ is the Lebesgue measure of the interval.

Additionally, if $\omega\gg 1$, we anticipate the smallest eigenvalues to have the asymptotic behaviour of the harmonic oscillator:

[TABLE]

With $\omega=25$, all three of these asymptotic régimes are on display in the left panel of Figure 1 when viewing the first $10,000$ eigenvalues of Eq. (40). To calculate them, we create a $16,667\times 16,667$ banded and symmetric-definite discretization, we compute the truncated eigenvalues, and we extract the lowest $60\%$, which we anticipate to be accurate based on the discounting factor of $\frac{\pi}{2}$ in the degree of Chebyshev/Legendre expansions required to resolve plane wavefunctions of the same frequency [9].

4.2 Generalized eigenvalue problems

Our Petrov–Galerkin approach also symmetrizes banded discretizations of generalized eigenvalue problems. The following sixth order generalized eigenvalue problem models the vibration of nonuniform rings [41]. In self-adjoint form222Formally self-adjoint problems need not be instantiated in self-adjoint form for the symmetry to transcend the discretization, but a more compact representation may be more efficient., the problem reads:

[TABLE]

where $f(x)=-4(r-1)x^{2}+4(r-1)x+1$, for a thickness parameter $0<r<\infty$, $\phi=f^{3}$, and the boundary conditions are $u=u^{\prime}=u^{\prime\prime\prime}=0$ at both $x=0$ and $x=1$.

With sixth order derivatives in the first operator in Eq. (47), this is by all accounts an ill-conditioned generalized eigenvalue problem. And as $r\searrow 0$, this ill-conditioning is compounded by the fact that both linear operators become singular as the smallest generalized eigenvalue tends to zero. This is explained by the fact that the variable coefficients:

[TABLE]

and $\phi$ have multiple roots at the midpoint of the computational domain when $r=0$, and some of these roots would split apart were $r$ permitted to cross the physical threshold and become negative.

For constant coefficients ($r=1$), Driscoll [3] accurately computes the first four eigenvalues and eigenfunctions in $64$-bit floating-point arithmetic, reporting the square roots of the eigenvalues as the vibration frequencies. The results are included here in Table 1.

Using our symmetric-definite banded discretization, we compute the generalized eigenvalues in 256-bit extended precision. The banded sparsity is crucial in extended precision and for ill-conditioned problems such as this one because every flop is so much more costly. Using the general Rayleigh quotient iteration [10, 11] for the pencil $(A,B)$:

[TABLE]

our extended precision results in Table 2 agree with every digit in Driscoll’s (correctly rounded) computations.

The right panel of Figure 1 plots the first $10$ generalized eigenvalues as functions of the thickness parameter $r$ computed by creating $1,600\times 1,600$ discretizations in extended precision and using the general Rayleigh quotient to iteratively refine approximations based on the spectral data at a slightly larger thickness parameter; that is, a homotopy marches across the figure from right to left. A banded $LDL^{\top}$ factorization is used to solve the linear systems in Eq. (49) in linear complexity in the discretization dimension. The discrete generalized eigenvalues are compared with those corresponding to $3,200\times 3,200$ discretizations and the pointwise relative difference on the left side of the figure is less than $3\times 10^{-8}$, sufficient for plotting purposes. We remark that some of the first $10$ generalized eigenvalues of the $800\times 800$ discretization agree with the $1,600\times 1,600$ only to two digits, a testament to the severe ill-conditioning of this problem. Based on our numerical experiments, we conjecture that the smallest generalized eigenvalue scales with the power law:

[TABLE]

4.3 Piecewise-defined bases

Eigenfunctions need not live in an unweighted Hilbert space. On the ray $D=[-1,\infty)$, we consider the Sturm–Liouville eigenvalue problem:

[TABLE]

with the piecewise-defined function:

[TABLE]

Formally, $\frac{1}{1+\left|{x}\right|}\left[-\mathcal{D}^{2}+V(x)\right]$ is a self-adjoint operator on $L^{2}([-1,\infty),(1+\left|{x}\right|){\rm\,d}x)$, and we could identify an orthonormal basis that spans this Hilbert space. However, division by $1+\left|{x}\right|$ is only approximately representable as a banded linear operator on weighted polynomial bases, resulting in an undue complication in the discretization. Instead, we view both sides of Eq. (51) as self-adjoint operators on $L^{2}([-1,\infty),{\rm\,d}x)$ and we utilize the following piecewise-defined bases as the starting point for a banded symmetric-definite discretization:

[TABLE]

where $\tilde{L}_{n}(x)\equiv L_{n}(x)$ are the orthonormalized Laguerre polynomials that happen to coincide with the standard normalization [24, §18.3], and {\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}_{I}(x) is the characteristic function on the set $I$.

Notice that by incorporating half the full Laguerre weight, $e^{-x}$, the orthonormal basis $\{\phi_{n}^{R}\}$ already satisfies the boundary condition at infinity. However, by interlacing the left and right bases:

[TABLE]

we must also enforce two additional continuity conditions at the junction, namely $u(0^{-})-u(0^{+})=u^{\prime}(0^{-})-u^{\prime}(0^{+})=0$, for a grand total of three infinite-dimensional linear algebraic conditions in $\mathcal{B}$.

In this problem, an exponentially small error would be retained if the domain were truncated to a reasonable size, such as $[-1,50]$, and two piecewise-defined scaled and shifted Legendre polynomial bases were elected to span $L^{2}([-1,50],{\rm\,d}x)$. However, recombinations of weighted Laguerre polynomial expansions allow one to recover high relative accuracy far along the tails of the eigenfunctions that may be utilized in sensitive subsequent calculations, such as the computation of conditional expectations.

Figure 2 shows the banded sparsity patterns in the symmetric-definite discretizations of all three problems in this section.

4.4 Spectral analysis of the model problem

In this section, we consider a spectral analysis of the negative Laplacian on $[-1,1]$ with Dirichlet boundary conditions, as discussed in § 2. Our spectral discretization results in the generalized symmetric-definite pencil $(D,M)$, where:

[TABLE]

and all other entries are zero. By the type of pencil, all generalized eigenvalues are real. Since $M$ is irreducible, all generalized eigenvalues are distinct. The remaining challenge is to obtain individual bounds.

A first idea is to employ Gerschgorin’s theorem, either on the pencil [12] or after congruence transformation to the standard eigenproblem $D^{-\frac{1}{2}}MD^{-\frac{1}{2}}$. Since $D$ is diagonal and positive, its square-root is real and agrees with the component-wise definition. However, most of the resulting Gerschgorin disks are not disjoint and so the approach fails to offer reasonable bounds. Furthermore, it is readily confirmed that diagonal similarity transformations cannot isolate each and every disk [42, 43].

We thus continue with a different approach based on Sylvester’s law of inertia, which has been explored and extended to pencils by Nakatsukasa and Noferini [44]. Given a shift $\sigma\in\mathbb{R}$, one may determine how many generalized eigenvalues of $(D,M)$ are larger than, equal to, or smaller than $\sigma$ by examining the inertia of the matrix $D-\sigma M$; they are equal. From Eq. (55), the left boundary of the $n^{\rm th}$ Gerschgorin disk of $D-\sigma M$ is:

[TABLE]

Since the term in parentheses is bounded above by $1$ provided that $n>0$, the left boundary satisfies:

[TABLE]

It follows that if the mild shift $\sigma_{n}=(n+1)(n+2)$ is chosen, then the trailing infinite section $(D-\sigma M)_{n+1:\infty,n+1:\infty}$ is positive-definite, and hence $\lambda_{n}(D,M)>(n+1)(n+2)$. For an upper bound, by the Courant–Fisher minimax theorem [1, Theorem 8.1.2], it suffices to consider the $(n+1)\times(n+1)$ principal finite sections of $D$ and $M$. Then:

[TABLE]

For estimates of the induced $2$-norms, we find the lower bound on the smallest Gerschgorin disk for $M_{0:n,0:n}$:

[TABLE]

and since $D$ is diagonal, $\left\|{D_{0:n,0:n}}\right\|_{2}=(n+1)(n+2)$ so that:

[TABLE]

Finally, we note that the inequality $2<\lambda_{0}(D,M)<8$ also holds since the $0^{\rm th}$ Gerschgorin disk for $(D,M)$ is isolated.

The analysis of the eigenvalues of the negative Laplacian are important for understanding the stability of the choice of spatial discretizations when marching forward in time, e.g. when solving the heat equation. It is also important as a model problem because all the eigenvalues are known analytically. For $n\times n$ pseudospectral discretizations in Chebyshev points, Gottlieb and Lustman [4] proved that the spectra are real, positive, and distinct, Charalambides and Waleffe [45] proved that the spectra of sequential discretizations interlace, and Weideman and Trefethen [9] proved an asymptotically constant lower-bound and quartic upper-bound for all of them. These results are impressive for a dense asymmetric matrix. But in the limit of a large number of collocation points, do they guarantee that the pseudospectral discretization has discrete spectrum that accumulates at infinity? Has continuous spectrum been ruled out? The preceding analysis of our symmetric-definite and banded discretization reveals independent bounds for each eigenvalue such that our discretization is faithful to the continuous operator all the way to infinity.

More generally, we argue that banded sparsity promotes a deeper, more penetrating, analysis that is equally applicable to the other problems in this section; please see Figure 2. While the analysis may be tedious, the sparsity in our discretization renders it in principle much simpler than many alternative (dense) discretizations.

Remark 4.1**.**

The diagonal entries, $d_{n}$, of $D$ in an $LDL^{\top}$ factorization of $D-\sigma M$ are given recursively by:

[TABLE]

Properties of the inertia of $D-\sigma M$ are obscured by the nonlinearity in this recurrence relation, but we should not preclude tighter bounds from being determined by this avenue.

5 Symmetric-(definite) banded eigensolvers

So far, our applications have only considered eigenvalue computations. From here to the end, we consider algorithms for the spectral decomposition of self-adjoint problems.

In this section and the next, we break with the notation already in use. In particular, regular capital letters represent finite-dimensional matrices, calligraphic letters represent infinite-dimensional matrices (rather than continuous operators), and indexing begins at $1,1$ rather than $0,0$. We also use shorthand notation of $a:b$ to represent all entries of an array from index $a$ to $b$ inclusive, and this is sometimes tensorized. Other notational shorthands attempt to keep the presentation light.

Proposition 5.1**.**

Let $A\in\mathbb{R}^{k\times k}$ be a symmetric banded matrix with bandwidth $b>0$ such that $A_{i,j}=0$ for $\left|{i-j}\right|>b$. The algorithms of [46, 47] compute the spectral decomposition:

[TABLE]

in $\mathcal{O}(k^{2}b)$ flops. The eigenvectors $Q$ are products of Givens rotations and the eigenvectors of the similar symmetric tridiagonal problem, and the dense matrices may be formed in $\mathcal{O}(k^{3})$ flops if required.

The algorithms reduce the symmetric banded problem to symmetric tridiagonal form and introduce orthogonal similarity transformations to preserve the bandwidth of the problem, chasing bulges beyond the bandwidth off the matrix. These algorithms are backward stable, and Kaufman improves the efficiency by demonstrating that multiple orthogonal similarity transformations may be vectorized. The symmetric tridiagonal spectral decomposition may be produced using one of many alternatives in $\mathcal{O}(k^{2})$ flops [1].

Proposition 5.2**.**

Let $A\in\mathbb{R}^{k\times k}$ and $B\in\mathbb{R}^{k\times k}$ be symmetric banded matrices both bandlimited by $b>0$, and let $B$ be positive-definite. The algorithms of [48, 49, 50] compute the generalized spectral decomposition:

[TABLE]

in $\mathcal{O}(k^{2}b)$ flops. In this case, the eigenvectors are $B$-orthogonal, satisfying $V^{\top}BV=I$.

In LAPACK, this algorithm uses the split-Cholesky factorization of $B=S^{\top}S$ to convert the symmetric-definite banded generalized eigenvalue problem into the regular symmetric banded eigenvalue problem:

[TABLE]

The matrix $X:=S^{-1}Q$ consists of the inverse of one of the split-Cholesky factors combined with orthogonal transformations that are designed to ensure that $C$ is banded. This improvement due to Wilkinson [49] is based on Crawford’s algorithm [48] which originally uses a Cholesky factorization of $B=LL^{\top}$.

Afterward, the symmetric banded spectral decomposition of $C$ follows from Proposition 5.1:

[TABLE]

In terms of $S$, $Q$, and $W$, the generalized eigenvectors are given by the product $V=S^{-1}QW$.

The stability of the algorithms is related to the norms of the split-Cholesky factors. Thus, ill-conditioning in $B$ may obscure the result.

6 An adaptive spectral decomposition

Several authors have recently investigated algorithms for infinite-dimensional numerical linear algebra, including an adaptive $QR$ decomposition by Olver and Townsend [6], algorithms for the spectra of finite-rank perturbations of Toeplitz symmetric tridiagonal Jacobi operators by Webb [51], and related algorithms investigated by Hansen [32] and Colbrook and Hansen [52]. In infinite dimensions, similarity and congruence transformations may begin “at infinity” and terminate at the base of the operator, or they may begin at the base, iteratively working on finite-sections of the operator. In our adaptive spectral decomposition below, we espouse the latter perspective.

6.1 An adaptive spectral decomposition

In this section, we describe an algorithm for the partial spectral decomposition of:

[TABLE]

with bandwidth $b$.

Algorithm 6.1**.**

Initialize $k:=2b$. Partition $\mathcal{L}$ into:

[TABLE]

In $\mathcal{O}(k^{2}b)$ operations, compute the symmetric banded spectral decomposition:

[TABLE]

and use the orthogonal matrix $Q$ to rewrite the partition as a similarity transformation on the operator:

[TABLE]

Define $M:=L_{b}Q$ and set $j$ to be the largest integer multiple of $b$ that satisfies the inequality:

[TABLE]

for some tolerance ${\tt tol}>0$. For eigenvalues only, ${\tt tol}=\epsilon$ and for both eigenvalues and eigenfunctions, ${\tt tol=\epsilon^{2}}$. If there is no integer $j$ that satisfies the inequality, then double $k$ and reiterate until deflation sets in.

After determining that $j$ eigenvalues are to be deflated, we still have the operator partitioned as above. Since $M$ is dense, the structure of the similar operator is now:

[TABLE]

We perturb the operator by setting $M_{1:b,1:j}\leftarrow 0$. Thus:

[TABLE]

It is now clear that the first $j$ eigenvalues of the perturbed operator may be deflated. Furthermore, the remaining part of the similar operator is of the same form:

[TABLE]

We now perform orthogonal similarity transformations to return the operator to symmetric banded form. Below, we describe the effect of applying Householder reflectors [1, §5.1] from the left. The analogous Householder reflectors are also applied from the right, and the figures show the effect of both as similarity transformations. To simplify the exposition, let $\ell:=k-j$.

Begin by applying Householder reflectors to lower triangularize the entries in rows $1:b$ and columns $\ell+1:\ell+b$, resulting in the similar operator:

[TABLE]

Black symbols represent the cumulative modification of nonzero submatrices. 2. 2.

Next, apply Householder reflectors to lower triangularize the entries in rows $1:2b$ and columns $\ell+1:\ell+b$:

[TABLE] 3. 3.

Before continuing to annihilate entries in the $\ell+1:\ell+b$ columns, apply Householder reflectors to lower triangularize the entries in rows $1:b$ and columns $b+1:2b$:

[TABLE] 4. 4.

Next, apply Householder reflectors to lower triangularize the entries in rows $b+1:3b$ and columns $\ell+1:\ell+b$:

[TABLE]

Note that by lower triangularizing in the previous step, we extend entries into the rows $1:b$ and columns $2b+1:3b$ as a lower triangular submatrix. 5. 5.

Apply Householder reflectors to lower triangularize the entries in rows $1:2b$ and columns $2b+1:3b$:

[TABLE]

And apply Householder reflectors to lower triangularize the entries in rows $1:b$ and columns $b+1:2b$:

[TABLE] 6. 6.

Continuing, we apply Householder reflectors to lower triangularize the entries in rows $2b+1:4b$ and columns $\ell+1:\ell+b$:

[TABLE] 7. 7.

At this point, we are ready to paraphrase the remaining orthogonal similarity transformations. The bulge that has been created in the banded submatrix is chased back up by Householder reflectors toward the top left corner to return to banded form:

[TABLE] 8. 8.

Finally, further sequences of Householder reflectors continue to annihilate entries in the $\ell+1:\ell+b$ columns, chasing bulges in the banded submatrix in between, until a we arrive at a similar banded operator:

[TABLE]

6.2 Complexity analysis

Based on the results in [46, 47], this algorithm requires:

$\mathcal{O}(k^{2}b)$ flops for the symmetric banded spectral decomposition $L_{k}=Q\Lambda Q^{\top}$; 2. 2.

$\mathcal{O}(kb^{2})$ flops for the banded matrix-matrix product $M=L_{b}Q$ ($L_{b}$ has at most $b(b+1)/2$ nontrivial entries); 3. 3.

$\mathcal{O}(jb)$ flops for computing the deflation criterion; and, 4. 4.

$\mathcal{O}(\ell b^{2})$ flops for folding in the symmetric spikes and $\mathcal{O}(\ell^{2}b)$ for chasing the symmetric bulges up and off the band.

To ensure quadratic complexity, we must ensure the number of deflated eigenvalues in each iteration is not too small relative to the dimensions of the deflation window. To make this precise, let $j_{1}$ be the number of deflated eigenvalues after one step and let $j_{2}$ be the number after another. The total cost is:

[TABLE]

where $C_{i}$, for $i=1:5$, are positive constants that depend on the implementation.

Compare this with the same algorithm on the larger window of dimensions $(k+j_{1})\times(k+j_{1})$. This window theoretically deflates the same number of eigenvalues as two consecutive steps of fixed dimensions $k\times k$:

[TABLE]

If we assume that $C_{1}\gg C_{2}+C_{3}+C_{4}+C_{5}$, then this implies that the spectral decomposition is more expensive than computing the banded matrix-matrix multiplication, deflation criterion, bulge chasing, and spike folding. With this assumption, for the two step cost to break even with the larger one step cost, the number of eigenvalues that must be deflated is:

[TABLE]

If this is not the case, then we double the dimensions of the deflation window.

6.3 Determining the deflation criterion

In aggressive early deflation, it is important to identify which eigenvalues have converged in the working precision, subject to rounding errors. Based on Eq. (69), partition the similar operator as:

[TABLE]

where $M_{j}=M_{1:\infty,1:j}$, $M_{k}=M_{1:\infty,j+1:k}$, $\Lambda_{j}=\Lambda_{1:j,1:j}$, $\Lambda_{k}=\Lambda_{j+1:k,j+1:k}$, and $\mathcal{L}_{\infty}=\mathcal{L}_{k+1:\infty}$ for brevity. To determine the deflation criterion, we appeal to the block determinant formulæ:

[TABLE]

provided the inverses exist. If we use the formulæ twice reflecting the partition in Eq. (86), then:

[TABLE]

The key to understanding Eq. (88) is the second determinant in the product:

[TABLE]

This determinant informs us how close the first $j$ eigenvalues are to deflation, for if $M_{j}\equiv 0$, then complete deflation is admissible. Numerically, we deflate when $M_{j}^{\top}(\mathcal{L}_{\infty}-\lambda I)^{-1}M_{j}$ is relatively small compared with $\Lambda_{j}-\lambda I$. Since the operator is self-adjoint, we are interested in the relative norms over the spectrum of $\Lambda_{j}$:

[TABLE]

for any induced $p$-norm. On the other hand, the analysis of:

[TABLE]

is simplified in the induced $2$-norm, for the norm of the inverse of the infinite-dimensional operator $\mathcal{L}_{\infty}-\lambda I$ may be rephrased in terms of the smallest eigenvalue of $(\mathcal{L}_{\infty}-\lambda I)^{-1}$ for nearest $\lambda$, that is, the separation distance between the spectra of $\Lambda_{j}$ and $\mathcal{L}_{\infty}$.

Using the above analysis, we derive the heuristic deflation criterion in Eq. (70) by assuming that the smallest eigenvalue of $\mathcal{L}_{\infty}$ is larger than the largest eigenvalue of $\Lambda$, in absolute value, and by using the equivalence of (finite-dimensional) Frobenius and induced $2$-norms. The Frobenius norm simplifies the computation, since no singular values need be computed, and since the square of the $2$-norm of the $j^{\rm th}$ column of $M$ may be added to $\left\|{M_{1:b,1:j-1}}\right\|_{F}^{2}$ recursively:

[TABLE]

6.4 An adaptive generalized spectral decomposition

As in the applications of § 4, our Petrov–Galerkin scheme for a regular eigenvalue problem may require reformulation as a generalized eigenvalue problem:

[TABLE]

We will assume that the $b$ is the maximum of the bandwidths of $\mathcal{A}$ and $\mathcal{B}$ and that $\mathcal{B}$ is also positive-definite.

Algorithm 6.2**.**

Initialize $k:=2b$. Partition $\mathcal{A}$ into:

[TABLE]

and partition $\mathcal{B}$ conformably:

[TABLE]

In $\mathcal{O}(k^{2}b)$ operations, compute the symmetric-definite banded generalized spectral decomposition:

[TABLE]

and use the eigenvectors $V$ to rewrite the partition as a congruence transformation on the operators:

[TABLE]

Define $M:=A_{b}V$ and $N:=B_{b}V$ and set $j$ to be the largest integer multiple of $b$ that satisfies the inequalities:

[TABLE]

for some tolerance ${\tt tol}>0$. For eigenvalues only, ${\tt tol}=\epsilon$ and for both eigenvalues and eigenfunctions, ${\tt tol=\epsilon^{2}}$. If there is no integer $j$ that satisfies the inequalities, then double $k$ and reiterate until deflation sets in.

After determining that $j$ eigenvalues are to be deflated, we have a symmetric-definite operator pencil, and the second operator has unit diagonal elements for the first $k$ entries. We perturb the operators by setting $M_{1:b,1:j}\leftarrow 0$ and $N_{1:b,1:j}\leftarrow 0$ in order to deflate the first $j$ generalized eigenvalues. Thus, the remaining parts of the congruent operators are of the same form:

[TABLE]

We now perform congruence transformations to return the operator pencil to symmetric-definite banded form.

The following symmetric positive-definite matrix has a sparse Cholesky factorization:

[TABLE]

Embedding the Cholesky factors in the identity operator:

[TABLE]

Similarly:

[TABLE]

This congruence transformation results in the operator pencil:

[TABLE]

As in Algorithm 6.1, black symbols represent the cumulative modification of nonzero submatrices. The first operator is still symmetric and the second operator also retains symmetry and positive-definiteness.

Finally, we are free to use the same orthogonal similarity transformations described in Algorithm 6.1; these similarity transformations have no effect on the second operator since the first $\ell$ diagonal entries are unit, guaranteeing our return to symmetric-definite banded form:

[TABLE]

6.5 Summary of the data structures

The adaptive spectral decomposition returns eigenfunctions as a product of orthogonal operators:

[TABLE]

The orthogonal operators $\mathcal{Q}_{i}$ consist of the orthogonal matrices that implement the similarity transformations embedded in identity operators.

If the first $n$ eigenfunctions are desired, then the product of the orthogonal operators up to and beyond convergence of the eigenvalues is formed, at a cost of $\mathcal{O}(n^{3})$. If only a subset of the eigenfunctions is desired, then it is possible to reconstruct them in lower complexity. If the operator is applied to a function, then this may be performed in $\mathcal{O}(n^{2})$ operations.

In the same spirit, the adaptive generalized spectral decomposition returns generalized eigenfunctions as a product of operators:

[TABLE]

In certain cases, the complexity of our adaptive spectral decomposition may be linear. Thus, if we only need to apply the operator (or some function of it) to a function, then the product form returned by the adaptive spectral decomposition is ideal. We have observed linear complexity for a bounded perturbation of an (unbounded) spectral operator defining the expansion basis. In this setting, the $n^{\rm th}$ eigenfunction requires at most $n+\mathcal{O}(\log\epsilon^{-1})$ degrees of freedom to resolve. This is reasonable from an approximation-theoretic point of view since the asymptotic distribution of roots of the sought eigenfunctions is the same as the spectral basis used to represent them. Precise bounds for the special cases of the Mathieu and spheroidal wave equations have been described by Volkmer [53].

7 The fractional Schrödinger equation

There are many applications that benefit from the spectral decomposition of a linear differential operator. Consider the time-dependent fractional Schrödinger equation [54, 55]:

[TABLE]

and where $\mathcal{H}_{\alpha}$ is a self-adjoint time-independent Hamiltonian and the Lévy index $0<\alpha\leq 2$ denotes the strength of the implied fractional derivative. The formal solution of this equation is:

[TABLE]

If the partial spectral decomposition of $\mathcal{H}_{\alpha}$ is available, then this problem no longer requires time-stepping. We consider $\mathcal{H}_{\alpha}$ to be a fractional Mathieu operator [56, 57] with periodic boundary conditions and a periodic potential on $[0,2\pi)$:

[TABLE]

The real number $q$ determines the relative strength of the potential compared with the differential operator. We use a real-valued orthonormal Fourier series to represent the solution:

[TABLE]

where the ordering is chosen so that the regularity of the solution is observed in the decay of the coefficients. When $q=0$, the eigenfunctions are pure Fourier modes and otherwise, the fractional Mathieu operator is a bounded perturbation of the (unbounded) fractional Laplacian [58, 59]. Using the differentiation and linearization properties of trigonometric functions, the fractional Mathieu operator is:

[TABLE]

Figure 3 shows the evolution of the fractional Schrödinger equation for various values of $\alpha$ and $q$. For an initial condition whose frequency is on the order of $q$, the particle densities of the classical Schrödinger equation barely notice the potential. On the other hand, for large deviations from a uniform potential, the fractional Schrödinger equation appears to focus the particle densities around the wells when $\alpha\ll 2$.

It is straightforward to confirm that the fractional Schrödinger operator is unitary; therefore, the $L^{2}$-norm of the solution is a constant of motion. While unitarity preserving time-stepping schemes are of great interest for such problems, we note that unitarity preservation is also a property of a partial spectral decomposition. In Figure 4, we monitor the numerical error in the $L^{2}$-norm of the Fourier coefficients starting with two $2$-normalized initial conditions: first, we start with $n=10,001$ coefficients that mimic an analytic function $u_{\pm i}=\exp(-36(i+1)/\lceil\frac{n}{2}\rceil)$; second, with pseudorandom standard normal coefficients $u_{\pm i}=\mathcal{N}(0,1)$, for $i=0,\ldots,\lfloor\frac{n}{2}\rfloor$. In both cases, the norm of the solution is very well preserved on long time scales. We also show the calculation time for the partial spectral decomposition and execution time of the eigenmatrix-vector product as a function of the truncation degree of the Fourier series, confirming the $\mathcal{O}(n)$ complexity predicted in § 6.5.

8 Conclusion & outlook

There are many avenues for future work. It is straightforward to consider symmetrizing eigenvalue problems from singular integral equations [30]. Another natural extension is to eigenvalue problems in multiple dimensions. This includes tensor-product spectral methods [60] and total-degree spectral methods on spheres, disks, and triangles [61]. In the former, operators must have a type of separable structure to permit fast linear algebra such as an eigen-analogue of the generalized Sylvester matrix equation, whereas in the latter, the logic of Algorithms 6.1 and 6.2 is the same but the blocks are now banded matrices themselves whose dimensions increase with the block indices.

We presented an algorithm for an adaptive spectral decomposition where approximations are made in the deflation criterion and in the use of floating-point arithmetic. The FMM has been adapted to provide asymptotically fast approximate application of the eigenvectors of symmetric arrowhead matrices [62, 63], where the approximate precision is a tuneable parameter. The FMM may in principle be utilized to accelerate the spectral decomposition of the block arrow matrices that arise in Eq. (71), though locating the spectrum may be more challenging in the blocked generalization.

Other variants of infinite-dimensional linear algebra could be pursued. These include an adaptive singular value decomposition and an adaptive Schur decomposition for non-self-adjoint problems. The symmetrizer may also be useful in the design and analysis of contours for spectral projectors in, e.g. the contFEAST algorithm [64].

Acknowledgments

We thank Yuji Nakatsukasa for a private communication on the spectral analysis of the model problem in § 4.4. We thank Andrew Horning and Marcus Webb for a careful reading of the manuscript. JLA is supported by “la Caixa” Banking Foundation, through the Postdoctoral Junior Leader Fellowship Programme. RMS is supported by the Natural Sciences and Engineering Research Council of Canada, through a Discovery Grant (RGPIN-2017-05514).