Search

In this article we present an algorithm to efficiently evaluate the exchange matrix in periodic systems when Gaussian basis set with pseudopotentials are used. The usual algorithm for evaluating exchange matrix scales cubically with the system size because one has to perform O(N2) fast Fourier transforms (FFT). Here we introduce an algorithm that retains the cubic scaling but reduces the prefactor significantly by eliminating the need to do FFTs during each exchange build. This is accomplished by representing the products of Gaussian basis function using a linear combination of an auxiliary basis the number of which scales linearly with the size of the system. We store the potential due to these auxiliary functions in memory which allows us to obtain the exchange matrix without the need to do FFT, albeit at the cost of additional memory requirement. Although the basic idea of using auxiliary functions is not new, our algorithm is cheaper due to a combination of three ingredients: (a) we use robust Pseudospectral method that allows us to use a relatively small number of auxiliary basis to obtain high accuracy (b) we use occ-RI exchange which eliminates the need to construct the full exchange matrix and (c) we use the (interpolative separable density fitting) ISDF algorithm to construct these auxiliary basis that are used in the robust pseudospectral method. The resulting algorithm is accurate and we note that the error in the final energy decreases exponentially rapidly with the number of auxiliary functions.

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. Importantly, the multivariate mixtures do not have to be a separated representation of a function. 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. 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 ) ∼ d 3 since we need to diagonalize a d × d matrix, whereas for separated representations p ( d ) ∼ d (there is no need for diagonalization). 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.

We introduce a new functional representation of probability density functions (PDFs) of non-negative random variables via a product of a monomial factor and linear combinations of decaying exponentials with complex exponents. This approximate representation of PDFs is obtained for any finite, user-selected accuracy. Using a fast algorithm involving Hankel matrices, we develop a general numerical method for computing the PDF of the sums, products, or quotients of any number of non-negative random variables yielding the result in the same type of functional representation. We present several examples to demonstrate the accuracy of the approach.

We present a new fast algorithm for computing the Boys function using a nonlinear approximation of the integrand via exponentials. The resulting algorithms evaluate the Boys function with real and complex valued arguments and are competitive with previously developed algorithms for the same purpose., Comment: 9 pages, two figures and two tables

By using Poisson's summation formula, we calculate periodic integrals over Gaussian basis functions by partitioning the lattice summations between the real and reciprocal space, where both sums converge exponentially fast with a large exponent. We demonstrate that the summation can be performed efficiently to calculate two-center Gaussian integrals over various kernels including overlap, kinetic, and Coulomb. The summation in real space is performed using an efficient flavor of the McMurchie-Davidson recurrence relation. The expressions for performing summation in the reciprocal space are also derived and implemented. The algorithm for reciprocal space summation allows us to reuse several terms and leads to a significant improvement in efficiency when highly contracted basis functions with large exponents are used. We find that the resulting algorithm is only between a factor of 5 and 15 slower than that for molecular integrals, indicating the very small number of terms needed in both the real and reciprocal space summations. An outline of the algorithm for calculating three-center Coulomb integrals is also provided.

We report the first fully numerical approach for relativistic quantum chemical calculations applicable to molecules. The approach uses an adaptive basis of multiwavelet functions to solve the full four-component Dirac-Coulomb equation to a user-specified accuracy. The accuracy of the code is demonstrated by comparison with ground state energy calculations of atoms performed in GRASP, and the applicability to molecules is shown via ground state calculations of some simple molecules, including water analogs up to H2Po. In the case of molecules, comparison is made with Gaussian basis set calculations in DIRAC.

We construct high-order derivative operators for smooth functions represented via discontinuous multiwavelet bases. The need for such operators arises in order to avoid artifacts when computing functionals involving high-order derivatives of solutions of integral equations. Previously high-order derivatives had to be formed by repeated application of a first-derivative operator that, while uniquely defined, has a spectral norm that grows quadratically with polynomial order and, hence, greatly amplifies numerical noise (truncation error) in the multiwavelet computation. The new constructions proceed via least-squares projection onto smooth bases and provide substantially improved numerical properties as well as permitting direct construction of high-order derivatives. We employ either b-splines or bandlimited exponentials as the intermediate smooth basis, with the former maintaining the concept of approximation order while the latter preserves the pure imaginary spectrum of the first-derivative operator and provides more direct control over the bandlimit and accuracy of computation. We demonstrate the properties of these new operators via several numerical tests as well as application to a problem in nuclear physics. Keywords: Multiwavelets, Multiresolution, Derivatives, Numerical, Discontinuous, Bandlimited exponentials

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 to approximate solutions within any user-selected accuracy. Our new method achieves adaptivity without using a multiresolution basis. Instead, as 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.

We introduce a new method for functional representation of oscillatory integrals within any user-supplied accuracy. Our approach is based on robust methods for nonlinear approximation of functions via exponentials. The complexity of evaluation of the resulting representations of the oscillatory integrals does not depend or depends only mildly on the size of the parameter responsible for the oscillatory behavior.

