[ANN] `BandlimitedOperators.jl`: fast actions of bandlimited kernel matrices

Package Announcements
1

Hi all—

Just a quick announcement here of some code I’ve finally cleaned up into a proper package. BandlimitedOperators.jl gives you fast methods for implicitly representing the action of matrices of the form

pts1 = [...] # Vector{Float64} or Vector{SVector{D,Float64}} for 1 \leq D \leq 3
pts2 = [...] # same
dense_matrix = [fun(xj - xk) for xj in pts1, xk in pts2]

where fun is some bandlimited function. Examples include powers of sinc, or (up to numerical precision at least) Gaussian functions. Here is a quick example where fun is a Gaussian function with a bandwidth parameter that means that the dense matrix isn’t so rank degenerate that a simple low-rank approximation would work:

using BandlimitedOperators

# (No need to sort, it just makes a prettier matrix plot of M when you are debugging)
pts1 = sort(rand(1000)).*10
pts2 = sort(rand(1200)).*10

const ALPHA = 1000.0 # shape parameter for Gaussian function
gauss_kernel(t) = exp(-ALPHA*abs2(t))
M = [gauss_kernel(xj-xk) for xj in pts1, xk in pts2] # the matrix to approximate

# The Gaussian function isn't actually bandlimited, but its FT is also a
# Gaussian function and so numerically speaking it may as well be. The BANDWIDTH
# value here is the cutoff at which the FT goes below 1e-18.
const BANDWIDTH = sqrt(-(ALPHA/(pi^2))*log(sqrt(ALPHA/pi)*1e-18))
gauss_ft(omega) = sqrt(pi/ALPHA)*exp(-abs2(pi*omega)/ALPHA)

# This object now has mul!, *, size, and the usual standard implicit operator methods.
fastM = FastBandlimited(pts1, pts2, gauss_ft, BANDWIDTH)

x      = randn(length(pts2))
Mx     = M*x
fastMx = fastM*x
@show maximum(abs, Mx - fastMx) # ~ 1e-12

But internally the action of fastM comes at the cost of just two type 3 NUFFTs, and so you can represent the action of some seriously large matrices to basically machine precision:

julia> bigpts = rand(100_000).*10;

julia> bigfastM = FastBandlimited(bigpts, bigpts, gauss_ft, BANDWIDTH);

julia> x = randn(100_000);

julia> using BenchmarkTools

julia> @btime $bigfastM*$x; # nice.
  13.128 ms (61 allocations: 6.18 MiB)

The key to this algorithm being fast is to sort of represent the matrix as a Fourier multiplier. That is fast, but it also means that you need to use a sufficiently dense quadrature rule to resolve the highest oscillation—which depends on your point locations. So while this code will be very fast in many cases (turning O(n^2) work into O(n \log n)) and some clever ways to keep that quadrature rule size down exist, it can also be thwarted. If you have the 2 \times 2 matrices that look like

mean_matrix(n) = [fun(xj-xk) for xj in (1, n^2), xk in (1, n^2)],

for example, you can bonk the good complexity here. I would love to enforce checks like this in the code and throw an error rather than hand somebody a struct that mul!s poorly, but it’s not entirely obvious to me how best to do that. I would love feedback on how to integrate checks that prevent people from footguns like that.

Also, I would guess that there are other packages that offer similar functionality and I’d love to add a little segment to the README to point users to similar functionality in other places. Maybe such matrix operations can also be accelerated with ApproxFun, for example, @dlfivefifty?

Anyways. Happy matvec-ing and please do shoot me a note if you have thoughts on ways to improve this package. I’m really trying to get a bit better about writing code that is polished and thoughtful enough that it makes sense for other people to use it. I’ve got some tests (although no docstrings at least, which I know is bad and it is a high todo list item). But I’m definitely still in the learning process of releasing code beyond the standard of “people coming from the paper can run these scripts”, if you know what I mean, and so any feedback on how to set this up pro-socially would be very welcome.

Cheers,
CJG

5 Likes