Lexing Ying

Xun Tang, Haoxuan Chen, Yuehaw Khoo et al. · stanford

We introduce Sketch Tomography, an efficient procedure for quantum state tomography based on the classical shadow protocol used for quantum observable estimations. The procedure applies to the case where the ground truth quantum state is a matrix product state (MPS). The density matrix of the ground truth state admits a tensor train ansatz as a result of the MPS assumption, and we estimate the tensor components of the ansatz through a series of observable estimations, thus outputting an approximation of the density matrix. The procedure is provably convergent with a sample complexity that scales quadratically in the system size. We conduct extensive numerical experiments to show that the procedure outputs an accurate approximation to the quantum state. For observable estimation tasks involving moderately large subsystems, we show that our procedure gives rise to a more accurate estimation than the classical shadow protocol. We also show that sketch tomography is more accurate in observable estimation than quantum states trained from the maximum likelihood estimation formulation.

Jack Poulson, Björn Engquist, Siwei Li et al.

A parallelization of a sweeping preconditioner for 3D Helmholtz equations without large cavities is introduced and benchmarked for several challenging velocity models. The setup and application costs of the sequential preconditioner are shown to be O(γ^2 N^{4/3}) and O(γ N log N), where γ(ω) denotes the modestly frequency-dependent number of grid points per Perfectly Matched Layer. Several computational and memory improvements are introduced relative to using black-box sparse-direct solvers for the auxiliary problems, and competitive runtimes and iteration counts are reported for high-frequency problems distributed over thousands of cores. Two open-source packages are released along with this paper: "Parallel Sweeping Preconditioner (PSP)" and the underlying distributed multifrontal solver, "Clique".

Jikai Jin, Yiping Lu, Jose Blanchet et al. · stanford

Learning mappings between infinite-dimensional function spaces has achieved empirical success in many disciplines of machine learning, including generative modeling, functional data analysis, causal inference, and multi-agent reinforcement learning. In this paper, we study the statistical limit of learning a Hilbert-Schmidt operator between two infinite-dimensional Sobolev reproducing kernel Hilbert spaces. We establish the information-theoretic lower bound in terms of the Sobolev Hilbert-Schmidt norm and show that a regularization that learns the spectral components below the bias contour and ignores the ones that are above the variance contour can achieve the optimal learning rate. At the same time, the spectral components between the bias and variance contours give us flexibility in designing computationally feasible machine learning algorithms. Based on this observation, we develop a multilevel kernel operator learning algorithm that is optimal when learning linear operators between infinite-dimensional function spaces.

Yinuo Ren, Hongli Zhao, Yuehaw Khoo et al. · stanford

We propose the tensorizing flow method for estimating high-dimensional probability density functions from the observed data. The method is based on tensor-train and flow-based generative modeling. Our method first efficiently constructs an approximate density in the tensor-train form via solving the tensor cores from a linear system based on the kernel density estimators of low-dimensional marginals. We then train a continuous-time flow model from this tensor-train density to the observed empirical distribution by performing a maximum likelihood estimation. The proposed method combines the optimization-less feature of the tensor-train with the flexibility of the flow-based generative models. Numerical results are included to demonstrate the performance of the proposed method.

Yiping Lu, Jose Blanchet, Lexing Ying · stanford

In this paper, we study the statistical limits in terms of Sobolev norms of gradient descent for solving inverse problem from randomly sampled noisy observations using a general class of objective functions. Our class of objective functions includes Sobolev training for kernel regression, Deep Ritz Methods (DRM), and Physics Informed Neural Networks (PINN) for solving elliptic partial differential equations (PDEs) as special cases. We consider a potentially infinite-dimensional parameterization of our model using a suitable Reproducing Kernel Hilbert Space and a continuous parameterization of problem hardness through the definition of kernel integral operators. We prove that gradient descent over this objective function can also achieve statistical optimality and the optimal number of passes over the data increases with sample size. Based on our theory, we explain an implicit acceleration of using a Sobolev norm as the objective function for training, inferring that the optimal number of epochs of DRM becomes larger than the number of PINN when both the data size and the hardness of tasks increase, although both DRM and PINN can achieve statistical optimality.

Emmanuel Candes, Laurent Demanet, Lexing Ying

