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 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.

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

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 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 @e, so that they integrate functions e^i^b^x, for all |b|=

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 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, 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 need to compute small con-eigenvalues and the associated con-eigenvectors of positive-definite Cauchy matrices naturally arises when constructing rational approximations with a (near) optimally small $L^{\infty}$ error. Specifically, given a rational function with $n$ poles in the unit disk, a rational approximation with $m\ll n$ poles in the unit disk may be obtained from the $m$th con-eigenvector of an $n\times n$ Cauchy matrix, where the associated con-eigenvalue $\lambda_{m}>0$ gives the approximation error in the $L^{\infty}$ norm. Unfortunately, standard algorithms do not accurately compute small con-eigenvalues (and the associated con-eigenvectors) and, in particular, yield few or no correct digits for con-eigenvalues smaller than the machine roundoff. We develop a fast and accurate algorithm for computing con-eigenvalues and con-eigenvectors of positive-definite Cauchy matrices, yielding even the tiniest con-eigenvalues with high relative accuracy. The algorithm computes the $m$th con-eigenval...

The cubed-sphere gravitational model is a modification of a base model, e.g., the spherical harmonic model, to allow for the fast evaluation of acceleration. The model consists of concentric spheres, each mapped to the surface of a cube and combined with an appropriate interpolation scheme. The paper presents a brief description of the cubed-sphere model and a comparison of it with the spherical harmonic model. The model was configured to achieve a desired accuracy so that dynamical tests, e.g., evaluation of the integration constant, closely approximate that of the spherical harmonic model. The new model closely approximates the spherical harmonic model, with propagated orbits deviating by a fraction of a millimeter at or above feasible Earth-centered altitudes.

We consider the problem of reconstructing a compactly supported function with singularities either from values of its Fourier transform available only in a bounded interval or from a limited number of its Fourier coefficients. Our results are based on several observations and algorithms in [G. Beylkin, L. Monzón, On approximation of functions by exponential sums, Appl. Comput. Harmon. Anal. 19 (1) (2005) 17–48]. We avoid both the Gibbs phenomenon and the use of windows or filtering by constructing approximations to the available Fourier data via a short sum of decaying exponentials. Using these exponentials, we extrapolate the Fourier data to the whole real line and, on taking the inverse Fourier transform, obtain an efficient rational representation in the spatial domain. An important feature of this rational representation is that the positions of its poles indicate location of singularities of the function. We consider these representations in the absence of noise and discuss the impact of adding white noise to the Fourier data. We also compare our results with those obtained by other techniques. As an example of application, we consider our approach in the context of the kernel polynomial method for estimating density of states (eigenvalues) of Hermitian operators. We briefly consider the related problem of approximation by rational functions and provide numerical examples using our approach.

We construct near-optimal quadratures for the sphere that are invariant under the icosahedral rotation group. These quadratures integrate all ( N +1) 2 linearly independent functions in a rotationally invariant subspace of maximal order and degree N . The nodes of these quadratures are nearly uniformly distributed, and the number of nodes is only marginally more than the optimal ( N +1) 2 /3 nodes. Using these quadratures, we discretize the reproducing kernel on a rotationally invariant subspace to construct an analogue of Lagrange interpolation on the sphere. This representation uses function values at the quadrature nodes. In addition, the representation yields an expansion that uses a single function centred and mostly concentrated at nodes of the quadrature, thus providing a much better localization than spherical harmonic expansions. We show that this representation may be localized even further. We also describe two algorithms of complexity for using these grids and representations. Finally, we note that our approach is also applicable to other discrete rotation groups.

