Gregory Beylkin

Gregory Beylkin, Kristian Sandberg

We use generalized Gaussian quadratures for exponentials to develop a new ODE solver. Nodes and weights of these quadratures are computed for a given bandlimit $c$ and user selected accuracy $ε$, so that they integrate functions $e^{ibx}$, for all $|b|\le c$, with accuracy $ε$. Nodes of these quadratures do not concentrate excessively near the end points of an interval as those of the standard, polynomial-based Gaussian quadratures. Due to this property, the usual implicit Runge Kutta (IRK) collocation method may be used with a large number of nodes, as long as the method chosen for solving the nonlinear system of equations converges. We show that the resulting ODE solver is symplectic and demonstrate (numerically) that it is A-stable. We use this solver, dubbed Band-limited Collocation (BLC-IRK), in the problem of orbit determination. Since BLC-IRK minimizes the number of nodes needed to obtain the solution, in this problem we achieve speed close to that of explicit multistep methods.

Gregory Beylkin, Martin J. Mohlenkamp, Fernando Pérez

The wavefunction for the multiparticle Schrödinger equation is a function of many variables and satisfies an antisymmetry condition, so it is natural to approximate it as a sum of Slater determinants. Many current methods do so, but they impose additional structural constraints on the determinants, such as orthogonality between orbitals or an excitation pattern. We present a method without any such constraints, by which we hope to obtain much more efficient expansions, and insight into the inherent structure of the wavefunction. We use an integral formulation of the problem, a Green's function iteration, and a fitting procedure based on the computational paradigm of separated representations. The core procedure is the construction and solution of a matrix-integral system derived from antisymmetric inner products involving the potential operators. We show how to construct and solve this system with computational complexity competitive with current methods.

Matthew J Reynolds, Gregory Beylkin, Alireza Doostan

We introduce a new, quadratically convergent algorithm for finding maximum absolute value entries of tensors represented in the canonical format. The computational complexity of the algorithm is linear in the dimension of the tensor. We show how to use this algorithm to find global maxima of non-convex multivariate functions in separated form. We demonstrate the performance of the new algorithms on several examples.

Gregory Beylkin, Lucas Monzon, Xinshuo Yang

We present a new adaptive method for electronic structure calculations based on novel fast algorithms for reduction of multivariate mixtures. In our calculations, spatial orbitals are maintained as Gaussian mixtures whose terms are selected in the process of solving equations. Using a fixed basis leads to the so-called "basis error" since orbitals may not lie entirely within the linear span of the basis. To avoid such an error, multiresolution bases are used in adaptive algorithms so that basis functions are selected from a fixed collection of functions, large enough as to approximate solutions within any user-selected accuracy. Our new method achieves adaptivity without using a multiresolution basis. Instead, as a part of an iteration to solve nonlinear equations, our algorithm selects the "best" subset of linearly independent terms of a Gaussian mixture from a collection that is much larger than any possible basis since the locations and shapes of the Gaussian terms are not fixed in advance. Approximating an orbital within a given accuracy, our algorithm yields significantly fewer terms than methods using multiresolution bases. We demonstrate our approach by solving the Hartree-Fock equations for two diatomic molecules, HeH+ and LiH, matching the accuracy previously obtained using multiwavelet bases.

Gregory Beylkin, Vani Cheruvu, Fernando Pérez

We present a fast, adaptive multiresolution algorithm for applying integral operators with a wide class of radially symmetric kernels in dimensions one, two and three. This algorithm is made efficient by the use of separated representations of the kernel. We discuss operators of the class $(-Δ+μ^{2}I)^{-α}$, where $μ\geq0$ and $0<α<3/2$, and illustrate the algorithm for the Poisson and Schrödinger equations in dimension three. The same algorithm may be used for all operators with radially symmetric kernels approximated as a weighted sum of Gaussians, making it applicable across multiple fields by reusing a single implementation. This fast algorithm provides controllable accuracy at a reasonable cost, comparable to that of the Fast Multipole Method (FMM). It differs from the FMM by the type of approximation used to represent kernels and has an advantage of being easily extendable to higher dimensions.

Gregory Beylkin, Lucas Monzon, Xinshuo Yang

