Approximation of integral operators by Green quadrature and nested cross approximation

We present a fast algorithm that constructs a data-sparse approximation of matrices arising in the context of integral equation methods for elliptic partial differential equations. The new algorithm uses Green's representation formula in combination with quadrature to obtain a first approximation of the kernel function and then applies nested cross approximation to obtain a more efficient representation. The resulting $\mathcal{H}^2$-matrix representation requires $\mathcal{O}(n k)$ units of storage for an $n\times n$ matrix, where $k$ depends on the prescribed accuracy

SMASH: Structured matrix approximation by separation and hierarchy

This paper presents an efficient method to perform Structured Matrix Approximation by Separation and Hierarchy (SMASH), when the original dense matrix is associated with a kernel function. Given points in a domain, a tree structure is first constructed based on an adaptive partitioning of the computational domain to facilitate subsequent approximation procedures. In contrast to existing schemes based on either analytic or purely algebraic approximations, SMASH takes advantage of both approaches and greatly improves the efficiency. The algorithm follows a bottom-up traversal of the tree and is able to perform the operations associated with each node on the same level in parallel. A strong rank-revealing factorization is applied to the initial analytic approximation in the separation regime so that a special structure is incorporated into the final nested bases. As a consequence, the storage is significantly reduced on one hand and a hierarchy of the original grid is constructed on the other hand. Due to this hierarchy, nested bases at upper levels can be computed in a similar way as the leaf level operations but on coarser grids. The main advantages of SMASH include its simplicity of implementation, its flexibility to construct various hierarchical rank structures and a low storage cost. Rigorous error analysis and complexity analysis are conducted to show that this scheme is fast and stable. The efficiency and robustness of SMASH are demonstrated through various test problems arising from integral equations, structured matrices, etc

Approximation of boundary element matrices using GPGPUs and nested cross approximation

The efficiency of boundary element methods depends crucially on the time required for setting up the stiffness matrix. The far-field part of the matrix can be approximated by compression schemes like the fast multipole method or $\mathcal{H}$-matrix techniques. The near-field part is typically approximated by special quadrature rules like the Sauter-Schwab technique that can handle the singular integrals appearing in the diagonal and near-diagonal matrix elements. Since computing one element of the matrix requires only a small amount of data but a fairly large number of operations, we propose to use general-purpose graphics processing units (GPGPUs) to handle vectorizable portions of the computation: near-field computations are ideally suited for vectorization and can therefore be handled very well by GPGPUs. Modern far-field compression schemes can be split into a small adaptive portion that exhibits divergent control flows, and should therefore be handled by the CPU, and a vectorizable portion that can again be sent to GPGPUs. We propose a hybrid algorithm that splits the computation into tasks for CPUs and GPGPUs. Our method presented in this article is able to reduce the setup time of boundary integral operators by a significant factor of 19-30 for both the Laplace and the Helmholtz equation in 3D when using two consumer GPGPUs compared to a quad-core CPU

O(N)O(N) Iterative and O(NlogN)O(NlogN) Fast Direct Volume Integral Equation Solvers with a Minimal-Rank H2{\cal H}^2-Representation for Large-Scale 33-D Electrodynamic Analysis

Linear complexity iterative and log-linear complexity direct solvers are developed for the volume integral equation (VIE) based general large-scale electrodynamic analysis. The dense VIE system matrix is first represented by a new cluster-based multilevel low-rank representation. In this representation, all the admissible blocks associated with a single cluster are grouped together and represented by a single low-rank block, whose rank is minimized based on prescribed accuracy. From such an initial representation, an efficient algorithm is developed to generate a minimal-rank ${\cal H}^2$-matrix representation. This representation facilitates faster computation, and ensures the same minimal rank's growth rate with electrical size as evaluated from singular value decomposition. Taking into account the rank's growth with electrical size, we develop linear-complexity ${\cal H}^2$-matrix-based storage and matrix-vector multiplication, and thereby an $O(N)$ iterative VIE solver regardless of electrical size. Moreover, we develop an $O(NlogN)$ matrix inversion, and hence a fast $O(NlogN)$ \emph{direct} VIE solver for large-scale electrodynamic analysis. Both theoretical analysis and numerical simulations of large-scale $1$-, $2$- and $3$-D structures on a single-core CPU, resulting in millions of unknowns, have demonstrated the low complexity and superior performance of the proposed VIE electrodynamic solvers. %The algorithms developed in this work are kernel-independent, and hence applicable to other IE operators as well.Comment: 13 pages, 15 figures. This paper was submitted to IEEE Trans. Microwave Theory Tech in June 201

Efficient arithmetic operations for rank-structured matrices based on hierarchical low-rank updates

