Qin Li

Thomas Y. Hou, Qin Li, Pengchuan Zhang

We propose a stochastic multiscale finite element method (StoMsFEM) to solve random elliptic partial differential equations with a high stochastic dimension. The key idea is to simultaneously upscale the stochastic solutions in the physical space for all random samples and explore the low stochastic dimensions of the stochastic solution within each local patch. We propose two effective methods to achieve this simultaneous local upscaling. The first method is a high order interpolation method in the stochastic space that explores the high regularity of the local upscaled quantities with respect to the random variables. The second method is a reduced-order method that explores the low rank property of the multiscale basis functions within each coarse grid patch. Our complexity analysis shows that compared with the standard FEM on a fine grid, the StoMsFEM can achieve computational saving in the order of $(H/h)^{d}/(\log(H/h))^k$, where $H/h$ is the ratio between the coarse and the fine gird sizes, $d$ is the physical dimension and $k$ is the local stochastic dimension. Several numerical examples are presented to demonstrate the accuracy and effectiveness of the proposed methods. In the high contrast example, we observe a factor of 2000 speed-up.

Ke Chen, Qin Li, Jianfeng Lu et al.

In the framework of generalized finite element methods for elliptic equations with rough coefficients, efficiency and accuracy of the numerical method depend critically on the use of appropriate basis functions. This work explores several random sampling strategies that construct approximations to the optimal set of basis functions of a given dimension, and proposes a quantitative criterion to analyze and compare these sampling strategies. Numerical evidence shows that the best results are achieved by two strategies, Random Gaussian and Smooth boundary sampling.

Yasser Abduallah, Vania K. Jordanova, Hao Liu et al.

Solar energetic particles (SEPs) are an essential source of space radiation, which are hazards for humans in space, spacecraft, and technology in general. In this paper we propose a deep learning method, specifically a bidirectional long short-term memory (biLSTM) network, to predict if an active region (AR) would produce an SEP event given that (i) the AR will produce an M- or X-class flare and a coronal mass ejection (CME) associated with the flare, or (ii) the AR will produce an M- or X-class flare regardless of whether or not the flare is associated with a CME. The data samples used in this study are collected from the Geostationary Operational Environmental Satellite's X-ray flare catalogs provided by the National Centers for Environmental Information. We select M- and X-class flares with identified ARs in the catalogs for the period between 2010 and 2021, and find the associations of flares, CMEs and SEPs in the Space Weather Database of Notifications, Knowledge, Information during the same period. Each data sample contains physical parameters collected from the Helioseismic and Magnetic Imager on board the Solar Dynamics Observatory. Experimental results based on different performance metrics demonstrate that the proposed biLSTM network is better than related machine learning algorithms for the two SEP prediction tasks studied here. We also discuss extensions of our approach for probabilistic forecasting and calibration with empirical evaluation.

Thomas Y. Hou, Qin Li, Pengchuan Zhang

Suppose that $A \in \mathbb{R}^{N \times N}$ is symmetric positive semidefinite with rank $K \le N$. Our goal is to decompose $A$ into $K$ rank-one matrices $\sum_{k=1}^K g_k g_k^T$ where the modes $\{g_{k}\}_{k=1}^K$ are required to be as sparse as possible. In contrast to eigen decomposition, these sparse modes are not required to be orthogonal. Such a problem arises in random field parametrization where $A$ is the covariance function and is intractable to solve in general. In this paper, we partition the indices from 1 to $N$ into several patches and propose to quantify the sparseness of a vector by the number of patches on which it is nonzero, which is called patch-wise sparseness. Our aim is to find the decomposition which minimizes the total patch-wise sparseness of the decomposed modes. We propose a domain-decomposition type method, called intrinsic sparse mode decomposition (ISMD), which follows the "local-modes-construction + patching-up" procedure. The key step in the ISMD is to construct local pieces of the intrinsic sparse modes by a joint diagonalization problem. Thereafter a pivoted Cholesky decomposition is utilized to glue these local pieces together. Optimal sparse decomposition, consistency with different domain decomposition and robustness to small perturbation are proved under the so called regular-sparse assumption (see Definition 1.2). We provide simulation results to show the efficiency and robustness of the ISMD. We also compare the ISMD to other existing methods, e.g., eigen decomposition, pivoted Cholesky decomposition and convex relaxation of sparse principal component analysis [25] and [40].

Haodi Jiang, Qin Li, Yan Xu et al.

