一周活动预告（12.5-12.11）

本次报告聚焦于在求解大规模复杂系统中常用的手段之一，数学优化的建模，求解与软件开发。聚焦于数学规划开发中的一些经典问题，以及在一些大规模复杂管理系统中常用的线性，整数，非线性建模的思路与解决手段。

在这个报告中，我们将介绍一个求解Stefan自由边界问题的基于位势理论的直角网格方法。Stefan问题的自由边界会随时间推进产生极大变形，依赖表面张力系数的值，还可能长出非常复杂的指状图案。用贴体非结构网格方法求解Stefan问题的计算开销很大。这个报告中将提到的方法首先把包含不规则自由边界的矩形计算区域划分成不贴体的直角网格，不要求直角网格线与自由边界匹配。在求解整个Stefan问题的过程中，直角网格固定不变，做到了网格开销极小。在直角网格上离散Stefan问题应用基于位势理论的（涉及边界积分和体积分）方法进行离散、修正和结合快速算法进行求解。我们还将在这个报告中演示一些数值模拟结果。

There are many important practical optimization problems whose feasible regions are not known to be nonempty or not, and optimizers of the objective function with the least constraint violation prefer to be found. A natural way for dealing with these problems is to extend the nonlinear optimization problem as the one optimizing the objective function over the set of points with the least constraint violation. This leads to the study of the shifted problem. This report focuses on the constrained convex optimization problem. The sufficient condition for the closedness of the set of feasible shifts is presented and the continuity properties of the optimal value function and the solution mapping for the shifted problem are studied. Properties of the conjugate dual of the shifted problem are discussed through the relations between the dual function and the optimal value function. The solvability of the dual of the optimization problem with the least constraint violation is investigated. It is shown that, if the least violated shift is in the domain of the subdifferential of the optimal value function, then this dual problem has an unbounded solution set. Under this condition, the optimality conditions for the problem with the least constraint violation are established in term of the augmented Lagrangian. It is shown that the augmented Lagrangian method has the properties that the sequence of shifts converges to the least violated shift and the sequence of multipliers is unbounded. Moreover, it is proved that the augmented Lagrangian method is able to find an approximate solution to the problem with the least constraint violation and it has linear rate of convergence under an error bound condition. The augmented Lagrangian method is applied to an illustrative convex second-order cone constrained optimization problem with least violation constraint and numerical results verify the theoretical results obtained in this paper. This is a joint work with Professor Yu-Hong Dai.

While deep learning has achieved great success in practice, its theory is still mysterious to us. Experience shows that when adopting the regularization layer (e.g., Batch Normalization, Layer Normalization), the performance of neural networks in practice will be much better. From a theoretical point of view, the regularization layers make the neural networks scale invariant, which is a helpful property for the understanding of the optimization, especially the learning rate. With the convenience of the scale invariant and effective learning rate, we can partly explain some phenomena in deep learning.

In this talk, we will briefly introduce the scale invariant architectures, as well as the effective learning rate it brings out. We will further explain how it helps explain several phenomena in deep learning. Some works based on the effective learning rate itself will also be discussed.

For many years now, mathematicians have been using computers to help us compute. But right now we seem to be at the beginning of a new era where mathematicians can use computers to help us to reason. Professor Buzzard will speak about how neural networks and interactive theorem provers are giving mathematicians new ways to use computers in their research.

The talk will be suitable for a general mathematical audience; no background knowledge of either neural networks or interactive theorem provers will be assumed.

Time discretization is an important issue for time-dependent partial differential equations (PDEs). For the k-th ( $\mathrm{z} \geq$ 2) order PDEs, the explicit method may suffer from a severe time step restriction $\Delta t=O(\Delta x k)$ for stability. Implicit methods are generally unconditionally stable, however, they are cumbersome for nonlinear equations, since a nonlinear algebraic system must be solved (e.g. by Newton iteration) at each time step. The implicit-explicit (IMEX) methods, which treat the stiffer terms implicitly and the rest of the terms explicitly, can not only alleviate time step constraint, but also reduce the difficulty of solving the algebraic system especially when the stiffer terms are linear. We have analyzed the stability of various large time-stepping IMEX schemes for the high order PDEs such as the convection-dispersion equation in conjunction with high order finite difference method for spatial discretization.

Furthermore, for the equations with nonlinear high derivative terms, the IMEX methods are still too expensive to use. A better alternative is to use the explicit-implicit-null (EIN) method. The EIN method does not need any nonlinear iterative solver, and the severe time step restriction for explicit methods can be removed. Coupled with the EIN time-marching method, we have discussed the stability of high order finite difference and the LDG schemes for solving the convection-dispersion and the biharmonic type equations, respectively. Our main contribution is to show rigorously that the resulting numerical schemes are stable for large time step if stabilization terms of appropriate size are chosen. Numerical experiments are given to assess accuracy and stability.

