Recently Olver and Townsend presented a fast spectral method that relies on bases of ultraspherical polynomials to give differentiation matrices that are almost banded. The almost-banded structure allowed them to develop efficient algorithms for solving ceratin discretized systems in linear time. We show that one can also design fast algorithms for standard spectral methods because the underlying matrices, though dense, have the same rank structure as those of Olver and Townsend

Aurentz, Jared L

Mathematical Institute Eprints Archive

FAST ALGORITHMS FOR SPECTRAL DIFFERENTIATION MATRICES ∗
JARED L. AURENTZ†
Abstract. Recently Olver and Townsend presented a fast spectral method that relies on bases of ultraspherical
polynomials to give differentiation matrices that are almost banded. The almost-banded structure allowed them to
develop efficient algorithms for solving ceratin discretized systems in linear time. We show that one can also design
fast algorithms for standard spectral methods because the underlying matrices, though dense, have the same rank
structure as those of Olver and Townsend.
Key words. spectral methods, rank structured matrices
AMS subject classifications. 65N35, 33C45, 65F05
1. Introduction. In [13] Olver and Townsend presented a fast spectral method for solv-
ing linear differential equations using bases of ultraspherical polynomials, which has sub-
sequently been exploited in Chebfun [6] and ApproxFun [12]. The “fast” in these methods
appears to come from the fact that the discretized systems involve almost-banded matrices,
making it possible to solve the systems in linear time via structured QR factorizations.12
More standard bases, such as Chebyshev, give dense matrices that appear to require more
expensive linear algebra. Here we show that, although standard spectral matrices are dense,
their underlying rank structure is equivalent mathematically to being almost-banded and the
systems of equations can still be solved in linear time.
2. Simple example. We will present the main idea by means of a particular example,
namely, given a continuous function f defined on [−1, 1], find u such that
u′ + u = f, u(1) = 0.
As is well established in [13], this model serves as a prototype for a larger class of problems.
Our first task is to rewrite this equation in terms of coefficients of orthogonal polynomials.
We will do this in two ways. The first method will construct the matrices using the tech-
niques in [13]. The second will construct equivalent matrices that depend only on Chebyshev
coefficients and make no conversions to other ultraspherical representations. When clarity is
needed we will subscript variables corresponding to the first method by ultra and variables
corresponding to the second method with cheb.
Our linear operator can be written as the sum of two basic operators: differentiation (D)
and identity (I). In the first method, the ultraspherical method of Olver and Townsend, D is
the matrix that maps the Chebyshev coefficients of u to the ultraspherical coefficients of u′.
∗This research was partially supported by the European Research Council under the European Union’s Seventh
Framework Programme (FP7/2007–2013)/ERC grant agreement no. 291068. The views expressed in this article are
not those of the ERC or the European Commission, and the European Union is not liable for any use that may be
made of the information contained here.
†Mathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Road, OX2 6GG Oxford,
UK; (aurentz@maths.ox.ac.uk).
1 This is true when the coefficients of the differential equation are low degree polynomials.
2 Almost-banded matrices also appear in the Chebyshev-based spectral methods in [3] and [9]. These techniques
rely on first reformulating the differential equation as an integral equation.
1
2 J. L. Aurentz
Let {ui}
∞
i=0 be the Chebyshev coefficients of u. The matrixDultra is defined as
Dultrau =


0 1
2
3
. . .




u0
u1
u2
u3
...


.
Similarly, one can write down the matrix Iultra that takes Chebyshev coefficients to ultras-
pherical coefficients (also known as the conversion operator),
Iultrau =
1
2


2 0 −1
1 0 −1
1 0
. . .




u0
u1
u2
u3
...


.
To account for the boundary condition we need an additional functional for computing u(1).
For Chebyshev coefficients this takes the form of a simple inner product,
u(1) =
[
1 1 1 1 · · ·
]


u0
u1
u2
u3
...


.
Combining all these pieces together we get the following system of equations:
(2.1)
1
2


2 2 2 2 · · ·
2 2 −1
1 4 −1
1 6
. . .




u0
u1
u2
u3
...


=


0
fˆ0
fˆ1
fˆ2
...


.
Here the {fˆi} are the ultraspherical coefficients of f computed by multiplying f by Iultra.
We will refer to this matrix as Lultra.
To facilitate our discussion, we focus on the nonzero pattern. Below is the almost-banded
pattern of the 10× 10 Lultra matrix,
(2.2)


× × × × × × × × × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× ×


.
Fast spectral methods 3
Before we start discussing equivalence, we need the matrices for method 2, the more
standard spectral method based on Chebyshev coefficients (see e.g. [11]). The matrices D
and I now take the following forms:
Dcheb =


0 1 3 5 7 · · ·
4 8 12 16
6 10 14
8 12 16
. . .
. . .


, Icheb =


1
1
1
. . .

 .
Notice how all the rows of Dcheb with even/odd indices are multiples of each other away
from the main diagonal. This hints at a special structure that can be used to describe Dcheb
using only O(n) parameters.
The entire system in a single matrix has the following form:
(2.3)


