Convergence of Stochastic Gradient Methods for Wide Two-Layer Physics-Informed Neural Networks

TL;DR

This paper proves that stochastic gradient descent converges linearly when training over-parameterized two-layer physics-informed neural networks for PDEs, providing theoretical guarantees for stochastic training methods.

Contribution

It extends previous convergence results from gradient descent to stochastic gradient descent for PINNs, handling dynamic randomness in training.

Findings

Establishes linear convergence of SGD for over-parameterized PINNs

Handles general activation functions in the analysis

Provides high probability convergence guarantees

Abstract

Physics informed neural networks (PINNs) represent a very popular class of neural solvers for partial differential equations. In practice, one often employs stochastic gradient descent type algorithms to train the neural network. Therefore, the convergence guarantee of stochastic gradient descent is of fundamental importance. In this work, we establish the linear convergence of stochastic gradient descent / flow in training over-parameterized two layer PINNs for a general class of activation functions in the sense of high probability. These results extend the existing result [18] in which gradient descent was analyzed. The challenge of the analysis lies in handling the dynamic randomness introduced by stochastic optimization methods. The key of the analysis lies in ensuring the positive definiteness of suitable Gram matrices during the training. The analysis sheds insight into the dynamics of the optimization process, and provides guarantees on the neural networks trained by stochastic algorithms.

Equations

