Using $\mathcal{H}^2$-Matrices

PetrovGalerkinNCA constructs an $\mathcal{H}^2$-matrix from selected entries of a kernel matrix. Dense blocks store near interactions, while nested bases and small coupling matrices represent the far field. The result implements Julia's matrix-vector product interface.

Point kernels provide a self-contained way to use and test the nested construction without a boundary-element space or quadrature rule. Define a callable kernel with an eltype, cluster its test and trial positions, and pass the kernel, positions, and block tree to PetrovGalerkinNCA.

The currently tested scope comprises Helmholtz point kernels and the EFIE single-layer operator through BEAST.jl. The EFIE workflow uses the same constructor and is shown under Application Examples. Double-layer operators are not currently part of the tested NestedCrossApproximation workflow.

Point-Kernel Matrices

Define a callable kernel and its matrix entry type. The two point sets below are separated so that the kernel has no singular entries.

using NestedCrossApproximation: PetrovGalerkinNCA, TopDown, storage
using AdaptiveCrossApproximation
using H2Trees
using LinearAlgebra
using ParallelKMeans
using Random
using StaticArrays

struct HelmholtzKernel
    k::Float64
end

Base.eltype(::HelmholtzKernel) = ComplexF64
function (kernel::HelmholtzKernel)(x, y)
    r = norm(x - y)
    return cis(kernel.k * r) / (4pi * r)
end

testpoints = [
    SVector(x, y, 0.0)
    for x in range(0, 1; length=32), y in range(0, 1; length=32)
][:]
trialpoints = testpoints .+ Ref(SVector(1.1, 0.0, 0.2))
kernel = HelmholtzKernel(2pi / 10)

The block tree combines one cluster tree for the test points with one for the trial points. For general use, a k-means tree is the recommended backend because its balanced clusters adapt to the geometry. Loading ParallelKMeans.jl enables this backend in H2Trees.jl.

testbuilder = KMeansTreeBuilder(;
    numberofclusters=2,
    minvalues=32,
    splitterkwargs=(; rng=MersenneTwister(1)),
)
trialbuilder = KMeansTreeBuilder(;
    numberofclusters=2,
    minvalues=32,
    splitterkwargs=(; rng=MersenneTwister(1)),
)
tree = BlockTree(
    KMeansTree(testpoints; builder=testbuilder),
    KMeansTree(trialpoints; builder=trialbuilder),
)

numberofclusters is the number of children created when a cluster is split, and minvalues limits the number of points in a leaf. Smaller leaves expose more compression opportunities; larger leaves reduce tree depth and create larger dense terminal blocks.

Each side of the matrix has its own nested basis. The testcompressor constructs the test basis and the trialcompressor constructs the trial basis. A direct first choice is TopDown with ACA:

tol = 1e-4
h2mat = PetrovGalerkinNCA(
    kernel,
    testpoints,
    trialpoints,
    tree;
    testcompressor=TopDown(; factorization=ACA(; tol=tol)),
    trialcompressor=TopDown(; factorization=ACA(; tol=tol)),
)
size(h2mat)

(1024, 1024)

This example has 1,024 points on each side and a nonempty compressed far field. storage reports the stored memory and its ratio to a dense matrix:

storage(h2mat)

storage: 0.0011943680000000001 GB
summary size: 0.001460345 GB
compression ratio: 0.07118988037109375

The factorization tolerance controls every local basis construction. The error of the complete matrix is usually of the same order, but is not bounded by tol because local errors can accumulate through the cluster tree.

Working with the Result

PetrovGalerkinNCA is a LinearMap and supports ordinary, transposed, and adjoint products. A dense matrix is useful for checking a small problem:

dense = [kernel(t, s) for t in testpoints, s in trialpoints]
Random.seed!(2)
x = randn(ComplexF64, length(trialpoints))
z = randn(ComplexF64, length(testpoints))

norm(h2mat * x - dense * x) / norm(dense * x)

0.00034903619423811284

norm(transpose(h2mat) * z - transpose(dense) * z) / norm(transpose(dense) * z)

0.00028565216851535045