We consider fast deterministic algorithms to identify the "best" linearly independent terms in multivariate mixtures and use them to compute, up to a user-selected accuracy, an equivalent representation with fewer terms. One algorithm employs a pivoted Cholesky decomposition of the Gram matrix constructed from the terms of the mixture to select what we call skeleton terms and the other uses orthogonalization for the same purpose. Importantly, the multivariate mixtures do not have to be a separated representation of a function. Both algorithms require $O(r^2 N + p(d) r N) $ operations, where $N$ is the initial number of terms in the multivariate mixture, $r$ is the number of selected linearly independent terms, and $p(d)$ is the cost of computing the inner product between two terms of a mixture in $d$ variables. For general Gaussian mixtures $p(d) \sim d^3$ since we need to diagonalize a $d\times d$ matrix, whereas for separated representations $p(d) \sim d$. Due to conditioning issues, the resulting accuracy is limited to about one half of the available significant digits for both algorithms. We also describe an alternative algorithm that is capable of achieving higher accuracy but is only applicable in low dimensions or to multivariate mixtures in separated form. We describe a number of initial applications of these algorithms to solve partial differential and integral equations and to address several problems in data science. For data science applications in high dimensions,we consider the kernel density estimation (KDE) approach for constructing a probability density function (PDF) of a cloud of points, a far-field kernel summation method and the construction of equivalent sources for non-oscillatory kernels (used in both, computational physics and data science) and, finally, show how to use the new algorithm to produce seeds for subdividing a cloud of points into groups.

Gregory Beylkin

We present an unconditionally stable algorithm for applying matrix transfer function of a linear time invariant system (LTI) in time domain. The state matrix of an LTI system used for modeling long range dependencies in state space models (SSMs) has eigenvalues close to $1$. The standard recursion defining LTI system becomes unstable if the $m\times m$ state matrix has just one eigenvalue with absolute value even slightly greater than 1. This may occur when approximating a state matrix by a structured matrix to reduce the cost of matrix-vector multiplication from $\mathcal{O}\left(m^{2}\right)$ to $\mathcal{O}\left(m\right)$ or $\mathcal{O}\left(m\log m\right).$ We introduce an unconditionally stable algorithm that uses an approximation of the rational transfer function in the z-domain by a matrix polynomial of degree $2^{N+1}-1$, where $N$ is chosen to achieve any user-selected accuracy. Using a cascade implementation in time domain, applying such transfer function to compute $L$ states requires no more than $2L$ matrix-vector multiplications (whereas the standard recursion requires $L$ matrix-vector multiplications). However, using unconditionally stable algorithm, it is not necessary to assure that an approximate state matrix has all eigenvalues with absolute values strictly less than 1 i.e., within the desired accuracy, the absolute value of some eigenvalues may possibly exceed $1$. Consequently, this algorithm allows one to use a wider variety of structured approximations to reduce the cost of matrix-vector multiplication and we briefly describe several of them to be used for this purpose.

Gregory Beylkin, Lucas Monzon, Ignas Satkauskas

We introduce a new approximate multiresolution analysis (MRA) using a single Gaussian as the scaling function, which we call Gaussian MRA (GMRA). As an initial application, we employ this new tool to accurately and efficiently compute the probability density function (PDF) of the product of independent random variables. In contrast with Monte-Carlo (MC) type methods (the only other universal approach known to address this problem), our method not only achieves accuracies beyond the reach of MC but also produces a PDF expressed as a Gaussian mixture, thus allowing for further efficient computations. We also show that an exact MRA corresponding to our GMRA can be constructed for a matching user-selected accuracy.

Matthew Reynolds, Alireza Doostan, Gregory Beylkin

This paper introduces a randomized variation of the alternating least squares (ALS) algorithm for rank reduction of canonical tensor formats. The aim is to address the potential numerical ill-conditioning of least squares matrices at each ALS iteration. The proposed algorithm, dubbed randomized ALS, mitigates large condition numbers via projections onto random tensors, a technique inspired by well-established randomized projection methods for solving overdetermined least squares problems in a matrix setting. A probabilistic bound on the condition numbers of the randomized ALS matrices is provided, demonstrating reductions relative to their standard counterparts. Additionally, results are provided that guarantee comparable accuracy of the randomized ALS solution at each iteration. The performance of the randomized algorithm is studied with three examples, including manufactured tensors and an elliptic PDE with random inputs. In particular, for the latter, tests illustrate not only improvements in condition numbers, but also improved accuracy of the iterative solver for the PDE solution represented in a canonical tensor format.

Robert J. Harrison, Gregory Beylkin, Florian A. Bischoff et al.

MADNESS (multiresolution adaptive numerical environment for scientific simulation) is a high-level software environment for solving integral and differential equations in many dimensions that uses adaptive and fast harmonic analysis methods with guaranteed precision based on multiresolution analysis and separated representations. Underpinning the numerical capabilities is a powerful petascale parallel programming environment that aims to increase both programmer productivity and code scalability. This paper describes the features and capabilities of MADNESS and briefly discusses some current applications in chemistry and several areas of physics.