Search

A mean-field game (MFG) seeks the Nash Equilibrium of a game involving a continuum of players, where the Nash Equilibrium corresponds to a fixed point of the best-response mapping. However, simple fixed-point iterations do not always guarantee convergence. Fictitious play is an iterative algorithm that leverages a best-response mapping combined with a weighted average. Through a thorough study of the best-response mapping, this paper develops a simple and unified convergence analysis, providing the first explicit convergence rate for the fictitious play algorithm in MFGs of general types, especially non-potential MFGs. We demonstrate that the convergence and rate can be controlled through the weighting parameter in the algorithm, with linear convergence achievable under a general assumption. Building on this analysis, we propose two strategies to accelerate fictitious play. The first uses a backtracking line search to optimize the weighting parameter, while the second employs a hierarchical grid strategy to enhance stability and computational efficiency. We demonstrate the effectiveness of these acceleration techniques and validate our convergence rate analysis with various numerical examples.

In this work, we propose a balanced multi-component and multi-layer neural network (MMNN) structure to approximate functions with complex features with both accuracy and efficiency in terms of degrees of freedom and computation cost. The main idea is motivated by a multi-component, each of which can be approximated effectively by a single-layer network, and multi-layer decomposition in a "divide-and-conquer" type of strategy to deal with a complex function. While an easy modification to fully connected neural networks (FCNNs) or multi-layer perceptrons (MLPs) through the introduction of balanced multi-component structures in the network, MMNNs achieve a significant reduction of training parameters, a much more efficient training process, and a much improved accuracy compared to FCNNs or MLPs. Extensive numerical experiments are presented to illustrate the effectiveness of MMNNs in approximating high oscillatory functions and its automatic adaptivity in capturing localized features., Comment: Our codes and implementation details are available at https://github.com/ShijunZhangMath/MMNN

This paper explores the expressive power of deep neural networks for a diverse range of activation functions. An activation function set $\mathscr{A}$ is defined to encompass the majority of commonly used activation functions, such as $\mathtt{ReLU}$, $\mathtt{LeakyReLU}$, $\mathtt{ReLU}^2$, $\mathtt{ELU}$, $\mathtt{CELU}$, $\mathtt{SELU}$, $\mathtt{Softplus}$, $\mathtt{GELU}$, $\mathtt{SiLU}$, $\mathtt{Swish}$, $\mathtt{Mish}$, $\mathtt{Sigmoid}$, $\mathtt{Tanh}$, $\mathtt{Arctan}$, $\mathtt{Softsign}$, $\mathtt{dSiLU}$, and $\mathtt{SRS}$. We demonstrate that for any activation function $\varrho\in \mathscr{A}$, a $\mathtt{ReLU}$ network of width $N$ and depth $L$ can be approximated to arbitrary precision by a $\varrho$-activated network of width $3N$ and depth $2L$ on any bounded set. This finding enables the extension of most approximation results achieved with $\mathtt{ReLU}$ networks to a wide variety of other activation functions, albeit with slightly increased constants. Significantly, we establish that the (width,$\,$depth) scaling factors can be further reduced from $(3,2)$ to $(1,1)$ if $\varrho$ falls within a specific subset of $\mathscr{A}$. This subset includes activation functions such as $\mathtt{ELU}$, $\mathtt{CELU}$, $\mathtt{SELU}$, $\mathtt{Softplus}$, $\mathtt{GELU}$, $\mathtt{SiLU}$, $\mathtt{Swish}$, and $\mathtt{Mish}$.

In this work, a comprehensive numerical study involving analysis and experiments shows why a two-layer neural network has difficulties handling high frequencies in approximation and learning when machine precision and computation cost are important factors in real practice. In particular, the following basic computational issues are investigated: (1) the minimal numerical error one can achieve given a finite machine precision, (2) the computation cost to achieve a given accuracy, and (3) stability with respect to perturbations. The key to the study is the conditioning of the representation and its learning dynamics. Explicit answers to the above questions with numerical verifications are presented.