This paper is concerned with the fast computation of Fourier integral operators of the general form $\int_{\R^d} e^{2\piıΦ(x,k)} f(k) d k$, where $k$ is a frequency variable, $Φ(x,k)$ is a phase function obeying a standard homogeneity condition, and $f$ is a given input. This is of interest for such fundamental computations are connected with the problem of finding numerical solutions to wave equations, and also frequently arise in many applications including reflection seismology, curvilinear tomography and others. In two dimensions, when the input and output are sampled on $N \times N$ Cartesian grids, a direct evaluation requires $O(N^4)$ operations, which is often times prohibitively expensive. This paper introduces a novel algorithm running in $O(N^2 \log N)$ time, i. e. with near-optimal computational complexity, and whose overall structure follows that of the butterfly algorithm [Michielssen and Boag, IEEE Trans Antennas Propagat 44 (1996), 1086-1093]. Underlying this algorithm is a mathematical insight concerning the restriction of the kernel $e^{2\piıΦ(x,k)}$ to subsets of the time and frequency domains. Whenever these subsets obey a simple geometric condition, the restricted kernel has approximately low-rank; we propose constructing such low-rank approximations using a special interpolation scheme, which prefactors the oscillatory component, interpolates the remaining nonoscillatory part and, lastly, remodulates the outcome. A byproduct of this scheme is that the whole algorithm is highly efficient in terms of memory requirement. Numerical results demonstrate the performance and illustrate the empirical properties of this algorithm.

Yuehaw Khoo, Jianfeng Lu, Lexing Ying

The curse of dimensionality is commonly encountered in numerical partial differential equations (PDE), especially when uncertainties have to be modeled into the equations as random coefficients. However, very often the variability of physical quantities derived from a PDE can be captured by a few features on the space of the coefficient fields. Based on such an observation, we propose using a neural-network (NN) based method to parameterize the physical quantity of interest as a function of input coefficients. The representability of such quantity using a neural-network can be justified by viewing the neural-network as performing time evolution to find the solutions to the PDE. We further demonstrate the simplicity and accuracy of the approach through notable examples of PDEs in engineering and physics.

Yiping Lu, Jiajin Li, Lexing Ying et al. · stanford

The optimal design of experiments typically involves solving an NP-hard combinatorial optimization problem. In this paper, we aim to develop a globally convergent and practically efficient optimization algorithm. Specifically, we consider a setting where the pre-treatment outcome data is available and the synthetic control estimator is invoked. The average treatment effect is estimated via the difference between the weighted average outcomes of the treated and control units, where the weights are learned from the observed data. {Under this setting, we surprisingly observed that the optimal experimental design problem could be reduced to a so-called \textit{phase synchronization} problem.} We solve this problem via a normalized variant of the generalized power method with spectral initialization. On the theoretical side, we establish the first global optimality guarantee for experiment design when pre-treatment data is sampled from certain data-generating processes. Empirically, we conduct extensive experiments to demonstrate the effectiveness of our method on both the US Bureau of Labor Statistics and the Abadie-Diemond-Hainmueller California Smoking Data. In terms of the root mean square error, our algorithm surpasses the random design by a large margin.

Björn Engquist, Lexing Ying

This paper introduces a new sweeping preconditioner for the iterative solution of the variable coefficient Helmholtz equation in two and three dimensions. The algorithms follow the general structure of constructing an approximate $LDL^t$ factorization by eliminating the unknowns layer by layer starting from an absorbing layer or boundary condition. The central idea of this paper is to approximate the Schur complement matrices of the factorization using moving perfectly matched layers (PMLs) introduced in the interior of the domain. Applying each Schur complement matrix is equivalent to solving a quasi-1D problem with a banded LU factorization in the 2D case and to solving a quasi-2D problem with a multifrontal method in the 3D case. The resulting preconditioner has linear application cost and the preconditioned iterative solver converges in a number of iterations that is essentially indefinite of the number of unknowns or the frequency. Numerical results are presented in both two and three dimensions to demonstrate the efficiency of this new preconditioner.

Björn Engquist, Lexing Ying

The paper introduces the sweeping preconditioner, which is highly efficient for iterative solutions of the variable coefficient Helmholtz equation including very high frequency problems. The first central idea of this novel approach is to construct an approximate factorization of the discretized Helmholtz equation by sweeping the domain layer by layer, starting from an absorbing layer or boundary condition. Given this specific order of factorization, the second central idea of this approach is to represent the intermediate matrices in the hierarchical matrix framework. In two dimensions, both the construction and the application of the preconditioners are of linear complexity. The GMRES solver with the resulting preconditioner converges in an amazingly small number of iterations, which is essentially independent of the number of unknowns. This approach is also extended to the three dimensional case with some success. Numerical results are provided in both two and three dimensions to demonstrate the efficiency of this new approach.

Lin Lin, Jianfeng Lu, Lexing Ying et al.

Kohn-Sham density functional theory is one of the most widely used electronic structure theories. In the pseudopotential framework, uniform discretization of the Kohn-Sham Hamiltonian generally results in a large number of basis functions per atom in order to resolve the rapid oscillations of the Kohn-Sham orbitals around the nuclei. Previous attempts to reduce the number of basis functions per atom include the usage of atomic orbitals and similar objects, but the atomic orbitals generally require fine tuning in order to reach high accuracy. We present a novel discretization scheme that adaptively and systematically builds the rapid oscillations of the Kohn-Sham orbitals around the nuclei as well as environmental effects into the basis functions. The resulting basis functions are localized in the real space, and are discontinuous in the global domain. The continuous Kohn-Sham orbitals and the electron density are evaluated from the discontinuous basis functions using the discontinuous Galerkin (DG) framework. Our method is implemented in parallel and the current implementation is able to handle systems with at least thousands of atoms. Numerical examples indicate that our method can reach very high accuracy (less than 1meV) with a very small number ($4\sim 40$) of basis functions per atom.

