Sivaram Ambikasaran

Ritesh Khan, Sivaram Ambikasaran

We present two new algebraic multilevel hierarchical matrix algorithms to perform fast matrix-vector product (MVP) for $N$-body problems in $d$ dimensions, namely efficient $\mathcal{H}^2_{*}$ (fully nested algorithm, i.e., $\mathcal{H}^2$ matrix-like algorithm) and $(\mathcal{H}^2 + \mathcal{H})_{*}$ (semi-nested algorithm, i.e., cross of $\mathcal{H}^2$ and $\mathcal{H}$ matrix-like algorithms). The efficient $\mathcal{H}^2_{*}$ and $(\mathcal{H}^2 + \mathcal{H})_{*}$ hierarchical representations are based on our recently introduced weak admissibility condition in higher dimensions, where the admissible clusters are the far-field and the vertex-sharing clusters. Due to the use of nested form of the bases, the proposed hierarchical matrix algorithms are more efficient than the non-nested algorithms ($\mathcal{H}$ matrix algorithms). We rely on purely algebraic low-rank approximation techniques (e.g., ACA and NCA) and develop both algorithms in a black-box fashion. Another noteworthy contribution of this article is that we perform a comparative study of the proposed algorithms with different algebraic (NCA or ACA-based compression) fast MVP algorithms in $2$D and $3$D. The fast algorithms are tested on various kernel matrices and applied to get fast iterative solutions of a dense linear system arising from the discretized integral equations and radial basis function interpolation. Notably, all the algorithms are developed in a similar fashion in $\texttt{C++}$ and tested within the same environment, allowing for meaningful comparisons. The numerical results demonstrate that the proposed algorithms are competitive to the NCA-based standard $\mathcal{H}^2$ matrix algorithm with respect to the memory and time. The C++ implementation of the proposed algorithms is available at https://github.com/riteshkhan/H2weak/.

Sivaram Ambikasaran, Krithika Narayanaswamy

Using accurate multi-component diffusion treatment in numerical combustion studies remains formidable due to the computational cost associated with solving for diffusion velocities. To obtain the diffusion velocities, for low density gases, one needs to solve the Stefan-Maxwell equations along with the zero diffusion flux criteria, which scales as $\mathcal{O}(N^3)$, when solved exactly. In this article, we propose an accurate, fast, direct and robust algorithm to compute multi-component diffusion velocities. To our knowledge, this is the first provably accurate algorithm (the solution can be obtained up to an arbitrary degree of precision) scaling at a computational complexity of $\mathcal{O}(N)$ in finite precision. The key idea involves leveraging the fact that the matrix of the reciprocal of the binary diffusivities, $V$, is low rank, with its rank being independent of the number of species involved. The low rank representation of matrix $V$ is computed in a fast manner at a computational complexity of $\mathcal{O}(N)$ and the Sherman-Morrison-Woodbury formula is used to solve for the diffusion velocities at a computational complexity of $\mathcal{O}(N)$. Rigorous proofs and numerical benchmarks illustrate the low rank property of the matrix $V$ and scaling of the algorithm.

Sivaram Ambikasaran

This article looks at a simple modification to speedup the conventional Gauss Elimination. The proposed modification speeds up the conventional Gauss Elimination by a factor of nearly 9/7 (in the asymptotic limit).

Sivaram Ambikasaran, Carlos Borges, Lise-Marie Imbert-Gerard et al.

We present a fast direct solver for two dimensional scattering problems, where an incident wave impinges on a penetrable medium with compact support. We represent the scattered field using a volume potential whose kernel is the outgoing Green's function for the exterior domain. Inserting this representation into the governing partial differential equation, we obtain an integral equation of the Lippmann-Schwinger type. The principal contribution here is the development of an automatically adaptive, high-order accurate discretization based on a quad tree data structure which provides rapid access to arbitrary elements of the discretized system matrix. This permits the straightforward application of state-of-the-art algorithms for constructing compressed versions of the solution operator. These solvers typically require $O(N^{3/2})$ work, where $N$ denotes the number of degrees of freedom. We demonstrate the performance of the method for a variety of problems in both the low and high frequency regimes.

Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard et al.

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the $n$-dimensional setting, however, it requires the inversion of an $n \times n$ covariance matrix, $C$, as well as the evaluation of its determinant, $\det(C)$. In many cases, such as regression using Gaussian processes, the covariance matrix is of the form $C = σ^2 I + K$, where $K$ is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix $C$ is typically dense, causing standard direct methods for inversion and determinant evaluation to require $\mathcal O(n^3)$ work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix $C$ can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an $\mathcal O (n\log^2 n) $ algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant $\det(C)$, permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining $K$. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.

Amirhossein Aminfar, Sivaram Ambikasaran, Eric Darve

This article presents a fast solver for the dense "frontal" matrices that arise from the multifrontal sparse elimination process of 3D elliptic PDEs. The solver relies on the fact that these matrices can be efficiently represented as a hierarchically off-diagonal low-rank (HODLR) matrix. To construct the low-rank approximation of the off-diagonal blocks, we propose a new pseudo-skeleton scheme, the boundary distance low-rank approximation, that picks rows and columns based on the location of their corresponding vertices in the sparse matrix graph. We compare this new low-rank approximation method to the adaptive cross approximation (ACA) algorithm and show that it achieves betters speedup specially for unstructured meshes. Using the HODLR direct solver as a preconditioner (with a low tolerance) to the GMRES iterative scheme, we can reach machine accuracy much faster than a conventional LU solver. Numerical benchmarks are provided for frontal matrices arising from 3D finite element problems corresponding to a wide range of applications.