# Compute double multipole expansion

**URL:** <https://discourse.julialang.org/t/compute-double-multipole-expansion/67059>\
**Category:** Numerics\
**Created:** [August 26, 2021, 8:56pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059 "2021-08-26T20:56:07Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 26, 2021, 8:56pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/1 "2021-08-26T20:56:07Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [August 29, 2021, 11:04am UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/2 "2021-08-29T11:04:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 30, 2021, 5:14pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/3 "2021-08-30T17:14:33Z")

</div>

> [@dlfivefifty](#):
>
> HierarchicalMatrices.jl

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?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [August 30, 2021, 5:15pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/4 "2021-08-30T17:15:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 30, 2021, 5:39pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/5 "2021-08-30T17:39:15Z")

</div>

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:

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [August 30, 2021, 5:40pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/6 "2021-08-30T17:40:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [August 30, 2021, 6:53pm UTC](https://discourse.julialang.org/t/compute-double-multipole-expansion/67059/7 "2021-08-30T18:53:06Z")

</div>

Just passing through with another option…

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