Victor Minden, Kenneth L. Ho, Anil Damle et al.

We introduce the strong recursive skeletonization factorization (RS-S), a new approximate matrix factorization based on recursive skeletonization for solving discretizations of linear integral equations associated with elliptic partial differential equations in two and three dimensions (and other matrices with similar hierarchical rank structure). Unlike previous skeletonization-based factorizations, RS-S uses a simple modification of skeletonization, strong skeletonization, which compresses only far-field interactions. This leads to an approximate factorization in the form of a product of many block unit-triangular matrices that may be used as a preconditioner or moderate-accuracy direct solver, with dramatically reduced rank growth. We further combine the strong skeletonization procedure with alternating near-field compression to obtain the hybrid recursive skeletonization factorization (RS-WS), a modification of RS-S that exhibits reduced storage cost in many settings. Under suitable rank assumptions both RS-S and RS-WS exhibit linear computational complexity, which we demonstrate with a number of numerical examples.

Lexing Ying

We introduce a fast algorithm for computing sparse Fourier transforms supported on smooth curves or surfaces. This problem appear naturally in several important problems in wave scattering and reflection seismology. The main observation is that the interaction between a frequency region and a spatial region is approximately low rank if the product of their radii are bounded by the maximum frequency. Based on this property, equivalent sources located at Cartesian grids are used to speed up the computation of the interaction between these two regions. The overall structure of our algorithm follows the recently-introduced butterfly algorithm. The computation is further accelerated by exploiting the tensor-product property of the Fourier kernel in two and three dimensions. The proposed algorithm is accurate and has an $O(N \log N)$ complexity. Finally, we present numerical results in both two and three dimensions.

Fei Liu, Lexing Ying

We introduce a new additive sweeping preconditioner for the Helmholtz equation based on the perfect matched layer (PML). This method divides the domain of interest into thin layers and proposes a new transmission condition between the subdomains where the emphasis is on the boundary values of the intermediate waves. This approach can be viewed as an effective approximation of an additive decomposition of the solution operator. When combined with the standard GMRES solver, the iteration number is essentially independent of the frequency. Several numerical examples are tested to show the efficiency of this new approach.

Yuwei Fan, Jordi Feliu-Faba, Lin Lin et al.

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.

Chao Ma, Daniel Kunin, Lei Wu et al.

A quadratic approximation of neural network loss landscapes has been extensively used to study the optimization process of these networks. Though, it usually holds in a very small neighborhood of the minimum, it cannot explain many phenomena observed during the optimization process. In this work, we study the structure of neural network loss functions and its implication on optimization in a region beyond the reach of a good quadratic approximation. Numerically, we observe that neural network loss functions possesses a multiscale structure, manifested in two ways: (1) in a neighborhood of minima, the loss mixes a continuum of scales and grows subquadratically, and (2) in a larger region, the loss shows several separate scales clearly. Using the subquadratic growth, we are able to explain the Edge of Stability phenomenon [5] observed for the gradient descent (GD) method. Using the separate scales, we explain the working mechanism of learning rate decay by simple examples. Finally, we study the origin of the multiscale structure and propose that the non-convexity of the models and the non-uniformity of training data is one of the causes. By constructing a two-layer neural network problem we show that training data with different magnitudes give rise to different scales of the loss function, producing subquadratic growth and multiple separate scales.

Laurent Demanet, Lexing Ying

This paper deals with efficient numerical representation and manipulation of differential and integral operators as symbols in phase-space, i.e., functions of space $x$ and frequency $ξ$. The symbol smoothness conditions obeyed by many operators in connection to smooth linear partial differential equations allow to write fast-converging, non-asymptotic expansions in adequate systems of rational Chebyshev functions or hierarchical splines. The classical results of closedness of such symbol classes under multiplication, inversion and taking the square root translate into practical iterative algorithms for realizing these operations directly in the proposed expansions. Because symbol-based numerical methods handle operators and not functions, their complexity depends on the desired resolution $N$ very weakly, typically only through $\log N$ factors. We present three applications to computational problems related to wave propagation: 1) preconditioning the Helmholtz equation, 2) decomposing wavefields into one-way components and 3) depth-stepping in reflection seismology.

Lexing Ying

