HMatrices.jl

A package for assembling and factoring hierarchical matrices

Install from the Pkg REPL:

pkg> add HMatrices

This package provides some functionality for assembling as well as for doing linear algebra with hierarchical matrices with a strong focus in applications arising in boundary integral equation methods.

For the purpose of illustration, let us consider an abstract matrix K with entry i,j given by the evaluation of some kernel function G on points X[i] and Y[j], where X and Y are vector of points (in 3D here); that is, K[i,j]=G(X[i],Y[j]). This object can be constructed as follows:

using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
# sample some points on a sphere
m = 100_000
X = Y = [Point3D(sin(θ)cos(ϕ),sin(θ)*sin(ϕ),cos(θ)) for (θ,ϕ) in zip(π*rand(m),2π*rand(m))]
function G(x,y) 
  d = norm(x-y) + 1e-8
  1/(4π*d)
end
K = KernelMatrix(G,X,Y)

where we took G to be the free-space Greens function of Laplace's equation in 3D (to avoid division-by-zero we added 1e-8 to the distance between points).

The object K corresponds to a dense matrix, so converting it to a matrix can be costly both in terms of memory and flops. Instead, we can construct an approximation to K as a hierarchical matrix using:

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

Tip: For a smaller problem size (say m=10_000), you may try

using Plots
plot(H)

to visualize the underlying block-structure. You should see something similar to the figure below:

Calling HMatrices.compression_ratio(H) reveals that storing a dense version of K would take roughly 25 times as much space (and probably would not fit in most laptops). We can now use H as an approximation to K for some linear algebra operations, such as:

x = rand(m)
y = H*x

To check that this is indeed an approximation, we can compare against the exact value at a given entry:

y[42] - sum(K[42,j]*x[j] for j in 1:m)
# about 2e-7

It is also possibly to factor H by calling e.g. lu(H;atol=1e-6) (this may take a few minutes on a reasonable machine for the 100_000 × 100_000 problem size and specified tolerance). The result is an LU factorization object with a hierarchical low-rank structure, and the factored object can be used both in a direct solver or as a preconditioner for H in an iterative solver.

For more information, see the documentation.

Below are some good references on hierarchical matrices and their application to boundary integral equations:

[1] Hackbusch, Wolfgang. Hierarchical matrices: algorithms and analysis. Vol. 49. Heidelberg: Springer, 2015.

[2] Bebendorf, Mario. Hierarchical matrices. Springer Berlin Heidelberg, 2008.

If you are interested in hierarchical matrices and Julia, check out also the following packages:

• HierarchicalMatrices.jl: a flexible framework for hierarchical matrices implementing an abstract infrastructure.
• KernelMatrices.jl: a library implementing the Hierarchically Off-Diagonal Low-Rank structure (HODLR).
• HSSMatrices.jl: an implementation of the Hierarchically semi-separable structure (HSS).