Obtaining high-quality magnetic and velocity fields through Stokes inversion is crucial in solar physics. In this paper, we present a new deep learning method, named Stacked Deep Neural Networks (SDNN), for inferring line-of-sight (LOS) velocities and Doppler widths from Stokes profiles collected by the Near InfraRed Imaging Spectropolarimeter (NIRIS) on the 1.6 m Goode Solar Telescope (GST) at the Big Bear Solar Observatory (BBSO). The training data of SDNN is prepared by a Milne-Eddington (ME) inversion code used by BBSO. We quantitatively assess SDNN, comparing its inversion results with those obtained by the ME inversion code and related machine learning (ML) algorithms such as multiple support vector regression, multilayer perceptrons and a pixel-level convolutional neural network. Major findings from our experimental study are summarized as follows. First, the SDNN-inferred LOS velocities are highly correlated to the ME-calculated ones with the Pearson product-moment correlation coefficient being close to 0.9 on average. Second, SDNN is faster, while producing smoother and cleaner LOS velocity and Doppler width maps, than the ME inversion code. Third, the maps produced by SDNN are closer to ME's maps than those from the related ML algorithms, demonstrating the better learning capability of SDNN than the ML algorithms. Finally, comparison between the inversion results of ME and SDNN based on GST/NIRIS and those from the Helioseismic and Magnetic Imager on board the Solar Dynamics Observatory in flare-prolific active region NOAA 12673 is presented. We also discuss extensions of SDNN for inferring vector magnetic fields with empirical evaluation.

Qin Li, Lorenzo Pareschi

We consider the development of exponential methods for the robust time discretization of space inhomogeneous Boltzmann equations in stiff regimes. Compared to the space homogeneous case, or more in general to the case of splitting based methods, studied in Dimarco Pareschi (SIAM J. Num. Anal. 2011) a major difficulty is that the local Maxwellian equilibrium state is not constant in a time step and thus needs a proper numerical treatment. We show how to derive asymptotic preserving (AP) schemes of arbitrary order and in particular using the Shu-Osher representation of Runge-Kutta methods we explore the monotonicity properties of such schemes, like strong stability preserving (SSP) and positivity preserving. Several numerical results confirm our analysis.

Matthew J. Colbrook, Qin Li, Ryan V. Raut et al.

Koopman operators linearize nonlinear dynamical systems, making their spectral information of crucial interest. Numerous algorithms have been developed to approximate these spectral properties, and Dynamic Mode Decomposition (DMD) stands out as the poster child of projection-based methods. Although the Koopman operator itself is linear, the fact that it acts in an infinite-dimensional space of observables poses challenges. These include spurious modes, essential spectra, and the verification of Koopman mode decompositions. While recent work has addressed these challenges for deterministic systems, there remains a notable gap in verified DMD methods for stochastic systems, where the Koopman operator measures the expectation of observables. We show that it is necessary to go beyond expectations to address these issues. By incorporating variance into the Koopman framework, we address these challenges. Through an additional DMD-type matrix, we approximate the sum of a squared residual and a variance term, each of which can be approximated individually using batched snapshot data. This allows verified computation of the spectral properties of stochastic Koopman operators, controlling the projection error. We also introduce the concept of variance-pseudospectra to gauge statistical coherency. Finally, we present a suite of convergence results for the spectral information of stochastic Koopman operators. Our study concludes with practical applications using both simulated and experimental data. In neural recordings from awake mice, we demonstrate how variance-pseudospectra can reveal physiologically significant information unavailable to standard expectation-based dynamical models.

Yewei Xu, Shi Chen, Qin Li

Does the use of auto-differentiation yield reasonable updates for deep neural networks (DNNs)? Specifically, when DNNs are designed to adhere to neural ODE architectures, can we trust the gradients provided by auto-differentiation? Through mathematical analysis and numerical evidence, we demonstrate that when neural networks employ high-order methods, such as Linear Multistep Methods (LMM) or Explicit Runge-Kutta Methods (ERK), to approximate the underlying ODE flows, brute-force auto-differentiation often introduces artificial oscillations in the gradients that prevent convergence. In the case of Leapfrog and 2-stage ERK, we propose simple post-processing techniques that effectively eliminates these oscillations, correct the gradient computation and thus returns the accurate updates.

Qin Li, Jianfeng Lu

We design a numerical scheme for transport equations with oscillatory periodic scattering coefficients. The scheme is asymptotic preserving in the diffusion limit as Knudsen number goes to zero. It also captures the homogenization limit as the length scale of the scattering coefficient goes to zero. The proposed method is based on the construction of multiscale finite element basis and a Galerkin projection based on the even-odd decomposition. The method is analyzed in the asymptotic regime, as well as validated numerically.

