Fast Direct Solver Research Articles

Fast Direct Solvers With Arbitrary Admissibility Using Generalized Source Integral Equations

Fast Direct Solvers With Arbitrary Admissibility Using Generalized Source Integral Equations

A High-Order Fast Direct Solver for Surface PDEs

A High-Order Fast Direct Solver for Surface PDEs

A fast semi-analytical meshless method in two-dimensions

This paper introduces a fast direct algorithm for the singular boundary method (SBM) in two-dimensional (2D) problems, utilizing the hierarchical off-diagonal low-rank (HODLR) matrix concept as the foundation of the fast direct solver. The HODLR matrix is constructed by hierarchically partitioning the coefficient matrix into blocks using a binary tree, with all off-diagonal blocks at each level being represented as low-rank factors. Furthermore, The Sherman-Morrison-Woodbury formula can be used to efficiently compute the inverse of the HODLR matrix. The numerical experiment results demonstrate that the new fast solver can significantly improve the computational efficiency of the SBM while maintaining the same level of precision.

A Generalized Source Integral Equation for Enhanced Compression in Three Dimensions

An integral equation formulation, using generalized directional sources, for three-dimensional scattering by impenetrable and essentially convex bodies, is presented. This formulation increases the rank-deficiency of moment matrix blocks representing interactions between non-overlapping source and observer regions. Thereby, it enables enhanced low-rank approximation based matrix compression and the development of corresponding fast direct solvers. The directional sources are constructed by augmenting the conventional basis functions with spherical absorbing three-dimensional shields, on which auxiliary source distributions are defined. The bottleneck arising from the need to integrate over these distributions upon computing the modified Green’s function is removed by using efficient non-uniform sampling and tabulation of the modified Green’s function or its components, in a region-dependent manner. The formulation is studied and its favorable compressibility is demonstrated, for two fundamental types of compression strategies. The non-uniform sampling approach is employed also to facilitate the rank-revealing analysis of large off-diagonal blocks.

Algebraic Inverse Fast Multipole Method: A fast direct solver that is better than HODLR based fast direct solver

This article presents a fast direct solver, termed Algebraic Inverse Fast Multipole Method (from now on abbreviated as AIFMM), for linear systems arising out of N-body problems. AIFMM relies on the following three main ideas: (i) Certain sub-blocks in the matrix corresponding to N-body problems can be efficiently represented as low-rank matrices; (ii) The low-rank sub-blocks in the above matrix are leveraged to construct an extended sparse linear system; (iii) While solving the extended sparse linear system, certain fill-ins that arise in the elimination phase are represented as low-rank matrices and are “redirected” through other variables maintaining zero fill-in sparsity. The main highlights of this article are the following: (i) Our method is completely algebraic (as opposed to the existing Inverse Fast Multipole Method [1–3], from now on abbreviated as IFMM). We rely on our new Nested Cross Approximation [4] (from now on abbreviated as NNCA) to represent the matrix arising out of N-body problems. (ii) A significant contribution is that the algorithm presented in this article is more efficient than the existing IFMMs. In the existing IFMMs, the fill-ins are compressed and redirected as and when they are created. Whereas in this article, we update the fill-ins first without affecting the computational complexity. We then compress and redirect them only once. (iii) Another noteworthy contribution of this article is that we provide a comparison of AIFMM with Hierarchical Off-Diagonal Low-Rank (from now on abbreviated as HODLR) based fast direct solver and NNCA powered GMRES based fast iterative solver. (iv) Additionally, AIFMM is also demonstrated as a preconditioner.

An explicitly-sparse representation for oscillatory kernels with wave atom-like functions

An explicitly-sparse representation for oscillatory kernels is presented in this work by developing a wave atom based method. Multilevel wave atom-like functions are constructed as the transform of the original nodal basis. Then the system matrix in a new non-standard form is derived with respect to the wave atom basis. The wave atom representation is explicitly-sparse in the sense that the system matrix is sparse and computed explicitly. Its sparsity is further enhanced via a-posteriori compression. Finally its log-linear computational complexity with controllable accuracy is demonstrated with numerical results. This explicitly-sparse representation is expected to lay ground to future work related to fast direct solvers and effective preconditioners for high frequency problems. The algorithm may also be viewed as the generalization of wavelet based methods to high frequency cases, and used as a new wideband fast algorithm for wave problems.

A fast direct solver for surface-based whole-head modeling of transcranial magnetic stimulation