This paper explores the expressive power of deep neural networks through the framework of function compositions. We demonstrate that the repeated compositions of a single fixed-size ReLU network exhibit surprising expressive power, despite the limited expressive capabilities of the individual network itself. Specifically, we prove by construction that $\mathcal{L}_2\circ \boldsymbol{g}^{\circ r}\circ \boldsymbol{\mathcal{L}}_1$ can approximate $1$-Lipschitz continuous functions on $[0,1]^d$ with an error $\mathcal{O}(r^{-1/d})$, where $\boldsymbol{g}$ is realized by a fixed-size ReLU network, $\boldsymbol{\mathcal{L}}_1$ and $\mathcal{L}_2$ are two affine linear maps matching the dimensions, and $\boldsymbol{g}^{\circ r}$ denotes the $r$-times composition of $\boldsymbol{g}$. Furthermore, we extend such a result to generic continuous functions on $[0,1]^d$ with the approximation error characterized by the modulus of continuity. Our results reveal that a continuous-depth network generated via a dynamical system has immense approximation power even if its dynamics function is time-independent and realized by a fixed-size ReLU network.

Linear evolution PDE $\partial_t u(x,t) = -\mathcal{L} u$, where $\mathcal{L}$ is a strongly elliptic operator independent of time, is studied as an example to show if one can superpose snapshots of a single (or a finite number of) solution(s) to construct an arbitrary solution. Our study shows that it depends on the growth rate of the eigenvalues, $\mu_n$, of $\mathcal{L}$ in terms of $n$. When the statement is true, a simple data-driven approach for model reduction and approximation of an arbitrary solution of a PDE without knowing the underlying PDE is designed. Numerical experiments are presented to corroborate our analysis.

In this work we study the problem about learning a partial differential equation (PDE) from its solution data. PDEs of various types are used as examples to illustrate how much the solution data can reveal the PDE operator depending on the underlying operator and initial data. A data driven and data adaptive approach based on local regression and global consistency is proposed for stable PDE identification. Numerical experiments are provided to verify our analysis and demonstrate the performance of the proposed algorithms., Comment: 44 pages

In this paper, we consider a Bayesian inverse problem modeled by elliptic partial differential equations (PDEs). Specifically, we propose a data-driven and model-based approach to accelerate the Hamiltonian Monte Carlo (HMC) method in solving large-scale Bayesian inverse problems. The key idea is to exploit (model-based) and construct (data-based) the intrinsic approximate low-dimensional structure of the underlying problem which consists of two components - a training component that computes a set of data-driven basis to achieve significant dimension reduction in the solution space, and a fast solving component that computes the solution and its derivatives for a newly sampled elliptic PDE with the constructed data-driven basis. Hence we achieve an effective data and model-based approach for the Bayesian inverse problem and overcome the typical computational bottleneck of HMC - repeated evaluation of the Hamiltonian involving the solution (and its derivatives) modeled by a complex system, a multiscale elliptic PDE in our case. We present numerical examples to demonstrate the accuracy and efficiency of the proposed method.

The photoacoustic tomography (PAT) is a hybrid modality that combines the optics and acoustics to obtain high resolution and high contrast imaging of heterogeneous media. In this work, our objective is to study the inverse problem in the quantitative step of PAT which aims to reconstruct the optical coefficients of the governing radiative transport equation from the ultrasound measurements. In our analysis, we take the simplified $P_N$ approximation of the radiative transport equation as the physical model and then show the uniqueness and stability for this modified inverse problem. Numerical simulations based on synthetic data are presented to validate our analysis., Comment: 32 pages