Qin Li, Xuefeng Liu

The paper develops an explicit a priori error estimate for finite element solution to nonhomogeneous Neumann problems. For this purpose, the hypercircle equation over finite element spaces is constructed and the explicit upper bound of the constant in the trace theorem is given. Numerical examples are shown in the final section, which implies the proposed error {estimate} has the convergence rate as $0.5$.

Kathrin Hellmuth, Christian Klingenberg, Qin Li

The quality of numerical reconstructions for unknown parameters in inverse problems depends fundamentally on the selection of experimental data. To ensure a robust reconstruction, it is crucial to select data that are sensitive to the parameters, a property typically characterized by the conditioning of the Fisher Information Matrix (FIM). In this work, we propose a general framework for an efficient down-sampling strategy that selects experimental setups that preserves the information content of the full-data FIM. Our approach leverages matrix sketching techniques from randomized numerical linear algebra to achieve a sensitivity-preserving approximation. The method involves drawing samples from a sensitivity-informed distribution, which we execute using gradient-free ensemble sampling methods to handle potentially non-smooth or discrete design spaces. Numerical experiments demonstrate the effectiveness of this framework in selecting optimal sensor locations for a Schroedinger potential reconstruction problem.

Borong Zhang, Leonardo Zepeda-Núñez, Qin Li

This paper introduces a novel deep neural network architecture for solving the inverse scattering problem in frequency domain with wide-band data, by directly approximating the inverse map, thus avoiding the expensive optimization loop of classical methods. The architecture is motivated by the filtered back-projection formula in the full aperture regime and with homogeneous background, and it leverages the underlying equivariance of the problem and compressibility of the integral operator. This drastically reduces the number of training parameters, and therefore the computational and sample complexity of the method. In particular, we obtain an architecture whose number of parameters scale sub-linearly with respect to the dimension of the inputs, while its inference complexity scales super-linearly but with very small constants. We provide several numerical tests that show that the current approach results in better reconstruction than optimization-based techniques such as full-waveform inversion, but at a fraction of the cost while being competitive with state-of-the-art machine learning methods.

Shi Chen, Qin Li, Oliver Tse et al.

The acceleration of gradient-based optimization methods is a subject of significant practical and theoretical importance, particularly within machine learning applications. While much attention has been directed towards optimizing within Euclidean space, the need to optimize over spaces of probability measures in machine learning motivates exploration of accelerated gradient methods in this context too. To this end, we introduce a Hamiltonian-flow approach analogous to momentum-based approaches in Euclidean space. We demonstrate that, in the continuous-time setting, algorithms based on this approach can achieve convergence rates of arbitrarily high order. We complement our findings with numerical examples.

Zhiyan Ding, Shi Chen, Qin Li et al.

Finding the optimal configuration of parameters in ResNet is a nonconvex minimization problem, but first-order methods nevertheless find the global optimum in the overparameterized regime. We study this phenomenon with mean-field analysis, by translating the training process of ResNet to a gradient-flow partial differential equation (PDE) and examining the convergence properties of this limiting process. The activation function is assumed to be $2$-homogeneous or partially $1$-homogeneous; the regularized ReLU satisfies the latter condition. We show that if the ResNet is sufficiently large, with depth and width depending algebraically on the accuracy and confidence levels, first-order optimization methods can find global minimizers that fit the training data.

Zhiyan Ding, Shi Chen, Qin Li et al.

Finding parameters in a deep neural network (NN) that fit training data is a nonconvex optimization problem, but a basic first-order optimization method (gradient descent) finds a global optimizer with perfect fit (zero-loss) in many practical situations. We examine this phenomenon for the case of Residual Neural Networks (ResNet) with smooth activation functions in a limiting regime in which both the number of layers (depth) and the number of weights in each layer (width) go to infinity. First, we use a mean-field-limit argument to prove that the gradient descent for parameter training becomes a gradient flow for a probability distribution that is characterized by a partial differential equation (PDE) in the large-NN limit. Next, we show that under certain assumptions, the solution to the PDE converges in the training time to a zero-loss solution. Together, these results suggest that the training of the ResNet gives a near-zero loss if the ResNet is large enough. We give estimates of the depth and width needed to reduce the loss below a given threshold, with high probability.

Zhiyan Ding, Qin Li

