`FastTransforms.jl`

allows the user to conveniently work with orthogonal polynomials with degrees well into the millions.

This package provides a Julia wrapper for the C library of the same name. Additionally, all three types of nonuniform fast Fourier transforms are available, as well as the Padua transform.

The build script, which works on macOS and Linux systems with x86_64 processors, downloads precompiled binaries of the latest version of FastTransforms. This library depends on `FFTW`

, `MPFR`

, and `OpenBLAS`

(on Linux). Installation may be as straightforward as:

```
pkg> add FastTransforms
julia> using FastTransforms, LinearAlgebra
```

The 26 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:25)`

. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:

```
julia> c = rand(8192);
julia> leg2cheb(c);
julia> cheb2leg(c);
julia> norm(cheb2leg(leg2cheb(c; normcheb=true); normcheb=true)-c)/norm(c)
1.1866591414786334e-14
```

The implementation separates pre-computation into an `FTPlan`

. This type is constructed with either `plan_leg2cheb`

or `plan_cheb2leg`

. Let's see how much faster it is if we pre-compute.

```
julia> p1 = plan_leg2cheb(c);
julia> p2 = plan_cheb2leg(c);
julia> @time leg2cheb(c);
0.433938 seconds (9 allocations: 64.641 KiB)
julia> @time p1*c;
0.005713 seconds (8 allocations: 64.594 KiB)
julia> @time cheb2leg(c);
0.423865 seconds (9 allocations: 64.641 KiB)
julia> @time p2*c;
0.005829 seconds (8 allocations: 64.594 KiB)
```

Furthermore, for orthogonal polynomial connection problems that are degree-preserving, we should expect to be able to apply the transforms in-place:

```
julia> lmul!(p1, c);
julia> lmul!(p2, c);
julia> ldiv!(p1, c);
julia> ldiv!(p2, c);
```

Let `F`

be an array of spherical harmonic expansion coefficients with columns arranged by increasing order in absolute value, alternating between negative and positive orders. Then `sph2fourier`

converts the representation into a bivariate Fourier series, and `fourier2sph`

converts it back. Once in a bivariate Fourier series on the sphere, `plan_sph_synthesis`

converts the coefficients to function samples on an equiangular grid that does not sample the poles, and `plan_sph_analysis`

converts them back.

```
julia> F = sphrandn(Float64, 1024, 2047); # convenience method
julia> P = plan_sph2fourier(F);
julia> PS = plan_sph_synthesis(F);
julia> PA = plan_sph_analysis(F);
julia> @time G = PS*(P*F);
0.090767 seconds (12 allocations: 31.985 MiB, 1.46% gc time)
julia> @time H = P\(PA*G);
0.092726 seconds (12 allocations: 31.985 MiB, 1.63% gc time)
julia> norm(F-H)/norm(F)
2.1541073345177038e-15
```

Due to the structure of the spherical harmonic connection problem, these transforms may also be performed in-place with `lmul!`

and `ldiv!`

.

The NUFFTs are implemented thanks to Alex Townsend:

`nufft1`

assumes uniform samples and noninteger frequencies;`nufft2`

assumes nonuniform samples and integer frequencies;`nufft3 ( = nufft)`

assumes nonuniform samples and noninteger frequencies;`inufft1`

inverts an`nufft1`

; and,`inufft2`

inverts an`nufft2`

.

Here is an example:

```
julia> n = 10^4;
julia> c = complex(rand(n));
julia> ω = collect(0:n-1) + rand(n);
julia> nufft1(c, ω, eps());
julia> p1 = plan_nufft1(ω, eps());
julia> @time p1*c;
0.002383 seconds (6 allocations: 156.484 KiB)
julia> x = (collect(0:n-1) + 3rand(n))/n;
julia> nufft2(c, x, eps());
julia> p2 = plan_nufft2(x, eps());
julia> @time p2*c;
0.001478 seconds (6 allocations: 156.484 KiB)
julia> nufft3(c, x, ω, eps());
julia> p3 = plan_nufft3(x, ω, eps());
julia> @time p3*c;
0.058999 seconds (6 allocations: 156.484 KiB)
```

The Padua transform and its inverse are 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`

.

```
julia> n = 200;
julia> N = div((n+1)*(n+2), 2);
julia> v = rand(N); # The length of v is the number of Padua points
julia> @time norm(ipaduatransform(paduatransform(v)) - v)/norm(v)
0.007373 seconds (543 allocations: 1.733 MiB)
3.925164683252905e-16
```

[1] B. Alpert and V. Rokhlin. A fast algorithm for the evaluation of Legendre expansions, *SIAM J. Sci. Stat. Comput.*, **12**:158—179, 1991.

[2] N. Hale and A. Townsend. A fast, simple, and stable Chebyshev—Legendre transform using an asymptotic formula, *SIAM J. Sci. Comput.*, **36**:A148—A167, 2014.

[3] J. Keiner. Computing with expansions in Gegenbauer polynomials, *SIAM J. Sci. Comput.*, **31**:2151—2171, 2009.

[4] D. Ruiz—Antolín and A. Townsend. A nonuniform fast Fourier transform based on low rank approximation, arXiv:1701.04492, 2017.

[5] R. M. Slevinsky. On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform, *IMA J. Numer. Anal.*, **38**:102—124, 2018.

[6] R. M. Slevinsky. Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.

[7] R. M. Slevinsky, Conquering the pre-computation in two-dimensional harmonic polynomial transforms, arXiv:1711.07866, 2017.

[8] A. Townsend, M. Webb, and S. Olver. Fast polynomial transforms based on Toeplitz and Hankel matrices, in press at *Math. Comp.*, 2017.

02/07/2016

7 days ago

318 commits