In this work, a simple and efficient dual iterative refinement (DIR) method is proposed for dense correspondence between two nearly isometric shapes. The key idea is to use dual information, such as spatial and spectral, or local and global features, in a complementary and effective way, and extract more accurate information from current iteration to use for the next iteration. In each DIR iteration, starting from current correspondence, a zoom-in process at each point is used to select well matched anchor pairs by a local mapping distortion criterion. These selected anchor pairs are then used to align spectral features (or other appropriate global features) whose dimension adaptively matches the capacity of the selected anchor pairs. Thanks to the effective combination of complementary information in a data-adaptive way, DIR is not only efficient but also robust to render accurate results within a few iterations. By choosing appropriate dual features, DIR has the flexibility to handle patch and partial matching as well. Extensive experiments on various data sets demonstrate the superiority of DIR over other state-of-the-art methods in terms of both accuracy and efficiency., Comment: 9 pages, 9 figures and 1 table

In this work, we introduce a novel local pairwise descriptor and then develop a simple, effective iterative method to solve the resulting quadratic assignment through sparsity control for shape correspondence between two approximate isometric surfaces. Our pairwise descriptor is based on the stiffness and mass matrix of finite element approximation of the Laplace-Beltrami differential operator, which is local in space, sparse to represent, and extremely easy to compute while containing global information. It allows us to deal with open surfaces, partial matching, and topological perturbations robustly. To solve the resulting quadratic assignment problem efficiently, the two key ideas of our iterative algorithm are: 1) select pairs with good (approximate) correspondence as anchor points, 2) solve a regularized quadratic assignment problem only in the neighborhood of selected anchor points through sparsity control. These two ingredients can improve and increase the number of anchor points quickly while reducing the computation cost in each quadratic assignment iteration significantly. With enough high-quality anchor points, one may use various pointwise global features with reference to these anchor points to further improve the dense shape correspondence. We use various experiments to show the efficiency, quality, and versatility of our method on large data sets, patches, and point clouds (without global meshes)., Comment: 8 pages, 6 figures. Compared to the version to be published in the 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Figure 1 has been changed to a more illustrative example and run time table 1 has been updated by our recently optimized code

In this work, we introduce a fast numerical algorithm to solve the time-dependent radiative transport equation (RTE). Our method uses the integral formulation of RTE and applies the treecode algorithm to reduce the computational complexity from O(M^{2+1/d} ) to O(M^{1+1/d} log M ), where M is the number of points in the physical domain. The error analysis is presented and numerical experiments are performed to validate our algorithm.

We prove the Marchenko-Pastur law for the eigenvalues of $p \times p$ sample covariance matrices in two new situations where the data does not have independent coordinates. In the first scenario - the block-independent model - the $p$ coordinates of the data are partitioned into blocks in such a way that the entries in different blocks are independent, but the entries from the same block may be dependent. In the second scenario - the random tensor model - the data is the homogeneous random tensor of order $d$, i.e. the coordinates of the data are all $\binom{n}{d}$ different products of $d$ variables chosen from a set of $n$ independent random variables. We show that Marchenko-Pastur law holds for the block-independent model as long as the size of the largest block is $o(p)$ and for the random tensor model as long as $d = o(n^{1/3})$. Our main technical tools are new concentration inequalities for quadratic forms in random variables with block-independent coordinates, and for random tensors., Comment: 23 pages, 4 figures

Phase retrieval with prior information can be cast as a nonsmooth and nonconvex optimization problem. We solve the problem by graph projection splitting (GPS), where the two proximity subproblems and the graph projection step can be solved efficiently. With slight modification, we also propose a robust graph projection splitting (RGPS) method to stabilize the iteration for noisy measurements. Contrary to intuition, RGPS outperforms GPS with fewer iterations to locate a satisfying solution even for noiseless case. Based on the connection between GPS and Douglas-Rachford iteration, under mild conditions on the sampling vectors, we analyze the fixed point sets and provide the local convergence of GPS and RGPS applied to noiseless phase retrieval without prior information. For noisy case, we provide the error bound of the reconstruction. Compared to other existing methods, thanks for the splitting approach, GPS and RGPS can efficiently solve phase retrieval with prior information regularization for general sampling vectors which are not necessarily isometric. For Gaussian phase retrieval, compared to existing gradient flow approaches, numerical results show that GPS and RGPS are much less sensitive to the initialization. Thus they markedly improve the phase transition in noiseless case and reconstruction in the presence of noise respectively. GPS shows sharpest phase transition among existing methods including RGPS, while it needs more iterations than RGPS when the number of measurement is large enough. RGPS outperforms GPS in terms of stability for noisy measurements. When applying RGPS to more general non-Gaussian measurements with prior information, such as support, sparsity and TV minimization, RGPS either outperforms state-of-the-art solvers or can be combined with state-of-the-art solvers to improve their reconstruction quality.