This report will introduce a numerical discretization method based on tensor decomposition for solving partial differential equations. Based on this idea, a tensor neural network and its corresponding machine learning algorithm are introduced. The neural grid is built based on the form of tensor product, which can directly integrate high-dimensional functions, and convert high-dimensional integrals with exponential complexity into polynomial tensor integrals without the help of Monte Carlo process. Next, we use tensor neural network to design machine learning algorithms for solving high-dimensional partial differential equations and eigenvalue problems, hoping to bring more degrees of freedom and operability to the solution of high-dimensional partial differential equations. The application of tensor discrete method and machine learning algorithm of tensor neural network in solving high-dimensional eigenvalue problems and multibody problems will be introduced in the report.

For many decades, physics-based PDEs have been commonly employed for modeling complex system responses, then traditional numerical methods were employed to solve the PDEs and provide predictions. However, when governing laws are unknown or when high degrees of heterogeneity present, these classical models may become inaccurate. In this talk we propose to use data-driven modeling which directly utilizes high-fidelity simulation and experimental measurements to learn the hidden physics and provide further predictions. In particular, we develop PDE-inspired neural operator architectures, to learn the mapping between loading conditions and the corresponding system response. By parameterizing the increment between layers as an integral operator, our neural operator can be seen as the analog of a time-dependent nonlocal equation, which captures the long-range dependencies in the feature space and is guaranteed to be resolution-independent. Moreover, when applying to (hidden) PDE solving tasks, our neural operator provides a universal approximator to a fixed point iterative procedure, and partial physical knowledge can be incorporated to further improve the model’s generalizability and transferability. As an application, we learn the material models directly from digital image correlation (DIC) displacement tracking measurements on a porcine tricuspid valve leaflet tissue, and show that the learnt model substantially outperforms conventional constitutive models.

In this talk, we present superconvergence properties of the local discontinuous Galerkin (LDG) methods for solving nonlinear convection-diffusion equations in one space dimension. The main technicality is an elaborate estimate to terms involving projection errors. By introducing a new projection and constructing some correction functions, we prove the (2k+1)th order superconvergence for the cell averages and the numerical flux in the discrete L2 norm with polynomials of degree k≥1, no matter whether the flow direction f'(u) changes or not. Superconvergence of order k +2 (k +1) is obtained for the LDG error (its derivative) at interior right (left) Radau points, and the convergence ord er for the error derivative at Radau points can be improved to k+2 when the direction of the flow doesn't change. Finally, a supercloseness result of order k+2 towards a special Gauss-Radau projection of the exact solution is shown. The superconvergence analysis can be extended to the generalized numerical fluxes and the mixed boundary conditions. All theoretical findings are confirmed by numerical experiments.

Adaptive computation is of great importance in numerical simulations. The ideas for adaptive computations can be dated back to adaptive finite element methods in 1970s. In this talk, we shall first review some recent development for adaptive method with applications. Then, we shall propose a deep adaptive sampling method for solving PDEs where deep neural networks are utilized to approximate the solutions. In particular, we propose the failure informed PINNs (FI-PINNs), which can adaptively refine the training set with the goal of reducing the failure probability. Compared to the neural network approximation obtained with uniformly distributed collocation points, the developed algorithms can significantly improve the accuracy, especially for low regularity and high-dimensional problems.

Two novel orthogonalization-free algorithms are proposed to solve extreme eigenvalue problems. One proposed algorithm modifies the multi-column gradient such that earlier columns are decoupled from later ones. The other proposed algorithm adopts a novel objective function to achieve the goal. Global convergence to eigenvectors instead of eigenspace is guaranteed almost surely for both algorithms. Locally, algorithms converge linearly with a convergence rate depending on eigengaps. Momentum acceleration, exact linesearch, and column locking are incorporated to accelerate both algorithms further and reduce their computational costs. We demonstrate the efficiency of both algorithms on matrices from full configuration interaction problems and compute their ground state and low-lying excited states.

There is rich in unconventional oil and gas resources such as shale/tight and fractured-cavity carbonate rocks in China, which are realistic fields for increasing reserves and production in the future. Compared with foreign unconventional oil and gas reservoirs, the unconventional oil and gas reservoirs have various types of pores in the storage and seepage space in China, including matrix pores, complex fractures, and dissolution pores, which are highly heterogeneous and span multiple scales. In the middle stage of the process, cracks, pores and dissolution cavities are prone to deformation and have strong stress sensitivity. In this talk we will mainly introduce the numerical simulation method, simulator development and application of multi-scale and heat-fluid-solid coupling problems in the development of shale/tight, fractured-cavity carbonate rock unconventional reservoirs.

