In physics, it is sometimes desirable to compute the so-called \emph{Density Of States} (DOS), also known as the \emph{spectral density}, of a real symmetric matrix $A$. The spectral density can be viewed as a probability density distribution that measures the likelihood of finding eigenvalues near some point on the real line. The most straightforward way to obtain this density is to compute all eigenvalues of $A$. But this approach is generally costly and wasteful, especially for matrices of large dimension. There exists alternative methods that allow us to estimate the spectral density function at much lower cost. The major computational cost of these methods is in multiplying $A$ with a number of vectors, which makes them appealing for large-scale problems where products of the matrix $A$ with arbitrary vectors are relatively inexpensive. This paper defines the problem of estimating the spectral density carefully, and discusses how to measure the accuracy of an approximate spectral density. It then surveys a few known methods for estimating the spectral density, and proposes some new variations of existing methods. All methods are discussed from a numerical linear algebra point of view.
In recent years, deep learning has led to impressive results in many fields. In this paper, we introduce a multi-scale artificial neural network for high-dimensional non-linear maps based on the idea of hierarchical nested bases in the fast multipole method and the $\mathcal{H}^2$-matrices. This approach allows us to efficiently approximate discretized nonlinear maps arising from partial differential equations or integral equations. It also naturally extends our recent work based on the generalization of hierarchical matrices [Fan et al. arXiv:1807.01883] but with a reduced number of parameters. In particular, the number of parameters of the neural network grows linearly with the dimension of the parameter space of the discretized PDE. We demonstrate the properties of the architecture by approximating the solution maps of non-linear Schr{ö}dinger equation, the radiative transfer equation, and the Kohn-Sham map.
Phonon calculations based on first principle electronic structure theory, such as the Kohn-Sham density functional theory, have wide applications in physics, chemistry and material science. The computational cost of first principle phonon calculations typically scales steeply as $\mathcal{O}(N_e^4)$, where $N_e$ is the number of electrons in the system. In this work, we develop a new method to reduce the computational complexity of computing the full dynamical matrix, and hence the phonon spectrum, to $\mathcal{O}(N_e^3)$. The key concept for achieving this is to compress the polarizability operator adaptively with respect to the perturbation of the potential due to the change of the atomic configuration. Such adaptively compressed polarizability operator (ACP) allows accurate computation of the phonon spectrum. The reduction of complexity only weakly depends on the size of the band gap, and our method is applicable to insulators as well as semiconductors with small band gaps. We demonstrate the effectiveness of our method using one-dimensional and two-dimensional model problems.
For a large Hermitian matrix $A\in \mathbb{C}^{N\times N}$, it is often the case that the only affordable operation is matrix-vector multiplication. In such case, randomized method is a powerful way to estimate the spectral density (or density of states) of $A$. However, randomized methods developed so far for estimating spectral densities only extract information from different random vectors independently, and the accuracy is therefore inherently limited to $\mathcal{O}(1/\sqrt{N_{v}})$ where $N_{v}$ is the number of random vectors. In this paper we demonstrate that the "$\mathcal{O}(1/\sqrt{N_{v}})$ barrier" can be overcome by taking advantage of the correlated information of random vectors when properly filtered by polynomials of $A$. Our method uses the fact that the estimation of the spectral density essentially requires the computation of the trace of a series of matrix functions that are numerically low rank. By repeatedly applying $A$ to the same set of random vectors and taking different linear combination of the results, we can sweep through the entire spectrum of $A$ by building such low rank decomposition at different parts of the spectrum. Under some assumptions, we demonstrate that a robust and efficient implementation of such spectrum sweeping method can compute the spectral density accurately with $\mathcal{O}(N^2)$ computational cost and $\mathcal{O}(N)$ memory cost. Numerical results indicate that the new method can significantly outperform existing randomized methods in terms of accuracy. As an application, we demonstrate a way to accurately compute a trace of a smooth matrix function, by carefully balancing the smoothness of the integrand and the regularized density of states using a deconvolution procedure.
We present the first systematic work for deriving a posteriori error estimates for general non-polynomial basis functions in an interior penalty discontinuous Galerkin (DG) formulation for solving eigenvalue problems associated with second order linear operators. Eigenvalue problems of such types play important roles in scientific and engineering applications, particularly in theoretical chemistry, solid state physics and material science. Based on the framework developed in [{\it L. Lin, B. Stamm, http://dx.doi.org/10.1051/m2an/2015069}] for second order PDEs, we develop residual type upper and lower bound error estimates for measuring the a posteriori error for eigenvalue problems. The main merit of our method is that the method is parameter-free, in the sense that all but one solution-dependent constants appearing in the upper and lower bound estimates are explicitly computable by solving local and independent eigenvalue problems, and the only non-computable constant can be reasonably approximated by a computable one without affecting the overall effectiveness of the estimates in practice. Compared to the PDE case, we find that a posteriori error estimators for eigenvalue problems must neglect certain terms, which involves explicitly the exact eigenvalues or eigenfunctions that are not accessible in numerical simulations. We define such terms carefully, and justify numerically that the neglected terms are indeed numerically high order terms compared to the computable estimators. Numerical results for a variety of problems in 1D and 2D demonstrate that both the upper bound and lower bound are effective for measuring the error of eigenvalues and eigenfunctions.
Given a sparse Hermitian matrix $A$ and a real number $μ$, we construct a set of sparse vectors, each approximately spanned only by eigenvectors of $A$ corresponding to eigenvalues near $μ$. This set of vectors spans the column space of a localized spectrum slicing (LSS) operator, and is called an LSS basis set. The sparsity of the LSS basis set is related to the decay properties of matrix Gaussian functions. We present a divide-and-conquer strategy with controllable error to construct the LSS basis set. This is a purely algebraic process using only submatrices of $A$, and can therefore be applied to general sparse Hermitian matrices. The LSS basis set leads to sparse projected matrices with reduced sizes, which allows the projected problems to be solved efficiently with techniques using sparse linear algebra. As an example, we demonstrate that the LSS basis set can be used to solve interior eigenvalue problems for a discretized second order partial differential operator in one-dimensional and two-dimensional domains, as well as for a matrix of general sparsity pattern.
In this paper, we propose a new Green's function embedding method called PEXSI-$Σ$ for describing complex systems within the Kohn-Sham density functional theory (KSDFT) framework, after revisiting the physics literature of Green's function embedding methods from a numerical linear algebra perspective. The PEXSI-$Σ$ method approximates the density matrix using a set of nearly optimally chosen Green's functions evaluated at complex frequencies. For each Green's function, the complex boundary conditions are described by a self energy matrix $Σ$ constructed from a physical reference Green's function, which can be computed relatively easily. In the linear regime, such treatment of the boundary condition can be numerically exact. The support of the $Σ$ matrix is restricted to degrees of freedom near the boundary of computational domain, and can be interpreted as a frequency dependent surface potential. This makes it possible to perform KSDFT calculations with $\mathcal{O}(N^2)$ computational complexity, where $N$ is the number of atoms within the computational domain. Green's function embedding methods are also naturally compatible with atomistic Green's function methods for relaxing the atomic configuration outside the computational domain. As a proof of concept, we demonstrate the accuracy of the PEXSI-$Σ$ method for graphene with divacancy and dislocation dipole type of defects using the DFTB+ software package.
Fine-grained spectral properties of quantum Hamiltonians, including both eigenvalues and their multiplicities, provide useful information for characterizing many-body quantum systems as well as for understanding phenomena such as topological order. Extracting such information with small additive error is $\#\textsf{BQP}$-complete in the worst case. In this work, we introduce QFAMES (Quantum Filtering and Analysis of Multiplicities in Eigenvalue Spectra), a quantum algorithm that efficiently identifies clusters of closely spaced dominant eigenvalues and determines their multiplicities under physically motivated assumptions, which allows us to bypass worst-case complexity barriers. QFAMES also enables the estimation of observable expectation values within targeted energy clusters, providing a powerful tool for studying quantum phase transitions and other physical properties. We validate the effectiveness of QFAMES through numerical demonstrations, including its applications to characterizing quantum phases in the transverse-field Ising model and estimating the ground-state degeneracy of a topologically ordered phase in the two-dimensional toric code model. We also generalize QFAMES to the setting of mixed initial states. Our approach offers rigorous theoretical guarantees and significant advantages over existing subspace-based quantum spectral analysis methods, particularly in terms of the sample complexity and the ability to resolve degeneracies.
Wannier functions provide a localized representation of spectral subspaces of periodic Hamiltonians, and play an important role for interpreting and accelerating Hartree-Fock and Kohn-Sham density functional theory calculations in quantum physics and chemistry. For systems with isolated band structure, the existence of exponentially localized Wannier functions and numerical algorithms for finding them are well studied. In contrast, for systems with entangled band structure, Wannier functions must be generalized to span a subspace larger than the spectral subspace of interest to achieve favorable spatial locality. In this setting, little is known about the theoretical properties of these Wannier functions, and few algorithms can find them robustly. We develop a variational formulation to compute these generalized maximally localized Wannier functions. When paired with an initial guess based on the selected columns of the density matrix (SCDM) method, our method can robustly find Wannier functions for systems with entangled band structure. We formulate the problem as a constrained nonlinear optimization problem, and show how the widely used disentanglement procedure can be interpreted as a splitting method to approximately solve this problem. We demonstrate the performance of our method using real materials including silicon, copper, and aluminum. To examine more precisely the localization properties of Wannier functions, we study the free electron gas in one and two dimensions, where we show that the maximally-localized Wannier functions only decay algebraically. We also explain using a one dimensional example how to modify them to obtain super-algebraic decay.
A fundamental problem in quantum physics is to encode functions that are completely anti-symmetric under permutations of identical particles. The Barron space consists of high-dimensional functions that can be parameterized by infinite neural networks with one hidden layer. By explicitly encoding the anti-symmetric structure, we prove that the anti-symmetric functions which belong to the Barron space can be efficiently approximated with sums of determinants. This yields a factorial improvement in complexity compared to the standard representation in the Barron space and provides a theoretical explanation for the effectiveness of determinant-based architectures in ab-initio quantum chemistry.
Nilin Abrahamsen, Zhiyan Ding, Gil Goldshlager et al.
We provide theoretical convergence bounds for the variational Monte Carlo (VMC) method as applied to optimize neural network wave functions for the electronic structure problem. We study both the energy minimization phase and the supervised pre-training phase that is commonly used prior to energy minimization. For the energy minimization phase, the standard algorithm is scale-invariant by design, and we provide a proof of convergence for this algorithm without modifications. The pre-training stage typically does not feature such scale-invariance. We propose using a scale-invariant loss for the pretraining phase and demonstrate empirically that it leads to faster pre-training.
Due to the ever growing amounts of data leveraged for machine learning and scientific computing, it is increasingly important to develop algorithms that sample only a small portion of the data at a time. In the case of linear least-squares, the randomized block Kaczmarz method (RBK) is an appealing example of such an algorithm, but its convergence is only understood under sampling distributions that require potentially prohibitively expensive preprocessing steps. To address this limitation, we analyze RBK when the data is sampled uniformly, showing that its iterates converge in a Monte Carlo sense to a $\textit{weighted}$ least-squares solution. Unfortunately, for general problems the condition number of the weight matrix and the variance of the iterates can become arbitrarily large. We control these issues by incorporating regularization into the RBK iterations, yielding the regularized algorithm ReBlocK. Numerical experiments including examples arising from natural gradient optimization demonstrate that ReBlocK can outperform both RBK and minibatch stochastic gradient descent for inconsistent problems with rapidly decaying singular values.
Neural network wavefunctions optimized using the variational Monte Carlo method have been shown to produce highly accurate results for the electronic structure of atoms and small molecules, but the high cost of optimizing such wavefunctions prevents their application to larger systems. We propose the Subsampled Projected-Increment Natural Gradient Descent (SPRING) optimizer to reduce this bottleneck. SPRING combines ideas from the recently introduced minimum-step stochastic reconfiguration optimizer (MinSR) and the classical randomized Kaczmarz method for solving linear least-squares problems. We demonstrate that SPRING outperforms both MinSR and the popular Kronecker-Factored Approximate Curvature method (KFAC) across a number of small atoms and molecules, given that the learning rates of all methods are optimally tuned. For example, on the oxygen atom, SPRING attains chemical accuracy after forty thousand training iterations, whereas both MinSR and KFAC fail to do so even after one hundred thousand iterations.
Variational quantum algorithms stand at the forefront of simulations on near-term and future fault-tolerant quantum devices. While most variational quantum algorithms involve only continuous optimization variables, the representational power of the variational ansatz can sometimes be significantly enhanced by adding certain discrete optimization variables, as is exemplified by the generalized quantum approximate optimization algorithm (QAOA). However, the hybrid discrete-continuous optimization problem in the generalized QAOA poses a challenge to the optimization. We propose a new algorithm called MCTS-QAOA, which combines a Monte Carlo tree search method with an improved natural policy gradient solver to optimize the discrete and continuous variables in the quantum circuit, respectively. We find that MCTS-QAOA has excellent noise-resilience properties and outperforms prior algorithms in challenging instances of the generalized QAOA.
The combination of neural networks and quantum Monte Carlo methods has arisen as a path forward for highly accurate electronic structure calculations. Previous proposals have combined equivariant neural network layers with an antisymmetric layer to satisfy the antisymmetry requirements of the electronic wavefunction. However, to date it is unclear if one can represent antisymmetric functions of physical interest, and it is difficult to measure the expressiveness of the antisymmetric layer. This work attempts to address this problem by introducing explicitly antisymmetrized universal neural network layers as a diagnostic tool. We first introduce a generic antisymmetric (GA) layer, which we use to replace the entire antisymmetric layer of the highly accurate ansatz known as the FermiNet. We demonstrate that the resulting FermiNet-GA architecture can yield effectively the exact ground state energy for small systems. We then consider a factorized antisymmetric (FA) layer which more directly generalizes the FermiNet by replacing products of determinants with products of antisymmetrized neural networks. Interestingly, the resulting FermiNet-FA architecture does not outperform the FermiNet. This suggests that the sum of products of antisymmetries is a key limiting aspect of the FermiNet architecture. To explore this further, we investigate a slight modification of the FermiNet called the full determinant mode, which replaces each product of determinants with a single combined determinant. The full single-determinant FermiNet closes a large part of the gap between the standard single-determinant FermiNet and FermiNet-GA. Surprisingly, on the nitrogen molecule at a dissociating bond length of 4.0 Bohr, the full single-determinant FermiNet can significantly outperform the standard 64-determinant FermiNet, yielding an energy within 0.4 kcal/mol of the best available computational benchmark.
Variational quantum eigensolvers have recently received increased attention, as they enable the use of quantum computing devices to find solutions to complex problems, such as the ground energy and ground state of strongly-correlated quantum many-body systems. In many applications, it is the optimization of both continuous and discrete parameters that poses a formidable challenge. Using reinforcement learning (RL), we present a hybrid policy gradient algorithm capable of simultaneously optimizing continuous and discrete degrees of freedom in an uncertainty-resilient way. The hybrid policy is modeled by a deep autoregressive neural network to capture causality. We employ the algorithm to prepare the ground state of the nonintegrable quantum Ising model in a unitary process, parametrized by a generalized quantum approximate optimization ansatz: the RL agent solves the discrete combinatorial problem of constructing the optimal sequences of unitaries out of a predefined set and, at the same time, it optimizes the continuous durations for which these unitaries are applied. We demonstrate the noise-robust features of the agent by considering three sources of uncertainty: classical and quantum measurement noise, and errors in the control unitary durations. Our work exhibits the beneficial synergy between reinforcement learning and quantum control.
The efficient treatment of long-range interactions for point clouds is a challenging problem in many scientific machine learning applications. To extract global information, one usually needs a large window size, a large number of layers, and/or a large number of channels. This can often significantly increase the computational cost. In this work, we present a novel neural network layer that directly incorporates long-range information for a point cloud. This layer, dubbed the long-range convolutional (LRC)-layer, leverages the convolutional theorem coupled with the non-uniform Fourier transform. In a nutshell, the LRC-layer mollifies the point cloud to an adequately sized regular grid, computes its Fourier transform, multiplies the result by a set of trainable Fourier multipliers, computes the inverse Fourier transform, and finally interpolates the result back to the point cloud. The resulting global all-to-all convolution operation can be performed in nearly-linear time asymptotically with respect to the number of input points. The LRC-layer is a particularly powerful tool when combined with local convolution as together they offer efficient and seamless treatment of both short and long range interactions. We showcase this framework by introducing a neural network architecture that combines LRC-layers with short-range convolutional layers to accurately learn the energy and force associated with a $N$-body potential. We also exploit the induced two-level decomposition and propose an efficient strategy to train the combined architecture with a reduced number of samples.
The quantum alternating operator ansatz (QAOA) is a prominent example of variational quantum algorithms. We propose a generalized QAOA called CD-QAOA, which is inspired by the counterdiabatic driving procedure, designed for quantum many-body systems and optimized using a reinforcement learning (RL) approach. The resulting hybrid control algorithm proves versatile in preparing the ground state of quantum-chaotic many-body spin chains by minimizing the energy. We show that using terms occurring in the adiabatic gauge potential as generators of additional control unitaries, it is possible to achieve fast high-fidelity many-body control away from the adiabatic regime. While each unitary retains the conventional QAOA-intrinsic continuous control degree of freedom such as the time duration, we consider the order of the multiple available unitaries appearing in the control sequence as an additional discrete optimization problem. Endowing the policy gradient algorithm with an autoregressive deep learning architecture to capture causality, we train the RL agent to construct optimal sequences of unitaries. The algorithm has no access to the quantum state, and we find that the protocol learned on small systems may generalize to larger systems. By scanning a range of protocol durations, we present numerical evidence for a finite quantum speed limit in the nonintegrable mixed-field spin-1/2 Ising and Lipkin-Meshkov-Glick models, and for the suitability to prepare ground states of the spin-1 Heisenberg chain in the long-range and topologically ordered parameter regimes. This work paves the way to incorporate recent success from deep learning for the purpose of quantum many-body control.
Jiefu Zhang, Leonardo Zepeda-Núñez, Yuan Yao et al.
The task of using machine learning to approximate the mapping $\mathbf{x}\mapsto\sum_{i=1}^d x_i^2$ with $x_i\in[-1,1]$ seems to be a trivial one. Given the knowledge of the separable structure of the function, one can design a sparse network to represent the function very accurately, or even exactly. When such structural information is not available, and we may only use a dense neural network, the optimization procedure to find the sparse network embedded in the dense network is similar to finding the needle in a haystack, using a given number of samples of the function. We demonstrate that the cost (measured by sample complexity) of finding the needle is directly related to the Barron norm of the function. While only a small number of samples is needed to train a sparse network, the dense network trained with the same number of samples exhibits large test loss and a large generalization gap. In order to control the size of the generalization gap, we find that the use of explicit regularization becomes increasingly more important as $d$ increases. The numerically observed sample complexity with explicit regularization scales as $\mathcal{O}(d^{2.5})$, which is in fact better than the theoretically predicted sample complexity that scales as $\mathcal{O}(d^{4})$. Without explicit regularization (also called implicit regularization), the numerically observed sample complexity is significantly higher and is close to $\mathcal{O}(d^{4.5})$.
The quantum approximate optimization algorithm (QAOA), as a hybrid quantum/classical algorithm, has received much interest recently. QAOA can also be viewed as a variational ansatz for quantum control. However, its direct application to emergent quantum technology encounters additional physical constraints: (i) the states of the quantum system are not observable; (ii) obtaining the derivatives of the objective function can be computationally expensive or even inaccessible in experiments, and (iii) the values of the objective function may be sensitive to various sources of uncertainty, as is the case for noisy intermediate-scale quantum (NISQ) devices. Taking such constraints into account, we show that policy-gradient-based reinforcement learning (RL) algorithms are well suited for optimizing the variational parameters of QAOA in a noise-robust fashion, opening up the way for developing RL techniques for continuous quantum control. This is advantageous to help mitigate and monitor the potentially unknown sources of errors in modern quantum simulators. We analyze the performance of the algorithm for quantum state transfer problems in single- and multi-qubit systems, subject to various sources of noise such as error terms in the Hamiltonian, or quantum uncertainty in the measurement process. We show that, in noisy setups, it is capable of outperforming state-of-the-art existing optimization algorithms.
Leonardo Zepeda-Núñez, Yixiao Chen, Jiefu Zhang et al.
The recently developed Deep Potential [Phys. Rev. Lett. 120, 143001, 2018] is a powerful method to represent general inter-atomic potentials using deep neural networks. The success of Deep Potential rests on the proper treatment of locality and symmetry properties of each component of the network. In this paper, we leverage its network structure to effectively represent the mapping from the atomic configuration to the electron density in Kohn-Sham density function theory (KS-DFT). By directly targeting at the self-consistent electron density, we demonstrate that the adapted network architecture, called the Deep Density, can effectively represent the electron density as the linear combination of contributions from many local clusters. The network is constructed to satisfy the translation, rotation, and permutation symmetries, and is designed to be transferable to different system sizes. We demonstrate that using a relatively small number of training snapshots, Deep Density achieves excellent performance for one-dimensional insulating and metallic systems, as well as systems with mixed insulating and metallic characters. We also demonstrate its performance for real three-dimensional systems, including small organic molecules, as well as extended systems such as water (up to $512$ molecules) and aluminum (up to $256$ atoms).
Quantum embedding theories are playing an increasingly important role in bridging different levels of approximation to the many body Schrödinger equation in physics, chemistry and materials science. In this paper, we present a linear algebra perspective of the recently developed projection based embedding theory (PET) [Manby et al, J. Chem. Theory Comput. 8, 2564, 2012], restricted to the context of Kohn-Sham density functional theory. By partitioning the global degrees of freedom into a `system' part and a `bath' part, and by choosing a proper projector from the bath, PET is an in principle exact formulation to confine the calculation to the system part only, and hence can be carried out with reduced computational cost. Viewed from the perspective of the domain decomposition method, one particularly interesting feature of PET is that it does not enforce a boundary condition explicitly, and remains applicable even when the discretized Hamiltonian matrix is dense, such as in the context of the planewave discretization. In practice, the accuracy of PET depends on the accuracy of the projector for the bath. Based on the linear algebra reformulation, we develop a first order perturbation correction to the projector from the bath to improve its accuracy. Numerical results for real chemical systems indicate that with a proper choice of reference system, the perturbatively corrected PET can be sufficiently accurate even when strong perturbation is applied to very small systems, such as the computation of the ground state energy of a SiH$_3$F molecule, using a SiH$_4$ molecule as the reference system.
Localized representation of the Kohn-Sham subspace plays an important role in quantum chemistry and materials science. The recently developed selected columns of the density matrix (SCDM) method [J. Chem. Theory Comput. 11, 1463, 2015] is a simple and robust procedure for finding a localized representation of a set of Kohn-Sham orbitals from an insulating system. The SCDM method allows the direct construction of a well conditioned (or even orthonormal) and localized basis for the Kohn-Sham subspace. The SCDM algorithm avoids the use of an optimization procedure and does not depend on any adjustable parameters. The most computationally expensive step of the SCDM method is a column pivoted QR factorization that identifies the important columns for constructing the localized basis set. In this paper, we develop a two stage approximate column selection strategy to find the important columns at much lower computational cost. We demonstrate the effectiveness of this process using the dissociation process of a BH$_3$NH$_3$ molecule, an alkane chain and a supercell with 256 water molecules. Numerical results for the large collection of water molecules show that two stage localization procedure can be more than 30 times faster than the original SCDM algorithm and compares favorably with the popular Wannier90 package.
This paper generalizes the parallel selected inversion algorithm called PSelInv to sparse non- symmetric matrices. We assume a general sparse matrix A has been decomposed as PAQ = LU on a distributed memory parallel machine, where L, U are lower and upper triangular matrices, and P, Q are permutation matrices, respectively. The PSelInv method computes selected elements of A-1. The selection is confined by the sparsity pattern of the matrix AT . Our algorithm does not assume any symmetry properties of A, and our parallel implementation is memory efficient, in the sense that the computed elements of A-T overwrites the sparse matrix L+U in situ. PSelInv involves a large number of collective data communication activities within different processor groups of various sizes. In order to minimize idle time and improve load balancing, tree-based asynchronous communication is used to coordinate all such collective communication. Numerical results demonstrate that PSelInv can scale efficiently to 6,400 cores for a variety of matrices.
Recently, we have proposed the adaptive local basis set for electronic structure calculations based on Kohn-Sham density functional theory in a pseudopotential framework. The adaptive local basis set is efficient and systematically improvable for total energy calculations. In this paper, we present the calculation of atomic forces, which can be used for a range of applications such as geometry optimization and molecular dynamics simulation. We demonstrate that, under mild assumptions, the computation of atomic forces can scale nearly linearly with the number of atoms in the system using the adaptive local basis set. We quantify the accuracy of the Hellmann-Feynman forces for a range of physical systems, benchmarked against converged planewave calculations, and find that the adaptive local basis set is efficient for both force and energy calculations, requiring at most a few tens of basis functions per atom to attain accuracy required in practice. Since the adaptive local basis set has implicit dependence on atomic positions, Pulay forces are in general nonzero. However, we find that the Pulay force is numerically small and systematically decreasing with increasing basis completeness, so that the Hellmann-Feynman force is sufficient for basis sizes of a few tens of basis functions per atom. We verify the accuracy of the computed forces in static calculations of quasi-1D and 3D disordered Si systems, vibration calculation of a quasi-1D Si system, and molecular dynamics calculations of H$_2$ and liquid Al-Si alloy systems, where we find excellent agreement with independent benchmark results in literature.
We develop a hierarchical matrix construction algorithm using matrix-vector multiplications, based on the randomized singular value decomposition of low-rank matrices. The algorithm uses $\mathcal{O}(\log n)$ applications of the matrix on structured random test vectors and $\mathcal{O}(n \log n)$ extra computational cost, where $n$ is the dimension of the unknown matrix. Numerical examples on constructing Green's functions for elliptic operators in two dimensions show efficiency and accuracy of the proposed algorithm.