Kernel matrices

While in the introduction we presented a somewhat general way to assemble an HMatrix, abstract matrices associated with an underlying kernel function are common enough in boundary integral equations that a special interface exists for facilitating their use.

The AbstractKernelMatrix interface is used to represent matrices K with i,j entry given by f(X[i],Y[j]), where X=rowelements(K) and Y=colelements(K). The row and columns elements may be as simple as points in $\mathbb{R}^d$ (as is the case for Nyström methods), but they can also be more complex objects such as triangles or basis functions. In such cases, it is required that center(X[i]) and center(Y[j]) return a point as an SVector, and that f(X[i],Y[j]) make sense.

A concrete implementation of AbstractKernelMatrix is provided by the KernelMatrix type. Creating the matrix associated with the Helmholtz free-space Greens function, for example, can be accomplished through:

using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
X = rand(Point3D,10_000)
Y = rand(Point3D,10_000)
const k = 2π
function G(x,y)
    d = norm(x-y)
    exp(im*k*d)/(4π*d)
end
K = KernelMatrix(G,X,Y)

10000×10000 KernelMatrix{typeof(Main.var"Main".G), Vector{StaticArraysCore.SVector{3, Float64}}, Vector{StaticArraysCore.SVector{3, Float64}}, ComplexF64}:
 -0.0453107-0.106765im      -0.110609-0.0751819im   …   -0.168167+0.0662665im
  -0.111832-0.0741135im     0.0372654-0.0906286im       0.0791568-0.0289281im
  0.0676824-0.0590757im      0.078971+0.00336069im      -0.168207+0.055641im
  -0.011236-0.107909im     -0.0594146-0.103454im        0.0547437-0.0760409im
  0.0789261+0.00358147im   -0.0412307-0.10741im         0.0273735+0.0614285im
  -0.139238+0.173598im      0.0698619-0.0552364im   …   -0.104981+0.226147im
  -0.158897-0.000804682im   -0.107067-0.0781464im       -0.139529-0.0418062im
  -0.139419-0.0419757im     0.0782592+0.00655391im     -0.0306436-0.108443im
  -0.163266+0.106986im      -0.168243+0.0622898im       -0.069906-0.0998659im
  0.0462629-0.0839558im     -0.148078-0.0269882im       0.0568523-0.0737583im
           ⋮                                        ⋱  
  0.0261996-0.0970027im     -0.114263-0.0719186im      -0.0604536-0.103142im
  0.0594282-0.0707599im     -0.118276-0.068076im       -0.0978604-0.0850243im
  -0.154278-0.0134622im      -0.16809+0.068532im        0.0220981-0.0989382im
  0.0592928-0.0709237im    -0.0840872-0.0933689im       0.0461978-0.0840098im
   0.026941-0.0966298im      0.330132+0.432549im    …   0.0568319-0.0737811im
 -0.0972792-0.0854213im    -0.0558381+0.276954im        -0.102151-0.0819603im
  0.0748691-0.0442438im     0.0163102-0.101323im        0.0688802-0.0570179im
 -0.0795635-0.095658im      -0.107801-0.0775474im       0.0530888-0.0777354im
  -0.166649+0.0865839im     -0.107644-0.0776759im       0.0119131-0.102882im

Compressing K is now as simple as:

H = assemble_hmatrix(K;rtol=1e-6)

HMatrix of ComplexF64 with range 1:10000 × 1:10000
	 number of nodes in tree: 14681
	 number of leaves: 11011 (2776 admissible + 8235 full)
	 min rank of sparse blocks : 20
	 max rank of sparse blocks : 86
	 min length of dense blocks : 450
	 max length of dense blocks : 2850
	 min number of elements per leaf: 450
	 max number of elements per leaf: 445550
	 depth of tree: 9
	 compression ratio: 2.465165

It is worth noting that several default choices are made during the compression above. See the Assembling an HMatrix section or the documentation of assemble_hmatrix for information on how to obtain a more granular control of the assembling stage.

As before, you can multiply H by a vector, or do an lu factorization of it.

Finally, here is a somewhat contrived example of how to use a KernelMatrix when the rowelements and colelements are not simply points (as required e.g. in a Galerkin discretization of boundary integral equations):

using HMatrices, LinearAlgebra, StaticArrays
# create a simple structure to represent a segment
struct Segment
    start::SVector{2,Float64}
    stop::SVector{2,Float64}
end
# extend the function center to work with segments
HMatrices.center(s::Segment) = 0.5*(s.start + s.stop)
# P1 mesh of a circle
npts = 10_000
nodes = [SVector(cos(s), sin(s)) for s in range(0,stop=2π,length=npts+1)]
segments = [Segment(nodes[i],nodes[i+1]) for i in 1:npts]
# Now define a kernel function that takes two segments (instead of two points) and returns a scalar
function G(target::Segment,source::Segment)
    x, y = HMatrices.center(target), HMatrices.center(source)
    d = norm(x-y)
    return -log(d + 1e-10)
end
K = KernelMatrix(G,segments,segments)
# compress the kernel matrix
H = assemble_hmatrix(K;rtol=1e-6)

HMatrix of Float64 with range 1:10000 × 1:10000
	 number of nodes in tree: 3341
	 number of leaves: 2506 (798 admissible + 1708 full)
	 min rank of sparse blocks : 4
	 max rank of sparse blocks : 7
	 min length of dense blocks : 625
	 max length of dense blocks : 4300
	 min number of elements per leaf: 625
	 max number of elements per leaf: 2775556
	 depth of tree: 9
	 compression ratio: 23.061837

Support for tensor kernels

For vector-valued partial differential equations such as Stokes or time-harmonic Maxwell's equation, the underlying integral operator has a kernel function which is a tensor. This package currently provides some limited support for these types of operators. The example below illustrates how to build an HMatrix representing a KernelMatrix corresponding to Stokes Greens function for points on a sphere:

using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
m = 5_000
X = Y = [Point3D(sin(θ)cos(ϕ),sin(θ)*sin(ϕ),cos(θ)) for (θ,ϕ) in zip(π*rand(m),2π*rand(m))]
const μ = 5
function G(x,y)
    r = x-y
    d = norm(r) + 1e-10
    1/(8π*μ) * (1/d*I + r*transpose(r)/d^3)
end
K = KernelMatrix(G,X,Y)
H = assemble_hmatrix(K;atol=1e-4)

HMatrix of StaticArraysCore.SMatrix{3, 3, Float64, 9} with range 1:5000 × 1:5000
	 number of nodes in tree: 4005
	 number of leaves: 3004 (812 admissible + 2192 full)
	 min rank of sparse blocks : 7
	 max rank of sparse blocks : 28
	 min length of dense blocks : 144
	 max length of dense blocks : 11368
	 min number of elements per leaf: 144
	 max number of elements per leaf: 103632
	 depth of tree: 10
	 compression ratio: 3.307420

You can now multiply H by a density σ, where σ is a Vector of SVector{3,Float64}

σ = rand(SVector{3,Float64},m)
y = H*σ
# test the output agains the exact value for a given `i`
i = 42
exact = sum(K[i,j]*σ[j] for j in 1:m)
@show y[i] - exact

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 4.3138861656188965e-6
 8.121132850646973e-6
 1.1026859283447266e-6