We study in this work an integral formulation for the radiative transfer equation (RTE) in anisotropic media with truncated approximation to the scattering phase function. The integral formulation consists of a coupled system of integral equations for the angular moments of the transport solution. We analyze the approximate separability of the kernel functions in these integral formulations, deriving asymptotic lower and upper bounds on the number of terms needed in a separable approximation of the kernel functions as the moment grows. Our analysis provides the mathematical understanding on when low-rank approximations to the discretized integral kernels can be used to develop fast numerical algorithms for the corresponding system of integral equations.

We propose a data-driven approach to solve multiscale elliptic PDEs with random coefficients based on the intrinsic low dimension structure of the underlying elliptic differential operators. Our method consists of offline and online stages. At the offline stage, a low dimension space and its basis are extracted from the data to achieve significant dimension reduction in the solution space. At the online stage, the extracted basis will be used to solve a new multiscale elliptic PDE efficiently. The existence of low dimension structure is established by showing the high separability of the underlying Green's functions. Different online construction methods are proposed depending on the problem setup. We provide error analysis based on the sampling error and the truncation threshold in building the data-driven basis. Finally, we present numerical examples to demonstrate the accuracy and efficiency of the proposed method.

In this work, we study the instability of an inverse problem of radiative transport equation with angularly averaged measurement near the diffusion limit, i.e. the normalized mean free path (the Knudsen number) $0 < \eps \ll 1$. It is well-known that there is a transition of stability from H\"{o}lder type to logarithmic type with $\eps\to 0$, the theory of this transition of stability is still an open problem. In this study, we show the transition of stability by establishing the balance of two different regimes depending on the relative sizes of $\eps$ and the perturbation in measurements. When $\eps$ is sufficiently small, we obtain exponential instability, which stands for the diffusive regime, and otherwise we obtain H\"{o}lder instability instead, which stands for the transport regime., Comment: 20 pages

We aim to find a solution $\bm{x}\in\mathbb{C}^n$ to a system of quadratic equations of the form $b_i=\lvert\bm{a}_i^*\bm{x}\rvert^2$, $i=1,2,\ldots,m$, e.g., the well-known NP-hard phase retrieval problem. As opposed to recently proposed state-of-the-art nonconvex methods, we revert to the semidefinite relaxation (SDR) PhaseLift convex formulation and propose a successive and incremental nonconvex optimization algorithm, termed as \texttt{IncrePR}, to indirectly minimize the resulting convex problem on the cone of positive semidefinite matrices. Our proposed method overcomes the excessive computational cost of typical SDP solvers as well as the need of a good initialization for typical nonconvex methods. For Gaussian measurements, which is usually needed for provable convergence of nonconvex methods, \texttt{IncrePR} with restart strategy outperforms state-of-the-art nonconvex solvers with a sharper phase transition of perfect recovery and typical convex solvers in terms of computational cost and storage. For more challenging structured (non-Gaussian) measurements often occurred in real applications, such as transmission matrix and oversampling Fourier transform, \texttt{IncrePR} with several restarts can be used to find a good initial guess. With further refinement by local nonconvex solvers, one can achieve a better solution than that obtained by applying nonconvex solvers directly when the number of measurements is relatively small. Extensive numerical tests are performed to demonstrate the effectiveness of the proposed method., Comment: 20 pages, 25 figures