1 1 1 1 1 1 1 1 1 · · ·
1 1 3 5 7 · · ·
1 4 8 12 16
1 6 10 14
1 8 12 16
. . .
. . .




u0
u1
u2
u3
u4
u5
u6
u7
u8
...


=


0
f0
f1
f2
f3
...


.
We will refer to this matrix as Lcheb.
Here we can see why using only Chebyshev coefficients looks expensive. The new sys-
tem is essentially dense in the upper-triangular part and appears to require O(n2) storage.
Below we illustrate the nonzero pattern of the 10× 10 Lcheb matrix,
(2.4)


× × × × × × × × × ×
× × × × × ×
× × × × ×
× × × × ×
× × × ×
× × × ×
× × ×
× × ×
× ×
× ×


.
3. Rank structured matrices. Now that we have our matrices, (2.1)-(2.2) and (2.3)-
(2.4) in hand, we can describe how these two have essentially the same rank structure. Rank
structured matrices are quite prevalent in applied mathematics and in the last few decades
considerable work has been done to develop efficient algorithms. The central idea is that even
dense matrices can be represented by very few parameters if certain submatrices have small
rank. Let’s look more closely at Lultra. Below we’ve highlighted a rectangular submatrix
4 J. L. Aurentz
whose lower-left corner intersects the main diagonal,


× × × × × × × × × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× ×


.
Notice that the highlighted matrix has rank at most 3. This remains true for any submatrix
whose lower-left corner intersects the main diagonal,


× × × × × × × × × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× ×


,


× × × × × × × × × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× × ×
× ×


.
Now we draw a similar box on the Lcheb matrix,


× × × × × × × × × ×
× × × × × ×
× × × × ×
× × × × ×
× × × ×
× × × ×
× × ×
× × ×
× ×
× ×
×


.
Fast spectral methods 5
This submatrix too has rank at most 3! In fact, any submatrix whose lower-left corner inter-
sects the main diagonal will have rank at most 3. To see this we decompose Lcheb into three
parts and highlight the same submatrix,
Lcheb =


1 1 1 1 1 1 1 1 1 1
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0


+


0 0 0 0 0 0 0 0 0 0
1 3 5 7 9
0
6 10 14 18
0
10 14 18
0
14 18
0
18


+


0 0 0 0 0 0 0 0 0 0
0
4 8 12 16
0
8 12 16
0
12 16
0
16
0


.
Since the highlighted submatrix in each of the 3 parts has rank at most 1, the highlighted
block of Lcheb has rank at most 3.
This type of rank structure is a generalization of semiseparability.3 The study of rank
structured matrices originates with Gantmacher and Kreı˘n [10]. For more recent develop-
ments see, for example, [4] or [15].
4. Fast solutions viaQR factorization. In [13] the discretized systems are solved using
QR factorizations. The authors use the almost-bandedness to show that both of the factors
Q and R in QR = Lultra have a special structure that makes it possible to compute Q and
perform the back solve in linear time. We now observe that it’s not the almost-bandedness
that matters, it’s really the rank structure. This means that it is just as easy to develop a linear
time QR based solver for the Lcheb system.
The details are a bit complicated and will be spelled out with numerical applications in
a future publication. Here we highlight the key points. The first thing to note is that both
Lultra and Lcheb are upper-Hessenberg. This implies that Qultra and Qcheb (QR = L) are
upper-Hessenberg. Consequently they and their inverses can be represented by products of n
elementary 2× 2 matrices and applied to a vector in linear time.
3Applicable definitions include both sequentially-semiseparable [5] and quasiseparable [7].
6 J. L. Aurentz
The R’s are much more interesting. In general they are dense,
Q∗ultraLultra = Rultra =


× × × × × × × × × · · ·
× × × × × × × ×
× × × × × × ×
× × × × × ×
× × × × ×
× × × ×
× × ×
. . .


,
but any submatrix whose lower left corner intersects the diagonal has rank at most 4,


× × × × × × × × × · · ·
× × × × × × × ×
× × × × × × ×
× × × × × ×
× × × × ×
× × × ×
× × ×
. . .


.
One can explain the rank structure by recalling that Lultra is rank structured and unitary
matrices preserve rank. Exactly the same reasoning can be applied to Lcheb, with the same
conclusion.
The back substitution, in both methods, involves an upper-triangular matrix with a def-
inite rank structure. Olver and Townsend exploited this structure to develop a linear time
algorithm for performing a back solve with Rultra. Their algorithm is an extension of a sim-
ilar method that was first developed by Chandrasekaran and Gu [2]. The same technique can
also be used for Rcheb.
4.1. Conditioning. It should be noted that the simple form of the Chebyshev differenti-
ation matrices presented here is subject to the same ill-conditioning as other Chebyshev-based
spectral methods [1, 8]. In [13] the authors construct a structure-preserving preconditioner to
ensure that the ultraspherical-based systems are well-conditioned. A similar technique can be
applied to the Chebyshev-based methods proposed in this note.
5. Other orthogonal families. Above, we considered a specific Chebyshev example.
We conclude this note by showing that orthogonal families more generally possess rank struc-
tured differentiation matrices. The matrices for these examples are derived from formulas
in [14].
The first example is the 10× 10 version of Lleg , the spectral discretization of u
′ + u =
Fast spectral methods 7
f, u(1) = 0 using Legendre coefficients,
Lleg =


