Chao Yang

Lin Lin, Yousef Saad, Chao Yang

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.

Ruipeng Li, Yuanzhe Xi, Eugene Vecharynski et al.

Polynomial filtering can provide a highly effective means of computing all eigenvalues of a real symmetric (or complex Hermitian) matrix that are located in a given interval, anywhere in the spectrum. This paper describes a technique for tackling this problem by combining a Thick-Restart version of the Lanczos algorithm with deflation (`locking') and a new type of polynomial filters obtained from a least-squares technique. The resulting algorithm can be utilized in a `spectrum-slicing' approach whereby a very large number of eigenvalues and associated eigenvectors of the matrix are computed by extracting eigenpairs located in different sub-intervals independently from one another.

R. Huang, J. Sun, C. Yang

Recently, a non-classical eigenvalue solver, called RIM, was proposed to compute (all) eigenvalues in a region on the complex plane. Without solving any eigenvalue problem, it tests if a region contains eigenvalues using an approximate spectral projection. Regions that contain eigenvalues are subdivided and tested recursively until eigenvalues are isolated with a specified precision. This makes RIM an eigensolver distinct from all existing methods. Furthermore, it requires no a priori spectral information. In this paper, we propose an improved version of {\bf RIM} for non-Hermitian eigenvalue problems. Using Cayley transformation and Arnoldi's method, the computation cost is reduced significantly. Effectiveness and efficiency of the new method are demonstrated by numerical examples and compared with 'eigs' in Matlab.

Meiyue Shao, Chao Yang

This is the user manual for the software package BSEPACK (Bethe--Salpeter Eigenvalue Solver Package).

Eugene Vecharynski, Chao Yang

We describe preconditioned iterative methods for estimating the number of eigenvalues of a Hermitian matrix within a given interval. Such estimation is useful in a number of applications.In particular, it can be used to develop an efficient spectrum-slicing strategy to compute many eigenpairs of a Hermitian matrix. Our method is based on the Lanczos- and Arnoldi-type of iterations. We show that with a properly defined preconditioner, only a few iterations may be needed to obtain a good estimate of the number of eigenvalues within a prescribed interval. We also demonstrate that the number of iterations required by the proposed preconditioned schemes is independent of the size and condition number of the matrix. The efficiency of the methods is illustrated on several problems arising from density functional theory based electronic structure calculations.

Chao Chao Yang, Tao Wang, Xiaoping Xie

This paper develops and analyses numerical approximation for linear-quadratic optimal control problem governed by elliptic interface equations. We adopt variational discretization concept to discretize optimal control problem, and apply an interface-unfitted finite element method due to [A. Hansbo and P. Hansbo. An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. Comput. Methods Appl. Mech. Engrg., 191(47-48): 5537-5552, 2002] to discretize corresponding state and adjoint equations, where piecewise cut basis functions around interface are enriched into standard conforming finite element space. Optimal error estimates in both $L^2$ norm and a mesh-dependent norm are derived for optimal state, co-state and control under different regularity assumptions. Numerical results verify the theoretical results.

Tao Wang, Chao Chao Yang, Xiaoping Xie

This paper analyzes two eXtended finite element methods (XFEMs) for linear quadratic optimal control problems governed by Poisson equation in non-convex domains. We follow the variational discretization concept to discretize the continuous problems, and apply an XFEM with a cut-off function and a classic XFEM with a fixed enrichment area to discretize the state and co-state equations. Optimal error estimates are derived for the state, co-state and control. Numerical results confirm our theoretical results.

Hardeep Bassi, Yuanran Zhu, Senwei Liang et al.

In this paper, we propose using LSTM-RNNs (Long Short-Term Memory-Recurrent Neural Networks) to learn and represent nonlinear integral operators that appear in nonlinear integro-differential equations (IDEs). The LSTM-RNN representation of the nonlinear integral operator allows us to turn a system of nonlinear integro-differential equations into a system of ordinary differential equations for which many efficient solvers are available. Furthermore, because the use of LSTM-RNN representation of the nonlinear integral operator in an IDE eliminates the need to perform a numerical integration in each numerical time evolution step, the overall temporal cost of the LSTM-RNN-based IDE solver can be reduced to $O(n_T)$ from $O(n_T^2)$ if a $n_T$-step trajectory is to be computed. We illustrate the efficiency and robustness of this LSTM-RNN-based numerical IDE solver with a model problem. Additionally, we highlight the generalizability of the learned integral operator by applying it to IDEs driven by different external forces. As a practical application, we show how this methodology can effectively solve the Dyson's equation for quantum many-body systems.

Hardeep Bassi, Yuanran Zhu, Erika Ye et al.

We present RELift (Restrict, Evolve, Lift), a two-phase learning framework that couples coarse-grid numerical solvers with neural operators to super-resolve and forecast fine-grid dynamics for time-dependent partial differential equations (PDEs). In Phase 1, RELift learns a super-resolution operator that maps the solution on a coarse grid to a fine grid. In Phase 2, this learned operator is composed with a coarse-grid numerical integrator to construct an effective fine-grid propagator for the governing equation. We benchmark RELift on three canonical two-dimensional PDEs of increasing dynamical complexity -- the heat equation, the wave equation, and the incompressible Navier--Stokes equations -- and we further demonstrate its performance on a kinetic electron dynamics case study via the 1D1V Vlasov--Poisson system. Across all examples, RELift delivers high-fidelity super-resolution (Phase 1) and accurate long-horizon rollouts (Phase 2), outperforming standard super-resolution and neural operator baselines in both field-level error metrics and physics-relevant diagnostics. Finally, we provide error analysis of the effective fine-grid propagator, characterizing how approximation errors accumulate over time and explaining the observed numerical stability of the RELift framework.

Senwei Liang, Aditya N. Singh, Yuanran Zhu et al.

We propose a reinforcement learning based method to identify important configurations that connect reactant and product states along chemical reaction paths. By shooting multiple trajectories from these configurations, we can generate an ensemble of configurations that concentrate on the transition path ensemble. This configuration ensemble can be effectively employed in a neural network-based partial differential equation solver to obtain an approximation solution of a restricted Backward Kolmogorov equation, even when the dimension of the problem is very high. The resulting solution, known as the committor function, encodes mechanistic information for the reaction and can in turn be used to evaluate reaction rates.

Aydin Buluc, Tamara G. Kolda, Stefan M. Wild et al.

Randomized algorithms have propelled advances in artificial intelligence and represent a foundational research area in advancing AI for Science. Future advancements in DOE Office of Science priority areas such as climate science, astrophysics, fusion, advanced materials, combustion, and quantum computing all require randomized algorithms for surmounting challenges of complexity, robustness, and scalability. This report summarizes the outcomes of that workshop, "Randomized Algorithms for Scientific Computing (RASC)," held virtually across four days in December 2020 and January 2021.

Gianina Alina Negoita, James P. Vary, Glenn R. Luecke et al.

Ab initio approaches in nuclear theory, such as the no-core shell model (NCSM), have been developed for approximately solving finite nuclei with realistic strong interactions. The NCSM and other approaches require an extrapolation of the results obtained in a finite basis space to the infinite basis space limit and assessment of the uncertainty of those extrapolations. Each observable requires a separate extrapolation and most observables have no proven extrapolation method. We propose a feed-forward artificial neural network (ANN) method as an extrapolation tool to obtain the ground state energy and the ground state point-proton root-mean-square (rms) radius along with their extrapolation uncertainties. The designed ANNs are sufficient to produce results for these two very different observables in $^6$Li from the ab initio NCSM results in small basis spaces that satisfy the following theoretical physics condition: independence of basis space parameters in the limit of extremely large matrices. Comparisons of the ANN results with other extrapolation methods are also provided.

Gaigong Zhang, Lin Lin, Wei Hu et al.

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.

Eugene Vecharynski, Chao Yang, Fei Xue

We introduce the Generalized Preconditioned Locally Harmonic Residual (GPLHR) method for solving standard and generalized non-Hermitian eigenproblems. The method is particularly useful for computing a subset of eigenvalues, and their eigen- or Schur vectors, closest to a given shift. The proposed method is based on block iterations and can take advantage of a preconditioner if it is available. It does not need to perform exact shift-and-invert transformation. Standard and generalized eigenproblems are handled in a unified framework. Our numerical experiments demonstrate that GPLHR is generally more robust and efficient than existing methods, especially if the available memory is limited.

Qingming Tang, Chao Yang, Jian Peng et al.

This paper considers the problem of estimating multiple related Gaussian graphical models from a $p$-dimensional dataset consisting of different classes. Our work is based upon the formulation of this problem as group graphical lasso. This paper proposes a novel hybrid covariance thresholding algorithm that can effectively identify zero entries in the precision matrices and split a large joint graphical lasso problem into small subproblems. Our hybrid covariance thresholding method is superior to existing uniform thresholding methods in that our method can split the precision matrix of each individual class using different partition schemes and thus split group graphical lasso into much smaller subproblems, each of which can be solved very fast. In addition, this paper establishes necessary and sufficient conditions for our hybrid covariance thresholding algorithm. The superior performance of our thresholding method is thoroughly analyzed and illustrated by a few experiments on simulated data and real gene expression data.