Hierarchical matrix

In numerical mathematics, hierarchical matrices (H-matrices)
are used as data-sparse approximations of non-sparse matrices.
While a sparse matrix of dimension formula_1 can be represented efficiently in formula_2 units of storage
by storing only its non-zero entries, a non-sparse matrix would require formula_3 units of storage, and using this type
of matrices for large problems would therefore be prohibitively expensive in terms of storage and computing time.
Hierarchical matrices provide an approximation requiring only formula_4 units of storage, where formula_5 is a
parameter controlling the accuracy of the approximation.
In typical applications, e.g., when discretizing integral equations
,
preconditioning the resulting systems of linear equations
or solving elliptic partial differential equations
a rank proportional to formula_6 with a small constant formula_7 is sufficient to ensure an
accuracy of formula_8.
Compared to many other data-sparse representations of non-sparse matrices, hierarchical matrices offer a major advantage:
the results of matrix arithmetic operations like matrix multiplication, factorization or inversion can be approximated
in formula_9 operations, where formula_10
Basic idea.
Hierarchical matrices rely on local low-rank approximations:
let formula_11 be index sets, and let formula_12 denote the matrix we have to approximate.
In many applications (see above), we can find subsets formula_13 such that formula_14
can be approximated by a rank-formula_5 matrix.
This approximation can be represented in factorized form formula_16 with factors
formula_17.
While the standard representation of the matrix formula_14 requires formula_19 units of storage,
the factorized representation requires only formula_20 units.
If formula_5 is not too large, the storage requirements are reduced significantly.
In order to approximate the entire matrix formula_22, it is split into a family of submatrices.
Large submatrices are stored in factorized representation, while small submatrices are stored in standard representation
in order to improve efficiency.
Low-rank matrices are closely related to degenerate expansions used in panel clustering and the fast multipole method
to approximate integral operators.
In this sense, hierarchical matrices can be considered the algebraic counterparts of these techniques.
Application to integral operators.
Hierarchical matrices are successfully used to treat integral equations, e.g., the single and double layer potential operators
appearing in the boundary element method.
A typical operator has the form
The Galerkin method leads to matrix entries of the form
where formula_25 and formula_26 are families of finite element basis functions.
If the kernel function formula_27 is sufficiently smooth, we can approximate it by polynomial interpolation to obtain
where formula_29 is the family of interpolation points and formula_30
is the corresponding family of Lagrange polynomials.
Replacing formula_27 by formula_32 yields an approximation
with the coefficients
If we choose formula_13 and use the same interpolation points for all formula_37, we obtain
formula_38.
Obviously, any other approximation separating the variables formula_39 and formula_40, e.g., the multipole expansion,
would also allow us to split the double integral into two single integrals and thus arrive at a similar factorized low-rank matrix.
Of particular interest are cross approximation techniques
that use only the entries of the original matrix formula_22 to construct a low-rank approximation.
Application to elliptic partial differential equations.
Since the solution operator of an elliptic partial differential equation can be expressed as an integral operator involving
Green's function, it is not surprising that the inverse of the stiffness matrix arising from the finite element method
and spectral method
can be approximated by a hierarchical matrix.
Green's function depends on the shape of the computational domain, therefore it is usually not known.
Nevertheless, approximate arithmetic operations can be employed to compute an approximate inverse without knowing the
function explicitly.
Surprisingly, it is possible to prove that the inverse can be approximated even if
the differential operator involves non-smooth coefficients and Green's function is therefore not smooth.
Arithmetic operations.
The most important innovation of the hierarchical matrix method is the development of efficient algorithms for performing
(approximate) matrix arithmetic operations on non-sparse matrices, e.g., to compute approximate inverses, LU decompositions
and solutions to matrix equations.
The central algorithm is the efficient matrix-matrix multiplication, i.e., the computation of formula_42
for hierarchical matrices formula_43 and a scalar factor formula_44.
The algorithm requires the submatrices of the hierarchical matrices to be organized in a block tree structure and takes
advantage of the properties of factorized low-rank matrices to compute the updated formula_45 in
formula_46 operations.
Taking advantage of the block structure, the inverse can be computed by using recursion to compute inverses and
Schur complements of diagonal blocks and combining both using the matrix-matrix multiplication.
In a similar way, the LU decomposition
can be constructed using only recursion and multiplication.
Both operations also require formula_46 operations.
H2-matrices.
In order to treat very large problems, the structure of hierarchical matrices can be improved:
H2-matrices
replace the general low-rank structure of the blocks by a hierarchical representation closely related to the
fast multipole method in order to reduce the storage complexity to formula_48.
In the context of boundary integral operators, replacing the fixed rank formula_5 by block-dependent ranks
leads to approximations that preserve the rate of convergence of the underlying boundary element method
at a complexity of formula_50
Arithmetic operations like multiplication, inversion, and Cholesky or LR factorization of H2-matrices
can be implemented based on two fundamental operations: the matrix-vector multiplication with submatrices
and the low-rank update of submatrices.
While the matrix-vector multiplication is straightforward, implementing efficient low-rank updates with
adaptively optimized cluster bases poses a significant challenge.
Software.
HLib is a C software library implementing the most important algorithms for hierarchical and formula_51-matrices.
AHMED is a C++ software library that can be downloaded for educational purposes.
HLIBpro is an implementation of the core hierarchical matrix algorithms for commercial applications.
H2Lib is an open source implementation of hierarchical matrix algorithms intended for research and teaching.
awesome-hierarchical-matrices is a repository containing a list of other H-Matrices implementations.
HierarchicalMatrices.jl is a Julia package implementing hierarchical matrices.