The formal representation of the quasi-periodic Helmholtz Green's function obtained by the method of images is only conditionally convergent and, thus, requires an appropriate summation convention for its evaluation. Instead of using this formal sum, we derive a candidate Green's function as a sum of two rapidly convergent series, one to be applied in the spatial domain and the other in the Fourier domain (as in Ewald's method). We prove that this representation of Green's function satisfies the Helmholtz equation with the quasi-periodic condition and, furthermore, leads to a fast algorithm for its application as an operator.We approximate the spatial series by a short sum of separable functions given by Gaussians in each variable. For the series in the Fourier domain, we exploit the exponential decay of its terms to truncate it. We use fast and accurate algorithms for convolving functions with this approximation of the quasi-periodic Green's function. The resulting method yields a fast solver for the Helmholtz equation with the quasi-periodic boundary condition. The algorithm is adaptive in the spatial domain and its performance does not significantly deteriorate when Green's function is applied to discontinuous functions or potentials with singularities. We also construct Helmholtz Green's functions with Dirichlet, Neumann or mixed boundary conditions on simple domains and use a modification of the fast algorithm for the quasi-periodic Green's function to apply them.The complexity, in dimensiond≥2, of these algorithms is(κd log κ+C(log ϵ−1)d), whereϵis the desired accuracy,κis proportional to the number of wavelengths contained in the computational domain andCis a constant. We illustrate our approach with examples.

For a finite but arbitrary precision, we construct efficient low separation rank representations for the Poisson kernel and for the projector on the divergence free functions in the dimension d = 3 . Our construction requires computing only one-dimensional integrals. We use scaling functions of multiwavelet bases, thus making these representations available for a variety of multiresolution algorithms. Besides having many applications, these two operators serve as examples of weakly singular and singular operators for which our approach is applicable. Our approach provides a practical implementation of separated representations of a class of weakly singular and singular operators in dimensions d ⩾ 2 .

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.

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 that are 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.

We develop a two-dimensional solver for the acoustic wave equation with spatially varying coefficients. In what is a new approach, we use a basis of approximate prolate spheroidal wavefunctions and construct derivative operators that incorporate boundary and interface conditions. Writing the wave equation as a first-order system, we evolve the equation in time using the matrix exponential. Computation of the matrix exponential requires efficient representation of operators in two dimensions and for this purpose we use short sums of one-dimensional operators. We also use a partitioned low-rank representation in one dimension to further speed up the algorithm. We demonstrate that the method significantly reduces numerical dispersion and computational time when compared with a fourth-order finite difference scheme in space and an explicit fourth-order Runge–Kutta solver in time.

Multiresolution techniques in multiwavelet bases, made practical in three and higher dimensions by separated representations, have enabled significant advances in the accuracy and manner of computation of molecular electronic structure. The mathematical and numerical techniques are described in the article by Fann. This paper summarizes the major accomplishments in computational chemistry which represent the first substantial application of most of these new ideas in three and higher dimensions. These include basis set limit computation with linear scaling for Hartree-Fock and Density Functional Theory with a wide variety of functionals including hybrid and asymptotically corrected forms. Current capabilities include energies, analytic derivatives, and excitation energies from linear response theory. Direct solution in 6-D of the two-particle wave equation has also been demonstrated. These methods are written using MADNESS which provides a high level of composition using functions and operators with guarantees of speed and precision.

We describe some recent mathematical results in constructing computational methods that lead to the development of fast and accurate multiresolution numerical methods for solving problems in computational chemistry (the so-called multiresolution quantum chemistry). Using low separation rank representations of functions and operators and representations in multiwavelet bases, we developed a multiscale solution method for integral and differential equations and integral transforms. The Poisson equation and the Schrodinger equation, the projector on the divergence free functions, provide important examples with a wide range of applications in computational chemistry, computational electromagnetic and fluid dynamics. We have implemented these ideas that include adaptive representations of operators and functions in the multiwavelet basis and low separation rank approximation of operators and functions. These methods have been implemented into a software package called Multiresolution Adaptive Numerical Evaluation for Scientific Simulation (MADNESS).

In a previous study [R. J. Harrison et al., J. Chem. Phys. (in press)] we reported an efficient, accurate multiresolution solver for the Kohn–Sham self-consisitent field (KS-SCF) method for general polyatomic molecules. This study presents an efficient numerical algorithm to evalute Hartree–Fock (HF) exchange in the multiresolution SCF method to solve the HF equations. The algorithm employs fast integral convolution with the Poission kernel in the nonstandard form, screening the sparse multiwavelet representation to compute results of the integral operator only where required by the nonlocal exchange operator. Localized molecular obitals are used to attain near linear scaling. Results for atoms and molecules demonstrate reliable precision and speed. Calculations for small water clusters demonstrate a total cost to compute the HF exchange potential for all n occ occpuied MOs scaling as O(n occ 1.5 ).

An efficient and accurate analytic gradient method is presented for Hartree–Fock and density functional calculations using multiresolution analysis in multiwavelet bases. The derivative is efficiently computed as an inner product between compressed forms of the density and the differentiated nuclear potential through the Hellmann–Feynman theorem. A smoothed nuclear potential is directly differentiated, and the smoothing parameter required for a given accuracy is empirically determined from calculations on six homonuclear diatomic molecules. The derivatives of N 2 molecule are shown using multiresolution calculation for various accuracies with comparison to correlation consistent Gaussian-type basis sets. The optimized geometries of several molecules are presented using Hartree–Fock and density functional theory. A highly precise Hartree–Fock optimization for the H 2 O molecule produced six digits for the geometric parameters.

We review some recent results on multiwavelet methods for solving integral and partial differential equations and present an efficient representation of operators using discontinuous multiwavelet bases, including the case for singular integral operators. Numerical calculus using these representations produces fast O(N) methods for multiscale solution of integral equations when combined with low separation rank methods. Using this formulation, we compute the Hilbert transform and solve the Poisson and SchrA¶dinger equations. For a fixed order of multiwavelets and for arbitrary but finite- precision computations, the computational complexity is O(N). The computational structures are similar to fast multipole methods but are more generic in yielding fast O(N) algorithm development.

When an algorithm in dimension one is extended to dimension d , in nearly every case its computational cost is taken to the power d . This fundamental difficulty is the single greatest impediment to solving many important problems and has been dubbed the curse of dimensionality . For numerical analysis in dimension d , we propose to use a representation for vectors and matrices that generalizes separation of variables while allowing controlled accuracy. Basic linear algebra operations can be performed in this representation using one-dimensional operations, thus bypassing the exponential scaling with respect to the dimension. Although not all operators and algorithms may be compatible with this representation, we believe that many of the most important ones are. We prove that the multiparticle Schrödinger operator, as well as the inverse Laplacian, can be represented very efficiently in this form. We give numerical evidence to support the conjecture that eigenfunctions inherit this property by computing the ground-state eigenfunction for a simplified Schrödinger operator with 30 particles. We conjecture and provide numerical evidence that functions of operators inherit this property, in which case numerical operator calculus in higher dimensions becomes feasible.

We introduce new families of Gaussian-type quadratures for weighted integrals of exponential functions and consider their applications to integration and interpolation of bandlimited functions. We use a generalization of a representation theorem due to Caratheodory to derive these quadratures. For each positive measure, the quadratures are parameterized by eigenvalues of the Toeplitz matrix constructed from the trigonometric moments of the measure. For a given accuracy ϵ, selecting an eigenvalue close to ϵ yields an approximate quadrature with that accuracy. To compute its weights and nodes, we present a new fast algorithm. These new quadratures can be used to approximate and integrate bandlimited functions, such as prolate spheroidal wave functions, and essentially bandlimited functions, such as Bessel functions. We also develop, for a given precision, an interpolating basis for bandlimited functions on an interval.

Singular and hypersingular operators are ubiquitous in problems of physics, and their use requires a careful numerical interpretation. Although analytical methods for their regularization have long been known, the classical approach does not provide numerical procedures for constructing or applying the regularized operator. We present a multiresolution definition of regularization for integral operators with convolutional kernels which are homogeneous or associated homogeneous functions. We show that our procedure yields the same operator as the classical definition. Moreover, due to the constructive nature of our definition, we provide concise numerical procedures for the construction and application of singular and hypersingular operators. As an application, we present an algorithm for fast computation of discrete sums and briefly discuss several other examples.

New compactly supported wavelets for which both the scaling and wavelet functions have a high number of vanishing moments are presented. Such wavelets are a generalization of the so-called coiflets and they are useful in applications where interpolation and linear phase are of importance. The new approach is to parameterize coiflets by the first moment of the scaling function. By allowing noninteger values for this parameter, the interpolation and linear phase properties of coiflets are optimized. Besides giving a new definition for coiflets, a new system for the filter coefficients is introduced. This system has a minimal set of defining equations and can be solved with algebraic or numerical methods. Examples are given of the various types of coiflets that can be obtained from such systems. The corresponding filter coefficients are listed and their properties are illustrated.

The multiresolution analysis (MRA) strategy for the reduction of a nonlinear differential equation is a procedure for constructing an equation directly for the coarse scale component of the solution. The MRA homogenization process is a method for building a simpler equation whose solution has the same coarse behavior as the solution to a more complex equation. We present two multiresolution reduction methods for nonlinear differential equations: a numerical procedure and an analytic method. We also discuss one possible appproach to the homogenization method.

In many practical problems coefficients of PDEs are changing across many spatial or temporal scales, whereas we might be interested in the behavior of the solution only on some relatively coarse scale. We approach the problem of capturing the influence of fine scales on the behavior of the solution on a coarse scale using the multiresolution strategy. Considering adjacent scales of a multiresolution analysis, we explicitly eliminate variables associated with the finer scale, which leaves us with a coarse-scale equation. We use the term reduction to designate a recursive application of this procedure over a finite number of scales. We present a multiresolution strategy for reduction of self-adjoint, strictly elliptic operators in one and two dimensions. It is known that the non-standard form for a wide class of operators has fast off-diagonal decay and the rate of decay is controlled by the number of vanishing moments of the wavelet. We prove that the reduction procedure preserves the rate of decay over any finite number of scales and therefore results in sparse matrices for computational purposes. Furthermore, the reduction procedure approximately preserves small eigenvalues of self-adjoint, strictly elliptic operators. We also introduce a modified reduction procedure which preserves the small eigenvalues with greater accuracy than the standard reduction procedure and obtain estimates for the perturbation of those eigenvalues. Finally, we discuss potential extensions of the reduction procedure to parabolic and hyperbolic problems.

This work develops fast and adaptive algorithms for numerically solving nonlinear partial differential equations of the formut= Lu+ Nf(u), where L and N are linear differential operators andf(u) is a nonlinear function. These equations are adaptively solved by projecting the solutionuand the operators L and N into a wavelet basis. Vanishing moments of the basis functions permit a sparse representation of the solution and operators. Using these sparse representations fast and adaptive algorithms that apply operators to functions and evaluate nonlinear functions, are developed for solving evolution equations. For a wavelet representation of the solutionuthat containsNssignificant coefficients, the algorithms update the solution usingO(Ns) operations. The approach is applied to a number of examples and numerical results are given.

We introduce tensor Interpolative Decomposition (tensor ID) for the reduction of the separation rank of Canonical Tensor Decompositions (CTDs). Tensor ID selects, for a user-defined accuracy \epsilon, a near optimal subset of terms of a CTD to represent the remaining terms via a linear combination of the selected terms. Tensor ID can be used as an alternative to or a step of the Alternating Least Squares (ALS) algorithm. In addition, we briefly discuss Q-factorization to reduce the size of components within an ALS iteration. Combined, tensor ID and Q-factorization lead to a new paradigm for the reduction of the separation rank of CTDs. In this context, we also discuss the spectral norm as a computational alternative to the Frobenius norm. We reduce the problem of finding tensor IDs to that of constructing Interpolative Decompositions of certain matrices. These matrices are generated via either randomized projection or randomized sampling of the given tensor. We provide cost estimates and several examples of the new approach to the reduction of separation rank.

We consider implementation of operators via filter banks in the framework of the multiresolution analysis. Our method is particularly efficient for convolution operators. Although our method of applying operators to functions may be used with any wavelet basis with a sufficient number of vanishing moments, we distinguish two particular settings, namely, orthogonal bases and the autocorrelation shell. We apply our method to evaluate the Hilbert transform of signals and derive a fast algorithm capable of achieving any given accuracy. We consider the case where the wavelet is the autocorrelation function of another wavelet associated with an orthonormal basis and where our method provides a fast algorithm for the computation of the modulus and the phase of signals. Moreover, the resulting wavelet may be viewed as being (approximately, but with any given accuracy) in the Hardy spaceH2(open face R).

We consider a simple approach for the fast evaluation of the Fourier transform of functions with singularities based on projecting such functions on a subspace of Multiresolution Analysis. We obtain an explicit approximation of the Fourier Transform of generalized functions and develop a fast algorithm based on its evaluation. In particular, we construct an algorithm for the Unequally Spaced Fast Fourier Transform and test its performance in one and two dimensions. The number of operations required by algorithms of this paper is O ( N · log N + N p · (− log ϵ)) in one dimension and O ( N 2 · log N + N p · (− log ϵ) 2 ) in two dimensions, where ϵ is the precision of computation, N is the number of computed frequencies and N p is the number of nodes. We also address the problem of using approximations of generalized functions for solving partial differential equation with singular coefficients or source terms.

Homogenization may be defined as an analysis in which we construct equations describing coarse-scale behavior of the solution while ignoring fine-scale detail. This work concerns a multiresolution strategy for homogenization of differential equations. As the first step towards a more general treatment of nonlinear ODEs and PDEs, we consider the homogenization via multiresolution analysis (MRA) of systems of linear ODEs with variable coefficients and forcing terms. We develop an efficient numerical approach which generates the coefficients of the homogenized equation. As one of the examples we treat wave propagation in a stratified medium.

Significant advances were made on all objectives of the research program. We have developed fast multiresolution methods for performing electronic structure calculations with emphasis on constructing efficient representations of functions and operators. We extended our approach to problems of scattering in solids, i.e. constructing fast algorithms for computing above the Fermi energy level. Part of the work was done in collaboration with Robert Harrison and George Fann at ORNL. Specific results (in part supported by this grant) are listed here and are described in greater detail. (1) We have implemented a fast algorithm to apply the Green's function for the free space (oscillatory) Helmholtz kernel. The algorithm maintains its speed and accuracy when the kernel is applied to functions with singularities. (2) We have developed a fast algorithm for applying periodic and quasi-periodic, oscillatory Green's functions and those with boundary conditions on simple domains. Importantly, the algorithm maintains its speed and accuracy when applied to functions with singularities. (3) We have developed a fast algorithm for obtaining and applying multiresolution representations of periodic and quasi-periodic Green's functions and Green's functions with boundary conditions on simple domains. (4) We have implemented modifications to improve the speed of adaptive multiresolution algorithms formore » applying operators which are represented via a Gaussian expansion. (5) We have constructed new nearly optimal quadratures for the sphere that are invariant under the icosahedral rotation group. (6) We obtained new results on approximation of functions by exponential sums and/or rational functions, one of the key methods that allows us to construct separated representations for Green's functions. (7) We developed a new fast and accurate reduction algorithm for obtaining optimal approximation of functions by exponential sums and/or their rational representations.« less

International audience; no abstract

This paper describes exact and explicit representations of the differential operators, ${{d^n } / {dx^n }}$, $n = 1,2, \cdots $, in orthonormal bases of compactly supported wavelets as well as the representations of the Hilbert transform and fractional derivatives. The method of computing these representations is directly applicable to multidimensional convolution operators.Also, sparse representations of shift operators in orthonormal bases of compactly supported wavelets are discussed and a fast algorithm requiring $O(N\log N)$ operations for computing the wavelet coefficients of all N circulant shifts of a vector of the length $N = 2^n $ is constructed. As an example of an application of this algorithm, it is shown that the storage requirements of the fast algorithm for applying the standard form of a pseudodifferential operator to a vector (see [G. Beylkin, R. R. Coifman, and V. Rokhlin, Comm. Pure. Appl. Math., 44 (1991), pp. 141–183]) may be reduced from $O(N)$ to $O(\log ^2 N)$ significant entries.

The wavefunction for the multiparticle Schr\"odinger 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., Comment: 30 pages

Interpolating scaling functions give a faithful representation of a localized charge distribution by its values on a grid. For such charge distributions, using a Fast Fourier method, we obtain highly accurate electrostatic potentials for free boundary conditions at the cost of O(N log N) operations, where N is the number of grid points. Thus, with our approach, free boundary conditions are treated as efficiently as the periodic conditions via plane wave methods., Comment: 5 pages, 3 figures

We describe a multiresolution solver for the all-electron local density approximation Kohn-Sham equations for general polyatomic molecules. The resulting solutions are obtained to a user-specified precision and the computational cost of applying all operators scales linearly with the number of parameters. The construction and use of separated forms for operators (here, the Green's functions for the Poisson and bound-state Helmholtz equations) enable practical computation in three and higher dimensions. Initial applications include the alkali-earth atoms down to strontium and the water and benzene molecules.

International audience; Nearly every numerical analysis algorithm has computational complexity that scales exponentially in the underlying physical dimension. The separated representation, introduced previously , allows many operations to be performed with scaling that is formally linear in the dimension. In this paper we further develop this representation by (i) discussing the variety of mechanisms that allow it to be surprisingly efficient; (ii) addressing the issue of conditioning; (iii) presenting algorithms for solving linear systems within this framework; and (iv) demonstrating methods for dealing with antisymmetric functions, as arise in the multiparticle Schrödinger equation in quantum mechanics. Numerical examples are given.