A hybrid data assimilation algorithm is developed for complex dynamical systems with partial observations. The method starts with applying a spectral decomposition to the entire spatiotemporal fields, followed by creating a machine learning model that builds a nonlinear map between the coefficients of observed and unobserved state variables for each spectral mode. A cheap low-order nonlinear stochastic parameterized extended Kalman filter (SPEKF) model is employed as the forecast model in the ensemble Kalman filter to deal with each mode associated with the observed variables. The resulting ensemble members are then fed into the machine learning model to create an ensemble of the corresponding unobserved variables. In addition to the ensemble spread, the training residual in the machine learning-induced nonlinear map is further incorporated into the state estimation that advances the quantification of the posterior uncertainty. The hybrid data assimilation algorithm is applied to a precipitation quasi-geostrophic (PQG) model, which includes the effects of water vapor, clouds, and rainfall beyond the classical two-layer QG model. The complicated nonlinearities in the PQG equations prevent traditional methods from building simple and accurate reduced-order forecast models. In contrast, the SPEKF model is skillful in recovering the intermittent observed states, and the machine learning model effectively estimates the chaotic unobserved signals. Utilizing the calibrated SPEKF and machine learning models under the situation of a moderate cloud fraction, the resulting hybrid data assimilation remains reasonably accurate when applied to other geophysical scenarios with nearly clear skies or relatively heavy rainfall, implying the robustness of the algorithm.

We consider the problem of testing whether a single coefficient is equal to zero in high-dimensional fixed-design linear models. In the high-dimensional setting where the dimension of covariates $p$ is allowed to be in the same order of magnitude as sample size $n$, to achieve finite-population validity, existing methods usually require strong distributional assumptions on the noise vector (such as Gaussian or rotationally invariant), which limits their applications in practice. In this paper, we propose a new method, called \emph{residual permutation test} (RPT), which is constructed by projecting the regression residuals onto the space orthogonal to the union of the column spaces of the original and permuted design matrices. RPT can be proved to achieve finite-population size validity under fixed design with just exchangeable noises, whenever $p < n / 2$. Moreover, RPT is shown to be asymptotically powerful for heavy tailed noises with bounded $(1+t)$-th order moment when the true coefficient is at least of order $n^{-t/(1+t)}$ for $t \in [0,1]$. We further proved that this signal size requirement is essentially optimal in the minimax sense. Numerical studies confirm that RPT performs well in a wide range of simulation settings with normal and heavy-tailed noise distributions.

Deep neural networks, as a powerful system to represent high dimensional complex functions, play a key role in deep learning. Convergence of deep neural networks is a fundamental issue in building the mathematical foundation for deep learning. We investigated the convergence of deep ReLU networks and deep convolutional neural networks in two recent papers, where only the Rectified Linear Unit (ReLU) activation function was studied. The important pooling strategy was not considered therein either. In this paper, we study the convergence of deep neural networks as the depth tends to infinity for general activation functions which cover most of commonly-used activation functions in artificial neural networks. Pooling will also be studied. Specifically, we adapt the linear method developed in our recent papers to prove that the major condition there is still sufficient for the neural networks defined by non-expansive activation functions to converge, despite their nonlinearity. For contractive activation functions such as the logistic sigmoid function, we establish a uniform and exponential convergence of the associated deep neural networks.

A linearly implicit renormalized lumped mass finite element method is considered for solving the equations describing heat flow of harmonic maps, of which the exact solution naturally satisfies the pointwise constraint $|\m|=1$. At every time level, the method first computes an auxiliary numerical solution by a linearly implicit lumped mass method and then renormalizes it at all finite element nodes before proceeding to the next time level. It is shown that such a renormalized finite element method has an error bound of $O(\tau+h^{r+1})$ for tensor-product finite elements of degree $r\ge1$. The proof of the error estimates is based on a geometric relation between the auxiliary and renormalized numerical solutions. The extension of the error analysis to triangular mesh is straightforward and discussed in the conclusion section.

In this talk, a thorough review of numerical considerations for fractional sub diffusion problems as well as its extension to generalized problems will be introduced first. Then, two class of numerical schemes for generalized problems are presented sequentially: (i) generalized interpolation based approach; (ii) generalized convolution quadrature based approach. After that, the convergence and stability of related numerical schemes are analyzed using the well-known energy method. In the end, some numerical results will be presented to demonstrate the accuracy of these schemes.