We present for the first time real-space, arbitrarily-accurate representations of the operators required for up to second-order Douglas-Kroll-Hess (DKH), a model for constructing quasi-relativistic electronic Hamiltonians. The approach can be extended to other operator-based quasi-relativistic models. The representations are in the form of sums of Gaussian functions with positive coefficients and thus enable efficient and numerically-accurate formulations using conventional Gaussian basis sets or other bases such as multiwavelets. The operators are demonstrated with application to hydrogen-like systems using the relativistic-kinematic and first-order DKH Hamiltonians.

The standard particle-in-cell algorithm suffers from grid heating. There exists a gridless alternative which bypasses the deposition step and calculates each Fourier mode of the charge density directly from the particle positions. We show that a gridless method can be computed efficiently through the use of an Unequally Spaced Fast Fourier Transform (USFFT) algorithm. After a spectral field solve, the forces on the particles are calculated via the inverse USFFT (a rapid solution of an approximate linear system). We provide one and two dimensional implementations of this algorithm with an asymptotic runtime of $O(N_p + N_m^D \log N_m^D)$ for each iteration, identical to the standard PIC algorithm (where $N_p$ is the number of particles and $N_m$ is the number of Fourier modes, and $D$ is the spatial dimensionality of the problem) We demonstrate superior energy conservation and reduced noise, as well as convergence of the energy conservation at small time steps., 17 pages, 12 figures

The standard particle-in-cell algorithm suffers from finite grid effects which break energy conservation, cause numerical dispersion, and create numerical instabilities. We present a gridless alternative, bypassing the deposition step and calculating each Fourier mode of the charge density directly from the particle positions. This can be done efficiently through the use of an Unequally Spaced Fast Fourier Transform (USFFT) algorithm1 2. After a spectral field solve, the forces on the particles are calculated via the inverse USFFT (a rapid solution of an approximate linear system). We provide a 1D implementation of this algorithm with an asymptotic runtime of O(Np + Nm log Nm) for each iteration, identical to the standard PIC algorithm (where Np is the number of particles and Nm is the number of Fourier modes). We demonstrate superior energy conservation and reduced noise, as well as convergence at small time steps.

A fully numerical method for the time-dependent Hartree-Fock and density functional theory (TD-HF/DFT) with the Tamm-Dancoff (TD) approximation is presented in a multiresolution analysis (MRA) approach. From a reformulation with effective use of the density matrix operator, we obtain a general form of the HF/DFT linear response equation in the first quantization formalism. It can be readily rewritten as an integral equation with the bound-state Helmholtz (BSH) kernel for the Green's function. The MRA implementation of the resultant equation permits excited state calculations without virtual orbitals. The integral equation is efficiently and adaptively solved using a numerical multiresolution solver with multiwavelet bases. Our implementation of the TD-HF/DFT methods is applied for calculating the excitation energies of H2, Be, N2, H2O, and C2H4 molecules. The numerical errors of the calculated excitation energies converge in proportion to the residuals of the equation in the molecular orbitals and response functions. The energies of the excited states at a variety of length scales ranging from short-range valence excitations to long-range Rydberg-type ones are consistently accurate. It is shown that the multiresolution calculations yield the correct exponential asymptotic tails for the response functions, whereas those computed with Gaussian basis functions are too diffuse or decay too rapidly. We introduce a simple asymptotic correction to the local spin-density approximation (LSDA) so that in the TDDFT calculations, the excited states are correctly bound.

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.

We introduce a new computationally efficient algorithm for constructing near optimal rational approximations of large (one-dimensional) data sets. In contrast to wavelet-type approximations, these new approximations are effectively shift invariant. We note that the complexity of current algorithms for computing near optimal rational approximations prevents their use for large data sets. In order to obtain a near optimal rational approximation of a large data set, we first construct its B-spline representation. Then, by using a new rational approximation of B-splines, we arrive at a suboptimal rational approximation of the data set. We then use a recently developed fast and accurate reduction algorithm for obtaining a near optimal rational approximation from a suboptimal one. Our approach requires first splitting the data into large segments, which may later be merged together, if needed. We also describe a fast algorithm for evaluating these rational approximations. In particular, this allows us to interpolate the original data to any grid. One of the practical applications of our algorithm is the compression of audio signals. To demonstrate the potential competitiveness of our approach, we construct a near optimal rational approximation of a piano recording.

We review the methods in Beylkin and Monzon (2002) [4] and Xiao et al. (2001) [24] for constructing quadratures for bandlimited exponentials and introduce a new algorithm for the same purpose. As in Beylkin and Monzon (2002) [4] , our approach also yields generalized Gaussian quadratures for exponentials integrated against a non-sign-definite weight function. In addition, we compute quadrature weights via l 2 and l ∞ minimization and compare the corresponding quadrature errors.

