Compute double multipole expansion

Specific Domains Numerics
1

Hello,

I am working on the calculation of the pressure induced at d locations \{\xi_1, \ldots, \xi_d\} due to a set of n singularities located at \{z_1, \ldots, z_n\} with strength \{ S_1, \ldots, S_n\}, where \xi_i, z_J, S_J are complex numbers.

I end up having to compute the following double sum:
\sum_J \sum_{K \neq J} \frac{S_J}{2\pi} \frac{\overline{S}_K}{2\pi} \frac{1}{\xi_i - z_J} \frac{1}{\overline{z}_J -\overline{z}_K} for all i = 1,\ldots, d

Let define the matrices C_{yx} \in \mathbb{C}^{d \times n} as [C_{yx}]_{i,j} = 1/(\xi_i - z_j), C_{xx} \in \mathbb{C}^{d \times n} as [C_{xx}]_{i,j} = (1-\delta_{i,j})/(z_i - z_j), and the vector S = [S_1, \ldots, S_n] \in \mathbb{C}^n.

We can compute this double sum for all the \xi_i's as the matrix product:
\frac{1}{4\pi^2}C_{yx} \; \text{Diagonal}(S) \; \overline{C_{xx}} \; \overline{S}.

This product can be efficiently computed by applying twice the fast multipole method. Do you know if there are “double fast multipole” methods that perform the double summation at once?

Also, which fast multipole Julia packages would you recommend?

1 Like
2

I don’t know of a fast multipolar pack but there’s HierarchicalMatrices.jl which may be of use

1 Like
3

Thank you @dlfivefifty for your answer! Can I use HierarchicalMatrices.jl for the Cauchy kernel k(x, y) = 1/(x - y) with x, y complex numbers?

4

It’s one of the examples https://github.com/JuliaMatrices/HierarchicalMatrices.jl/blob/master/examples/Kernel.jl

1 Like
5

I saw it but I wasn’t sure if complex numbers are handled as well.

For real entries, we can define the kernel matrix as follows:

using HierarchicalMatrices

@inline cauchykernel(x::T, y::T) where T = inv(x-y)
cauchykernel(X::Vector{T}, Y::Vector{T}) where T = T[cauchykernel(x, y) for x in X, y in Y]

N = 1000
x = Float64[i*(i+1) for i = N:-1:1]
y = Float64[(i+1/2)*(i+3/2) for i = N:-1:1]
xmin, xmax = extrema(x)
ymin, ymax = extrema(y)
K = KernelMatrix(cauchykernel, x, y, xmax, xmin, ymax, ymin)

But the signature of KernelMatrix is

KernelMatrix(::Function, ::Vector{T}, ::Vector{T}, ::T, ::T, ::T, ::T) where T. There is no notion of order for complex numbers, I don’t know what entries to put for xmax, xmin, ymax, ymin in this case.

N = 1000
x = randn(ComplexF64,N)
y = randn(ComplexF64,N)

K = KernelMatrix(cauchykernel, x, y, xmax, xmin, ymax, ymin)

Also, I ran your example script, but I don’t see any speed-up in the calculation of Ab with HierarchicalMatrices.jl. Isn’t HierarchicalMatrices using some kind of multipole expansion of the kernel to accelerate matrix vector products?

6

It does a low rank approximation which is essentially the same thing

1 Like
7

Just passing through with another option…

using KernelMatrices, KernelMatrices.HODLR, BenchmarkTools

cauchykernel(x,y) = inv(x-y)

const N = 10000
const x = Float64[i*(i+1) for i in N:-1:1]
const y = Float64[(i+1/2)*(i+3/2) for i in N:-1:1]

const K     = KernelMatrix(x, y, cauchykernel)
const trueK = KernelMatrices.full(K)
const H     = HODLR.RKernelHODLR(K, 1e-12)

const testx = randn(N)

print("Assembly time to 1e-12 ACA blocks:")
@btime HODLR.RKernelHODLR($K, 1e-12) # 69 ms, a ton of allocs :-(

print("Accuracy for random input vector: ")
println(norm(H*testx - trueK*testx)) # 3.5e-13

println("Matvec speeds:")
print("Exact:")
@btime trueK*testx # 16.2 ms, 2 alloc
print("HODLR:")
@btime H*testx # 1.12 ms, 984 allocs :-(
println()

HierarchicalMatrices.jl totally smokes KernelMatrices.jl in terms of code quality, which you can see in the allocs. But last I checked, KernelMatrices.jl is actually still neck and neck with it in terms of speed and resulting accuracy.

You’ll note that this calls the RKernelHODLR function and not the constructor in the readme. This is an alternative recursive struct I kind of quickly threw together at one point because I realized that the package was otherwise entirely focused on symmetric matrices. That code has gotten a lot less optimization attention than some other parts. But as you can see, when N gets big enough it’s still pretty helpful.

If you cranked up the N I’m sure the HierarchicalMatrices.jl example would also be much faster too, by the way. BLAS 2 operations are pretty darn fast, so you need reasonably large matrices before the overhead of fancier methods like this becomes worth it, in my experience. You’ll also probably only benefit from this if you’re going to apply the same matrix to many inputs.

Sorry, one more edit: the efficiency of these representations depends a lot on how you order the points. KernelMatrices.jl really doesn’t handle that for you, and assumes that you’re putting points in in the order that you want used for the matrix construction. For 2D and 3D points, I think probably the best off-the-shelf general-purpose ordering I’ve seen is using the Hilbert sorting method in GeometricalPredicates.jl. But just a mention to take care of that, because I’ve had several experiences of users telling me that it’s totally broken and not working because they passed the spatial locations in completely unsorted.

2 Likes