The Lippmann-Schwinger equation is an integral equation formulation for acoustic and electromagnetic scattering from an inhomogeneous media and quantum scattering from a localized potential. We present the sparsifying preconditioner for accelerating the iterative solution of the Lippmann-Schwinger equation. This new preconditioner transforms the discretized Lippmann-Schwinger equation into sparse form and leverages the efficient sparse linear algebra algorithms for computing an approximate inverse. This preconditioner is efficient and easy to implement. When combined with standard iterative methods, it results in almost frequency-independent iteration counts. We provide 2D and 3D numerical results to demonstrate the effectiveness of this new preconditioner.

Fei Liu, Lexing Ying

This paper presents an efficient preconditioner for the Lippmann-Schwinger equation that combines the ideas of the sparsifying and the sweeping preconditioners. Following first the idea of the sparsifying preconditioner, this new preconditioner starts by transforming the dense linear system of the Lippmann-Schwinger equation into a nearly sparse system. The key novelty is a newly designed perfectly matched layer (PML) stencil for the boundary degrees of freedoms. The resulting sparse system gives rise to fairly accurate solutions and hence can be viewed as an accurate discretization of the Helmholtz equation. This new PML stencil also paves the way for applying the moving PML sweeping preconditioner to invert the resulting sparse system approximately. When combined with the standard GMRES solver, this new preconditioner for the Lippmann-Schwinger equation takes only a few iterations to converge for both 2D and 3D problems, where the iteration numbers are almost independent of the frequency. To the best of our knowledge, this is the first method that achieves near-linear cost to solve the 3D Lippmann-Schwinger equation in high frequency cases.

Lexing Ying, Sergey Fomel

We introduce two efficient algorithms for computing the partial Fourier transforms in one and two dimensions. Our study is motivated by the wave extrapolation procedure in reflection seismology. In both algorithms, the main idea is to decompose the summation domain of into simpler components in a multiscale way. Existing fast algorithms are then applied to each component to obtain optimal complexity. The algorithm in 1D is exact and takes $O(N\log^2 N)$ steps. Our solution in 2D is an approximate but accurate algorithm that takes $O(N^2 \log^2 N)$ steps. In both cases, the complexities are almost linear in terms of the degree of freedom. We provide numerical results on several test examples.

Yuwei Fan, Jing An, Lexing Ying

We investigate integral formulations and fast algorithms for the steady-state radiative transfer equation with isotropic and anisotropic scattering. When the scattering term is a smooth convolution on the unit sphere, a model reduction step in the angular domain using the Fourier transformation in 2D and the spherical harmonic transformation in 3D significantly reduces the number of degrees of freedoms. The resulting Fourier coefficients or spherical harmonic coefficients satisfy a Fredholm integral equation of the second kind. We study the uniqueness of the equation and proved an a priori estimate. For a homogeneous medium, the integral equation can be solved efficiently using the FFT and iterative methods. For an inhomogeneous medium, the recursive skeletonization factorization method is applied instead. Numerical simulations demonstrate the efficiency of the proposed algorithms in both homogeneous and inhomogeneous cases and for both transport and diffusion regimes.

Lin Lin, Ze Xu, Lexing Ying

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.

Laurent Demanet, Lexing Ying

This paper presents a numerical compression strategy for the boundary integral equation of acoustic scattering in two dimensions. These equations have oscillatory kernels that we represent in a basis of wave atoms, and compress by thresholding the small coefficients to zero. This phenomenon was perhaps first observed in 1993 by Bradie, Coifman, and Grossman, in the context of local Fourier bases \cite{BCG}. Their results have since then been extended in various ways. The purpose of this paper is to bridge a theoretical gap and prove that a well-chosen fixed expansion, the nonstandard wave atom form, provides a compression of the acoustic single and double layer potentials with wave number $k$ as $O(k)$-by-$O(k)$ matrices with $O(k^{1+1/\infty})$ nonnegligible entries, with a constant that depends on the relative $\ell_2$ accuracy $\eps$ in an acceptable way. The argument assumes smooth, separated, and not necessarily convex scatterers in two dimensions. The essential features of wave atoms that enable to write this result as a theorem is a sharp time-frequency localization that wavelet packets do not obey, and a parabolic scaling wavelength $\sim$ (essential diameter)${}^2$. Numerical experiments support the estimate and show that this wave atom representation may be of interest for applications where the same scattering problem needs to be solved for many boundary conditions, for example, the computation of radar cross sections.

Fei Liu, Lexing Ying

This paper introduces the recursive sweeping preconditioner for the numerical solution of the Helmholtz equation in 3D. This is based on the earlier work of the sweeping preconditioner with the moving perfectly matched layers (PMLs). The key idea is to apply the sweeping preconditioner recursively to the quasi-2D auxiliary problems introduced in the 3D sweeping preconditioner. Compared to the non-recursive 3D sweeping preconditioner, the setup cost of this new approach drops from $O(N^{4/3})$ to $O(N)$, the application cost per iteration drops from $O(N\log N)$ to $O(N)$, and the iteration count only increases mildly when combined with the standard GMRES solver. Several numerical examples are tested and the results are compared with the non-recursive sweeping preconditioner to demonstrate the efficiency of the new approach.