Many matrices appearing in numerical methods for partial differential equations and integral equations are rank-structured, i.e., they contain submatrices that can be approximated by matrices of low rank. A relatively general class of rank-structured matrices are $\mathcal{H}^2$-matrices: they can reach the optimal order of complexity, but are still general enough for a large number of practical applications. We consider algorithms for performing algebraic operations with $\mathcal{H}^2$-matrices, i.e., for approximating the matrix product, inverse or factorizations in almost linear complexity. The new approach is based on local low-rank updates that can be performed in linear complexity. These updates can be combined with a recursive procedure to approximate the product of two $\mathcal{H}^2$-matrices, and these products can be used to approximate the matrix inverse and the LR or Cholesky factorization. Numerical experiments indicate that the new method leads to preconditioners that require $\mathcal{O}(n)$ units of storage, can be evaluated in $\mathcal{O}(n)$ operations, and take $\mathcal{O}(n \log n)$ operations to set up

Adaptive compression of large vectors

Numerical algorithms for elliptic partial differential equations frequently employ error estimators and adaptive mesh refinement strategies in order to reduce the computational cost. We can extend these techniques to general vectors by splitting the vectors into a hierarchically organized partition of subsets and using appropriate bases to represent the corresponding parts of the vectors. This leads to the concept of \emph{hierarchical vectors}. A hierarchical vector with $m$ subsets and bases of rank $k$ requires $mk$ units of storage, and typical operations like the evaluation of norms and inner products or linear updates can be carried out in $\mathcal{O}(mk^2)$ operations. Using an auxiliary basis, the product of a hierarchical vector and an $\mathcal{H}^2$-matrix can also be computed in $\mathcal{O}(mk^2)$ operations, and if the result admits an approximation with $\widetilde m$ subsets in the original basis, this approximation can be obtained in $\mathcal{O}((m+\widetilde m)k^2)$ operations. Since it is possible to compute the corresponding approximation error exactly, sophisticated error control strategies can be used to ensure the optimal compression. Possible applications of hierarchical vectors include the approximation of eigenvectors and the solution of time-dependent problems with moving local irregularities

Optimal adaptivity for non-symmetric FEM/BEM coupling

We develop a framework which allows us to prove the essential general quasi-orthogonality for the non-symmetric Johnson-Nedelec finite element/boundary element coupling. General quasi-orthogonality was first proposed in [Axioms of Adaptivity, 2014] as a necessary ingredient of optimality proofs and is the major difficulty on the way to prove rate optimal convergence of adaptive algorithms for many strongly non-symmetric problems. The proof exploits a new connection between the general quasi-orthogonality and LU-factorization of infinite matrices. We then derive that a standard adaptive algorithm for the Johnson-Nedelec coupling converges with optimal rates. The developed techniques are fairly general and can most likely be applied to other problems like Stokes equation

Iterative representing set selection for nested cross approximation

A new fast algebraic method for obtaining an $\mathcal{H}^2$-approximation of a matrix from its entries is presented. The main idea behind the method is based on the nested representation and the maximum-volume principle to select submatrices in low-rank matrices. A special iterative approach for the computation of so-called representing sets is established. The main advantage of the method is that it uses only the hierarchical partitioning of the matrix and does not require special "proxy surfaces" to be selected in advance. The numerical experiments for the electrostatic problem and for the boundary integral operator confirm the effectiveness and robustness of the approach. The complexity is linear in the matrix size and polynomial in the ranks. The algorithm is implemented as an open-source Python package that is available online.Comment: Numer. Linear Algebra Appl. 201

Sparse Grid Discontinuous Galerkin Methods for High-Dimensional Elliptic Equations

This paper constitutes our initial effort in developing sparse grid discontinuous Galerkin (DG) methods for high-dimensional partial differential equations (PDEs). Over the past few decades, DG methods have gained popularity in many applications due to their distinctive features. However, they are often deemed too costly because of the large number of degrees of freedom of the approximation space, which are the main bottleneck for simulations in high dimensions. In this paper, we develop sparse grid DG methods for elliptic equations with the aim of breaking the \emph{curse of dimensionality}. Using a hierarchical basis representation, we construct a sparse finite element approximation space, reducing the degrees of freedom from the standard {$O(h^{-d})$ to $O(h^{-1}|\log_2 h|^{d-1})$} for $d$-dimensional problems, where $h$ is the uniform mesh size in each dimension. Our method, based on the interior penalty (IP) DG framework, can achieve accuracy of $O(h^{k}|\log_2 h|^{d-1})$ in the energy norm, where $k$ is the degree of polynomials used. Error estimates are provided and confirmed by numerical tests in multi-dimensions

Approximation of the high-frequency Helmholtz kernel by nested directional interpolation

We present and analyze an approximation scheme for a class of highly oscillatory kernel functions, taking the 2D and 3D Helmholtz kernels as examples. The scheme is based on polynomial interpolation combined with suitable pre- and postmultiplication by plane waves. It is shown to converge exponentially in the polynomial degree and supports multilevel approximation techniques. Our convergence analysis may be employed to establish exponential convergence of certain classes of fast methods for discretizations of the Helmholtz integral operator that feature polylogarithmic-linear complexity