$$\{\begin{array}{rl}\Delta u& =f,\text{mbox}in\Omega ,\\ u& =g,\text{mbox}on\partial \Omega ,\end{array}$$

$$\mathcal{L}(\mathbf{w},\mathbf{a})=\frac{1}{2}{\int }_{\Omega }{\left(\Delta [\varphi (\mathbf{x};\theta )]-f(\mathbf{x})\right)}^{2}dx+\frac{\gamma }{2}{\int }_{\partial \Omega }{\left(\varphi (\mathbf{y};\theta )-g(\mathbf{y})\right)}^{2}dy,$$

$$\varphi (\mathbf{x};\mathbf{w},\mathbf{a})$$

$$|{\sigma }^{(k)}(z_1)-{\sigma }^{(k)}(z_2)|\le C_M|z_1-z_2|,k=0,1,2,3.$$

$$C_M\le C(1+M^{\ell }),\forall M>0.$$

$$|{\sigma }^{(k)}(z)|\le |{\sigma }^{(k)}(z)-{\sigma }^{(k)}(0)|+|{\sigma }^{(k)}(0)|\le C_M(|z|+1).$$

$$L(\mathbf{w},\mathbf{a}):=\sum _{p=1}^{n_1}\frac{|\Omega |}{2n_1}{\left(\sum _{i=1}^d\frac{{\partial }^2\varphi }{\partial x_i^2}({\mathbf{x}}_p;\mathbf{w},\mathbf{a})-f({\mathbf{x}}_p)\right)}^2+\sum _{q=1}^{n_2}\frac{\gamma |\partial \Omega |}{2n_2}{\left(\varphi ({\mathbf{y}}_q;\mathbf{w},\mathbf{a})-g({\mathbf{y}}_q)\right)}^2,$$

$$L(\mathbf{w},\mathbf{a}):=\sum _{p=1}^{n_1}\frac{1}{2n_1}{\left(\sum _{i=1}^d\frac{{\partial }^2\varphi }{\partial x_i^2}({\mathbf{x}}_p;\mathbf{w},\mathbf{a})-f({\mathbf{x}}_p)\right)}^2+\sum _{q=1}^{n_2}\frac{\gamma }{2n_2}{\left(\varphi ({\mathbf{y}}_q;\mathbf{w},\mathbf{a})-g({\mathbf{y}}_q)\right)}^2.$$

$$s_p(\mathbf{w},\mathbf{a})=\frac{1}{\sqrt{n_1}}\left(\sum _{i=1}^d\frac{{\partial }^2\varphi }{\partial x_i^2}({\mathbf{x}}_p;\mathbf{w},\mathbf{a})-f({\mathbf{x}}_p)\right)\text{mbox}and{h}_q(\mathbf{w},\mathbf{a})=\sqrt{\frac{\gamma }{n_2}}(\varphi ({\mathbf{y}}_q;\mathbf{w},\mathbf{a})-g({\mathbf{y}}_q)),$$

$$\mathbf{s}(\mathbf{w},\mathbf{a})=[s_1(\mathbf{w},\mathbf{a})\cdots s_{n_1}(\mathbf{w},\mathbf{a}){]}^{\top }\in {\mathbb{R}}^{n_1}\text{mbox}and\mathbf{h}(\mathbf{w},\mathbf{a})=[h_1(\mathbf{w},\mathbf{a})\cdots h_{n_2}(\mathbf{w},\mathbf{a}){]}^{\top }\in {\mathbb{R}}^{n_2}.$$

$$L(\mathbf{w},\mathbf{a})=\frac{1}{2}(\Vert \mathbf{s}(\mathbf{w},\mathbf{a}){\Vert }_2^2+\Vert \mathbf{h}(\mathbf{w},\mathbf{a}){\Vert }_2^2).$$

$$\begin{array}{r}\frac{d\mathbf{w}(t)}{dt}=-\frac{\partial L(t)}{\partial \mathbf{w}}=-\sum _{p=1}^{n_1}{s}_p(t)\cdot \frac{\partial s_p(t)}{\partial \mathbf{w}}-\sum _{q=1}^{n_2}{h}_q(t)\cdot \frac{\partial h_q(t)}{\partial \mathbf{w}},\end{array}$$

$$\begin{array}{r}\frac{d\mathbf{a}(t)}{dt}=-\frac{\partial L(t)}{\partial \mathbf{a}}=-\sum _{p=1}^{n_1}{s}_p(t)\cdot \frac{\partial s_p(t)}{\partial \mathbf{a}}-\sum _{q=1}^{n_2}{h}_q(t)\cdot \frac{\partial h_q(t)}{\partial \mathbf{a}}.\end{array}$$

$$\frac{d{s}_p(t)}{dt}=⟨\frac{\partial s_p(t)}{\partial \mathbf{w}},\frac{d\mathbf{w}(t)}{dt}⟩+⟨\frac{\partial s_p(t)}{\partial \mathbf{a}},\frac{d\mathbf{a}(t)}{dt}⟩=⟨\frac{\partial s_p(t)}{\partial \theta },\frac{d\theta (t)}{dt}⟩.$$

$$\frac{d}{dt}\left[\begin{array}{c}\mathbf{s}(t)\\ \mathbf{h}(t)\end{array}\right]=-\left({\mathbf{G}}_{\mathbf{w}}(t)+{\mathbf{G}}_{\mathbf{a}}(t)\right)\left[\begin{array}{c}\mathbf{s}(t)\\ \mathbf{h}(t)\end{array}\right],$$

$${\mathbf{G}}_{\mathbf{w}}(\mathbf{w},\mathbf{a})={\mathbf{D}}_{\mathbf{w}}^{\top }{\mathbf{D}}_{\mathbf{w}},{\mathbf{D}}_{\mathbf{w}}=\left[\begin{array}{cccccc}\frac{\partial s_1(\mathbf{w},\mathbf{a})}{\partial \mathbf{w}}& \cdots & \frac{\partial s_{n_1}(\mathbf{w},\mathbf{a})}{\partial \mathbf{w}}& \frac{\partial h_1(\mathbf{w},\mathbf{a})}{\partial \mathbf{w}}& \cdots & \frac{\partial h_{n_2}(\mathbf{w},\mathbf{a})}{\partial \mathbf{w}}\end{array}\right],$$

$${\mathbf{G}}_{\mathbf{a}}(\mathbf{w},\mathbf{a})={\mathbf{D}}_{\mathbf{a}}^{\top }{\mathbf{D}}_{\mathbf{a}},{\mathbf{D}}_{\mathbf{a}}=\left[\begin{array}{cccccc}\frac{\partial s_1(\mathbf{w},\mathbf{a})}{\partial \mathbf{a}}& \cdots & \frac{\partial s_{n_1}(\mathbf{w},\mathbf{a})}{\partial \mathbf{a}}& \frac{\partial h_1(\mathbf{w},\mathbf{a})}{\partial \mathbf{a}}& \cdots & \frac{\partial h_{n_2}(\mathbf{w},\mathbf{a})}{\partial \mathbf{a}}\end{array}\right].$$

$${\mathbf{G}}_{\mathbf{w}}^{\infty }:={\mathbb{E}}_{\theta (0)}\left[{\mathbf{G}}_{\mathbf{w}}(\theta (0))\right]\text{mbox}and{\mathbf{G}}_{\mathbf{a}}^{\infty }:={\mathbb{E}}_{\theta (0)}\left[{\mathbf{G}}_{\mathbf{a}}(\theta (0))\right].$$

$${\lambda }_{\mathbf{w}}:={\lambda }_{\min }({\mathbf{G}}_{\mathbf{w}}^{\infty })>0\text{mbox}and{\lambda }_{\mathbf{a}}:={\lambda }_{\min }({\mathbf{G}}_{\mathbf{a}}^{\infty })>0.$$

$$⟨\frac{\partial s_j(\theta )}{\partial \mathbf{w}},\frac{\partial s_k(\theta )}{\partial \mathbf{w}}⟩=\frac{1}{m}\sum _{r=1}^m⟨\mathbf{v}({\mathbf{w}}_r,a_r;{\mathbf{x}}_j),\mathbf{v}({\mathbf{w}}_r,a_r;{\mathbf{x}}_k)⟩=:\frac{1}{m}\sum _{r=1}^{m}V({\mathbf{w}}_r,a_r;{\mathbf{x}}_j,{\mathbf{x}}_k),$$

$$\mathbf{v}({\mathbf{w}}_r,a_r;{\mathbf{x}}_j)=\frac{1}{\sqrt{n_1}}\sum _{i=1}^d\left(a_r{\sigma }^{\prime \prime \prime }({\mathbf{w}}_r^{\top }{\tilde{\mathbf{x}}}_j)w_{ri}^2{\tilde{\mathbf{x}}}_j+2a_r{\sigma }^{\prime \prime }({\mathbf{w}}_r^{\top }{\tilde{\mathbf{x}}}_j)w_{ri}{\mathbf{e}}_i\right).$$

$${\mathbf{G}}_{\mathbf{w}}^{\infty }[j,k]={\mathbb{E}}_{\theta }\left[\frac{1}{m}\sum _{r=1}^{m}V({\mathbf{w}}_r,a_r;{\mathbf{x}}_j,{\mathbf{x}}_k)\right]=\frac{1}{m}\sum _{r=1}^m{\mathbb{E}}_{{\mathbf{w}}_r,a_r}\left[V({\mathbf{w}}_r,a_r;{\mathbf{x}}_j,{\mathbf{x}}_k)\right].$$

$${\mathbf{G}}_{\theta }(\mathbf{w},\mathbf{a})={\mathbf{G}}_{\mathbf{a}}(\mathbf{w},\mathbf{a})+{\mathbf{G}}_{\mathbf{w}}(\mathbf{w},\mathbf{a})={\mathbf{D}}_{\theta }^{\top }{\mathbf{D}}_{\theta },$$

$${\mathbf{D}}_{\theta }=\left[\begin{array}{cccccc}\frac{\partial s_1(\mathbf{w},\mathbf{a})}{\partial \theta }& \cdots & \frac{\partial s_{n_1}(\mathbf{w},\mathbf{a})}{\partial \theta }& \frac{\partial h_1(\mathbf{w},\mathbf{a})}{\partial \theta }& \cdots & \frac{\partial h_{n_2}(\mathbf{w},\mathbf{a})}{\partial \theta }\end{array}\right],$$

$$\mathrm{SGD}:$$

$$\mathrm{SGF}:$$

$$\Vert {\mathbf{w}}_r(0){\Vert }_2\le \sqrt{2(d+1)\log \left(\frac{2m(d+1)}{\delta }\right)},\forall r\in [m].$$

$$E(M):=\bigcap _{r=1}^m\left\{\Vert {\mathbf{w}}_r(0){\Vert }_2\le M\right\}.$$

$$B:=1+\sqrt{2(d+1)\log \left(\frac{2m(d+1)}{\delta }\right)}.$$

$$m\gtrsim \log \left(\frac{n_1+n_2}{\delta }\right),$$

$$L(0)\lesssim C_B^2{d}^2{\log }^2\left(\frac{n_1+n_2}{\delta }\right).$$

Full text

Key words: physics informed neural network, stochastic gradient descent, convergence, neural tangent kernel

1. Introduction

Partial differential equations (PDEs) represent one very popular and flexible class of mathematical models in nearly all disciplines in natural science and engineering, and their solutions is of critical importance. One major hurdle to solve PDEs is the notorious curse of dimensionality, i.e., the computational complexity grows exponentially with the problem dimensionality. Recently, deep learning using deep neural networks (DNNs) has emerged as a powerful tool for solving PDEs, and it has received a lot of attention due to its tremendous potentials to break the curse [16, 19, 20]. Various deep learning methods have been proposed to solve PDEs [37, 16, 47, 10]. We refer interested readers to the reviews [3, 40] for details.

Among various existing neural PDE solvers, the most popular method is physics-informed neural network (PINN) due to Raissi et al [37], which can be traced back to the seminal work [28] in the 1990s. PINN is based on the principle of PDE residual minimization, i.e., the loss is a suitable weighted combination of the residuals of the PDE (in the domain) and boundary condition (on the boundary), commonly measured in the standard $L^{2}$ norms, and employs DNNs as the ansatz space to approximate the PDE solutions. By the construction, PINN directly integrates the physical knowledge encoded in the differential equation into the neural networks. In practice, one learns neural network parameters by training the empirical losses using suitable optimization algorithms, including gradient descent, stochastic gradient descent, Adam [26], and limited-memory BFGS [6]. The method enjoys a number of distinct features, e.g., the ease of implementation, flexibility of the equation type (elliptic, parabolic and hyperbolic), and strong empirical performance for many PDEs problems. Indeed, it has shown impressive empirical performance across a wide range of challenging direct and inverse problems associated with PDEs, e.g., Navier-Stokes equations [24, 17], Hamilton-Jacobi-Bellman equation [39], and various PDE inverse problems [48, 23, 8]. These results clearly show its significant potentials in diverse scenarios, and thus it has received immense attention within the scientific and engineering computing community [12].

Despite the empirical successes of PINNs, the theoretical understanding of PINNs is still largely in its infancy. Several works have investigated the generalization error, and derived error bounds between the DNN approximation and exact solution in terms of DNN architecture parameters and the number of sampling points etc. under suitable a priori regularity assumption on the exact solution [34, 13, 21]. The analysis is mostly conducted under the assumption of zero or small error in the optimization procedure. However, due to the nonlinearity of the activation function in the neural network, the empirical loss is highly nonconvex with respect to the DNN parameters, and thus it is very challenging to find a global minimizer to the loss. In practice, the choice of a suitable optimization algorithm is crucial for PINNs to achieve satisfactory numerical results [4, 41]. The mathematical study on the convergence of the optimizers for PINN training is of fundamental importance.

There are several works analyzing the optimization process within the framework of neural tangent kernel (NTK) [22]. The key observation of the NTK is that it is essentially deterministic at initialization and then keeps invariant during training process when the width of the network tends to infinity. Based on this important observation, Du et al [15] proved that gradient descent (GD) can find a global minimum of the least-squares loss involving over-parameterized two-layer neural networks for the regression task, and later the authors generalized the analysis to deep neural networks (DNNs) [14]. Allen-Zhu et al [1] investigated both the GD and stochastic gradient descent (SGD) for over-parameterized DNN, convolutional neural networks (CNN), and residual neural networks (ResNet). All these works consider mean squared loss for the standard regression task.

When compared with the standard regression task, the convergence analysis of PINNs is more involved: the PINN loss includes both interior and boundary type terms, and involves the partial derivatives of the neural network. The first breakthrough on the convergence of the optimization algorithms for training PINNs for the standard parabolic equation is due to Gao, Gu and Ng [18], who proved that GD can reach a global minimizer of over-parameterized two-layer PINNs with the $\mathrm{ReLU}^{3}$ activation. (The ReLU3 activation ensures the well definedness of the PINN loss.) The key part of the analysis, compared to the classical regression and classification problems, is to handle the intricate interaction among the trainable parameters caused by the physics-informed loss, and to establish the positive definiteness of the associated Gram matrices under the ReLU3 activation. This analysis was recently refined by Xu et al [45] using a different error decomposition. See also [32] for an analysis of gradient flow under some technical assumptions. More recently, Xu et al [46] studied the implicit gradient descent (IGD) for training PINNs, and prove that for smooth globally Lipschitz non-polynomial activation $\sigma$, IGD converges to a globally optimal solution at a linear rate. Moreover, they observed that the learning rate can be chosen more flexibly because of the unconditional stability of IGD. All the above-mentioned works focus on the over-parameterized regime, where the number of DNN parameters far exceeds the number of sampling points. In the under-parameterized regime, Niessen and Muller [36] very recently provided the optimization guarantee of projected gradient descent for two-layer PINNs with finite width, and derived an overall error bound based on an existing approximation estimate. However, all these important works focus on deterministic algorithms, and the analysis does not extend directly to stochastic methods.

Stochastic gradient descent (SGD) and its variants are simple and powerful methods for training PINNs in practice and have been widely adopted. They are prevalent in training neural PDE solvers because of their low computation cost and excellent generalization performance. To the best of our knowledge, the convergence analysis of SGD for training PINNs remains missing, although the previous works offer the analysis for regression problems [1, 7]. Due to the intricate loss of PINNs (involving various partial derivatives) and the dynamic randomness induced by SGM, it is more challenging to investigate the trajectory of the training dynamics. Indeed, the convergence analysis requires the bound on all parameters but the stochasticity of the iterations allow estimating the parameters in a high probability sense. Moreover, the number of trainable parameters will increase as the network gets wider, which requires new techniques to control them uniformly.

In view of the wide adoption SGD for training PINNs in practice, it is important to analyze the convergence of SGD for training PINNs. This motivates the convergence analysis of stochastic gradient methods (SGMs) for over-parameterized two-layer PINNs, including stochastic gradient descent (SGD) and stochastic gradient flow (SGF), i.e., the optimization problem with the objective given by (3) and trained by SGM. This represents the main contribution of this work, i.e., the linear convergence of SGMs to a global minimum for training PINN loss in the high probability sense. In the analysis, we assume that the nonlinear activation function is locally Lipschitz, which is milder than that in existing works [18, 46]. First, we prove that with high probability, the relevant Gram matrices stay positive definite during the training dynamics under some mild assumptions. Second, we track the parameters and the loss at each step and prove that if the Gram matrices are positive definite, the SGD algorithm can find a zero minimum in the sense of expectation for the optimization problem. Third, we establish and study the stochastic differential equation satisfied by the parameters, and prove that the linear convergence result is also valid for the SGF. In the course of the convergence analysis, we employ a general concentration inequality to deal with the quantities related to initialization. For SGD, we analyze each iterate and control the parameters and the loss using the bound from each step. For SGF, we first bound the logarithm of the loss using properties of Ito integral and then bound the parameters.

The rest of the paper is organized as follows. In section 2, we describe the problem setting, and give some preliminary estimates. In section 3, we provide the convergence analysis of SGD. In section 4, we consider the continuous time model. The proofs are given in the appendix. Throughout, we use $m$ to denote the width of the first hidden layer of the PINNs, and both the words ‘over-parameterized’ and ‘wide’ mean that $m$ approaches to infinity. Moreover, we also use $[m]$ to denote the set $\{1,2,\ldots,m\}$. Given two quantities $a$ and $b$, the symbol $a\gtrsim b$ means that there is an absolute constant $C$, independent of the network width $m$, the numbers $n_{1},n_{2}$ of sampling points and the probability tolerance $\delta$, such that $a\geq C\cdot b$, and the symbols ‘$\lesssim$’ and ‘$\approx$’ have the similar meaning.

2. Preliminary

In this section, we describe the problem formulation and give preliminary estimates.

2.1. Physics informed neural networks

Let $\Omega\subset\mathbb{R}^{d}$ be an open bounded domain with a boundary $\partial\Omega$. Consider the Poisson equation with the Dirichlet boundary condition:

[TABLE]

where $\bm{x}=(x_{1},x_{2},\cdots,x_{d})\in\mathbb{R}^{d}$, and $\Delta u=\sum_{i=1}^{d}\frac{\partial^{2}u}{\partial x_{i}^{2}}$ is the Laplace operator. The analysis below can be also extended to other linear second-order parabolic or hyperbolic PDEs with the related computation being slightly more involved. Moreover, one can also consider Neumann and Robin boundary conditions similarly.

To approximate the PDE solution $u$ using a neural network $\phi(\bm{x};\theta)$ (with the network parameters $\theta$), PINN constructs a physics-informed expected loss based on the principle of PDE residual minimization:

[TABLE]

where the hyperparameter $\gamma$ balances the interior and boundary losses. This loss reflects how well the neural network function $\phi(\bm{x};\theta)$ satisfies the differential equation and the boundary condition, henceforth the name physics-informed neural network.

In this work, we employ the standard fully connected two-layer neural networks:

[TABLE]

where $\sigma:\mathbb{R}\to\mathbb{R}$ is a nonlinear activation function and $\bm{w}=[\bm{w}_{1}^{\top}\,\cdots\,\bm{w}_{m}^{\top}]^{\top}\in\mathbb{R}^{m(d+1)}$, $\bm{a}=[a_{1}\,\cdots\,a_{m}]^{\top}\in\mathbb{R}^{m}$, $\bm{w}_{r}=[{w}_{r1}\,{w}_{r2}\cdots{w}_{rd}\,{w}_{r,d+1}]^{\top}\in\mathbb{R}^{d+1}$, $\widetilde{\bm{x}}=[\bm{x}^{\top}\,1]^{\top}\in\mathbb{R}^{d+1}$. The bias terms ${w}_{r,d+1}$ are absorbed into the vectors $\bm{w}_{r}$ and $\widetilde{\bm{x}}\in\mathbb{R}^{d+1}$ is the extended point of $\bm{x}$. Let the vector $\theta:=[\bm{w}^{\top}\,\bm{a}^{\top}]^{\top}\in\mathbb{R}^{m(d+2)}$ denote all the parameters in the neural network. We make the following assumption about the initialization and the activation function $\sigma$.

Assumption 2.1**.**

• (i)

For all $r\in[m]$, $\bm{w}_{r}(0)$ is sampled from the standard Gaussian distribution and $a_{r}(0)$ is sampled from Rademacher distribution, i.e., $\bm{w}_{r}(0)\sim\mathcal{N}(\bm{0},\bm{I}_{d+1})$ and $a_{r}(0)\sim\mathrm{Unif}\{-1,1\}$.

• (ii)

The activation function $\sigma$ is piecewise continuously differentiable up to third order, and has locally Lipschitz property, i.e., for all $|z_{1}|,|z_{2}|\leq M$, there exists $C_{M}>0$ such that

[TABLE]

The constant $C_{M}$ grows at most polynomially with $M$: there exists $C>0$ and $\ell>0$ such that

[TABLE]

For all $M>0$, the Lipschitz constant $C_{M}\geq 1$ and satisfies $|\sigma^{(k)}(0)|\leq C_{M}$ for $k=0,1,2,3$. Then we have for all $|z|\leq M$

[TABLE]

Remark 2.1**.**

For the initialization of $a_{r}(0)$, we only use the property of zero mean and boundedness, and one can also use other initialization schemes. The assumption on $\sigma$ is quite generous and covers many commonly used functions. Indeed, the RePU function with power not lower than three and smooth activations, e.g., hyperbolic tangent and softplus, satisfy the assumption. In practice, smooth activation functions are widely used in neural PDE solvers. In comparison, the works [18, 45] investigates PINN with ReLU3 activation, whereas the work [46] focuses on nonpolynomial analytic and three-times globally Lipschitz activations.

In practice, the population loss $\mathcal{L}(\bm{w},\bm{a})$ has to be further discretized since the involved integrals cannot be evaluated in closed form. This is commonly achieved using Monte Carlo methods, which is especially attractive in the high-dimensional case. Specifically, let the training samples $\{\bm{x}_{p}\}_{p=1}^{n_{1}}$ and $\{{\bm{y}}_{q}\}_{q=1}^{n_{2}}$ be drawn independently and identically distributed (i.i.d.) from the uniform distributions $\mathcal{U}(\Omega)$ and $\mathcal{U}(\partial\Omega)$ on the interior domain $\Omega$ and boundary $\partial\Omega$, respectively. Then the empirical loss $L(\bm{w},\bm{a})$ reads

[TABLE]

where $|\Omega|$ and $|\partial\Omega|$ denote their respective Lebesgue measures. To simplify the notation, we absorb $|\Omega|$ and $|\partial\Omega|$ into $\gamma$, and investigate the empirical loss

[TABLE]

The gradient of the loss $L(\bm{w},\bm{a})$ can be evaluated efficiently using automatic differentiation, and thus standard first-order algorithms such as SGD can be readily implemented.

Below we assume the samples $\|\bm{x}_{p}\|_{2},\|{\bm{y}}_{q}\|_{2}\leq 1$, so that the extended samples satisfy $\|\widetilde{\bm{x}}_{p}\|_{2},\|\widetilde{\bm{y}}_{q}\|_{2}\leq\sqrt{2}$. Throughout, we use the following notation. We denote the loss in the domain $\Omega$ and on the boundary $\partial\Omega$ by $s_{p}(\bm{w},\bm{a})$ and $h_{q}(\bm{w},\bm{a})$, respectively:

[TABLE]

and accordingly the following loss vectors

[TABLE]

Then the empirical loss $L$ can be written as

[TABLE]

The notation $t$ denotes the iterate index (for SGD) or time (for SGF). We use the following notation interchangeably: $\theta(t)=\theta_{t}=(\bm{w}_{t},\bm{a}_{t})=(\bm{w}(t),\bm{a}(t))$ and $L(t)=L(\theta_{t})=L(\bm{w}_{t},\bm{a}_{t})$. Moreover, when dealing with the randomness of the initialization, we may abuse the notations $\theta:=[\bm{w}^{\top}\,\bm{a}^{\top}]^{\top}$ and $\theta(0):=[\bm{w}(0)^{\top}\,\bm{a}(0)^{\top}]^{\top}$.

2.2. The idea of convergence analysis

We first use gradient flow to illustrate the idea of convergence for the gradient method. In gradient flow, the parameters satisfy the following ODE:

[TABLE]

By the chain rule, the domain loss $s_{p}(t)$ satisfies

[TABLE]

Consequently, the loss vector satisfies

[TABLE]

with

[TABLE]

[TABLE]

Now, if the smallest eigenvalues of the two Gram matrices $\bm{G}_{\bm{w}}(t)$ and $\bm{G}_{\bm{a}}(t)$ have positive lower bounds during the training, then the loss $L$ tends to [math] when $t$ tends to infinity. This can be achieved if the ‘infinite’ Gram matrices are positive definite and the training is in the ‘lazy training’ regime [11]. For the former condition, we make the following assumption, where $\mathbb{E}_{\theta(0)}$ denotes taking expectation with respect to the distribution of $\theta(0)$. These quantities essentially govern the training dynamics.

Assumption 2.2**.**

Let $\bm{G}_{\bm{w}}^{\infty}$ and $\bm{G}_{\bm{a}}^{\infty}\in\mathbb{R}^{(n_{1}+n_{2})\times(n_{1}+n_{2})}$ be the expectation of the initial Gram matrices, i.e.,

[TABLE]

The infinite Gram matrices $\bm{G}_{\bm{w}}^{\infty}$ and $\bm{G}_{\bm{a}}^{\infty}$ are positive definite:

[TABLE]

Remark 2.2**.**

Since the losses $s_{p}(\theta)$ and $h_{q}(\theta)$ depend on the value of the neural network $\phi(\bm{x};\bm{w},\bm{a})$ at the sample points, the Gram matrices $\bm{G}_{\bm{w}}^{\infty}$ and $\bm{G}_{\bm{a}}^{\infty}$ are related to the neural network and sampling points, while its size depends on the number of samples. However, they are independent of $m$ even though $\phi(\bm{x};\bm{w},\bm{a})$ depends on $m$. Indeed, for $1\leq j,k,\leq n_{1}$, the analysis in Appendix A gives

[TABLE]

with

[TABLE]

The $(j,k)$-th entry of $\bm{G}_{\bm{w}}^{\infty}$ is

[TABLE]

Since all the entries $(\bm{w}_{r},a_{r})$ are i.i.d., the right hand side reduces to $\mathbb{E}_{\bm{w}_{r},a_{r}}\left[V(\bm{w}_{r},a_{r};\bm{x}_{j},\bm{x}_{k})\right]$. Thus the smallest eigenvalue $\lambda_{\bm{w}}$ and $\lambda_{\bm{a}}$ of $\bm{G}_{\bm{w}}^{\infty}$ and $\bm{G}_{\bm{a}}^{\infty}$ are fixed once the sampling points are chosen. See [35] for the tight bounds on the smallest eigenvalue of for deep ReLU nets.

The proof of Assumption 2.2 often proceeds as follows [18, 46]. First suppose no two samples of $\{\bm{x}_{p}\}_{p=1}^{n_{1}}\cup\{{\bm{y}}_{q}\}_{q=1}^{n_{2}}$ are identical, which implies that no two samples of $\{\widetilde{\bm{x}}_{p}\}_{p=1}^{n_{1}}\cup\{\widetilde{\bm{y}}_{q}\}_{q=1}^{n_{2}}$ are parallel. Then one proves that the columns of $\bm{D}_{\bm{w}}(\theta(0))$ and $\bm{D}_{\bm{a}}(\theta(0))$ are linearly independent in the Hilbert space $\mathcal{H}$ (which consists of the measurable functions from $\mathbb{R}^{m(d+2)}$ to $\mathbb{R}^{m(d+1)}$ or $\mathbb{R}^{m}$ such that $\mathbb{E}_{\theta(0)}\left[\|f({\theta(0)})\|_{2}^{2}\right]<\infty$). We also use the notation

[TABLE]

with

[TABLE]

and $\lambda_{\theta}=\lambda_{\bm{w}}+\lambda_{\bm{a}}$. Note that by Weyl’s inequality, $\lambda_{\min}(\bm{G}_{\theta}^{\infty})=\lambda_{\min}(\bm{G}_{\bm{w}}^{\infty}+\bm{G}_{\bm{a}}^{\infty})\geq\lambda_{\bm{w}}+\lambda_{\bm{a}}=\lambda_{\theta}$.

Then we have the following result on the convergence of the SGD and SGF (informal).

Theorem 2.1** (informal).**

Let Assumptions 2.1 and 2.2 hold, and use SGD / SGF for training the loss $L(\bm{w,a})$. If the width $m$ of the neural network is large enough and the step size $\eta$ is small enough, then with high probability, the expected value of the loss $L(t)$ decays exponentially

[TABLE]

2.3. Static randomness

First we deal with the quantities related to random initialization, the so-called static randomness compared with the stochastic dynamics in the algorithm, which arises from the random selection of the data at each iteration. The proofs can be found in Appendix A.

Lemma 2.1**.**

For any $\delta\in(0,1)$, with probability at least $1-\delta$, there holds

[TABLE]

Given any $M>0$, let $E(M)$ be the event that $\|\bm{w}_{r}(0)\|_{2}\leq M$ for all $r\in[m]$, namely,

[TABLE]

Define the quantity

[TABLE]

Then the event in Lemma 2.1 is $E(B-1)$. This lemma indicates that with high probability, all the $\bm{w}_{r}(0)$ are bounded by $B-1$. We conduct the analysis within the event $E(B-1)$ below. The choice of $B$ is to ensure that all the weights $\|\bm{w}_{r}(t)\|_{2}$ are bounded by $B$ during the iteration, which can be realized by the initial bound $B-1$ in the ‘lazy training’ regime.

Remark 2.3**.**

The quantity $B$ and thus the Lipschitz constant $C_{B}$ are related to the width $m$ of neural network, which only grows with $\log m$ polynomially and can be dominated by the term linear in $m$.

Next we analyze the initial loss $L(0)$, Gram matrices and the continuity of Gram matrices with respect to the NN parameters $\theta$. The main tool in the analysis is concentration inequality, and the proof is similar under the local Lipschitz assumption.

Lemma 2.2**.**

Within the event $E(B-1)$, if

[TABLE]

then with probability at least $1-\delta$, we have

[TABLE]

Lemma 2.3**.**

Within the event $E(B-1)$, if

[TABLE]

then with probability at least $1-\delta$, we have

[TABLE]

Lemma 2.3 shows that the spectral norm of the difference between initial and infinite Gram matrices is small for wide NNs, which ensures the positive definiteness of the initial Gram matrices. This agrees with the result in [22], which asserts that the initial Gram matrices are deterministic and positive when $m$ tends to infinity and the loss is convex in the function space. Note that the differential operator in PDE (1) is linear, and thus the loss $L(\theta)$ is still convex with respect to the neural network $\phi$. The next lemma shows the continuity of the Gram matrices, i.e., $\bm{G_{w}}$ and $\bm{G_{a}}$ are stable under small perturbations. These two are the key points for the positive definiteness of Gram matrices during training.

Lemma 2.4**.**

Within the event $E(B-1)$, if for some $0<R_{\bm{w}},R_{\bm{a}}\leq 1$ we have

[TABLE]

then the corresponding Gram matrix $\bm{G}_{\bm{w}}(\bm{w},\bm{a})$ and $\bm{G}_{\bm{a}}(\bm{w},\bm{a})$ satisfies

[TABLE]

In particular, if $\lambda_{\min}(\bm{G}_{\bm{w}}(0))\geq\frac{3}{4}\lambda_{\bm{w}},\lambda_{\min}(\bm{G}_{\bm{a}}(0))\geq\frac{3}{4}\lambda_{\bm{a}}$, and the constants $R_{\bm{w}}$ and $R_{\bm{a}}$ satisfy

[TABLE]

then $\lambda_{\min}(\bm{G}_{\bm{w}}(\bm{w},\bm{a}))\geq\frac{1}{2}\lambda_{\bm{w}}$ and $\lambda_{\min}(\bm{G}_{\bm{a}}(\bm{w},\bm{a}))\geq\frac{1}{2}\lambda_{\bm{a}}.$

3. Convergence Analysis of SGD

Consider the mini-batch SGD, i.e., the the parameters are updated using a subset of the sampling points:

[TABLE]

where the indices $\mathcal{I}\subset[n_{1}]$ are sampled uniformly among $[n_{1}]$ and $\mathcal{J}\subset[n_{2}]$ are sampled uniformly among $[n_{2}]$, and $I=|\mathcal{I}|,J=|\mathcal{J}|$. Then the SGD iteration is given by

[TABLE]

Note that both $\widetilde{L}(\bm{w},\bm{a})$ and its gradient are unbiased:

[TABLE]

Next we handle the dynamic randomness, and we always assume the initial Gram matrices have their respective positive lower bounds, which happens with high probability by Lemma 2.3 if $m$ satisfies the inequality therein. The proofs of the lemmas can be found in Appendix B. Let $\mathcal{F}_{t}$ be the known information after $t$ iterations (i.e., the filtration generated by the random initialization and the random iteration index selection up to the $t$th iteration) and let $R_{\bm{w}}$ and $R_{\bm{a}}$ be equal to the quantities given in equation (11). Let $T$ be the stopping time such that whenever $t<T$, we have

[TABLE]

while at least one of above inequalities fails for $t=T$. Then within $t<T$, by Lemma 2.4, we have $\lambda_{\min}(\bm{G}_{\bm{w}}(t))\geq\frac{1}{2}\lambda_{\bm{w}}$, $\lambda_{\min}(\bm{G}_{\bm{a}}(t))\geq\frac{1}{2}\lambda_{\bm{a}}$ and then $\lambda_{\min}(\bm{G}_{\theta}(t))\geq\frac{1}{2}\lambda_{\bm{w}}+\frac{1}{2}\lambda_{\bm{a}}=\frac{1}{2}\lambda_{\theta}$.

Remark 3.1**.**

Note that the last inequality in (15) is not sharp since we expect that the loss $L(t)$ is decreasing and $m$ will tend to infinity. This is a technical tool used to overcome the dynamic randomness in stochastic algorithms. Meanwhile, it holds with high probability as the training proceeds. It indeed requires larger $m$ to overcome randomness for the convergence when compared with gradient descent; see the proof of Theorem 3.1.

During the iteration, we need to control the distance of $\bm{w}_{r}(t+1)$ from $\bm{w}_{r}(0)$ for all $r\in[m]$. We will control them uniformly and thus let:

[TABLE]

This also applies to other quantities and the counterpart related to $a_{r},r\in[m]$.

Lemma 3.1**.**

Within the event $E(B-1)$ and $t<T$, we have

[TABLE]

The above estimates hold for all deterministic time $t$. We add the terms ’within the stopping time $t<T$’ in the description to emphasize that such bounds could only be derived under relevant conditions. Next we estimate the expectation of the loss by mathematical induction within the stopping time. Note that $L(t)\leq\sqrt{m}L(0)$ is square integrable for $t<T$, and the characteristic function $\mathbf{1}_{T>t}$ is $\mathcal{F}_{t}$-measurable while $\mathbf{1}_{T>t+1}$ may be not.

Lemma 3.2**.**

Within the event $E(B-1)$, if $\eta\lesssim\frac{\lambda_{\theta}}{C_{B}^{4}B^{12}}$ and $m\gtrsim\frac{L(0)^{2}}{C_{B}^{4}B^{12}}$, then within $t+1<T$, we have

[TABLE]

Remark 3.2**.**

We discuss within $t+1<T$ since we need the bounds of $\bm{w}_{r}(t+1),a_{r}(t+1),L(t)$ in (15) and the square integrability of $L(t+1)\cdot\mathbf{1}_{T>t+1}$, which means that the estimate holds even if we do not have the bound in (15) for $L(t+1)$, and this subtlety will be used in the proof of Theorem 3.1.

Within $t<T$, using Lemma 3.2, we can estimate the expectation of $L(t)$ by mathematical induction:

[TABLE]

Therefore, within the stopping time, the loss can be controlled in the expectation sense.

Corollary 3.1**.**

Within the event $E(B-1)$, if $\eta\lesssim\frac{\lambda_{\theta}}{C_{B}^{4}B^{12}}$ and $m\gtrsim\frac{L(0)^{2}}{C_{B}^{4}B^{12}}$, then within $t<T$, there holds

[TABLE]

We next use the above estimates to bound the distance of parameters $\bm{w}_{*}(T)$ and $a_{*}(T)$ at time $T$ from the initialization.

Lemma 3.3**.**

Within the event $E(B-1)$, if $\eta\lesssim\frac{\lambda_{\theta}}{C_{B}^{4}B^{12}}$ and $m\gtrsim\frac{L(0)^{2}}{C_{B}^{4}B^{12}}$, then at the stopping time $T$, we have

[TABLE]

Finally, we combine the preceding lemmas and deduce that that $T=\infty$ with high probability, and then the estimate for the loss holds during the whole training process.

Theorem 3.1**.**

Let $\lambda_{\theta}=\lambda_{\bm{w}}+\lambda_{\bm{a}}$. For $\delta,\widetilde{\delta}\in(0,1)$, if $m$ and $\eta$ satisfy

[TABLE]

where

[TABLE]

and $C_{B}$ grows at most polynomially in $B$, then with probability at least $1-2\delta$ over the initialization, we have $\mathbb{P}\left(T=\infty\right)\geq(1-\widetilde{\delta})^{2}$, and then, we have for all $t\in\mathbb{N}$,

[TABLE]

Remark 3.3**.**

Note that $B$ is of order $\sqrt{\log m}$ and $C_{B}$ has at most polynomial growth of $B$. Thus they can be dominated by $m$ and the condition for $m$ in (16) is reasonable.

4. Continuous Time Model

Now, consider the stochastic gradient flow. The analysis requires the local Lipschitz property of the activation function $\sigma$ up to order four in Assumption 2.1 in order to apply Ito’s formula and estimate the Hessian of the loss. For the continuous time model, consider the following SDE

[TABLE]

or equivalently, the parameters $\theta$ satisfy

[TABLE]

where $W_{t}\in\mathbb{R}^{m(d+2)}$ is the standard Wiener process, $\eta$ is the stepsize and $\sigma(\theta_{t})\in\mathbb{R}^{m(d+2)\times m(d+2)}$ satisfies

[TABLE]

See [29] for more details about the derivation of the continuous time model. For the component $\bm{w}_{r}$ of $\theta$, we have

[TABLE]

where $\xi_{\bm{w}_{r}}$ denotes the corresponding part in the noise vector $\int_{0}^{t}\sigma(\theta_{t})dW_{t}$. There is a similar formula for the component $a_{r}$. Similar to the case of SGD, we define the stopping time $T$ by the infimum $t>0$ such that at least one of the following inequalities fails, i.e.,

[TABLE]

Throughout we assume that the SDE (17) for $\theta$ has a unique maximal local solution up to its (random) blow-up time, which is obviously behind the stopping time. The following lemmas show the existence of the SDE solution; See the explanation before Theorem 4.1. For the rigorous treatment of the deterministic case based on the PL property of $\nabla L$ discussed below, we refer to [9, Theorem 1.2].

We now prove the loss function is $\lambda_{\theta}$-Polyak-Lojasiewicz within stopping time; see [30, 9, 2, 31] for more about convergence related to Polyak-Lojasiewicz property. As in previous section, we always assume the event in Lemma 2.3 happens. First, by the chain rule

[TABLE]

For $t<T$, we have $\lambda_{\min}(\bm{G}_{\theta}(t))\geq\frac{\lambda_{\theta}}{2}$. Thus

[TABLE]

Note that loss function $L(\theta_{t})$ is nonnegative, we obtain the $\lambda_{\theta}$-Polyak-Lojasiewicz inequality

[TABLE]

Then using Ito’s formula and $\lambda_{\theta}$-Polyak-Lojasiewicz property, we get an estimate on the logarithm of the loss $L(\theta_{t})$. Note that we could discuss the problem only under the circumstance such that the loss $L(\theta)$ does not attain zero in finite time, otherwise we are done in the opposite case. The proofs can be found in Appendix C.

Lemma 4.1**.**

Within the stopping time, we have

[TABLE]

where $M_{t}$ and $\langle M\rangle_{t}$ are the local martingale and its quadratic variation defined respectively by

[TABLE]

We first analyze the terms related to the local martingale $M_{t}$. Let

[TABLE]

Then we have $dM_{t}=v(t)dW_{t}$ and the exponential martingale $\mathcal{E}_{t}:=\exp\left(M_{t}-\frac{1}{2}\langle M\rangle_{t}\right)$ of $M_{t}$ takes the form

[TABLE]

Using multi-dimensional Ito’s formula, we have

[TABLE]

which implies it is a supermartingale, i.e., $\mathbb{E}\left[\mathcal{E}_{t}\mid\mathcal{F}_{s}\right]\leq\mathcal{E}_{s}$ for $s\leq t$. We refer interested readers to [25, Theorem 1.2 and Example 1.2] for further details.

Lemma 4.2**.**

Within the stopping time, the stochastic process $\mathcal{E}_{t}=\exp\left(M_{t}-\frac{1}{2}\langle M\rangle_{t}\right)$ is a supermartingale. In particular, $\mathbb{E}[\mathcal{E}_{t}]\leq\mathbb{E}[\mathcal{E}_{0}]=1$.

We next deal with the trace term in Lemma 4.1.

Lemma 4.3**.**

Within the event $E(B-1)$ and $t<T$, we have

[TABLE]

Now we can get a decay of the loss, and then the distance of parameters from origin within the stopping time.

Lemma 4.4**.**

Within the event $E(B-1)$, if $\eta\lesssim\lambda_{\theta}\left(C_{B}^{3}B^{9}\left(C_{B}B^{3}+\sqrt{\sqrt{m}L(0)}\right)\right)^{-1}$, then within $t<T$, we have

[TABLE]

Lemma 4.5**.**

Within the event $E(B-1)$, if $\eta\lesssim\lambda_{\theta}\left(C_{B}^{3}B^{9}\left(C_{B}B^{3}+\sqrt{\sqrt{m}L(0)}\right)\right)^{-1}$, then within $t<T$, we have

[TABLE]

The following theorem is obtained by combining the preceding lemmas like Theorem 3.1. Specifically, first, let $t$ approaches the stopping time from left, then using Lemmas 4.4, 4.5 and (conditional) Markov inequality to obtain that $T=\infty$ with high probability. Conditioned on $T=\infty$, we get the exponential decay of the loss by Lemma 4.4 again. We omit the details of the proof.

Theorem 4.1**.**

Let $\lambda_{\theta}=\lambda_{\bm{w}}+\lambda_{\bm{a}}$. For $\delta,\widetilde{\delta}\in(0,1)$, if $m$ and $\eta$ satisfy

[TABLE]

where

[TABLE]

and $C_{B}$ grows at most polynomially in $B$, then with probability at least $1-2\delta$ over the initialization, we have $\mathbb{P}\left(T=\infty\right)\geq(1-\widetilde{\delta})^{2}$, and then, we have for all $t\in\mathbb{N}$,

[TABLE]

5. Discussion

In this work, we have established that both stochastic gradient descent and stochastic gradient flow can find global minima of the empirical PINN loss, under very generous assumptions on the activation function $\sigma$ so long as the strict positivity of the infinite Grammian holds (cf. Assumption 2.2). When compared with the standard gradient descent, stochastic algorithms requires less computational effort per iterations but require more neurons to guarantee the convergence. There are several avenues deserving further investigations. First, the analysis can only handle linear PDEs since the positive definiteness of Gram matrices may fail in the nonlinear case [5]. It remains an interesting open question to develop the convergence analysis for nonlinear PDEs. Second, the smallest eigenvalues of Gram matrices play a central role in the convergence analysis. However, the precise characterization of the smallest eigenvalues is still largely missing. Empirically, they are deeply related to the network architecture, activation functions and the locations of the sampling points. We refer interested readers to the works [43, 44] for the spectral bias of PINNs. A thorough theoretical analysis of the spectral behavior of the Gram matrices for PINNs is an important open problem.

Acknowledgments

Bangti Jin is supported by Hong Kong RGC General Research Fund (Project 14306423 and Project 14306824), and a start-up fund from The Chinese University of Hong Kong. Longjun Wu is supported by International Joint PhD Supervision Scheme of The Chinese University of Hong Kong.

Appendix A Proofs for section 2

A.1. Preliminary estimates

First we estimate the derivatives of the loss and then give the concentration inequality. The neural network $\phi(\bm{x};\bm{w},\bm{a})$ and its derivatives with respect to the input are given by

[TABLE]

Then the interior loss $s_{p}(\bm{w},\bm{a})$ and its derivatives with respect to the neural network parameters $\bm{w}_{r},a_{r},r\in[m],$ are given by

[TABLE]

where $\{\bm{e}_{1},\bm{e}_{2},\ldots,\bm{e}_{d}\}$ denotes the standard Cartesian orthonormal basis of $\mathbb{R}^{d}$. Similarly, the boundary loss $h_{q}(\bm{w},\bm{a})$ and its derivatives are given by

[TABLE]

The derivatives of the loss $L(\bm{w},\bm{a})$ can also be computed by the chain rule

[TABLE]

If $\bm{w}_{r}$ is bounded by $B\geq 1$ and $a_{r}$ is bounded by $2$, then by Assumption 2.1 and the Cauchy-Schwarz inequality, we have the following bounds on the derivatives of the interior loss $s_{p}(\bm{w},\bm{a})$:

[TABLE]

The counterparts for the boundary loss $h_{q}(\bm{w},\bm{a})$ are given by

[TABLE]

Note that the expressions of $s_{p}$ and its derivatives involve $w_{ri}^{2}$, which are not necessarily bounded, and thus one cannot apply Hoeffding’s inequality to bound the relevant quantities. Instead we recall the concentration inequality for sub-Weibull random variables defined below, which is a generalization of Hoeffding’s inequality.

Definition A.1**.**

Let $\psi_{\alpha}:=\exp(x^{\alpha})-1$, with $\alpha>0$. For a real-valued random variable $X$, its $\psi_{\alpha}$-Orlicz norm is defined by

[TABLE]

The random variable $X$ is said to be sub-Weibull$(\alpha)$ if $\|X\|_{\psi_{\alpha}}<\infty$. For sub-Weibull$(\alpha)$ variable $X$, we have for any $t\geq 0$,

[TABLE]

Moreover, $X$ is sub-Weibull$(\alpha)$ if it satisfies this inequality.

Then the following properties hold for sub-Weibull random variables [42]:

• •

First, bounded random variables are sub-Weibull$(\alpha)$ for all $\alpha>0$;

• •

Second, if $X$ is sub-Weibull$(\alpha)$, then $X^{k}$ is sub-Weibull$\left(\frac{\alpha}{k}\right)$ for $k>0$;

• •

Third, $\|\cdot\|_{\psi_{\alpha}}$ is a norm for $\alpha\geq 1$ and a quasi-norm for $0<\alpha<1$. In particular, if both $X$ and $Y$ are sub-Weibull$(\alpha)$, then so is their sum $X+Y$, and

[TABLE]

Lemma A.1** ([27, Theorem 3.1]).**

If $X_{1},\ldots,X_{n}$ are independent mean zero random variables with $\|X_{i}\|_{\psi_{\alpha}}<\infty$ for all $1\leq i\leq n$ and some $\alpha>0$, then for any vector $a=(a_{1},\ldots,a_{n})\in\mathbb{R}^{n}$, the following bound holds:

[TABLE]

where $b=(a_{1}\|X_{1}\|_{\psi_{\alpha}},\ldots,a_{n}\|X_{i}\|_{\psi_{\alpha}})\in\mathbb{R}^{n}$,

[TABLE]

and for $\beta(\alpha)=\infty$ when $\alpha\leq 1$ and $\beta(\alpha)=\alpha/(\alpha-1)$ when $\alpha>1$,

[TABLE]

and the quantity $L_{n}^{*}(\alpha)=L_{n}(\alpha)C(\alpha)\|b\|_{2}/\|b\|_{\beta(\alpha)}$.

A.2. Proof of Lemma 2.1

Since $\bm{w}_{r}\sim\mathcal{N}(\bm{0},\bm{I}_{d+1})$ for each $r\in[m]$, it suffices to control each of its entries. For any $M>0$, if the length of $\bm{w}_{r}$ is not less than $M$, then there exists at least one of its components whose length is not less than $\frac{M}{\sqrt{d+1}}$. By the standard concentration inequality for sub-Gaussian random variable $X\sim\mathrm{subG}(\sigma^{2})$ [38, Lemma 1.3]

[TABLE]

and noting for Gaussian variable $X\sim\mathcal{N}(0,\sigma^{2})$ satisfies $X\sim\mathrm{subG}(\sigma^{2})$, we get

[TABLE]

The lemma is obtained by setting the right hand side to $\delta/m$.

A.3. Proof of Lemma 2.2

For each interior loss $s_{p}(\bm{w},\bm{a})$, by the Cauchy-Schwarz inequality, we have

[TABLE]

Define the independent random variables

[TABLE]

Then within the event $E(B-1)$, by the local Lipschitz property of $\sigma^{\prime\prime}$, we have

[TABLE]

Note that $\|\bm{w}_{r}\|_{2}^{3}=\left(\sum_{i=1}^{d+1}w_{ri}^{2}\right)^{3/2}\lesssim\sum_{i=1}^{d+1}w_{ri}^{3}$ with each $w_{ri}\sim\mathcal{N}(0,1)$. Thus we have for $\alpha=\frac{2}{3}$,

[TABLE]

from which we derive

[TABLE]

We then have, by Lemma A.1, with probability at least $1-2\delta$,

[TABLE]

By discarding the constant term $f(\bm{x}_{p})^{2}$, we have

[TABLE]

Likewise, for the boundary loss $h_{q}(\bm{w},\bm{a})$, we have

[TABLE]

With $\alpha=2$, by Lemma A.1 we have with probability at least $1-2\delta$,

[TABLE]

By combining the preceding estimates, we have with probability at least $1-\delta$,

[TABLE]

Under the given condition on $m$, the second term is smaller than the first term, and the conclusion follows.

A.4. Proof of Lemma 2.3

Since the spectral norm of a matrix is bounded by its Frobenius norm, it suffices to bound each entry of the difference matrix

[TABLE]

There are three forms of its entries:

[TABLE]

The case $n_{1}+1\leq j\leq n_{1}+n_{2},\;1\leq k\leq n_{1}$ is the same as second case. We analyze the three cases of $\Delta\bm{G}_{\bm{w}}^{\infty}[j,k]$ separately. For any fixed pair $1\leq j\leq n_{1},1\leq k\leq n_{1}$, we write the $(j,k)$-th entry of $\Delta\bm{G}_{\bm{w}}^{\infty}$ as

[TABLE]

where $\{X_{r}\}_{r=1}^{m}$ are independent random variables. Within the event $E(B-1)$, $X_{r}$ satisfies

[TABLE]

Similar to Lemma 2.2, for $\alpha=1/3$, we have

[TABLE]

Thus we have with probability at least $1-2\delta$,

[TABLE]

Similarly for $1\leq j\leq n_{1},n_{1}+1\leq k\leq n_{1}+n_{2}$ and $n_{1}+1\leq j\leq n_{1}+n_{2},1\leq k\leq n_{1}$ we have $\alpha=1/2$, with probability at least $1-2\delta$,

[TABLE]

For $n_{1}+1\leq j,k\leq n_{1}+n_{2}$ we have $\alpha=1$, with probability at least $1-2\delta$,

[TABLE]

Square and sum over these elements, we have with probability at least $1-\delta$

[TABLE]

In order to bound $\|\Delta\bm{G}_{\bm{w}}^{\infty}\|_{F}$ by $\lambda_{\bm{w}}$, we need the right side smaller that $\lambda_{\bm{w}}^{2}$, from which we get

[TABLE]

Similarly, using the bound on the derivative with respect to $a_{r}$ given in (22), to bound $\|\Delta\bm{G}_{\bm{a}}^{\infty}\|_{F}$ by $\lambda_{\bm{a}}$, we need

[TABLE]

A.5. Proof of Lemma 2.4

Within the event $E(B-1)$, we have for all $r\in[m]$,

[TABLE]

The derivatives of the interior loss $s_{j}$ have been computed earlier, and we recall that

[TABLE]

Using the Lipschitz property of the activation function $\sigma$ and the polynomial function we get

[TABLE]

where the last step follows from the following inequalities

[TABLE]

The difference for the boundary loss $h_{j}$ is

[TABLE]

The entry of the difference matrix $\Delta\bm{G}_{\bm{w}}:=\bm{G}_{\bm{w}}(\bm{w},\bm{a})-\bm{G}_{\bm{w}}(\bm{w}(0),\bm{a}(0))$ involves three forms. For any fixed pair $1\leq j\leq n_{1},1\leq k\leq n_{1}$, using the estimate (21), we have

[TABLE]

Meanwhile, for $1\leq j\leq n_{1},n_{1}+1\leq k\leq n_{1}+n_{2}$ and $n_{1}+1\leq j\leq n_{1}+n_{2},1\leq k\leq n_{1}$, we have

[TABLE]

For $n_{1}+1\leq j,k\leq n_{1}+n_{2}$, we have

[TABLE]

Therefore, the following estimate holds

[TABLE]

Similarly, we have

[TABLE]

from which we deduce

[TABLE]

Combining the preceding estimates completes the proof of the lemma.

Appendix B Proofs for section 3

B.1. Proof of Lemma 3.1

According to the iteration process, for all $r\in[m]$,

[TABLE]

Within the event $t<T$, we take conditional expectation

[TABLE]

using the Cauchy-Schwarz inequality in the third line and the estimates (21) and (22) in the fourth line. The inequality holds uniformly for all $r\in[m]$, and thus the desired estimate follows. Similarly, we have for all $r\in[m]$,

[TABLE]

B.2. Proof of Lemma 3.2

The proof of the lemma requires the following result.

Lemma B.1**.**

Let $(\Omega,\mathcal{A},\mathbb{P})$ be a probability space and $\mathcal{F}\subset\mathcal{A}$ be a sub-$\sigma$-algebra. Let $X$ be a random variable with $\mathbb{E}\left[X^{2}\right]<\infty$. If $X=0\;a.e.$ on $A^{c}$ for some $A\in\mathcal{F}$, then the conditional expectation $Y:=\mathbb{E}\left[X\mid\mathcal{F}\right]$ also vanishes $a.e.$ on $A^{c}$.

Proof.

Let $\widetilde{Y}=Y\mathbf{1}_{A}$. Since $A\in\mathcal{F}$ and $Y$ is $\mathcal{F}$-measurable by the definition of conditional expectation, $\widetilde{Y}$ is $\mathcal{F}$-measurable. By the property of conditional expectation,

[TABLE]

Meanwhile, we have

[TABLE]

Combining these identities yields $Y^{2}\cdot\mathbf{1}_{A^{c}}\equiv 0$ almost everywhere. ∎

In the following discussions, we always discuss within the event $E(B-1)$ and $t+1<T$. We may omit the function $\mathbf{1}_{T>t+1}$ or $\mathbf{1}_{T>t}$ whenever there is no ambiguity. For $\tau\in[0,1]$, we define the linearly interpolated weight

[TABLE]

Note that $\|\bm{w}_{r}(t+1)\|_{2}\leq B,\|\bm{w}_{r}(t)\|_{2}\leq B$ and $\|a_{r}(t+1)\|_{r}\leq 2,\|a_{r}(t)\|_{2}\leq 2$. Then, by the convexity of norms, we also have

[TABLE]

for all $\tau\in[0,1]$, where $\bm{w}_{r}(t+\tau)$ and $a_{r}(t+\tau)$ are the corresponding component of $\theta_{t+\tau}$.

First, $L(t+1)\cdot\mathbf{1}_{T>t+1}$ is square integrable since $L(t+1)\leq\sqrt{m}L(0)$. By Lemma B.1, we have

[TABLE]

Second, consider the difference between interior loss $s_{p}$ at two successive steps:

[TABLE]

We take conditional expectation directly and obtain

[TABLE]

Meanwhile, using the argument in Lemma 2.4, we have

[TABLE]

Similarly,

[TABLE]

By taking conditional expectation and summing over $r\in[m]$, we get

[TABLE]

where we have used the estimates from proof of Lemma 3.1

[TABLE]

Likewise, for the boundary loss $h_{q}$, we have

[TABLE]

and the estimates

[TABLE]

Combining these estimates gives

[TABLE]

where for some constant $c_{1}>0$, the residual satisfies

[TABLE]

Moreover, recall that

[TABLE]

Using the local Lipschitz property of the activation function $\sigma$ and polynomial function, we get

[TABLE]

Squaring both sides and then taking the expectation yield

[TABLE]

A similar estimate holds for the boundary loss $h_{q}$. Then we obtain for some $c_{2}>0$,

[TABLE]

In summary, we obtain

[TABLE]

Using the inequality $\sqrt{L(t)}\leq\sqrt{\sqrt{m}L(0)}$ for the last term, we get

[TABLE]

If $\eta$ is small enough and $m$ is moderately large such that

[TABLE]

then we have

[TABLE]

Finally, under the above condition, we have

[TABLE]

B.3. Proof of Lemma 3.3

By the triangle inequality, we have

[TABLE]

Similarly, we have

[TABLE]

B.4. Proof of Theorem 3.1

First, with probability at least $1-\delta$, the event $E(B-1)$ occurs, and both the events in Lemma 2.2 and Lemma 2.3 happen with probability at least $1-\delta$ if

[TABLE]

Conditioned on these events, by Lemma 3.3 and Markov inequality, if

[TABLE]

where $L_{0}$ is the upper bound of the initial loss $L(0)$ in Lemma 2.2, we have

[TABLE]

with the quantities $R_{\bm{w}}$ and $R_{\bm{a}}$ defined in (11). Thus the inequalities for all $\bm{w}_{r}$ and $a_{r}$ in (15) hold with probability at least $1-\widetilde{\delta}$ if

[TABLE]

Conditioned on this event $\mathcal{E}$, all the $s_{p}(T),h_{q}(T)$ are bounded and thus $L(T)$ is square integrable. Then by Corollary 3.1 and the remark before it, we have

[TABLE]

and then by conditional Markov inequality

[TABLE]

Thus $L(T)\leq\sqrt{m}L(0)$ happens with probability at least $1-\widetilde{\delta}$ if $m\gtrsim\widetilde{\delta}^{-2}$. Note that if all the inequalities in (15) are satisfied at step $T$, then we have $T=\infty$ according to the definition of the stopping time. Therefore, if $\eta$ and $m$ satisfy the above conditions, we have that with probability at least $1-2\delta$, there holds

[TABLE]

In the case $T=\infty$, the second assertion is implied by Corollary 3.1.

Remark B.1**.**

If we just get the following estimate

[TABLE]

then by Markov inequality

[TABLE]

Then the right hand side has to be made smaller than $\delta/m$ in order to ensure that $\bm{w}_{r}(T)$ stays in the small ball of $\bm{w}_{r}(0)$ for all $r\in[m]$ with high probability, which however is impossible. Thus, we have to control $\|\bm{w}_{r}(T)-\bm{w}_{r}(0)\|_{2}$ and $\|a_{r}(T)-a_{r}(0)\|_{2}$ uniformly. The method is to control the distance using only the loss $L(t)$ and uniform bound of $\bm{w}_{r}$ within the stopping time.

Appendix C Proofs for section 4

C.1. Proof of Lemma 4.1

We may drop the variable $\theta_{t}$. The proof employs the following multi-dimensional Ito’s formula: for an Ito process $X_{t}$ defined by $dX_{t}=b(w,t)dt+\sigma(w,t)dW_{t}$, the process $Y_{t}=f(X_{t})$ satisfies

[TABLE]

We apply this result to $f=\log L$. Direct computation gives

[TABLE]

Note that within the stopping time, the loss $L(\theta)$ is $\lambda_{\theta}$-Polyak-Lojasiewicz. Using multi-dimensional Ito’s formula for the function $\log L(\theta_{t})$ gives

[TABLE]

with the indefinite integral $M_{t}$ and its quadratic variation $\langle M\rangle_{t}$ defined in (20). The integrability of $M_{t}$ is followed by the inequalities

[TABLE]

from which we obtain that $M_{t}$ is a local martingale.

C.2. Proof of Lemma 4.3

Since the matrix $\sigma(\theta_{t})\sigma(\theta_{t})^{\top}$ is positive semidefinite and $\mathrm{Hess}(L(\theta_{t}))$ is symmetric, by von Neumann’s trace inequality [33], we have

[TABLE]

Within the stopping time $t<T$, the estimates (21) and (22) hold. Consequently, by the proof in last lemma,

[TABLE]

With $\otimes$ being the Kronecker product, the Hessian of the loss $L(\theta)$ is given by

[TABLE]

For two vectors $\bm{u}$ and $\bm{v}$, we have $\|\bm{u}\otimes\bm{v}\|_{2}=\|\bm{u}\|_{2}\|\bm{v}\|_{2}$. Therefore, within $t<T$, we get

[TABLE]

Similarly, we have

[TABLE]

For all $r\in[m]$, we have

[TABLE]

and

[TABLE]

These second derivatives of $L$ are bounded:

[TABLE]

For $\alpha,\beta\in[m]$ such that $\alpha\neq\beta$, we have

[TABLE]

[TABLE]

These second derivatives of $L$ are also bounded:

[TABLE]

In summary, we have

[TABLE]

so that

[TABLE]

C.3. Proof of Lemma 4.4

Within $t<T$, by Lemma 4.3, we have for some $c_{1}$,

[TABLE]

By taking exponential, we get

[TABLE]

The desired conclusion follows by taking expectation of both sides on the inequality.

C.4. Proof of Lemma 4.5

Recall that the components of $\theta_{t}$ satisfy

[TABLE]

where $\xi_{\bm{w}_{r}}$ denotes the corresponding part in the noise vector $\int_{0}^{t}\sigma(\theta_{t})dW_{t}.$ Consequently,

[TABLE]

Meanwhile, we have

[TABLE]

By the property of Ito integral,

[TABLE]

Within $t<T$, by the Cauchy-Schwarz inequality and Lemma 4.4, we get

[TABLE]

In the same manner, we deduce

[TABLE]

This completes the proof of the lemma.