References
Base.Matrix
— MethodMatrix(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.
HMatrices.AbstractCompressor
— Typeabstract type AbstractCompressor
Types used to compress matrices.
HMatrices.AbstractKernelMatrix
— Typeabstract 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
.
HMatrices.AbstractSplitter
— Typeabstract type AbstractSplitter
An AbstractSplitter
is used to split a ClusterTree
. The interface requires the following methods:
should_split(clt,splitter)
: return aBool
determining if theClusterTree
should be further dividedsplit!(clt,splitter)
: perform the splitting of theClusterTree
handling the necessary data sorting.
See GeometricSplitter
for an example of an implementation.
HMatrices.CardinalitySplitter
— Typestruct 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
HMatrices.ClusterTree
— TypeClusterTree(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.
HMatrices.ClusterTree
— Typemutable struct ClusterTree{N,T}
Tree structure used to cluster poitns of type SVector{N,T}
into HyperRectangle
s.
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 ofloc2glob
permutation.children::Vector{ClusterTree{N,T}}
parent::ClusterTree{N,T}
HMatrices.DHMatrix
— Typemutable 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.
HMatrices.DHMatrix
— MethodDHMatrix{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.
HMatrices.GeometricSplitter
— Typestruct GeometricSplitter <: AbstractSplitter
Used to split a ClusterTree
in half along the largest axis. The children boxes are shrank to tighly fit the data.
HMatrices.HMatrix
— Typemutable struct HMatrix{R,T} <: AbstractMatrix{T}
A hierarchial matrix constructed from a rowtree
and coltree
of type R
and holding elements of type T
.
HMatrices.HMatrix
— MethodHMatrix{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.
HMatrices.HyperRectangle
— Typestruct HyperRectangle{N,T}
Axis-aligned hyperrectangle in N
dimensions given by low_corner::SVector{N,T}
and high_corner::SVector{N,T}
.
HMatrices.KernelMatrix
— TypeKernelMatrix{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
HMatrices.MulLinearOp
— Typestruct 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.
HMatrices.PartialACA
— Typestruct 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.
HMatrices.PermutedMatrix
— TypePermutedMatrix{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
.
HMatrices.PrincipalComponentSplitter
— Typestruct PrincipalComponentSplitter <: AbstractSplitter
HMatrices.RemoteHMatrix
— Typestruct RemoteHMatrix{S,T}
A light wrapper for a Future
storing an HMatrix
.
HMatrices.RkMatrix
— Typemutable 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.
HMatrices.StrongAdmissibilityStd
— Typestruct 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)
HMatrices.TSVD
— Typestruct 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.
HMatrices.VectorOfVectors
— Typestruct 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!
HMatrices.WeakAdmissibilityStd
— Typestruct WeakAdmissibilityStd
Two blocks are admissible under this condition if the distance
between them is positive.
Base.split
— Methodsplit(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.
HMatrices._aca_partial
— Function_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|.
HMatrices._aca_partial_pivot
— Method_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.
HMatrices._assemble_cpu!
— Method_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.
HMatrices._assemble_hmat_distributed
— Method_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
.
HMatrices._assemble_threads!
— Method_assemble_threads!(hmat::HMatrix,K,comp)
Like _assemble_cpu!
, but uses threads to assemble the leaves. Note that the threads are spanwned using Threads.@spawn
, which means they are spawned on the same worker as the caller.
HMatrices._build_block_structure!
— Method_build_block_structure!(adm_fun,current_node)
Recursive constructor for HMatrix
block structure. Should not be called directly.
HMatrices._hgemv_recursive!
— Method_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)
.
HMatrices._update_frob_norm
— Method_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.
HMatrices.assemble_hmatrix
— Methodassemble_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.
HMatrices.assemble_hmatrix
— Methodassemble_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)
.
HMatrices.binary_split!
— Methodbinary_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.
HMatrices.center
— Methodcenter(Ω)
Center of the smallest ball containing Ω
.
HMatrices.compress!
— Methodcompress!(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))
.
HMatrices.compress!
— Methodcompress!(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.
HMatrices.compression_ratio
— Methodcompression_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.
HMatrices.container
— Methodcontainer(clt::ClusterTree)
Return the object enclosing all the elements of the clt
.
HMatrices.depth
— Functiondepth(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).
HMatrices.diameter
— Methoddiameter(Ω)
Largest distance between x
and y
for x,y ∈ Ω
.
HMatrices.distance
— Methoddistance(X::ClusterTree, Y::ClusterTree)
Distance between the containers of X
and Y
.
HMatrices.distance
— Methoddistance(Ω1,Ω2)
Minimal Euclidean distance between a point x ∈ Ω1
and y ∈ Ω2
.
HMatrices.elements
— Methodelements(clt::ClusterTree)
Iterable list of the elements inside clt
.
HMatrices.filter_tree
— Functionfilter_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
.
HMatrices.filter_tree!
— Functionfilter_tree!(filter,nodes,tree,[isterminal=true])
Like filter_tree
, but appends results to nodes
.
HMatrices.getblock!
— Methodgetblock!(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.
HMatrices.getcol!
— Methodgetcol!(col, M, j)
Return the j
-th column of M
in col
.
HMatrices.getcol!
— Methodgetcol!(col,M::AbstractMatrix,j)
Fill the entries of col
with column j
of M
.
HMatrices.getcol
— Methodgetcol(M, j)
Return the j
-th column of M
.
HMatrices.glob2loc
— Methodglob2loc(clt::ClusterTree)
The inverse of loc2glob
.
HMatrices.hmul!
— Methodhmul!(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.
HMatrices.index_range
— Methodindex_range(clt::ClusterTree)
Indices of elements in root_elements(clt)
which lie inside clt
.
HMatrices.isclean
— Methodisclean(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.
HMatrices.leaves
— Methodleaves(tree)
Return a vector containing all the leaf nodes of tree
.
HMatrices.loc2glob
— Methodloc2glob(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.
HMatrices.newcol!
— Methodnewcol!(A::VectorOfVectors)
Append a new (unitialized) column to A
, and return a view of it.
HMatrices.nodes
— Methodleaves(tree)
Return a vector containing all the nodes of tree
.
HMatrices.parentnode
— Methodparentnode(clt::ClusterTree)
The node's parent. If t
is a root, then parent(t)==t
.
HMatrices.radius
— Methodradius(Ω)
Half the diameter
.
HMatrices.reset!
— Methodreset!(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.
HMatrices.root_elements
— Methodroot_elements(clt::ClusterTree)
The elements contained in the root of the tree to which clt
belongs.
HMatrices.should_split
— Functionshould_split(clt::ClusterTree, depth, splitter::AbstractSplitter)
Determine whether or not a ClusterTree
should be further divided.
HMatrices.split!
— Functionsplit!(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
.
HMatrices.use_global_index
— Methoduse_global_index()::Bool
Default choice of whether operations will use the global indexing system throughout the package.
HMatrices.use_threads
— Methoduse_threads()::Bool
Default choice of whether threads will be used throughout the package.
LinearAlgebra.cholesky!
— Methodcholesky!(M::HMatrix,comp)
Hierarhical cholesky facotrization of M
, using comp
to generate the compressed blocks during the multiplication routines.
LinearAlgebra.cholesky!
— Methodcholesky!(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.
LinearAlgebra.cholesky
— Methodcholesky(M::HMatrix,args...;kwargs...)
Hierarchical cholesky factorization. See cholesky!
for the available options.
LinearAlgebra.lu!
— Methodlu!(M::HMatrix,comp)
Hierarhical LU facotrization of M
, using comp
to generate the compressed blocks during the multiplication routines.
LinearAlgebra.lu!
— Methodlu!(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
— MethodLinearAlgebra.lu(M::HMatrix,args...;kwargs...)
Hierarchical LU factorization. See lu!
for the available options.
LinearAlgebra.mul!
— Functionmul!(y::AbstractVector,H::HMatrix,x::AbstractVector,a,b[;global_index,threads])
Perform y <-- H*x*a + y*b
in place.