References

Matrix(H::HMatrix;global_index=true)

Convert H to a Matrix. If global_index=true (the default), the entries are given in the global indexing system (see HMatrix for more information); otherwise the local indexing system induced by the row and columns trees are used.

abstract type AbstractCompressor

Types used to compress matrices.

abstract type AbstractKernelMatrix{T} <: AbstractMatrix{T}

Interface for abstract matrices represented through a kernel function f, target elements X, and source elements Y. The matrix entry i,j is given by f(X[i],Y[j]). Concrete subtypes should implement at least

`Base.getindex(K::AbstractKernelMatrix,i::Int,j::Int)`

If a more efficient implementation of getindex(K,I::UnitRange,I::UnitRange), getindex(K,I::UnitRange,j::Int) and getindex(adjoint(K),I::UnitRange,j::Int) is available (e.g. with SIMD vectorization), implementing such methods can improve the speed of assembling an HMatrix.

abstract type AbstractSplitter

An AbstractSplitter is used to split a ClusterTree. The interface requires the following methods:

• should_split(clt,splitter) : return a Bool determining if the ClusterTree should be further divided
• split!(clt,splitter) : perform the splitting of the ClusterTree handling the necessary data sorting.

See GeometricSplitter for an example of an implementation.

struct CardinalitySplitter <: AbstractSplitter

Used to split a ClusterTree along the largest dimension if length(tree)>nmax. The split is performed so the data is evenly distributed amongst all children.

See also: AbstractSplitter

ClusterTree(elements,splitter;[copy_elements=true, threads=false])

Construct a ClusterTree from the given elements using the splitting strategy encoded in splitter. If copy_elements is set to false, the elements argument are directly stored in the ClusterTree and are permuted during the tree construction.

mutable struct ClusterTree{N,T}

Tree structure used to cluster poitns of type SVector{N,T} into HyperRectangles.

Fields:

• _elements::Vector{SVector{N,T}} : vector containing the sorted elements.
• container::HyperRectangle{N,T} : container for the elements in the current node.
• index_range::UnitRange{Int} : indices of elements contained in the current node.
• loc2glob::Vector{Int} : permutation from the local indexing system to the original (global) indexing system used as input in the construction of the tree.
• glob2loc::Vector{Int} : inverse of loc2glob permutation.
• children::Vector{ClusterTree{N,T}}
• parent::ClusterTree{N,T}

mutable struct DHMatrix{R,T} <: AbstractMatrix{T}

Concrete type representing a hierarchical matrix with data distributed amongst various workers. Its structure is very similar to HMatrix, except that the leaves store a RemoteHMatrix object.

The data on the leaves of a DHMatrix may live on a different worker, so calling fetch on them should be avoided whenever possible.

DHMatrix{T}(rowtree,coltree;partition_strategy=:distribute_columns)

Construct the block structure of a distributed hierarchical matrix covering rowtree and coltree. Returns a DHMatrix with leaves that are empty.

The partition_strategy keyword argument determines how to partition the blocks for distributed computing. Currently, the only available options is distribute_columns, which will partition the columns of the underlying matrix into floor(log2(nw)) parts, where nw is the number of workers available.

struct GeometricSplitter <: AbstractSplitter

Used to split a ClusterTree in half along the largest axis. The children boxes are shrank to tighly fit the data.

mutable struct HMatrix{R,T} <: AbstractMatrix{T}

A hierarchial matrix constructed from a rowtree and coltree of type R and holding elements of type T.

HMatrix{T}(rowtree,coltree,adm)

Construct an empty HMatrix with rowtree and coltree using the admissibility condition adm. This function builds the skeleton for the hierarchical matrix, but does not compute data field in the blocks. See assemble_hmatrix for assembling a hierarhical matrix.

struct HyperRectangle{N,T}

Axis-aligned hyperrectangle in N dimensions given by low_corner::SVector{N,T} and high_corner::SVector{N,T}.

KernelMatrix{Tf,Tx,Ty,T} <:: AbstractKernelMatrix{T}

Generic kernel matrix representing a kernel function acting on two sets of elements. If K is a KernelMatrix, then K[i,j] = f(X[i],Y[j]) where f::Tf=kernel(K), X::Tx=rowelements(K) and Y::Ty=colelements(K).

