Content deleted Content added
m Open access bot: add arxiv identifier to citation with #oabot. |
m Moving Category:Matrices to Category:Matrices (mathematics) per Wikipedia:Categories for discussion/Speedy |
||
| (18 intermediate revisions by 10 users not shown) | |||
Line 1:
{{Short description|Approximation method}}
In [[numerical mathematics]], '''hierarchical matrices (H-matrices)'''<ref name="
are used as data-sparse approximations of non-sparse matrices. While a [[sparse matrix]] of dimension <math>n</math> can be represented efficiently in <math>O(n)</math> units of storage by storing only its non-zero entries, a non-sparse matrix would require <math>O(n^2)</math> 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 <math>O(n k\,\log(n))</math> units of storage, where <math>k</math> is a parameter controlling the accuracy of the approximation. In typical applications, e.g., when discretizing integral equations,<ref name="MB08">{{cite book|last=Bebendorf|first=Mario|date=2008|title=Hierarchical matrices: A means to efficiently solve elliptic boundary value problems|publisher=Springer}}</ref><ref name="HAKH00">{{cite journal|last1=Hackbusch|first1=Wolfgang|last2=Khoromskij|first2=Boris N.|date=2000|title=A sparse H-Matrix Arithmetic. Part II: Application to Multi-Dimensional Problems|journal=Computing|volume=64|pages=21–47|doi=10.1007/PL00021408 }}</ref><ref name="MB00">{{cite journal|last=Bebendorf|first=Mario|title=Approximation of boundary element matrices|date=2000|journal=Numer. Math.|volume=86|issue=4|pages=565–589|doi=10.1007/pl00005410|s2cid=206858339 }}</ref><ref name="BERJ03">{{cite journal|last1=Bebendorf|first1=Mario|last2=Rjasanow|first2=Sergej|date=2003|title=Adaptive low-rank approximation of collocation matrices|journal=Computing|volume=70|pages=1–24|doi=10.1007/s00607-002-1469-6|citeseerx=10.1.1.133.182|s2cid=16501661 }}</ref><ref name="BOGR05">{{cite journal|last1=Börm|first1=Steffen|last2=Grasedyck|first2=Lars|date=2005|title=Hybrid cross approximation of integral operators|journal=Numer. Math.|volume=101|issue=2|pages=221–249|doi=10.1007/s00211-005-0618-1|citeseerx=10.1.1.330.8950|s2cid=263882011 }}</ref><ref name="BOCH16">{{cite journal|last1=Börm|first1=Steffen|last2=Christophersen|first2=Sven|date=2016|title=Approximation of integral operators by Green quadrature and nested cross approximation|journal=Numer. Math.|volume=133|issue=3|pages=409–442|doi=10.1007/s00211-015-0757-y|arxiv=1404.2234|s2cid=253745725 }}</ref>
preconditioning the resulting systems of linear equations,<ref name="FAMEPR16">{{cite journal|
or solving elliptic partial differential equations,<ref name="BEHA03">{{cite journal|last1=Bebendorf|first1=Mario|last2=Hackbusch|first2=Wolfgang|date=2003|title=Existence of H-matrix approximants to the inverse FE-matrix of elliptic operators with <math>L^\infty</math>-coefficients|journal=Numer. Math.|volume=95|pages=1–28|doi=10.1007/s00211-002-0445-6|s2cid=263876883 }}</ref><ref name="BO10">{{cite journal|last=Börm|first=Steffen|date=2010|title=Approximation of solution operators of elliptic partial differential equations by H- and H<sup>2</sup>-matrices|journal=Numer. Math.|volume=115|issue=2|pages=165–193|doi=10.1007/s00211-009-0278-7|s2cid=7737211 }}</ref><ref name ="FAMEPR13">{{cite journal|last1=Faustmann|first1=Markus|last2=Melenk|first2=J. Markus|last3=Praetorius|first3=Dirk|date=2015|title=H-matrix approximability of the inverses of FEM matrices|journal=Numer. Math.|volume=131|issue=4|pages=615–642|doi=10.1007/s00211-015-0706-9|arxiv=1308.0499|s2cid=2619823 }}</ref><ref name ="SWX16">{{cite journal|last1=Shen|first1=Jie|last2=Wang|first2=Yingwei|last3=Xia|first3=Jianlin|date=2016|title=Fast structured direct spectral methods for differential equations with variable coefficients|journal= SIAM Journal on Scientific Computing|volume=38|issue=1|pages=A28–A54|doi=10.1137/140986815}}</ref> a rank proportional to <math>\log(1/\epsilon)^\gamma</math> with a small constant <math>\gamma</math> is sufficient to ensure an accuracy of <math>\epsilon</math>. 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 <math>O(n k^\alpha\,\log(n)^\beta)</math> operations, where <math>\alpha,\beta\in\{1,2,3\}.</math><ref name="GRHA02"/>
▲,<ref name="FAMEPR16">{{cite journal|last=Faustmann|first=Markus|last2=Melenk|first2=J. Markus|last3=Praetorius|first3=Dirk|date=2016|title=Existence of H-matrix approximants to the inverses of BEM matrices: The simple-layer operator|journal=Math. Comp.|volume=85|pages=119–152|url=https://dx.doi.org/10.1090/mcom/2990|doi=10.1090/mcom/2990}}</ref>
== Basic idea ==
Line 32 ⟶ 10:
let <math>I,J</math> be index sets, and let <math>G\in{\mathbb R}^{I\times J}</math> denote the matrix we have to approximate.
In many applications (see above), we can find subsets <math>t\subseteq I,s\subseteq J</math> such that <math>G|_{t\times s}</math>
can be approximated by a rank-<math>k</math> matrix. This approximation can be represented in factorized form <math>G|_{t\times s}\approx A B^*</math> with factors
<math>A\in{\mathbb R}^{t\times k},B\in{\mathbb R}^{s\times k}</math>.
While the standard representation of the matrix <math>G|_{t\times s}</math> requires <math>O((\#t)(\#s))</math> units of storage,
the factorized representation requires only <math>O(k(\#t+\#s))</math> units. If <math>k</math> is not too large, the storage requirements are reduced significantly.
In order to approximate the entire matrix <math>G</math>, 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
: <math>{\mathcal G}[u](x) = \int_\Omega \kappa(x,y) u(y) \,dy.</math>
Line 84 ⟶ 56:
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<ref name="MB00"/><ref name="BERJ03"/><ref name="TY00">{{cite journal|last=Tyrtyshnikov|first=Eugene|date=2000|title=Incomplete cross approximation in the mosaic-skeleton method|journal=Computing|volume=64|issue=4|pages=367–380|doi=10.1007/s006070070031|citeseerx=10.1.1.100.6153|s2cid=15850058 }}</ref>▼
▲<ref name="TY00">{{cite journal|last=Tyrtyshnikov|first=Eugene|date=2000|title=Incomplete cross approximation in the mosaic-skeleton method|journal=Computing|volume=64|pages=367–380|doi=10.1007/s006070070031}}</ref>
that use only the entries of the original matrix <math>G</math> to construct a [[low rank approximation|low-rank approximation]].
Line 93 ⟶ 62:
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<ref name="BEHA03"/><ref name="BO10"/><ref name="FAMEPR13"/><ref name="SWX16"/> 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 decomposition]]s and solutions to matrix equations.
The central algorithm is the efficient matrix-matrix multiplication, i.e., the computation of <math>Z = Z + \alpha X Y</math>
for hierarchical matrices <math>X,Y,Z</math> and a scalar factor <math>\alpha</math>.
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 <math>Z</math> in
<math>O(n k^2\,\log(n)^2)</math> operations.
Taking advantage of the block structure, the inverse can be computed by using recursion to compute inverses and [[Schur complement]]s of diagonal blocks and combining both using the matrix-matrix multiplication. In a similar way, the [[LU decomposition]]<ref name="BE07">{{cite journal|last=Bebendorf|first=Mario|date=2007|title=Why finite element discretizations can be factored by triangular hierarchical matrices|journal=SIAM J. Numer. Anal.|volume=45|issue=4|pages=1472–1494|doi=10.1137/060669747}}</ref><ref name="GRKRBO09">{{cite journal|last1=Grasedyck|first1=Lars|last2=Kriemann|first2=Ronald|last3=Le Borne|first3=Sabine|date=2009|title=Domain decomposition based H-LU preconditioning|journal=Numer. Math.|volume=112|issue=4|pages=565–600|doi=10.1007/s00211-009-0218-6|doi-access=free}}</ref>
can be constructed using only recursion and multiplication.
Both operations also require <math>O(n k^2\,\log(n)^2)</math> operations.
Line 123 ⟶ 83:
== H<sup>2</sup>-matrices ==
In order to treat very large problems, the structure of hierarchical matrices can be improved:
H<sup>2</sup>-matrices<ref name="HAKHSA02">{{cite book|last1=Hackbusch|first1=Wolfgang|last2=Khoromskij|first2=Boris N.|last3=Sauter|first3=Stefan|title=Lectures on Applied Mathematics |chapter=On H 2-Matrices |date=2002|pages=9–29|doi=10.1007/978-3-642-59709-1_2|isbn=978-3-642-64094-0}}</ref><ref name="BO10b">{{cite book|last=Börm|first=Steffen|date=2010|title=Efficient Numerical Methods for Non-local Operators: H<sup>2</sup>-Matrix Compression, Algorithms and Analysis|publisher=EMS Tracts in Mathematics|url=http://www.ems-ph.org/books/book.php?proj_nr=125|isbn=9783037190913}}</ref>▼
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 <math>O(n k)</math>.▼
▲<ref name="BO10b">{{cite book|last=Börm|first=Steffen|date=2010|title=Efficient Numerical Methods for Non-local Operators: H<sup>2</sup>-Matrix Compression, Algorithms and Analysis|publisher=EMS Tracts in Mathematics|url=http://www.ems-ph.org/books/book.php?proj_nr=125}}</ref>
▲replace the general low-rank structure of the blocks by a hierarchical representation closely related to the
In the context of boundary integral operators, replacing the fixed rank <math>k</math> by block-dependent ranks leads to approximations that preserve the rate of convergence of the underlying boundary element method at a complexity of <math>O(n).</math><ref name="SA00">{{cite journal|last=Sauter|first=Stefan|date=2000|title=Variable order panel clustering|journal=Computing|volume=64|issue=3|pages=223–261|doi=10.1007/s006070050045|s2cid=36813444 }}</ref><ref name="BOSA05">{{cite journal|last1=Börm|first1=Steffen|last2=Sauter|first2=Stefan|date=2005|title=BEM with linear complexity for the classical boundary integral operators|journal=Math. Comp.|volume=74|issue=251|pages=1139–1177|doi=10.1090/s0025-5718-04-01733-8|doi-access=free}}</ref>
Arithmetic operations like multiplication, inversion, and Cholesky or LR factorization of H<sup>2</sup>-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.<ref name="HARE14">{{cite journal|
▲adaptively optimized cluster bases poses a significant challenge.<ref name="HARE14">{{cite journal|last=Börm|first=Steffen|last2=Reimer|first2=Knut|date=2015|title=Efficient arithmetic operations for rank-structured matrices based on hierarchical low-rank updates|journal=Comp. Vis. Sci.|volume=16|issue=6|pages=247–258|url=http://www.arxiv.org/abs/1402.5056}}</ref>
== Literature ==
Line 147 ⟶ 99:
[http://www.hlib.org HLib] is a C software library implementing the most important algorithms for hierarchical and <math>{\mathcal H}^2</math>-matrices.
[https://www.wr.uni-bayreuth.de/en/software/ahmed AHMED] is a C++ software library that can be downloaded for educational purposes.
[http://www.hlibpro.com HLIBpro] is an implementation of the core hierarchical matrix algorithms for commercial applications.
Line 155 ⟶ 107:
[https://github.com/gchavez2/awesome-hierarchical-matrices awesome-hierarchical-matrices] is a repository containing a list of other H-Matrices implementations.
[https://github.com/JuliaMatrices/HierarchicalMatrices.jl HierarchicalMatrices.jl] is a Julia package implementing hierarchical matrices.
[[Category:Matrices]]▼
▲[[Category:Matrices (mathematics)]]
| |||