Random fields are commonly used for modeling of spatially (or timely) dependent stochastic processes. In this study, we provide a characterization of the intrinsic complexity of a random field in terms of its second order statistics, e.g., the covariance function, based on the Karhumen-Lo\'{e}ve expansion. We then show scaling laws for the intrinsic complexity of a random field in terms of the correlation length as it goes to 0. In the discrete setting, it becomes approximate embeddings of a set of random vectors. We provide a precise scaling law when the random vectors have independent and identically distributed entires using random matrix theory as well as when the random vectors has a specific covariance structure., Comment: 26 pages

We present a hybrid imaging method for a challenging travel time tomography problem which includes both unknown medium and unknown scatterers in a bounded domain. The goal is to recover both the medium and the boundary of the scatterers from the scattering relation data on the domain boundary. Our method is composed of three steps: 1) preprocess the data to classify them into three different categories of measurements corresponding to non-broken rays, broken-once rays, and others, respectively, 2) use the the non-broken ray data and an effective data-driven layer stripping strategy--an optimization based iterative imaging method--to recover the medium velocity outside the convex hull of the scatterers, and 3) use selected broken-once ray data to recover the boundary of the scatterers--a direct imaging method. By numerical tests, we show that our hybrid method can recover both the unknown medium and the not-too-concave scatterers efficiently and robustly., Comment: 28 pages

We propose a hybrid approach to solve the high-frequency Helmholtz equation with point source terms in smooth heterogeneous media. The method is based on the ray-based finite element method (ray-FEM), whose original version can not handle the singularity close to point sources accurately. This pitfall is addressed by combining the ray-FEM, which is used to compute the smooth far-field of the solution accurately, with a high-order asymptotic expansion close to the point source, which is used to properly capture the singularity of the solution in the near-field. The method requires a fixed number of grid points per wavelength to accurately represent the wave field with an asymptotic convergence rate of $\mathcal{O}(\omega^{-1/2})$, where $\omega$ is the frequency parameter in the Helmholtz equation. In addition, a fast sweeping-type preconditioner is used to solve the resulting linear system. We present numerical examples in 2D to show both accuracy and efficiency of our method as the frequency increases. In particular, we provide numerical evidence of the convergence rate, and we show empirically that the overall complexity is $\mathcal{O}(\omega^2)$ up to a poly-logarithmic factor.

We present a ray-based finite element method (ray-FEM) by learning basis adaptive to the underlying high-frequency Helmholtz equation in smooth media. Based on the geometric optics ansatz of the wave field, we learn local dominant ray directions by probing the medium using low-frequency waves with the same source. Once local ray directions are extracted, they are incorporated into the finite element basis to solve the high-frequency Helmholtz equation. This process can be continued to further improve approximations for both local ray directions and the high frequency wave field iteratively. The method requires a fixed number of grid points per wavelength to represent the wave field and achieves an asymptotic convergence as the frequency $\omega\rightarrow \infty$ without the pollution effect. A fast solver is developed for the resulting linear system with an empirical complexity $\mathcal{O}(\omega^d)$ up to a poly-logarithmic factor. Numerical examples in 2D are presented to corroborate the claims.

We propose a new and simple discretization, named the Modified Virtual Grid Difference (MVGD), for numerical approximation of the Laplace-Beltrami (LB) operator on manifolds sampled by point clouds. The key observation is that both the manifold and a function defined on it can both be parametrized in a local Cartesian coordinate system and approximated using least squares. Based on the above observation, we first introduce a local virtual grid with a scale adapted to the sampling density centered at each point. Then we propose a modified finite difference scheme on the virtual grid to discretize the LB operator. Instead of using the local least squares values on all virtual grid points like the typical finite difference method, we use the function value explicitly at the grid located at the center (coincided with the data point). The new discretization provides more diagonal dominance to the resulting linear system and improves its conditioning. We show that the linear system can be robustly, efficiently and accurately solved by existing fast solver such as the Algebraic Multigrid (AMG) method. We will present numerical tests and comparison with other exiting methods to demonstrate the effectiveness and the performance of the proposed approach.