The classical Langevin Monte Carlo method looks for samples from a target distribution by descending the samples along the gradient of the target distribution. The method enjoys a fast convergence rate. However, the numerical cost is sometimes high because each iteration requires the computation of a gradient. One approach to eliminate the gradient computation is to employ the concept of ``ensemble." A large number of particles are evolved together so the neighboring particles provide gradient information to each other. In this article, we discuss two algorithms that integrate the ensemble feature into LMC and the associated properties. In particular, we find that if one directly surrogates the gradient using the ensemble approximation, the algorithm, termed Ensemble Langevin Monte Carlo, is unstable due to a high variance term. If the gradients are replaced by the ensemble approximations only in a constrained manner, to protect from the unstable points, the algorithm, termed Constrained Ensemble Langevin Monte Carlo, resembles the classical LMC up to an ensemble error but removes most of the gradient computation.

Zhiyan Ding, Qin Li, Jianfeng Lu et al.

The Underdamped Langevin Monte Carlo (ULMC) is a popular Markov chain Monte Carlo sampling method. It requires the computation of the full gradient of the log-density at each iteration, an expensive operation if the dimension of the problem is high. We propose a sampling method called Random Coordinate ULMC (RC-ULMC), which selects a single coordinate at each iteration to be updated and leaves the other coordinates untouched. We investigate the computational complexity of RC-ULMC and compare it with the classical ULMC for strongly log-concave probability distributions. We show that RC-ULMC is always cheaper than the classical ULMC, with a significant cost reduction when the problem is highly skewed and high dimensional. Our complexity bound for RC-ULMC is also tight in terms of dimension dependence.

Zhiyan Ding, Qin Li, Jianfeng Lu et al.

Langevin Monte Carlo (LMC) is a popular Markov chain Monte Carlo sampling method. One drawback is that it requires the computation of the full gradient at each iteration, an expensive operation if the dimension of the problem is high. We propose a new sampling method: Random Coordinate LMC (RC-LMC). At each iteration, a single coordinate is randomly selected to be updated by a multiple of the partial derivative along this direction plus noise, and all other coordinates remain untouched. We investigate the total complexity of RC-LMC and compare it with the classical LMC for log-concave probability distributions. When the gradient of the log-density is Lipschitz, RC-LMC is less expensive than the classical LMC if the log-density is highly skewed for high dimensional problems, and when both the gradient and the Hessian of the log-density are Lipschitz, RC-LMC is always cheaper than the classical LMC, by a factor proportional to the square root of the problem dimension. In the latter case, our estimate of complexity is sharp with respect to the dimension.

Zhiyan Ding, Qin Li

Langevin Monte Carlo (LMC) is a popular Bayesian sampling method. For the log-concave distribution function, the method converges exponentially fast, up to a controllable discretization error. However, the method requires the evaluation of a full gradient in each iteration, and for a problem on $\mathbb{R}^d$, this amounts to $d$ times partial derivative evaluations per iteration. The cost is high when $d\gg1$. In this paper, we investigate how to enhance computational efficiency through the application of RCD (random coordinate descent) on LMC. There are two sides of the theory: 1 By blindly applying RCD to LMC, one surrogates the full gradient by a randomly selected directional derivative per iteration. Although the cost is reduced per iteration, the total number of iteration is increased to achieve a preset error tolerance. Ultimately there is no computational gain; 2 We then incorporate variance reduction techniques, such as SAGA (stochastic average gradient) and SVRG (stochastic variance reduced gradient), into RCD-LMC. It will be proved that the cost is reduced compared with the classical LMC, and in the underdamped case, convergence is achieved with the same number of iterations, while each iteration requires merely one-directional derivative. This means we obtain the best possible computational cost in the underdamped-LMC framework.

Zhiyan Ding, Qin Li

Sampling from a log-concave distribution function is one core problem that has wide applications in Bayesian statistics and machine learning. While most gradient free methods have slow convergence rate, the Langevin Monte Carlo (LMC) that provides fast convergence requires the computation of gradients. In practice one uses finite-differencing approximations as surrogates, and the method is expensive in high-dimensions. A natural strategy to reduce computational cost in each iteration is to utilize random gradient approximations, such as random coordinate descent (RCD) or simultaneous perturbation stochastic approximation (SPSA). We show by a counter-example that blindly applying RCD does not achieve the goal in the most general setting. The high variance induced by the randomness means a larger number of iterations are needed, and this balances out the saving in each iteration. We then introduce a new variance reduction approach, termed Randomized Coordinates Averaging Descent (RCAD), and incorporate it with both overdamped and underdamped LMC. The methods are termed RCAD-O-LMC and RCAD-U-LMC respectively. The methods still sit in the random gradient approximation framework, and thus the computational cost in each iteration is low. However, by employing RCAD, the variance is reduced, so the methods converge within the same number of iterations as the classical overdamped and underdamped LMC. This leads to a computational saving overall.