Examples

X = rand(SVector{2,Float64},2)
Y = rand(SVector{2,Float64},2)
K = KernelMatrix(X,Y) do x,y
    sum(x+y)
end

struct MulLinearOp{R,T} <: AbstractMatrix{T}

Abstract matrix representing the following linear operator:

    L = R + P + a * ∑ᵢ Aᵢ * Bᵢ

where R and P are of type RkMatrix{T}, Aᵢ,Bᵢ are of type HMatrix{R,T} and a is scalar multiplier. Calling compressor(L) produces a low-rank approximation of L, where compressor is an AbstractCompressor.

Note: this structure is used to group the operations required when multiplying hierarchical matrices so that they can later be executed in a way that minimizes recompression of intermediate computations.

struct PartialACA

Adaptive cross approximation algorithm with partial pivoting. This structure can be used to generate an RkMatrix from a matrix-like object M as follows:

using LinearAlgebra
rtol = 1e-6
comp = PartialACA(;rtol)
A = rand(10,2)
B = rand(10,2)
M = A*adjoint(B) # a low-rank matrix
R = comp(M, axes(M)...) # compress the entire matrix `M`
norm(Matrix(R) - M) < rtol*norm(M) # true

# output

true

Because it uses partial pivoting, the linear operator does not have to be evaluated at every i,j. This is usually much faster than TSVD, but due to the pivoting strategy the algorithm may fail in special cases, even when the underlying linear operator is of low rank.

PermutedMatrix{K,T} <: AbstractMatrix{T}

Structured used to reprensent the permutation of a matrix-like object. The original matrix is stored in the data::K field, and the permutations are stored in rowperm and colperm.

struct PrincipalComponentSplitter <: AbstractSplitter

struct RemoteHMatrix{S,T}

A light wrapper for a Future storing an HMatrix.

mutable struct RkMatrix{T}

Representation of a rank r matrix M in outer product format M = A*adjoint(B) where A has size m × r and B has size n × r.

The internal representation stores A and B, but R.Bt or R.At can be used to get the respective adjoints.

struct StrongAdmissibilityStd

Two blocks are admissible under this condition if the minimum of their diameter is smaller than eta times the distance between them, where eta::Float64 is a parameter.

Usage:

adm = StrongAdmissibilityStd(;eta=2.0)
adm(Xnode,Ynode)

struct TSVD

Compression algorithm based on a posteriori truncation of an SVD. This is the optimal approximation in Frobenius norm; however, it also tends to be very expensive and thus should be used mostly for "small" matrices.

struct VectorOfVectors{T}

A simple structure which behaves as a Vector{Vector{T}} but stores the entries in a contiguous data::Vector{T} field. All vectors in the VectorOfVectors are assumed to be of size m, and there are k of them, meaning this structure can be used to represent a m × k matrix.

Similar to a vector-of-vectors, calling A[i] returns a view to the i-th column.

See also: newcol!

struct WeakAdmissibilityStd

Two blocks are admissible under this condition if the distance between them is positive.

split(rec::HyperRectangle,[axis]::Int,[place])

Split a hyperrectangle in two along the axis direction at the position place. Returns a tuple with the two resulting hyperrectangles.

When no place is given, defaults to splitting in the middle of the axis.

When no axis and no place is given, defaults to splitting along the largest axis.

_aca_partial(K,irange,jrange,atol,rmax,rtol,istart=1)