Yiping Lu, Wenlong Ji, Zachary Izzo et al. · stanford

Although overparameterized models have shown their success on many machine learning tasks, the accuracy could drop on the testing distribution that is different from the training one. This accuracy drop still limits applying machine learning in the wild. At the same time, importance weighting, a traditional technique to handle distribution shifts, has been demonstrated to have less or even no effect on overparameterized models both empirically and theoretically. In this paper, we propose importance tempering to improve the decision boundary and achieve consistently better results for overparameterized models. Theoretically, we justify that the selection of group temperature can be different under label shift and spurious correlation setting. At the same time, we also prove that properly selected temperatures can extricate the minority collapse for imbalanced classification. Empirically, we achieve state-of-the-art results on worst group classification tasks using importance tempering.

Lexing Ying

We introduce a new coarse-graining algorithm, tensor network skeletonization, for the numerical computation of tensor networks. This approach utilizes a structure-preserving skeletonization procedure to remove short-range correlations effectively at every scale. This approach is first presented in the setting of 2D statistical Ising model and is then extended to higher dimensional tensor networks and disordered systems. When applied to the Euclidean path integral formulation, this approach also gives rise to new efficient representations of the ground states for 1D and 2D quantum Ising models.

Yingzhou Li, Haizhao Yang, Lexing Ying

This paper presents an efficient multiscale butterfly algorithm for computing Fourier integral operators (FIOs) of the form $(\mathcal{L} f)(x) = \int_{R^d}a(x,ξ) e^{2πıΦ(x,ξ)}\hat{f}(ξ) dξ$, where $Φ(x,ξ)$ is a phase function, $a(x,ξ)$ is an amplitude function, and $f(x)$ is a given input. The frequency domain is hierarchically decomposed into a union of Cartesian coronas. The integral kernel $a(x,ξ) e^{2πıΦ(x,ξ)}$ in each corona satisfies a special low-rank property that enables the application of a butterfly algorithm on the Cartesian phase-space grid. This leads to an algorithm with quasi-linear operation complexity and linear memory complexity. Different from previous butterfly methods for the FIOs, this new approach is simple and reduces the computational cost by avoiding extra coordinate transformations. Numerical examples in two and three dimensions are provided to demonstrate the practical advantages of the new algorithm.

Lin Lin, Jianfeng Lu, Lexing Ying et al.

We develop a technique for generating a set of optimized local basis functions to solve models in the Kohn-Sham density functional theory for both insulating and metallic systems. The optimized local basis functions are obtained by solving a minimization problem in an admissible set determined by a large number of primitive basis functions. Using the optimized local basis set, the electron energy and the atomic force can be calculated accurately with a small number of basis functions. The Pulay force is systematically controlled and is not required to be calculated, which makes the optimized local basis set an ideal tool for ab initio molecular dynamics and structure optimization. We also propose a preconditioned Newton-GMRES method to obtain the optimized local basis functions in practice. The optimized local basis set is able to achieve high accuracy with a small number of basis functions per atom when applied to a one dimensional model problem.

Victor Minden, Anil Damle, Kenneth L. Ho et al.

We present a method for updating certain hierarchical factorizations for solving linear integral equations with elliptic kernels. In particular, given a factorization corresponding to some initial geometry or material parameters, we can locally perturb the geometry or coefficients and update the initial factorization to reflect this change with asymptotic complexity that is polylogarithmic in the total number of unknowns and linear in the number of perturbed unknowns. We apply our method to the recursive skeletonization factorization and hierarchical interpolative factorization and demonstrate scaling results for a number of different 2D problem setups.

Jianfeng Lu, Lexing Ying

We propose an efficient algorithm for density fitting of Bloch waves for Hamiltonian operators with periodic potential. The algorithm is based on column selection and random Fourier projection of the orbital functions. The computational cost of the algorithm scales as $\mathcal{O}\bigl(N_{\text{grid}} N^2 + N_{\text{grid}} NK \log (NK)\bigr)$, where $N_{\text{grid}}$ is number of spatial grid points, $K$ is the number of sampling $k$-points in first Brillouin zone, and $N$ is the number of bands under consideration. We validate the algorithm by numerical examples in both two and three dimensions.

Xun Tang, Yoonhaeng Hur, Yuehaw Khoo et al.

In this paper, we present a density estimation framework based on tree tensor-network states. The proposed method consists of determining the tree topology with Chow-Liu algorithm, and obtaining a linear system of equations that defines the tensor-network components via sketching techniques. Novel choices of sketch functions are developed in order to consider graphical models that contain loops. Sample complexity guarantees are provided and further corroborated by numerical experiments.

Fei Liu, Lexing Ying