Zhiyan Ding, Yiding Chen, Qin Li et al.

Stochastic Gradient Descent (SGD) plays a central role in modern machine learning. While there is extensive work on providing error upper bound for SGD, not much is known about SGD error lower bound. In this paper, we study the convergence of constant step-size SGD. We provide error lower bound of SGD for potentially non-convex objective functions with Lipschitz gradients. To our knowledge, this is the first analysis for SGD error lower bound without the strong convexity assumption. We use experiments to illustrate our theoretical results.

Eric Chung, Yalchin Efendiev, Yanbo Li et al.

The Boltzmann equation, as a model equation in statistical mechanics, is used to describe the statistical behavior of a large number of particles driven by the same physics laws. Depending on the media and the particles to be modeled, the equation has slightly different forms. In this article, we investigate a model Boltzmann equation with highly oscillatory media in the small Knudsen number regime, and study the numerical behavior of the Generalized Multi-scale Finite Element Method (GMsFEM) in the fluid regime when high oscillation in the media presents. The Generalized Multi-scale Finite Element Method (GMsFEM) is a general approach to numerically treat equations with multi-scale structures. The method is divided into the offline and online steps. In the offline step, basis functions are prepared from a snapshot space via a well-designed generalized eigenvalue problem (GEP), and these basis functions are then utilized to patch up for a solution through DG formulation in the online step to incorporate specific boundary and source information. We prove the wellposedness of the method on the Boltzmann equation, and show that the GEP formulation provides a set of optimal basis functions that achieve spectral convergence. Such convergence is independent of the oscillation in the media, or the smallness of the Knudsen number, making it one of the few methods that simultaneously achieve numerical homogenization and asymptotic preserving properties across all scales of oscillations and the Knudsen number.

Ke Chen, Qin Li, Li Wang

We consider the inverse problem of reconstructing the scattering and absorption coefficients using boundary measurements for a time dependent radiative transfer equation (RTE). As the measurement is mostly polluted by errors, both experimental and computational, an important question is to quantify how the error is amplified in the process of reconstruction. In the forward setting, the solution to the RTE behaves differently in different regimes, and the stability of the inverse problem vary accordingly. In particular, we consider two scalings in this paper. The first one concerns with a diffusive scaling whose macroscopic limit is a diffusion equation. In this case, we showed, following the similar approach as in [Chen, Li and Wang, arXiv:1703.00097], that the stability degrades when the limit is taken. The second one considers a highly forward peaked scattering, wherein the scattering operator is approximated by a Fokker-Planck operator as a limit. In this case, we showed that a fully recover of the scattering coefficient is less possible in the limit, whereas obtaining a rescaled version of the scattering coefficient becomes more practice friendly.

Qin Li, Ruiwen Shu, Li Wang

The inverse radiative transfer problem finds broad applications in medical imaging, atmospheric science, astronomy, and many other areas. This problem intends to recover the optical properties, denoted as absorption and scattering coefficient of the media, through the source-measurement pairs. A typical computational approach is to form the inverse problem as a PDE-constraint optimization, with the minimizer being the to-be-recovered coefficients. The method is tested to be efficient in practice, but lacks analytical justification: there is no guarantee of the existence or uniqueness of the minimizer, and the error is hard to quantify. In this paper, we provide a different algorithm by levering the ideas from singular decomposition analysis. Our approach is to decompose the measurements into three components, two out of which encode the information of the two coefficients respectively. We then split the optimization problem into two subproblems and use those two components to recover the absorption and scattering coefficients separately. In this regard, we prove the well-posedness of the new optimization, and the error could be quantified with better precision. In the end, we incorporate the diffusive scaling and show that the error is harder to control in the diffusive limit.

Qin Li, Jianfeng Lu, Weiran Sun

We study half-space linear kinetic equations with general boundary conditions that consist of both given incoming data and various type of reflections, extending our previous work [LLS14] on half-space equations with incoming boundary conditions. As in [LLS14], the main technique is a damping adding-removing procedure. We establish the well-posedness of linear (or linearized) half-space equations with general boundary conditions and quasi-optimality of the numerical scheme. The numerical method is validated by examples including a two-species transport equation, a multi-frequency transport equation, and the linearized BGK equation in 2D velocity space.