1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1
1 3 3 3 3
1 5 5 5 5
1 7 7 7
1 9 9 9
1 11 11
1 13 13
1 15
1 17


.
Just as with Chebyshev, this submatrix has rank at most 3. A similar structure also exists for
trigonometric polynomials as illustrated by our second example. Below is the 9 × 9 version
of Lfour, the spectral discretization of u
′ = f, u(0) = 0 with
∫
f = 0, using Fourier
coefficients,
Lfour =


1 1 1 1 1 1 1 1 1
1i
−1i
2i
−2i
3i
−3i
4i
−4i


.
Here the rank is at most 2. The low rank structure can even be observed in spectral meth-
ods for unbounded domains. Our last example is the 10 × 10 version of Llag, the spectral
discretization of u′ + u = f, u(0) = 0 using Laguerre coefficients,
Llag =


1 1 1 1 1 1 1 1 1 1
1 −1 −1 −1 −1 −1 −1 −1 −1 −1
1 −1 −1 −1 −1 −1 −1 −1 −1
1 −1 −1 −1 −1 −1 −1 −1
1 −1 −1 −1 −1 −1 −1
1 −1 −1 −1 −1 −1
1 −1 −1 −1 −1
1 −1 −1 −1
1 −1 −1
1 −1


.
Here the rank is again, bounded by 2.
Conclusions. It has been shown that spectral differentiation matrices possess structure
that can be exploited to develop fast methods for a variety of orthogonal families. This ob-
servation opens the door to new classes of spectral methods for solving ordinary and partial
differential equations.
Acknowledgements. The author would like to thank the editor and referee for their
comments which helped improve the quality of this note, and for their incredibly fast response
time. The author would also like to thank Nick Trefethen for his help with both content and
aesthetics, and Raf Vandebril for his help clarifying various technical details.
8 J. L. Aurentz
REFERENCES
[1] J. P. BOYD, Chebyshev and Fourier Spectral Methods, Courier Corporation, 2001.
[2] S. CHANDRASEKARAN AND M. GU, Fast and stable algorithms for banded plus semiseparable systems of
linear equations, SIAM J. Matrix Anal. Appl., 25 (2003), pp. 373–384.
[3] C. W. CLENSHAW, The numerical solution of linear differential equations in chebyshev series, in Math. Proc.
Cambridge Philos. Soc., vol. 53, Cambridge Univ. Press, 1957, pp. 134–149.
[4] S. DELVAUX, Rank structured matrices, PhD thesis, KU Leuven, 2007.
[5] P. DEWILDE AND A.-J. VAN DER VEEN, Time-varying Systems and Computations, Springer, 1998.
[6] T. A. DRISCOLL, N. HALE, AND L. N. TREFETHEN,Chebfun Guide, Pafnuty Publications, Oxford, 2014.
[7] Y. EIDELMAN AND I. GOHBERG,On a new class of structured matrices, Integr. Equat. Oper. Th., 34 (1999),
pp. 293–324.
[8] B. FORNBERG, A Practical Guide to Pseudospectral Methods, Cambridge Univ. Press, 1998.
[9] L. FOX AND I. B. PARKER, Chebyshev polynomials in numerical analysis, Oxford Univ. Press, 1968.
[10] F. R. GANTMACHER AND M. G. KREI˘N, Sur les matrices oscillatoires et complètement non négatives,
Compos. Math., 4 (1937), pp. 445–476.
[11] A. KARAGEORGHIS, A note on the Chebyshev coefficients of the general order derivative of an infinitely
differentiable function, J. Comput. Appl. Math., 21 (1988), pp. 129–132.
[12] S. OLVER, G. GORETKIN, R. M. SLEVINSKY, AND A. TOWNSEND, ApproxFun.jl.
https://github.com/ApproxFun/ApproxFun.jl, 2015.
[13] S. OLVER AND A. TOWNSEND,A fast and well-conditioned spectral method, SIAMRev., 55 (2013), pp. 462–
489.
[14] J. SHEN, T. TANG, AND L.-L. WANG, Spectral Methods: Algorithms, Analysis and Applications, Springer,
2011.
[15] R. VANDEBRIL, M. VAN BAREL, AND N. MASTRONARDI,Matrix Computations and Semiseparable Ma-
trices: Linear Systems, vol. 1, Johns Hopkins Univ. Press, 2007.


Fast algorithms for spectral differentiation matrices

Abstract

Similar works

Full text

Available Versions

Mathematical Institute Eprints Archive