This paper presents the sparsifying preconditioner for the time-harmonic Maxwell's equations in the integral formulation. Following the work on sparsifying preconditioner for the Lippmann-Schwinger equation, this paper generalizes that approach from the scalar wave case to the vector case. The key idea is to construct a sparse approximation to the dense system by minimizing the non-local interactions in the integral equation, which allows for applying sparse linear solvers to reduce the computational cost. When combined with the standard GMRES solver, the number of preconditioned iterations remains small and essentially independent of the frequency. This suggests that, when the sparsifying preconditioner is adopted, solving the dense integral system can be done as efficiently as solving the sparse system from PDE discretization.

Yinuo Ren, Tesi Xiao, Tanmay Gangwani et al. · stanford

Multi-objective optimization (MOO) aims to optimize multiple, possibly conflicting objectives with widespread applications. We introduce a novel interacting particle method for MOO inspired by molecular dynamics simulations. Our approach combines overdamped Langevin and birth-death dynamics, incorporating a "dominance potential" to steer particles toward global Pareto optimality. In contrast to previous methods, our method is able to relocate dominated particles, making it particularly adept at managing Pareto fronts of complicated geometries. Our method is also theoretically grounded as a Wasserstein-Fisher-Rao gradient flow with convergence guarantees. Extensive experiments confirm that our approach outperforms state-of-the-art methods on challenging synthetic and real-world datasets.

Lexing Ying

The boundary integral method is an efficient approach for solving time-harmonic obstacle scattering problems by a bounded scatterer. This paper presents the directional preconditioner for the iterative solution of linear systems of the boundary integral method. This new preconditioner builds a data-sparse approximation of the integral operator, transforms it into a sparse linear system, and computes an approximate inverse with efficient sparse and hierarchical linear algebra algorithms. This preconditioner is efficient and results in small and almost frequency-independent iteration counts when combined with standard iterative solvers. Numerical results are provided to demonstrate the effectiveness of the new preconditioner.

Lexing Ying

This paper introduces the sparsifying preconditioner for the pseudospectral approximation of highly indefinite systems on periodic structures, which include the frequency-domain response problems of the Helmholtz equation and the Schrödinger equation as examples. This approach transforms the dense system of the pseudospectral discretization approximately into an sparse system via an equivalent integral reformulation and a specially-designed sparsifying operator. The resulting sparse system is then solved efficiently with sparse linear algebra algorithms and serves as a reasonably accurate preconditioner. When combined with standard iterative methods, this new preconditioner results in small iteration counts. Numerical results are provided for the Helmholtz equation and the Schrödinger in both 2D and 3D to demonstrate the effectiveness of this new preconditioner.

Lexing Ying

The boundary integral method is an efficient approach for solving time-harmonic acoustic obstacle scattering problems. The main computational task is the evaluation of an oscillatory boundary integral at each discretization point of the boundary. This paper presents a new fast algorithm for this task in two dimensions. This algorithm is built on top of directional low-rank approximations of the scattering kernel and uses oscillatory Chebyshev interpolation and local FFTs to achieve quasi-linear complexity. The algorithm is simple, fast, and kernel-independent. Numerical results are provided to demonstrate the effectiveness of the proposed algorithm.

Jianfeng Lu, Lexing Ying

We develop a robust and efficient method for soliton calculations for nonlinear Schrödinger equations. The method is based on the recently developed sparsifying preconditioner combined with Newton's iterative method. The performance of the method is demonstrated by numerical examples of gap solitons in the context of nonlinear optics.

Björn Engquist, Lexing Ying

This paper introduces a directional multiscale algorithm for the two dimensional $N$-body problem of the Helmholtz kernel with applications to high frequency scattering. The algorithm follows the approach in [Engquist and Ying, SIAM Journal on Scientific Computing, 29 (4), 2007] where the three dimensional case was studied. The main observation is that, for two regions that follow a directional parabolic geometric configuration, the interaction between the points in these two regions through the Helmholtz kernel is approximately low rank. We propose an improved randomized procedure for generating the low rank representations. Based on these representations, we organize the computation of the far field interaction in a multidirectional and multiscale way to achieve maximum efficiency. The proposed algorithm is accurate and has the optimal $O(N\log N)$ complexity for problems from two dimensional scattering applications. We present numerical results for several test examples to illustrate the algorithm and its application to two dimensional high frequency scattering problems.

Xin Ye, Jianlin Xia, Lexing Ying

It has been known in potential theory that, for some kernels matrices corresponding to well-separated point sets, fast analytical low-rank approximation can be achieved via the use of proxy points. This proxy point method gives a surprisingly convenient way of explicitly writing out approximate basis matrices for a kernel matrix. However, this elegant strategy is rarely known or used in the numerical linear algebra community. It still needs clear algebraic understanding of the theoretical background. Moreover, rigorous quantifications of the approximation errors and reliable criteria for the selection of the proxy points are still missing. In this work, we use contour integration to clearly justify the idea in terms of a class of important kernels. We further provide comprehensive accuracy analysis for the analytical compression and show how to choose nearly optimal proxy points. The analytical compression is then combined with fast rank-revealing factorizations to get compact low-rank approximations and also to select certain representative points. We provide the error bounds for the resulting overall low-rank approximation. This work thus gives a fast and reliable strategy for compressing those kernel matrices. Furthermore, it provides an intuitive way of understanding the proxy point method and bridges the gap between this useful analytical strategy and practical low-rank approximations. Some numerical examples help to further illustrate the ideas.