We solve viscous Burgersʼ equation using a fast and accurate algorithm—referred to here as the reduction algorithm—for computing near optimal rational approximations. Given a proper rational function with n poles, the reduction algorithm computes (for a desired L ∞ -approximation error) a rational approximation of the same form, but with a (near) optimally small number m ≪ n of poles. Although it is well known that (nonlinear) optimal rational approximations are much more efficient than linear representations of functions via a fixed basis (e.g. wavelets), their use in numerical computations has been limited by a lack of efficient, robust, and accurate algorithms. The reduction algorithm presented here computes reliably (near) optimal rational approximations with high accuracy (e.g., ≈ 10 − 14 ) and a complexity that is essentially linear in the number n of original poles. A key tool is a recently developed algorithm for computing small con-eigenvalues of Cauchy matrices with high relative accuracy, an impossible task for standard algorithms without extended precision. Using the reduction algorithm, we develop a numerical calculus for rational representations of functions. Indeed, while operations such as multiplication and convolution increase the number of poles in the representation, we use the reduction algorithm to maintain an optimally small number of poles. To demonstrate the efficiency, robustness, and accuracy of our approach, we solve Burgersʼ equation with small viscosity ν. It is well known that its solutions exhibit moving transition regions of width O ( ν ) , so that this equation provides a stringent test for adaptive PDE solvers. We show that optimal rational approximations capture the solutions with high accuracy using a small number of poles. In particular, we solve the equation with local accuracy ϵ = 10 − 9 for viscosity as small as ν = 10 − 5 .

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., Comment: 13 pages, 4 figures

We describe a systematic method for designing highly accurate and efficient infinite impulse response (IIR) and finite impulse response (FIR) filters given their specifications. In our approach, we first meet the specifications by constructing an IIR filter with, possibly, a large number of poles. We then construct, for any given accuracy, an optimal IIR version of such filter (with a minimal number of poles). Finally, also for any given accuracy, we convert the IIR filter to an efficient FIR filter cascade (either serial or parallel). Since in this FIR approximation the non-causal part of the IIR filter only introduces an additional delay (as a function of the desired accuracy), our IIR construction does not have to enforce causality. Thus, we obtain a simple method for constructing linear phase filters if the specifications so require. All of these procedures are accomplished via robust, fast algorithms. We provide several illustrative examples of our method.

We develop a multiresolution representation of a class of integral operators satisfying boundary conditions on simple domains in order to construct fast algorithms for their application. We also elucidate some delicate theoretical issues related to the construction of periodic Greenʼs functions for Poissonʼs equation. By applying the method of images to the non-standard form of the free space operator, we obtain lattice sums that converge absolutely on all scales, except possibly on the coarsest scale. On the coarsest scale the lattice sums may be only conditionally convergent and, thus, allow for some freedom in their definition. We use the limit of square partial sums as a definition of the limit and obtain a systematic, simple approach to the construction (in any dimension) of periodized operators with sparse non-standard forms. We illustrate the results on several examples in dimensions one and three: the Hilbert transform, the projector on divergence free functions, the non-oscillatory Helmholtz Greenʼs function and the Poisson operator. Remarkably, the limit of square partial sums yields a periodic Poisson Greenʼs function which is not a convolution. Using a short sum of decaying Gaussians to approximate periodic Greenʼs functions, we arrive at fast algorithms for their application. We further show that the results obtained for operators with periodic boundary conditions extend to operators with Dirichlet, Neumann, or mixed boundary conditions.

The cubed-sphere model provides rapid evaluation of the gravity field for more efficient orbit propagation. This paper characterizes the improved computational efficiency of sequential orbit determination, specifically the extended and unscented Kalman filters, when using this new model instead of the common spherical harmonic model. To use the new gravity model with the extended Kalman filter, capabilities to represent the Jacobian of the gravity acceleration are added to numerically integrate the state transition matrix. Filter tests consider improvements for several simulated satellite scenarios with several combinations of measurements provided for estimation. Since cubed-sphere models of higher degree require only a slight change in computation time, orbit propagation and determination systems may now use this model to improve fidelity without any significant change in cost. Using the cubed-sphere model reduces the computational burden of the orbit determination process, with larger benefits found for high-degree filter models. Differences in the estimated trajectories when using the disparate gravity models remain several orders of magnitude less than the absolute filter error for the cases examined.

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.

A new model, dubbed the MRQSphere, provides a multiresolution representation of the gravity field designed for its estimation. The multiresolution representation uses an approximation via Gaussians of the solution of the Laplace’s equation in the exterior of a sphere. Also, instead of the spherical harmonics, variations in the angular variables are modeled by a set of functions constructed using quadratures for the sphere invariant under the icosahedral group. When combined, these tools specify the spatial resolution of the gravity field as a function of altitude and required accuracy. We define this model, and apply it to representing and estimating the gravity field of the asteroid 433 Eros. We verified that a MRQSphere model derived directly from the true spherical harmonics gravity model satisfies the user defined precision. We also use the MRQSphere model to estimate the gravity field of Eros for a simulated satellite mission, yielding a solution with accuracy only limited by measurement errors and their spatial distribution.