FastTransforms.jl
The aim of this package is to provide new classes of fast transforms with low pre-computation. One approach is based on the use of asymptotic formulae to relate the transforms to a small number of fast Fourier transforms. Another approach is based on a Toeplitz-dot-Hankel decomposition of the matrix of connection coefficients. Both new classes of fast transforms do not require large pre-computation for fast execution and they are designed to work on expansions of functions with any degree of regularity.
The Chebyshev—Jacobi transform and its inverse are implemented. This allows the fast conversion of Chebyshev expansion coefficients to Jacobi expansion coefficients and back.
julia> Pkg.add("FastTransforms")
julia> using FastTransforms
julia> c = rand(10001);
julia> @time norm(icjt(cjt(c,0.1,-0.2),0.1,-0.2)-c,Inf)
0.258390 seconds (431 allocations: 6.278 MB)
1.4830359162942841e-12
julia> p1 = plan_cjt(c,0.1,-0.2);
julia> p2 = plan_icjt(c,0.1,-0.2);
julia> @time norm(p2*(p1*c)-c,Inf)
0.244842 seconds (17 allocations: 469.344 KB)
1.4830359162942841e-12
The design and implementation is analogous to FFTW: there is a type ChebyshevJacobiPlan
that stores pre-planned optimized DCT-I and DST-I plans, recurrence coefficients,
and temporary arrays to allow the execution of either the cjt or the icjt allocation-free.
This type is constructed with either plan_cjt or plan_icjt. Composition of transforms
allows the Jacobi—Jacobi transform, computed via jjt. The remainder in Hahn’s asymptotic expansion
is valid for the half-open square (α,β) ∈ (-1/2,1/2]^2. Therefore, the fast transform works best
when the parameters are inside. If the parameters (α,β) are not exceptionally beyond the square,
then increment/decrement operators are used with linear complexity (and linear conditioning) in the degree.
The Padua transform and its inverse are also implemented thanks to
Michael Clarke. These are optimized methods
designed for computing the bivariate Chebyshev coefficients by interpolating a
bivariate function at the Padua points on [-1,1]^2.
References:
-
N. Hale and A. Townsend. A fast, simple, and stable Chebyshev—Legendre transform using and asymptotic formula, SIAM J. Sci. Comput., 36:A148—A167, 2014.
-
R. M. Slevinsky. On the use of Hahn’s asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform, arXiv:1602.02618, 2016.
-
A. Townsend, M. Webb, and S. Olver. Fast polynomial transforms based on Toeplitz and Hankel matrices, arXiv:1604.07486, 2016.