Using Hierarchical Matrices

An HMatrix stores near-field interactions as dense blocks and compresses admissible far-field interactions with a low-rank compressor. It supports matrix-vector products, transpose and adjoint products, sums and products of H-matrices, LU factorization and solves, and conversion to a dense Matrix.

Basic Usage

For BEAST operators, AdaptiveCrossApproximation.assemble constructs the cluster tree and quadrature data automatically:

using AdaptiveCrossApproximation
using BEAST
using CompScienceMeshes
using H2Trees

mesh = CompScienceMeshes.meshsphere(1.0, 0.3)
space = raviartthomas(mesh)
operator = Maxwell3D.singlelayer(; wavenumber=1.0)

hmat = AdaptiveCrossApproximation.assemble(
    operator,
    space,
    space;
    tol=1e-4,
    maxrank=40,
)

The result contains the directly assembled near field and the ACA-compressed far field. The operator-specific default compressor is used unless compressor is provided explicitly. For supported BEAST double-layer operators, this is the robust ACA configuration combining maximum-value and random-sampling pivoting.

BEAST local operators such as Identity and NCross are assembled directly and do not produce an HMatrix, because their sparse local interactions do not benefit from hierarchical low-rank compression.

Point-Kernel Matrices

Point kernels provide a self-contained way to experiment with hierarchical matrices without a boundary-element space or quadrature rule. Define a callable kernel with an eltype, cluster its points, and pass both to HMatrix:

using AdaptiveCrossApproximation
using H2Trees
using LinearAlgebra
using ParallelKMeans
using StaticArrays

struct SoftenedKernel{T}
    radius::T
end

Base.eltype(::SoftenedKernel{T}) where {T} = T
(kernel::SoftenedKernel)(x, y) =
    inv(sqrt(sum(abs2, x - y) + kernel.radius^2))

points = [
    SVector(cos(t), sin(t), 0.2sin(2t))
    for t in range(0, 2pi; length=129)[1:(end - 1)]
]
builder = KMeansTreeBuilder(; numberofclusters=2, minvalues=16)
cluster = KMeansTree(points; builder=builder)
tree = BlockTree(cluster, cluster)
kernel = SoftenedKernel(0.05)

hmat = HMatrix(
    kernel,
    points,
    points,
    tree;
    tol=1e-4,
    maxrank=30,
)

dense = [kernel(x, y) for x in points, y in points]
norm(Matrix(hmat) - dense) / norm(dense)

The vector-space constructor creates an AdaptiveCrossApproximation.PointMatrix, which evaluates only the entries requested during near-field assembly and ACA compression. This makes point kernels useful both as toy problems and as a lightweight interface for kernel methods on point clouds.

Configuring Assembly

Compression and Admissibility

The most important compression controls are:

hmat = AdaptiveCrossApproximation.assemble(
    operator,
    space,
    space;
    tol=1e-4,
    maxrank=40,
    isnear=AdaptiveCrossApproximation.isnear(1.0),
)

tol controls each far-block compressor and maxrank places a hard limit on its rank. Reaching maxrank can indicate that the requested tolerance was not attained. The admissibility parameter controls which cluster pairs enter the far field; for the H2Trees backend, increasing it classifies more interactions as far field.

Pass compressor=ACA(...) to override the operator-specific default. In particular, avoid replacing the robust default for a double-layer operator unless the alternative has been validated for that problem.

Cluster-Tree Backend

For general use, the k-means backend is the recommended choice. Its balanced, geometry-adapted clusters usually provide better block partitions, particularly for large or nonuniformly distributed spaces. For BEAST assembly, load ParallelKMeans together with H2Trees to select it automatically:

using ParallelKMeans

hmat = AdaptiveCrossApproximation.assemble(
    operator,
    space,
    space;
    treekwargs=(; numberofclusters=2, minvalues=100),
    tol=1e-4,
    maxrank=40,
)

treekwargs is forwarded only when the tree is constructed automatically. Larger leaf sizes reduce tree depth and construction overhead but create larger terminal blocks. Smaller leaves expose more compression opportunities at the cost of more blocks and tree traversal. The TwoNTree backend remains a useful lightweight choice for small examples and regularly distributed spaces when ParallelKMeans is not loaded.

Supplying an Existing Tree

When a compatible block tree already exists, pass it to the explicit constructor:

hmat = HMatrix(
    operator,
    testspace,
    trialspace,
    tree;
    tol=1e-4,
    maxrank=40,
    compressor=ACA(; tol=1e-4),
    isnear=AdaptiveCrossApproximation.isnear(1.0),
    spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),
)

This entry point is useful when several operators share the same cluster structure or when the application controls tree construction itself.

Symmetric Assembly on One Space

For a transpose-symmetric operator, use the single-space call:

hmat = AdaptiveCrossApproximation.assemble(operator, space; tol=1e-4)

With an existing block tree, use the corresponding single-space constructor:

hmat = HMatrix(operator, space, tree; tol=1e-4)

Automatic assembly constructs one shared cluster tree. The explicit constructor uses only the supplied tree's test cluster tree for both axes. Compression, admissibility, and quadrature options follow general assembly.

These calls assume that the operator and any custom isnear predicate are transpose-symmetric; symmetry is not checked. For complex matrices the requirement is A == transpose(A), rather than Hermitian symmetry. The two-space calls assemble(operator, space, space) and HMatrix(operator, space, space, tree) continue to use general assembly, even when both arguments are the same object.

Each unordered terminal interaction is assembled or compressed once. The opposite orientation references the same physical block and applies a transpose, without conjugation. Diagonal near blocks are assembled and stored in full, then symmetrized by copying their upper triangles into their lower triangles. Both space-ordering modes are supported. Such matrices support products with vectors; hierarchical arithmetic requires general (full) assembly.

Near- and Far-Field Quadrature

The quadstrat keyword follows the convention of BEAST.assemble:

quadrature = BEAST.defaultquadstrat(operator, space, space)
hmat = AdaptiveCrossApproximation.assemble(
    operator,
    space,
    space;
    quadstrat=quadrature,
    tol=1e-4,
)

The strategy is used directly for dense near-field blocks. For BEAST's default Wilton–Sauter strategy, tofarquadstrat extracts the nonsingular far-field rules because singular quadrature is unnecessary for well-separated interactions. Use nearquadstrat and farquadstrat when the two parts require explicit, different strategies.

Space Ordering and Parallel Assembly

The default PreserveSpaceOrder() leaves the test and trial spaces unchanged and stores their tree indices with the blocks. PermuteSpaceInPlace() instead reorders the spaces to follow the tree layout in general two-space assembly:

hmat = AdaptiveCrossApproximation.assemble(
    operator,
    space,
    space;
    spaceordering=AdaptiveCrossApproximation.PermuteSpaceInPlace(),
)

The latter mutates the supplied spaces and should only be selected when the caller expects this reordered layout. Assembly runs on ntasks parallel tasks (keyword, default Threads.nthreads()). Kernel backends may implement localkernelmatrix to create independent block-local assembly caches while sharing the global kernel matrix between threads.

Working with the Result

HMatrix behaves like a matrix in products with vectors and matrices:

y = hmat * x
mul!(y, hmat, x)
mul!(y, transpose(hmat), xt, α, β)      # y = α*transpose(hmat)*xt + β*y
Y = hmat' * X                           # several right-hand sides at once

A real matrix also multiplies complex vectors. Products run on Threads.nthreads() parallel tasks. For repeated products, see Workspaces.

Converting to a dense matrix is useful for validation of small problems but defeats the storage advantage for large ones:

A = Matrix(hmat)

The near- and far-field contributions can be inspected separately:

near = AdaptiveCrossApproximation.nearmatrix(hmat)
far = AdaptiveCrossApproximation.farmatrix(hmat)
stored_gigabytes = AdaptiveCrossApproximation.storage(hmat)

The assembled object has concrete type HMatrixCore; HMatrix remains its public compatibility name. The numerical payload is available as hmat.near::Vector{Matrix{T}} and hmat.far::Vector{LowRankMatrix{T}}. hmat.csr is the canonical row-cluster ordering of terminal blocks. hmat.csc stores only column-cluster traversal metadata and integer references back into the CSR entries; it does not duplicate the matrices. The retained hmat.tree supplies the cluster hierarchy and index sets, and hmat.tol the assembly tolerance.

Workspaces

An assembled HMatrix holds data only: the cluster tree, the block partition (csr, csc), the dense and low-rank blocks and the tolerance. The struct is immutable, and its tree and partition never change. Operations read the matrix. Only the in-place operations overwrite its blocks, and they keep the tree and the partition: the argument of lu! and recompress!, and the destination of mul!, axpy!, axpby!, rmul! and lmul!. The matrix holds no buffers or other state of an operation, so any number of tasks can read it at the same time.

