# Compute double multipole expansion

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

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

1 Like

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?

Itâ€™s one of the examples HierarchicalMatrices.jl/Kernel.jl at master Â· JuliaMatrices/HierarchicalMatrices.jl Â· GitHub

1 Like

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?

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

1 Like

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