When modeling transcranial magnetic stimulation (TMS) in the brain, a fast and accurate electric field solver can support interactive neuronavigation tasks as well as comprehensive biophysical modeling. We formulate, test, and disseminate a direct (i.e., non-iterative) TMS solver that can accurately determine global TMS fields for any coil type everywhere in a high-resolution MRI-based surface model with ~ 200,000 or more arbitrarily selected observation points within approximately 5 s, with the solution time itself of 3 s. The solver is based on the boundary element fast multipole method (BEM-FMM), which incorporates the latest mathematical advancement in the theory of fast multipole methods—an FMM-based LU decomposition. This decomposition is specific to the head model and needs to be computed only once per subject. Moreover, the solver offers unlimited spatial numerical resolution. Despite the fast execution times, the present direct solution is numerically accurate for the default model resolution. In contrast, the widely used brain modeling software SimNIBS employs a first-order finite element method that necessitates additional mesh refinement, resulting in increased computational cost. However, excellent agreement between the two methods is observed for various practical test cases following mesh refinement, including a biophysical modeling task. The method can be readily applied to a wide range of TMS analyses involving multiple coil positions and orientations, including image-guided neuronavigation. It can even accommodate continuous variations in coil geometry, such as flexible H-type TMS coils. The FMM-LU direct solver is freely available to academic users.

Robust fast direct integral equation solver for three-dimensional doubly periodic scattering problems with a large number of layers

Frequency-domain wave scattering problems that arise in acoustics and electromagnetism can be often described by the Helmholtz equation. This work presents a boundary integral equation method for the Helmholtz equation in 3-D multilayered media with many doubly periodic smooth layers. Compared with conventional quasi-periodic Green's function method, the new method is robust at all scattering parameters. A periodizing scheme is used to decompose the solution into near- and far-field contributions. The near-field contribution uses the free-space Green's function in an integral equation on the interface in the unit cell and its immediate eight neighbors; the far-field contribution uses proxy point sources that enclose the unit cell. A specialized high-order quadrature is developed to discretize the underlying surface integral operators to keep the number of unknowns per layer small. We achieve overall linear computational complexity in the number of layers by reducing the linear system into block tridiagonal form and then solving the system directly via block LU decomposition. The new solver is capable of handling a 100-interface structure with 961.3k unknowns to 10−5 accuracy in less than 2 hours on a desktop workstation.

Effective elastic properties of three-dimensional multiple crack problems with the isogeometric boundary element parallel fast direct solver

We propose an isogeometric boundary element (IGABEM) parallel fast direct solver based on OpenMP to solve three-dimensional multiple crack problems. In this paper, the orderly and randomly distributed cracks are studied to observe the influence of the multiple crack interaction on the effective elastic properties. This paper designs the parallel implementation of the matrix low rank approximation and the matrix assembly for the hierarchical off-diagonal low rank (HODLR) matrix. In the matrix low rank approximation, we propose the m-n tree to improve the accelerated hybrid algorithm to balance the CPU performance and memory usage. The numerical results show that the effective elastic properties can respond to the strong interaction to a certain extent, but the maximum stress intensity factors (SIFs) of multiple cracks still need to be calculated to diagnose the possible crack initiation points. On the other hand, the numerical results also show that the IGABEM parallel fast direct solver can greatly improve the computational efficiency on the premise of ensuring the computational accuracy.

Fast Direct Localization for Millimeter Wave MIMO Systems via Deep ADMM Unfolding

Massive arrays deployed in millimeter-wave systems enable high angular resolution performance, which in turn facilitates sub-meter localization services. Albeit suboptimal, up to now the most popular localization approach has been based on a so-called two-step procedure, where triangulation is applied upon aggregation of the angle-of-arrival (AoA) measurements from the collaborative base stations. This is mainly due to the prohibitive computational cost of the existing direct localization approaches in large-scale systems. To address this issue, we propose a deep unfolding based fast direct localization solver. First, the direct localization is formulated as a joint <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$l_{1}$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$l_{2,1}$ </tex-math></inline-formula> norm sparse recovery problem, which is then solved by using alternating direction method of multipliers (ADMM). Next, we develop a deep ADMM unfolding network (DAUN) to learn the ADMM parameter settings from the training data and a position refinement algorithm is proposed for DAUN. Finally, simulation results showcase the superiority of the proposed DAUN over the baseline solvers in terms of better localization accuracy, faster convergence and significantly lower computational complexity.