Internal function implementing the adaptive cross-approximation algorithm with partial pivoting. The returned R::RkMatrix provides an approximation to K[irange,jrange] which has either rank is expected to satisfy|M - R| < max(atol,rtol*|M|)`, but this inequality may fail to hold due to the various errors involved in estimating the error and |M|.

_aca_partial_pivot(v,I)

Find in the valid set I the index of the element x ∈ v maximizing its smallest singular value. This is equivalent to minimizing the spectral norm of the inverse of x.

When x is a scalar, this is simply the element with largest absolute value.

This general implementation should work for both scalar as well as tensor-valued kernels; see (https://www.sciencedirect.com/science/article/pii/S0021999117306721)[https://www.sciencedirect.com/science/article/pii/S0021999117306721] for more details.

_assemble_cpu!(hmat::HMatrix,K,comp)

Assemble data on the leaves of hmat. The admissible leaves are compressed using the compressor comp. This function assumes the structure of hmat has already been intialized, and therefore should not be called directly. See HMatrix information on constructors.

_assemble_hmat_distributed(K,rtree,ctree;adm=StrongAdmissibilityStd(),comp=PartialACA();global_index=true,threads=false)

Internal methods called after the DHMatrix structure has been initialized in order to construct the HMatrix on each of the leaves of the DHMatrix.

_assemble_threads!(hmat::HMatrix,K,comp)

Like _assemble_cpu!, but uses threads to assemble the leaves. Note that the threads are spawned using Threads.@spawn, which means they are spawned on the same worker as the caller.

_build_block_structure!(adm_fun,current_node)

Recursive constructor for HMatrix block structure. Should not be called directly.

_hgemv_recursive!(C,A,B,offset)

Internal function used to compute C[I] <-- C[I] + A*B[J] where I = rowrange(A) - offset[1] and J = rowrange(B) - offset[2].

The offset argument is used on the caller side to signal if the original hierarchical matrix had a pivot other than (1,1).

_update_frob_norm(acc,A,B)

Given the Frobenius norm of Rₖ = A[1:end-1]*adjoint(B[1:end-1]) in acc, compute the Frobenius norm of Rₖ₊₁ = A*adjoint(B) efficiently.

assemble_hmatrix(K::AbstractKernelMatrix[; atol, rank, rtol, kwargs...])

Construct an approximation of K as an HMatrix using the partial ACA algorithm for the low rank blocks. The atol, rank, and rtol optional arguments are passed to the PartialACA constructor, and the remaining keyword arguments are forwarded to the main assemble_hmatrix function.

assemble_hmatrix([T,], K, rowtree, coltree;
    adm=StrongAdmissibilityStd(),
    comp=PartialACA(),
    threads=true,
    distributed=false,
    global_index=true)

Main routine for assembling a hierarchical matrix. The argument K represents the matrix to be approximated, rowtree and coltree are tree structure partitioning the row and column indices, respectively, adm can be called on a node of rowtree and a node of coltree to determine if the block is compressible, and comp is a function/functor which can compress admissible blocks.

It is assumed that K supports getindex(K,i,j), and that comp can be called as comp(K,irange::UnitRange,jrange::UnitRange) to produce a compressed version of K[irange,jrange] in the form of an RkMatrix.

The type paramter T is used to specify the type of the entries of the matrix, by default is inferred from K using eltype(K).

binary_split!(cluster::ClusterTree,predicate)

Split a ClusterTree into two, sorting all elements in the process according to predicate. cluster is assigned as parent to each children.

Each point is sorted according to whether f(x) returns true (point sorted on the "left" node) or false (point sorted on the "right" node). At the end a minimal HyperRectangle containing all left/right points is created.

center(Ω)

Center of the smallest ball containing Ω.

compress!(M::RkMatrix,tsvd::TSVD)

Recompress the matrix R using a truncated svd of R. The implementation uses the qr-svd strategy to efficiently compute svd(R) when rank(R) ≪ min(size(R)).

compress!(M::Matrix,tsvd::TSVD)

Recompress the matrix M using a truncated svd and output an RkMatrix. The data in M is invalidated in the process.

compression_ratio(H::HTypes)

The ratio of the uncompressed size of H to its compressed size. A compression_ratio of 10 means it would have taken 10 times more memory to store H as a dense matrix.

container(clt::ClusterTree)

Return the object enclosing all the elements of the clt.

depth(tree,acc=0)

Recursive function to compute the depth of node in a a tree-like structure.

Overload this function if your structure has a more efficient way to compute depth (e.g. if it stores it in a field).

diameter(Ω)

Largest distance between x and y for x,y ∈ Ω.

distance(X::ClusterTree, Y::ClusterTree)

Distance between the containers of X and Y.

distance(Ω1,Ω2)

Minimal Euclidean distance between a point x ∈ Ω1 and y ∈ Ω2.

elements(clt::ClusterTree)

Iterable list of the elements inside clt.

filter_tree(f,tree,isterminal=true)

Return a vector containing all the nodes of tree such that f(node)==true. The argument isterminal can be used to control whether to continue the search on children of nodes for which f(node)==true.

filter_tree!(filter,nodes,tree,[isterminal=true])

Like filter_tree, but appends results to nodes.

getblock!(block,K,irange,jrange)

Fill block with K[i,j] for i ∈ irange, j ∈ jrange, where block is of size length(irange) × length(jrange).

A default implementation exists which relies on getindex(K,i,j), but this method can be overloaded for better performance if e.g. a vectorized way of computing a block is available.

getcol!(col, M, j)

Return the j-th column of M in col.

getcol!(col,M::AbstractMatrix,j)

Fill the entries of col with column j of M.

getcol(M, j)

Return the j-th column of M.

glob2loc(clt::ClusterTree)

The inverse of loc2glob.

hmul!(C::HMatrix,A::HMatrix,B::HMatrix,a,b,compressor)

Similar to mul! : compute C <-- A*B*a + C*b, where A,B,C are hierarchical matrices and compressor is a function/functor used in the intermediate stages of the multiplication to avoid growring the rank of admissible blocks after addition is performed.

index_range(clt::ClusterTree)

Indices of elements in root_elements(clt) which lie inside clt.

isclean(H::HMatrix)

Return true if all leaves of H have data, and if the leaves are the only nodes containing data. This is the normal state of an ℋ-matrix, but during intermediate stages of a computation data may be associated with non-leaf nodes for convenience.

leaves(tree)

Return a vector containing all the leaf nodes of tree.

loc2glob(clt::ClusterTree)

The permutation from the (local) indexing system of the elements of the clt to the (global) indexes used upon the construction of the tree.

newcol!(A::VectorOfVectors)

Append a new (unitialized) column to A, and return a view of it.

leaves(tree)

Return a vector containing all the nodes of tree.

parentnode(clt::ClusterTree)

The node's parent. If t is a root, then parent(t)==t.

radius(Ω)

Half the diameter.

reset!(A::VectorOfVectors)

Set the number of columns of A to zero, and the number of rows to zero, but does not resize! the underlying data vector.

root_elements(clt::ClusterTree)

The elements contained in the root of the tree to which clt belongs.

should_split(clt::ClusterTree, depth, splitter::AbstractSplitter)

Determine whether or not a ClusterTree should be further divided.

split!(clt::ClusterTree,splitter::AbstractSplitter)

Divide clt using the strategy implemented by splitter. This function is reponsible of assigning the children and parent fields, as well as of permuting the data of clt.

use_global_index()::Bool

Default choice of whether operations will use the global indexing system throughout the package.

use_threads()::Bool

Default choice of whether threads will be used throughout the package.

cholesky!(M::HMatrix,comp)

Hierarhical cholesky facotrization of M, using comp to generate the compressed blocks during the multiplication routines.

cholesky!(M::HMatrix;atol=0,rank=typemax(Int),rtol=atol>0 ||
rank<typemax(Int) ? 0 : sqrt(eps(Float64)))

Hierarhical cholesky facotrization of M, using the PartialACA(;atol,rtol;rank) compressor.

cholesky(M::HMatrix,args...;kwargs...)

Hierarchical cholesky factorization. See cholesky! for the available options.

lu!(M::HMatrix,comp)

Hierarhical LU facotrization of M, using comp to generate the compressed blocks during the multiplication routines.

lu!(M::HMatrix;atol=0,rank=typemax(Int),rtol=atol>0 ||
rank<typemax(Int) ? 0 : sqrt(eps(Float64)))

Hierarhical LU facotrization of M, using the PartialACA(;atol,rtol;rank) compressor.

LinearAlgebra.lu(M::HMatrix,args...;kwargs...)

Hierarchical LU factorization. See lu! for the available options.

mul!(y::AbstractVector,H::HMatrix,x::AbstractVector,a,b[;global_index,threads])

Perform y <-- H*x*a + y*b in place.