norm(h2mat' * z - dense' * z) / norm(dense' * z)

0.00035818987254865575

The compressed far field and the dense near field can be applied separately:

using NestedCrossApproximation: farmatrix, nearmatrix

far = farmatrix(h2mat)
near = nearmatrix(h2mat)
norm(far * x + near * x - h2mat * x) / norm(h2mat * x)

5.506565825509466e-16

Concurrent Products

h2mat consists of a read-only core, which holds the near field and the nested far-field representation, and a workspace with the buffers of the matrix-vector product. The workspace is the LinearMap, so the buffers are allocated once and reused by every product. One workspace runs one product at a time; a second product on the same workspace throws a ConcurrentWorkspaceError. To multiply with the same matrix from several tasks, create one workspace per task over the shared core:

using NestedCrossApproximation: H2Workspace, core
using OhMyThreads: SerialScheduler

workspaces = [H2Workspace(core(h2mat); scheduler=SerialScheduler()) for _ in 1:2]
inputs = [randn(ComplexF64, length(trialpoints)) for _ in workspaces]
tasks = [Threads.@spawn ws * xi for (ws, xi) in zip(workspaces, inputs)]

maximum(norm(fetch(t) - h2mat * xi) / norm(h2mat * xi) for (t, xi) in zip(tasks, inputs))

0.0

The core alone is assembled with PetrovGalerkinNCACore, which takes the arguments of PetrovGalerkinNCA.

Choosing a Compressor

A compressor combines the tree traversal with a low-rank factorization from AdaptiveCrossApproximation.jl.

BottomUp with IACA

BottomUp compresses a cluster against its complete far field and restricts the candidates of every parent to pivots inherited from its children. Consequently, tree-aware geometric pivoting is the natural choice. The pivoting on the test side uses the trial tree, and vice versa:

using NestedCrossApproximation: BottomUp

convergence = FNormExtrapolator(tol)
testpivoting = TreeMimicryPivoting(
    testpoints, trialpoints, H2Trees.trialtree(tree)
)
trialpivoting = TreeMimicryPivoting(
    trialpoints, testpoints, H2Trees.testtree(tree)
)

testcompressor = BottomUp(;
    factorization=IACA(MaximumValue(), testpivoting, convergence),
)
trialcompressor = BottomUp(;
    factorization=IACA(trialpivoting, MaximumValue(), convergence),
)

MaximumValue() selects a pivot from an evaluated residual row or column. TreeMimicryPivoting selects the complementary pivot geometrically, without first evaluating the corresponding complete residual vector. This is the incomplete adaptive cross approximation described in [4].

Other Configurations

TopDown can also use IACA. Its far field is a set of points rather than tree nodes, so it requires MimicryPivoting:

testcompressor = TopDown(;
    factorization=IACA(
        MaximumValue(),
        MimicryPivoting(testpoints, trialpoints),
        convergence,
    ),
)
trialcompressor = TopDown(;
    factorization=IACA(
        MimicryPivoting(trialpoints, testpoints),
        MaximumValue(),
        convergence,
    ),
)

Both traversal directions accept ordinary ACA as well:

bottomup_aca = BottomUp(; factorization=ACA(; tol=tol))

BatchedBottomUp first selects bounded far-field subsets, assembles the resulting small blocks together, and then applies a dense fully pivoted ACA to every block. It is intended for kernel backends whose main advantage is batched evaluation.

Far-Field Representors

A representor reduces the complete far field of a cluster before ACA or IACA is applied. This can substantially reduce the number of matrix entries evaluated by TopDown or BottomUp. The resulting factorization tests convergence only on the representative block, so the reduction can also increase the approximation error.

Three strategies are available:

The representor is passed to each compressor independently. RandomRep needs no geometrical data:

using NestedCrossApproximation: RandomRep

testcompressor = TopDown(;
    factorization=ACA(; tol=tol),
    representor=RandomRep(),
)
trialcompressor = TopDown(;
    factorization=ACA(; tol=tol),
    representor=RandomRep(),
)

ChebyshevRep receives the positions on the far side of the corresponding compressor. The test basis is compressed against the trial far field, while the trial basis is compressed against the test far field:

using NestedCrossApproximation: ChebyshevRep

testcompressor = BottomUp(;
    factorization=ACA(; tol=tol),
    representor=ChebyshevRep(trialpoints),
)
trialcompressor = BottomUp(;
    factorization=ACA(; tol=tol),
    representor=ChebyshevRep(testpoints),
)

MimicryRep wraps a TreeMimicryPivoting constructed with the tree on the far side. It works only with BottomUp, because that compressor retains the cluster-node representation of the inherited far field:

using NestedCrossApproximation: MimicryRep

testrepresentor = MimicryRep(
    TreeMimicryPivoting(testpoints, trialpoints, H2Trees.trialtree(tree)),
)
trialrepresentor = MimicryRep(
    TreeMimicryPivoting(trialpoints, testpoints, H2Trees.testtree(tree)),
)

testcompressor = BottomUp(;
    factorization=ACA(; tol=tol),
    representor=testrepresentor,
)
trialcompressor = BottomUp(;
    factorization=ACA(; tol=tol),
    representor=trialrepresentor,
)

BatchedBottomUp has no separate representor keyword. Its first phase already performs a geometry-based reduction of every far field through its pivotselector.

Admissibility and Alternative Tree Backend

By default, a pair of clusters $t$ and $s$ enters the far field when

\[\max\{\operatorname{diam}(t),\operatorname{diam}(s)\} \leq \eta\,\operatorname{dist}(t,s),\]

with $\eta=1$. Pass isnear(η) to change the separation parameter:

h2mat = PetrovGalerkinNCA(
    kernel,
    testpoints,
    trialpoints,
    tree;
    isnear=isnear(0.5),
    testcompressor,
    trialcompressor,
)

A smaller η demands a larger separation and therefore stores more interactions in the dense near field. For electrically large clusters, use the mixed-frequency condition described in Using Directional $\mathcal{H}^2$-Matrices.

The basic example uses the recommended geometry-adaptive k-means tree. A TwoNTree can also be used. It provides a lightweight quadtree or octree for regular point sets and small examples:

twonbuilder = TwoNTreeBuilder(; minvalues=32)
tree = BlockTree(
    testpoints,
    trialpoints;
    builder=BlockTreeBuilder(; test=twonbuilder, trial=twonbuilder),
)

Assembly Options

The principal keywords of PetrovGalerkinNCA are:

Reaching maxrank can indicate that the requested local tolerance was not attained. Assembly and matrix-vector products run on ntasks parallel tasks; ntasks=1 runs on the calling task. The scheduling itself is chosen internally.

For the nested representation and the construction of its bases, transfer matrices, and coupling matrices, see the $\mathcal{H}^2$-matrix theory.