Everything else an operation needs is in a workspace. Its constructor does the work that depends only on the structure of the matrices:

• it analyses the trees and the block partition;
• it plans the block operations and the order in which they may run;
• it distributes them over ntasks parallel tasks;
• it allocates the buffers.

The operation itself then only computes.

The plain forms hmat * x, hmat \ b, lu(hmat), F \ b, A + B, A * B and recompress!(hmat) build the workspace they need, use it once and discard it. For a single operation, this is all you need. Construct a workspace yourself in the following cases:

• The same operation repeats. An iterative solver multiplies with the same matrix in every iteration, and a simulation may solve with one factorization for many right-hand sides. With a workspace, the planning and the allocation happen once instead of in every call.
• Several matrices share one structure. The matrices of a frequency sweep on one geometry have the same trees and partition and differ only in their ranks. One workspace for the LU or a product serves all of them.
• Several tasks work at the same time. A workspace serves one operation at a time. Give every concurrent task its own workspace; the matrix is shared.
• The number of tasks matters. ntasks is a keyword of the workspace constructors; the plain forms use Threads.nthreads().

There are two kinds of workspaces:

• Bound workspaces. MatVecWorkspace(H) and SolveWorkspace(F) hold their matrix or factorization and are used in place of it. A copy for another task, MatVecWorkspace(A) or SolveWorkspace(S), shares the matrix and the plan and has its own buffers.

  A = AdaptiveCrossApproximation.MatVecWorkspace(hmat; ntasks=8)
  y = A * x
  mul!(y, A, x)                   # no numerical allocations once warm
  Y = A' * X

  Products with hmat itself build a temporary MatVecWorkspace on every call.

• Structural workspaces. The workspaces of the LU, products, sums and recompression are built from the matrices of an operation and passed as its last argument. The operation binds the workspace to its arguments. They need the same cluster trees and block partition as the matrices the workspace was built for; otherwise the operation throws.

All workspace types live in the AdaptiveCrossApproximation namespace. For a frequency sweep, for example:

ws = AdaptiveCrossApproximation.LUWorkspace(H₁)
F₁ = lu!(H₁, ws)
F₂ = lu!(H₂, ws)        # same partition: no new planning

The LUWorkspace holds the elimination schedule, every triangular-solve and Schur-update plan and their buffers. Factorizations do not depend on the workspace after lu! returns.

Hierarchical Arithmetic

Sums, products, recompression and the LU factorization of HMatrix objects use Julia's standard functions. Every result block is truncated to the tolerance tol chosen at assembly; for operands with different tolerances the larger one is used. Truncation is in the spectral norm, relative to each block. The theory section explains the algorithms.

Arithmetic requires full block storage: matrices from symmetric assembly support products with vectors only.

A Complete Example

The following point-kernel example factorizes a matrix and solves a system with it. The regularized single-layer kernel gives a well-conditioned system; in BEM applications, hmat comes from assemble.

using AdaptiveCrossApproximation
using H2Trees
using LinearAlgebra
using ParallelKMeans
using StaticArrays

BLAS.set_num_threads(1)

struct RegularizedKernel{T}
    delta::T
end

Base.eltype(::RegularizedKernel{T}) where {T} = T
(kernel::RegularizedKernel)(x, y) = inv(4pi * (norm(x - y) + kernel.delta))

points = [
    SVector(sinpi(θ) * cospi(φ), sinpi(θ) * sinpi(φ), cospi(θ)) for
    (θ, φ) in zip(range(0, 1; length=2000), range(0, 80; length=2000))
]
cluster = KMeansTree(points; builder=KMeansTreeBuilder(; numberofclusters=2, minvalues=32))
tree = BlockTree(cluster, cluster)
kernel = RegularizedKernel(0.01)

hmat = HMatrix(kernel, points, points, tree; tol=1e-8, maxrank=100)

b = rand(length(points))
x = hmat \ b                                         # hmat is unchanged

dense = [kernel(x, y) for x in points, y in points]
norm(dense * x - b) / norm(b)                        # ≈ 6e-10

Solving Systems

x = hmat \ b            # factorize a copy of hmat and solve