This paper presents a fast iterative solver for Lippmann-Schwinger equation for high-frequency waves scattered by a smooth medium with a compactly supported inhomogeneity. The solver is based on the sparsifying preconditioner and a domain decomposition approach similar to the method of polarized traces. The iterative solver has two levels, the outer level in which a sparsifying preconditioner for the Lippmann-Schwinger equation is constructed, and the inner level, in which the resulting sparsified system is solved fast using an iterative solver preconditioned with a bi-directional matrix-free variant of the method of polarized traces. The complexity of the construction and application of the preconditioner is $\mathcal{O}(N)$ and $\mathcal{O}(N\log{N})$ respectively, where $N$ is the number of degrees of freedom. Numerical experiments in 2D indicate that the number of iterations in both levels depends weakly on the frequency resulting in method with an overall $\mathcal{O}(N\log{N})$ complexity., Comment: small modifications

Traditionally, the field of computational Bayesian statistics has been divided into two main subfields: variational methods and Markov chain Monte Carlo (MCMC). In recent years, however, several methods have been proposed based on combining variational Bayesian inference and MCMC simulation in order to improve their overall accuracy and computational efficiency. This marriage of fast evaluation and flexible approximation provides a promising means of designing scalable Bayesian inference methods. In this paper, we explore the possibility of incorporating variational approximation into a state-of-the-art MCMC method, Hamiltonian Monte Carlo (HMC), to reduce the required gradient computation in the simulation of Hamiltonian flow, which is the bottleneck for many applications of HMC in big data problems. To this end, we use a {\it free-form} approximation induced by a fast and flexible surrogate function based on single-hidden layer feedforward neural networks. The surrogate provides sufficiently accurate approximation while allowing for fast exploration of parameter space, resulting in an efficient approximate inference algorithm. We demonstrate the advantages of our method on both synthetic and real data problems.

For big data analysis, high computational cost for Bayesian methods often limits their applications in practice. In recent years, there have been many attempts to improve computational efficiency of Bayesian inference. Here we propose an efficient and scalable computational technique for a state-of-the-art Markov Chain Monte Carlo (MCMC) methods, namely, Hamiltonian Monte Carlo (HMC). The key idea is to explore and exploit the structure and regularity in parameter space for the underlying probabilistic model to construct an effective approximation of its geometric properties. To this end, we build a surrogate function to approximate the target distribution using properly chosen random bases and an efficient optimization process. The resulting method provides a flexible, scalable, and efficient sampling algorithm, which converges to the correct target distribution. We show that by choosing the basis functions and optimization process differently, our method can be related to other approaches for the construction of surrogate functions such as generalized additive models or Gaussian process models. Experiments based on simulated and real data show that our approach leads to substantially more efficient sampling algorithms compared to existing state-of-the art methods.

Markov Chain Monte Carlo (MCMC) algorithms play an important role in statistical inference problems dealing with intractable probability distributions. Recently, many MCMC algorithms such as Hamiltonian Monte Carlo (HMC) and Riemannian Manifold HMC have been proposed to provide distant proposals with high acceptance rate. These algorithms, however, tend to be computationally intensive which could limit their usefulness, especially for big data problems due to repetitive evaluations of functions and statistical quantities that depend on the data. This issue occurs in many statistic computing problems. In this paper, we propose a novel strategy that exploits smoothness (regularity) of parameter space to improve computational efficiency of MCMC algorithms. When evaluation of functions or statistical quantities are needed at a point in parameter space, interpolation from precomputed values or previous computed values is used. More specifically, we focus on Hamiltonian Monte Carlo (HMC) algorithms that use geometric information for faster exploration of probability distributions. Our proposed method is based on precomputing the required geometric information on a set of grids before running sampling information at nearby grids at each iteration of HMC. Sparse grid interpolation method is used for high dimensional problems. Tests on computational examples are shown to illustrate the advantages of our method.