Detecting Resonance of Radio-Frequency Cavities Using Fast Direct Integral Equation Solvers and Augmented Bayesian Optimization

Detecting Resonance of Radio-Frequency Cavities Using Fast Direct Integral Equation Solvers and Augmented Bayesian Optimization

An Efficient Neural-Network and Finite-Difference Hybrid Method for Elliptic Interface Problems with Applications

A new and efficient neural-network and finite-difference hybrid method is developed for solving Poisson equation in a regular domain with jump discontinuities on embedded irregular interfaces. Since the solution has low regularity across the interface, when applying finite difference discretization to this problem, an additional treatment accounting for the jump discontinuities must be employed. Here, we aim to elevate such an extra effort to ease our implementation by machine learning methodology. The key idea is to decompose the solution into singular and regular parts. The neural network learning machinery incorporating the given jump conditions finds the singular solution, while the standard five-point Laplacian discretization is used to obtain the regular solution with associated boundary conditions. Regardless of the interface geometry, these two tasks only require supervised learning for function approximation and a fast direct solver for Poisson equation, making the hybrid method easy to implement and efficient. The two- and three-dimensional numerical results show that the present hybrid method preserves second-order accuracy for the solution and its derivatives, and it is comparable with the traditional immersed interface method in the literature. As an application, we solve the Stokes equations with singular forces to demonstrate the robustness of the present method.

Rapid calculation of large-scale acoustic scattering from complex targets by a dual-level fast direct solver

Rapid calculation of large-scale acoustic scattering from complex targets by a dual-level fast direct solver

A fast direct boundary element method for 3D acoustic problems based on hierarchical matrices

The boundary element method (BEM) for acoustic problems is a numerical method based on solving the discretized boundary integral equation (BIE) corresponding to the Helmholtz equation. A fast direct BEM for 3D acoustic problems is proposed in this paper, which is more suitable for broadband acoustic simulation of complex structures, such as in the design and analysis of acoustic metamaterials. The main idea of the fast direct solver is based on the hierarchical off-diagonal low-rank (HODLR) matrix, randomized interpolative decomposition and fast matrix inversion formula. Several numerical examples in solving both interior and exterior acoustic problems are presented in this paper, including radiation and scattering problems with distributed and complex structures. The numerical results show that the same level of accuracy and higher computational efficiency can be achieved by using this fast direct BEM compared with the conventional direct BEM.The second and third authors (Liu and Ye) would like to dedicate this paper to Professor Subrata Mukherjee who had enormous influences on their academic lives and research work on the BEM.

A Modified HODLR Solver Based on Higher Order Basis Functions for Solving Electromagnetic Scattering Problems

In this letter, we present a novel fast direct solver based on a modified hierarchically off-diagonal low-rank (M-HODLR) matrix format with higher order basis functions for analyzing electromagnetic scattering by electrically large objects. The M-HODLR matrix compresses the far-group submatrices by an extended admissibility condition, which enhances the efficiency and accuracy of the original HODLR. To further alleviate the burden of the fast direct solution of the conventional surface integral equation (SIE), the higher order hierarchical vector basis (HOHVB) functions defined over a curved surface are employed. To account for HOHVB, an efficient assembly scheme based on patch coupling is developed to accelerate the adaptive across approximation implementation in M-HODLR. The proposed method can significantly reduce the computational time and memory cost of a conventional SIE solver. Several numerical examples are presented to validate the accuracy and efficiency of the proposed method in solving electrically large objects. Owing to the characteristics of HOHVB, different patch sizes or basis functions with different orders can be chosen according to the local geometry of the object. Therefore, the proposed method is very appealing for large- and multiscale simulations.

Overlapping Domain Decomposition Preconditioner for Integral Equations

The discretization of certain integral equations, e.g., the first-kind Fredholm equation of Laplace's equation, leads to symmetric positive-definite linear systems, where the coefficient matrix is dense and often ill-conditioned. We introduce a new preconditioner based on a novel overlapping domain decomposition that can be combined efficiently with fast direct solvers. Empirically, we observe that the condition number of the preconditioned system is $O(1)$, independent of the problem size. Our domain decomposition is designed so that we can construct approximate factorizations of the subproblems efficiently. In particular, we apply the recursive skeletonization algorithm to subproblems associated with every subdomain. We present numerical results on problem sizes up to $16\,384^2$ in two dimensions (2D) and $256^3$ in three dimensions (3D), which were solved in less than 16 hours and three hours, respectively, on an Intel Xeon Platinum 8280M.

