GitHub

IFGF.jl

A Julia implementation of the Interpolated Factored Green Function (IFGF) Method on adaptive trees

Overview

This package implements an algorithm for approximating the forward map (i.e. the matrix-vector product) of certain kernel matrices, where kernel matrix will be used to refer to any matrix $A \in \mathbb{F}^{m\times n}$ whose entries can be computed from the following three components:

  • a vector of target elements $X = \left\{ \boldsymbol{x}_i \right\}_{i=1}^m$
  • a vector of source elements $Y = \left\{ \boldsymbol{y}_i \right\}_{i=1}^n$
  • a kernel function $K(x,y) : X \times Y \to \mathbb{F}$, where $\mathbb{F}$ is the return type of the kernel (and of the underlying matrix)

The entries of a kernel matrix are given by the explicit formula $A_{i,j} = K(\boldsymbol{x}_i,\boldsymbol{y}_j)$. Typically, $X$ and $Y$ are vectors of points in $\mathbb{R}^d$, and $K$ is a function mapping elements of $X$ and $Y$ into $\mathbb{R}$ or $\mathbb{C}$.

Provided $K$ has sufficient structure, an approximation of the product $Ax$ can be computed in linear or log-linear complexity by various well known algorithms such as the Fast Multipole Method or hierarchical matrices; this package implements another acceleration algorithm termed Interpolated Factored Green Function (IFGF) method.

IFGF.jl comes loaded with a few predefined kernels arising in mathematical physics; check out the predefined kernels section for more details. Otherwise, setting up your own kernel is straightforward, as shown next.

Simple example

To illustrate how to set up the IFGFOp, let $X$ and $Y$ be a set of random points on the unit cube, and take $K$ to be the Helmholtz Green function in three dimensions:

\[ K(\boldsymbol{x},\boldsymbol{y}) = \frac{e^{ik|\boldsymbol{x}-\boldsymbol{y}|}}{4\pi |\boldsymbol{x}-\boldsymbol{y}|}\]

Setting up the aforementioned kernel matrix can be done as follows:

using IFGF, LinearAlgebra, StaticArrays
import IFGF: wavenumber

# random points on a cube
const Point3D = SVector{3,Float64}
m = 100_000
X = Y = rand(Point3D,m)

# define a the kernel matrix
struct HelmholtzMatrix <: AbstractMatrix{ComplexF64}
    X::Vector{Point3D}
    Y::Vector{Point3D}
    k::Float64
end

# indicate that this is an ocillatory kernel with wavenumber `k`
wavenumber(A::HelmholtzMatrix) = A.k

# functor interface
function (K::HelmholtzMatrix)(x,y)
    k = wavenumber(K)
    d = norm(x-y)
    return (d!=0) * (exp(im*k*d)*inv(4*pi*d)) # skip x == y case
end

# abstract matrix interface
Base.size(::HelmholtzMatrix) = length(X), length(Y)
Base.getindex(A::HelmholtzMatrix,i::Int,j::Int) = A(A.X[i],A.Y[j])

# create the abstract matrix
k = 4π
A = HelmholtzMatrix(X,Y,k)
100000×100000 Main.var"Main".HelmholtzMatrix:
        0.0-0.0im          -0.319991+0.241123im   …    0.122104-0.137472im
  -0.319991+0.241123im           0.0-0.0im             0.116231+0.0858827im
 -0.0903114+0.0696573im    0.0556752+0.122581im      -0.0468605-0.0830986im
 -0.0153609-0.214921im     -0.183652-0.177895im       0.0825934-0.175653im
  0.0149771-0.0883162im    0.0565182-0.0639255im       0.156611+0.0127257im
  0.0787316+0.00686576im   0.0727485-0.0396898im  …  -0.0883295-0.0484428im
  -0.164238-0.188781im     -0.108167-0.209809im        0.104013-0.157602im
  0.0394206-0.0778096im    -0.100152-0.0255628im      0.0792438-0.0161379im
 -0.0724591-0.0665069im   -0.0132802-0.0911958im       0.147614-0.093847im
  -0.106219+0.00115193im   -0.106707+0.0110962im      0.0453872+0.125251im
           ⋮                                      ⋱  
  -0.175487-0.182724im      0.108718-0.152884im       -0.327314+0.210196im
  0.0650109+0.0398307im    0.0551536-0.0653244im      0.0691201-0.0469226im
  -0.105072-0.00731373im   -0.052421+0.10841im       -0.0916828+0.0673264im
  -0.026642+0.0634326im    0.0340724+0.0648345im      0.0694673-0.0463002im
 -0.0622265+0.0182138im   -0.0435138+0.0514577im  …   0.0767115-0.028724im
  -0.116203-0.207668im    -0.0526689-0.217195im       -0.269688-0.0973112im
 0.00624209-0.0901617im   -0.0432379-0.0846217im     -0.0431978+0.113777im
   0.075009+0.0209778im    0.0243798-0.0852832im      0.0742709-0.0360174im
   0.122104-0.137472im      0.116231+0.0858827im            0.0-0.0im

Note that A is a lazy $100000 \times 100000$ matrix which computes its entries on demand; as such, it has a light memory footprint. Trying to instantiate the underlying Matrix{ComplexF64} would most likely result in you running out of memory.

To build an approximation of A, simply do:

L = assemble_ifgf(A,X,Y; tol = 1e-4)
IFGFOp of ComplexF64 with range 1:100000 × 1:100000
|-- interpolation points per cone: (8, 8, 8)
|-- maximal meshsize per cone:     (1.57, 1.57, 1.57)
|-- number of nodes:               585
|-- number of leaves:              512
|-- number of active cones:        4452

Check the documentation of assemble_ifgf for more information on the available options.

You can now use L in lieu of A to approximate the matrix vector product, as illustrated below:

x     = randn(ComplexF64,m)
y     = L*x
100000-element Vector{ComplexF64}:
  19.679756692626356 - 27.475718636498918im
  15.790380080351788 + 2.8523393769358236im
  16.518111089884066 + 53.517249809681154im
 -15.111061766664669 + 50.454464689316325im
   41.69729629713238 + 3.9214497660929615im
 -3.3730864383547683 - 82.09521537776644im
  28.985195189388364 - 23.194087914403667im
 -13.638895461433673 - 15.435563196775515im
   58.93785966068788 + 6.771709224627459im
  -84.02565099297551 - 7.3117677043074725im
                     ⋮
   38.85131856936423 + 23.6555155954715im
  -70.44113771746518 + 22.331067171004804im
 -1.7261415759992262 + 31.368614903416052im
 -39.249929839001865 + 33.40263665839471im
   80.68422961182615 - 54.85077621703655im
  58.396468505181076 + 32.88707189942526im
   1.021853395208542 + 39.277753722245905im
   -93.2468311879816 - 42.69551722288213im
  26.111567375708102 - 11.254631797923249im

Note that L, while not as compact as A, still has a relatively light memory footprint:

gb = Base.summarysize(L) / 1e9
print("Size of L: $gb gigabytes")
Size of L: 0.072876012 gigabytes

The quality of the approximation may be verified by computing the exact value at a few randomly selected set of indices:

I     = rand(1:m,100)
exact = [sum(A[i,j]*x[j] for j in 1:m) for i in I] # matrix-vector product
er    = norm(y[I]-exact) / norm(exact) # relative error
print("approximate relative error: $er")
approximate relative error: 1.83766146631823e-5

To keep things simple, we neglected some important optimizations that can be performed when defining your own kernel. See the custom kernels for more information.