## FastTransforms.jl

:rocket: Julia package for orthogonal polynomial transforms :snowboarder:
Popularity
228 Stars
Updated Last
7 Months Ago
Started In
February 2016

# FastTransforms.jl

`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.

## Installation

Installation, which uses BinaryBuilder for all of Julia's supported platforms (in particular Sandybridge Intel processors and beyond), may be as straightforward as:

```pkg> add FastTransforms

julia> using FastTransforms, LinearAlgebra
```

## Fast orthogonal polynomial transforms

The orthogonal polynomial transforms are listed in `FastTransforms.Transforms` or `FastTransforms.kind2string.(instances(FastTransforms.Transforms))`. 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:

### The Chebyshev--Legendre transform

```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);
```

### The spherical harmonic transform

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!`.

See also FastSphericalHarmonics.jl for a simpler interface to the spherical harmonic transforms defined in this package.

## Nonuniform fast Fourier transforms

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

0.007373 seconds (543 allocations: 1.733 MiB)
3.925164683252905e-16
```

# References:

 D. Ruiz—Antolín and A. Townsend. A nonuniform fast Fourier transform based on low rank approximation, SIAM J. Sci. Comput., 40:A529–A547, 2018.

 S. Olver, R. M. Slevinsky, and A. Townsend. Fast algorithms using orthogonal polynomials, Acta Numerica, 29:573—699, 2020.

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

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

### Required Packages

View all packages