Quantitative bounds on impedance-to-impedance operators with applications to fast direct solvers for PDEs

We prove quantitative norm bounds for a family of operators involving impedance boundary conditions on convex, polygonal domains. A robust numerical construction of Helmholtz scattering solutions in variable media via the Dirichlet-to-Neumann operator involves a decomposition of the domain into a sequence of rectangles of varying scales and constructing impedance-to-impedance boundary operators on each subdomain. Our estimates in particular ensure the invertibility, with quantitative bounds in the frequency, of the merge operators required to reconstruct the original Dirichlet-to-Neumann operator in terms of these impedance-to-impedance operators of the sub-domains. A key step in our proof is to obtain Neumann and Dirichlet boundary trace estimates on solutions of the impedance problem, which are of independent interest. In addition to the variable media setting, we also construct bounds for similar merge operators in the obstacle scattering problem.

An Inhomogeneous Plane-Wave Based Single-Level Fast Direct Solver for the Scattering Analysis of Extremely Large Antenna Arrays

A semianalytic single-level fast direct solver (FDS) is presented for the method-of-moments (MoM) analysis of large planar antenna arrays with up to few thousands of complex elements on current desktop computers. Assuming a partition of the impedance matrix, the proposed scheme factorizes the near-field blocks via a singular value decomposition (SVD)-based low-rank approximation and the intermediate-field blocks via an analytical inhomogeneous plane-wave (IPW) expansion. The near- and intermediate-field bases are then unified to leverage the Woodbury formula for the inversion. The solutions for the excitation of each array ports are obtained at once with <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$O(N^{15/7} \ln ^{3} N)$ </tex-math></inline-formula> computational cost and <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$O(N^{11/7} \ln N)$ </tex-math></inline-formula> memory for array problems with <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$N$ </tex-math></inline-formula> basis functions. The complexities are traced precisely and the error is controlled by means of analytical formulas for the required number of IPWs. The FDS is also combined with a macrobasis function (MBF) approach to deal with antenna geometries containing many subwavelength details. This article includes simulations of regular arrays of about <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$169^{2}~3$ </tex-math></inline-formula> -D dipoles and an irregular array composed of seven stations of the low-frequency square kilometer array (SKA) telescope comprising each 256 log-periodic (LP) antennas. The proposed FDS allows the analysis of the coupling effect both within any one station but also between the stations of the telescope, a size and complexity of problem which has not been possible with existing published techniques.

A fast direct solver for non-intrusive reduced order modeling of vibroacoustic problems

This paper proposes a fast model order reduction method for large-scale three-dimensional vibro-acoustic coupling problems, which proceeds in three perspectives. Firstly, part of the coupling matrices is neglected and another part of the coupling matrices is treated as the partial right-hand side of the Schur complement of the corresponding subsystem in order to implement model reduction quickly. Secondly, the proposed fast method does not require the use of the underlying complex-valued and dense boundary element full-order model for the construction of a global orthonormal basis in the offline phase. Rather, its sparse approximation is exploited, which benefits from the consistent amplitude decay property of the kernel functions. This way it can be regarded as a non-intrusive method. Lastly, in addition to structural loads, acoustic loads, such as the load from incident waves impinging on the flexible structure from the outside, can be handled as well. Application of this solution method to the frequency sweep analyses of coupled structural acoustic systems, i.e. simple sphere and complex geometrical model submerged in an infinite fluid with the number of degrees of freedom ranging from hundreds of thousands up to millions, shows the improved computational efficiency and reduced memory requirement as compared to the brute force approach based on the direct evaluation of original high-dimensional models and also the previous acceleration scheme.

A fast direct solver for integral equations on locally refined boundary discretizations and its application to multiphase flow simulations

In transient simulations of particulate Stokes flow, to accurately capture the interaction between the constituent particles and the confining wall, the discretization of the wall often needs to be locally refined in the region approached by the particles. Consequently, standard fast direct solvers lose their efficiency since the linear system changes at each time step. This manuscript presents a new computational approach that avoids this issue by pre-constructing a fast direct solver for the wall ahead of time, computing a low-rank factorization to capture the changes due to the refinement, and solving the problem on the refined discretization via a Woodbury formula. Numerical results illustrate the efficiency of the solver in accelerating particulate Stokes simulations.