Fei Liu, Lexing Ying

This paper introduces the localized sparsifying preconditioner for the pseudospectral approximations of indefinite systems on periodic structures. The work is built on top of the recently proposed sparsifying preconditioner with two major modifications. First, the local potential information is utilized to improve the accuracy of the preconditioner. Second, an FFT based method to compute the local stencil is proposed to reduce the setup time of the algorithm. Numerical results show that the iteration number of this improved method grows only mildly as the problem size grows, which implies that solving pseudospectral approximation systems is computationally as efficient as solving sparse systems, up to a mildly growing factor.

Yuhua Zhu, Zachary Izzo, Lexing Ying

This paper revisits the bandit problem in the Bayesian setting. The Bayesian approach formulates the bandit problem as an optimization problem, and the goal is to find the optimal policy which minimizes the Bayesian regret. One of the main challenges facing the Bayesian approach is that computation of the optimal policy is often intractable, especially when the length of the problem horizon or the number of arms is large. In this paper, we first show that under a suitable rescaling, the Bayesian bandit problem converges toward a continuous Hamilton-Jacobi-Bellman (HJB) equation. The optimal policy for the limiting HJB equation can be explicitly obtained for several common bandit problems, and we give numerical methods to solve the HJB equation when an explicit solution is not available. Based on these results, we propose an approximate Bayes-optimal policy for solving Bayesian bandit problems with large horizons. Our method has the added benefit that its computational cost does not increase as the horizon increases.

Samarth Gupta, Daniel N. Hill, Lexing Ying et al.

In most applications of model-based Markov decision processes, the parameters for the unknown underlying model are often estimated from the empirical data. Due to noise, the policy learnedfrom the estimated model is often far from the optimal policy of the underlying model. When applied to the environment of the underlying model, the learned policy results in suboptimal performance, thus calling for solutions with better generalization performance. In this work we take a Bayesian perspective and regularize the objective function of the Markov decision process with prior information in order to obtain more robust policies. Two approaches are proposed, one based on $L^1$ regularization and the other on relative entropic regularization. We evaluate our proposed algorithms on synthetic simulations and on real-world search logs of a large scale online shopping store. Our results demonstrate the robustness of regularized MDP policies against the noise present in the models.

Chao Ma, Lexing Ying

In this paper, we show that structures similar to self-attention are natural to learn many sequence-to-sequence problems from the perspective of symmetry. Inspired by language processing applications, we study the orthogonal equivariance of seq2seq functions with knowledge, which are functions taking two inputs -- an input sequence and a ``knowledge'' -- and outputting another sequence. The knowledge consists of a set of vectors in the same embedding space as the input sequence, containing the information of the language used to process the input sequence. We show that orthogonal equivariance in the embedding space is natural for seq2seq functions with knowledge, and under such equivariance the function must take the form close to the self-attention. This shows that network structures similar to self-attention are the right structures to represent the target function of many seq2seq problems. The representation can be further refined if a ``finite information principle'' is considered, or a permutation equivariance holds for the elements of the input sequence.

Zhenning Cai, Yuwei Fan, Lexing Ying

We propose an entropic Fourier method for the numerical discretization of the Boltzmann collision operator. The method, which is obtained by modifying a Fourier Galerkin method to match the form of the discrete velocity method, can be viewed both as a discrete velocity method and as a Fourier method. As a discrete velocity method, it preserves the positivity of the solution and satisfies a discrete version of the H-theorem. As a Fourier method, it allows one to readily apply the FFT-based fast algorithms. A second-order convergence rate is validated by numerical experiments

Lexing Ying

This note considers the unstructured sparse recovery problems in a general form. Examples include rational approximation, spectral function estimation, Fourier inversion, Laplace inversion, and sparse deconvolution. The main challenges are the noise in the sample values and the unstructured nature of the sample locations. This note proposes the eigenmatrix, a data-driven construction with desired approximate eigenvalues and eigenvectors. The eigenmatrix offers a new way for these sparse recovery problems. Numerical results are provided to demonstrate the efficiency of the proposed method.

Hongrui Chen, Lexing Ying