F = lu(hmat)            # keep the factorization for several solves; hmat is unchanged
x = F \ b               # vector or matrix of right-hand sides
ldiv!(F, b)             # in place
x = transpose(F) \ b
x = F' \ b

lu!(hmat)               # factorize in place: hmat holds the factors afterwards

The factorization follows P*A ≈ L*U, where P permutes rows only inside each dense diagonal leaf. lu! overwrites the blocks of the matrix with the factors and warns once. An exactly zero pivot in a leaf raises a LeafLUFailure; with check=false, issuccess(F) is false instead.

Each solve with F allocates its buffers. A SolveWorkspace keeps them and acts like the factorization: S = SolveWorkspace(F), then S \ b or ldiv!(x, S, b). The factorization is read-only during solves: tasks solving concurrently with one factorization each use their own workspace, SolveWorkspace(S).

Preconditioning

With a coarse tolerance, the LU factorization is cheap and usually an effective preconditioner for an iterative solver. For example, with Krylov.jl:

using Krylov

coarse = HMatrix(kernel, points, points, tree; tol=1e-2, maxrank=100)
P = AdaptiveCrossApproximation.SolveWorkspace(lu(coarse))
A = AdaptiveCrossApproximation.MatVecWorkspace(hmat)

x, stats = Krylov.gmres(A, b; M=P, ldiv=true, rtol=1e-8)

For the example above, this reduces the GMRES iterations from 58 to 4. The preconditioner needs the same row and column count as hmat, but not the same partition.

Sums and Products

C = A + B
C = A - 2B
axpy!(α, A, C)          # C = α*A + C
axpby!(α, A, β, C)      # C = α*A + β*C
rmul!(C, α)             # exact scaling

C = A * B
mul!(C, A, B, α, β)     # C = α*A*B + β*C

Sums require the same cluster trees and block partition. Products require matching cluster trees (the column tree of A is the row tree of B), while the block partitions may differ; A * B has the partition of A. In-place operations overwrite only their destination, which must not be a factor of the product. Scaling is exact and does not change the ranks.

Recompression

recompress!(hmat)

ACA does not produce optimal ranks. recompress! truncates every low-rank block to the smallest rank meeting hmat.tol, without changing the partition. In the example above, it reduces the storage by about 10%. It overwrites hmat and warns once; AdaptiveCrossApproximation.nnz counts the stored scalars before and after.

Accuracy and Diagnostics

The tolerance bounds the error of every single truncation, not the error of a product or factorization: the errors of many truncations accumulate. Check the residual of a solve, or the error against a dense reference for small problems, before relying on a tolerance. A factorization reports its truncations:

F.converged          # every truncation met the tolerance
F.rank_limited       # some truncation was capped by a rank limit
F.max_block_error    # largest discarded spectral error of a block
F.flushes            # number of truncations

Threading

Every parallel operation takes the number of tasks as the keyword ntasks, with default Threads.nthreads(); ntasks=1 runs on the calling task. The scheduling itself is chosen internally for each operation:

hmat = HMatrix(operator, testspace, trialspace, tree; ntasks=8)
A = AdaptiveCrossApproximation.MatVecWorkspace(hmat; ntasks=8)
AdaptiveCrossApproximation.RecompressionWorkspace(hmat; ntasks=8)
AdaptiveCrossApproximation.LUWorkspace(hmat; ntasks=8)
AdaptiveCrossApproximation.ProductWorkspace(C, A, B; ntasks=8)
AdaptiveCrossApproximation.AdditionWorkspace(C, A; ntasks=8)

Assembly, matrix-vector products, sums and recompression consist of independent pieces of uneven cost; each task takes the next piece as soon as it is done. The block operations of the LU factorization and of products depend on each other. Their workspace records once which blocks every operation reads and writes; the tasks then take the earliest operation whose inputs are final. The factors and products are identical for every ntasks. The forms without a workspace (hmat * x, lu(hmat), A + B, A * B, recompress!(hmat)) use Threads.nthreads() tasks.

Start Julia with several threads and use a single BLAS thread (BLAS.set_num_threads(1)): the block operations are small, and threaded BLAS calls slow them down.

example/pointcloud.jl constructs and directly solves a Coulomb-kernel H-matrix on a point cloud with 10,000 points.

For block partitioning, admissibility, storage complexity, and the assembly model, see the H-matrix theory.