Diffusion models have achieved huge empirical success in data generation tasks. Recently, some efforts have been made to adapt the framework of diffusion models to discrete state space, providing a more natural approach for modeling intrinsically discrete data, such as language and graphs. This is achieved by formulating both the forward noising process and the corresponding reversed process as Continuous Time Markov Chains (CTMCs). In this paper, we investigate the theoretical properties of the discrete diffusion model. Specifically, we introduce an algorithm leveraging the uniformization of continuous Markov chains, implementing transitions on random time points. Under reasonable assumptions on the learning of the discrete score function, we derive Total Variation distance and KL divergence guarantees for sampling from any distribution on a hypercube. Our results align with state-of-the-art achievements for diffusion models in $\mathbb{R}^d$ and further underscore the advantages of discrete diffusion models in comparison to the $\mathbb{R}^d$ setting.

Yinuo Ren, Haoxuan Chen, Yuchen Zhu et al. · gatech, stanford

Discrete diffusion models have emerged as a powerful generative modeling framework for discrete data with successful applications spanning from text generation to image synthesis. However, their deployment faces challenges due to the high dimensionality of the state space, necessitating the development of efficient inference algorithms. Current inference approaches mainly fall into two categories: exact simulation and approximate methods such as $τ$-leaping. While exact methods suffer from unpredictable inference time and redundant function evaluations, $τ$-leaping is limited by its first-order accuracy. In this work, we advance the latter category by tailoring the first extension of high-order numerical inference schemes to discrete diffusion models, enabling larger step sizes while reducing error. We rigorously analyze the proposed schemes and establish the second-order accuracy of the $θ$-trapezoidal method in KL divergence. Empirical evaluations on GPT-2 level text and ImageNet-level image generation tasks demonstrate that our method achieves superior sample quality compared to existing approaches under equivalent computational constraints.

Rajat Vadiraj Dwaraknath, Lexing Ying

Metric graphs are structures obtained by associating edges in a standard graph with segments of the real line and gluing these segments at the vertices of the graph. The resulting structure has a natural metric that allows for the study of differential operators and stochastic processes on the graph. Brownian motions in these domains have been extensively studied theoretically using their generators. However, less work has been done on practical algorithms for simulating these processes. We introduce the first algorithm for simulating Brownian motions on metric graphs through a timestep splitting Euler-Maruyama-based discretization of their corresponding stochastic differential equation. By applying this scheme to Langevin diffusions on metric graphs, we also obtain the first algorithm for sampling on metric graphs. We provide theoretical guarantees on the number of timestep splittings required for the algorithm to converge to the underlying stochastic process. We also show that the exit probabilities of the simulated particle converge to the vertex-edge jump probabilities of the underlying stochastic differential equation as the timestep goes to zero. Finally, since this method is highly parallelizable, we provide fast, memory-aware implementations of our algorithm in the form of custom CUDA kernels that are up to ~8000x faster than a GPU implementation using PyTorch on simple star metric graphs. Beyond simple star graphs, we benchmark our algorithm on a real cortical vascular network extracted from a DuMuX tissue-perfusion model for tracer transport. Our algorithm is able to run stable simulations with timesteps significantly larger than the stable limit of the finite volume method used in DuMuX while also achieving speedups of up to ~1500x.

Haoxuan Chen, Yinuo Ren, Martin Renqiang Min et al. · stanford

Diffusion models (DMs) have proven to be effective in modeling high-dimensional distributions, leading to their widespread adoption for representing complex priors in Bayesian inverse problems (BIPs). However, current DM-based posterior sampling methods proposed for solving common BIPs rely on heuristic approximations to the generative process. To exploit the generative capability of DMs and avoid the usage of such approximations, we propose an ensemble-based algorithm that performs posterior sampling without the use of heuristic approximations. Our algorithm is motivated by existing works that combine DM-based methods with the sequential Monte Carlo (SMC) method. By examining how the prior evolves through the diffusion process encoded by the pre-trained score function, we derive a modified partial differential equation (PDE) governing the evolution of the corresponding posterior distribution. This PDE includes a modified diffusion term and a reweighting term, which can be simulated via stochastic weighted particle methods. Theoretically, we prove that the error between the true posterior distribution can be bounded in terms of the training error of the pre-trained score function and the number of particles in the ensemble. Empirically, we validate our algorithm on several inverse problems in imaging to show that our method gives more accurate reconstructions compared to existing DM-based methods.

Neeraj Gangwar, Rishabh Deshmukh, Michael Shavlovsky et al.

As model sizes continue to grow, parameter-efficient fine-tuning has emerged as a powerful alternative to full fine-tuning. While LoRA is widely adopted among these methods, recent research has explored vector-based adaptation methods due to their extreme parameter efficiency. However, these methods typically require substantially higher ranks than LoRA to match its performance, leading to increased training costs. This work introduces GiVA, a gradient-based initialization strategy for vector-based adaptation. It achieves training times comparable to LoRA and maintains the extreme parameter efficiency of vector-based adaptation. We evaluate GiVA across diverse benchmarks, including natural language understanding, natural language generation, and image classification. Experiments show that our approach consistently outperforms or achieves performance competitive with existing vector-based adaptation methods and LoRA while reducing rank requirements by a factor of eight ($8\times$).