Optimization and Learning With Nonlocal Calculus

Our primary contributions are the theory and algorithmic applications of nonlocal differential operators, mainly nonlocal first and second order (gradient and Hessian) operators to perform optimization tasks. Specifically, we develop nonlocal analogs of gradient (including stochastic) descent, as well as Newton’s method. We establish convergence results for these methods that show how the optima obtained from nonlocal optimization converge to the standard gradient/Hessian based optima, as the nonlocal interaction effects become less pronounced. While our definition of nonlocal optimization methods is valid for irregular problems, our analysis will assume the required regularity of the objective functions in order to compare the local and nonlocal approaches. As a stylized numerical application of our nonlocal learning framework, we apply our nonlocal gradient descent procedure to the problem of parameter estimation on a non-differentiable “translation" manifold as studied via multiscale smoothing in [51] and [50]. In our framework, this problem can be solved without resorting to a multiscale analysis of the manifold; instead, the nonlocal gradient operators we use can be viewed as providing the required smoothing. While we have emphasized the use of nonlocal operators to non-differentiable optimization as an alternative to traditional subgradient methods ([43]), the scope of the theory is much wider. Indeed, one may use nonlocal methods to model long range interactions ([37]), phase transitions ([15]) etc.

Nonlocal models have, in recent years, been fruitfully deployed in a wide variety of application areas: in image processing ([25]) nonlocal operators are used for total variation (TV) inpainting and also image denoising ([11, 33, 32, 52]) with the nonlocal means framework. Computational/engineering mechanics has profited greatly in recent years with a plethora of nonlocal models. Plasticity theory ([8]) and fracture mechanics ([27, 21]) are two prominent examples. Indeed, the field of peridynamics (see [46, 45, 39]) is a prime consumer of nonlocal analysis. Nonlocal approaches have also featured in crowd dynamics [14], and cell biology ([1]). A nonlocal theory of vector calculus has been developed along with several numerical implementations (see [20, 4, 18]). The nonlocal approach has also found its way into deep learning ([49]). On the theoretical side, nonlocal and fractional differential operators have been studied in the context of nonlocal diffusion equations ([5] and references therein) as well as purely in the context of nonlocal gradient/Hessian operators ([35, 30, 10, 9, 29]). The use of fractional calculus via Caputo derivatives and Riemann-Liouville fractional integrals for gradient descent (GD) and backpropagation has been considered in [54, 38, 12, 42, 53]. Our work uses instead the related notion of nonlocal calculus from [35, 30, 10, 18] as its basis for developing not only a nonlocal (stochastic) GD, but also a nonlocal Newton’s method, and is backed by rigorous convergence theory. In addition, our nonlocal stochastic GD approach interprets, in a very natural way, the nonlocal interaction kernel as a conditional probability density. The multiple definitions of fractional and nonlocal calculi are related (see for e.g. [16]).

In all that follows, $\Omega\subset\mathbb{R}^{D}$ will denote a simply connected bounded domain (i.e., a bounded, connected, simply connected, open subset of $\mathbb{R}^{D}$) with a smooth, compact boundary $\partial\Omega$ where $D$ is a fixed positive integer. $\Omega$ will be endowed with the subspace topology of $\mathbb{R}^{D}$. We will consistently denote by $x,y,z$ points in $\Omega$. Subscripts of vectorial (resp. tensorial) quantities will denote the corresponding components. For e.g., $x_{i}$ or $A_{i,j}$ will denote the $i^{\text{th}}$ component of the vector $x$ or $(i,j)^{\text{th}}$ entry of matrix $A$ respectively for $i,j=1,\ldots,D$. For vectors $x,y\in\mathbb{R}^{D}$ the notation $x\otimes y\in\mathbb{R}^{D\times D}$ denotes the outer product of $x$ and $y$. The letters $i,j,k,n,m,N$ will denote positive integers while $t,s,p,\epsilon,\delta,M,K$ denote arbitrary real quantities. By $B_{\epsilon}(x)$, we mean an open ball centered at $x$ with radius $\epsilon>0$. Furthermore, $u,v,f,g,\rho:\Omega\rightarrow\mathbb{R}$ will in general denote real-valued functions with domain $\Omega$. All integrals will be in the sense of Lebesgue, and a.e. means almost every or almost everywhere with respect to Lebesgue measure. Given an arbitrary collection of subsets $\mathcal{S}$ of $\mathbb{R}^{D}$, the $\sigma$-algebra generated by $\mathcal{S}$ will be denoted by $\sigma(\mathcal{S})$.

Given $u:\Omega\rightarrow\mathbb{R}$, we define the support $\text{supp}(u)\subset\Omega$ as

where the overbar denotes closure in $\mathbb{R}^{D}$. A multi-index $\bar{n}=(n_{1},\ldots,n_{D})$ is an ordered $D$-tuple of non-negative integers $n_{i},i=1,\ldots D$, and we let $|\bar{n}|:=\sum_{i=1}^{D}n_{i}$. For a suitably regular function $u$, we denote by

If we consider all multi-indices $\bar{n}$ with $|\bar{n}|=1$, we can identify $\frac{\partial^{\bar{n}}u}{\partial x^{\bar{n}}}$ with the usual gradient $\nabla u$ of $u$. In this case, we denote the $i^{\text{th}}$ component $\frac{\partial u}{\partial x_{i}}$ of $\nabla u$ as $u_{x_{i}}$ for $i=1,\ldots D$. Likewise, the Hessian of $u:=Hu$ at a point $x\in\Omega$ is the $D\times D$ matrix $(Hu(x))_{i,j}=\frac{\partial^{2}u}{\partial x_{i}\partial x_{j}}(x)$. We shall use the following function spaces:

We let $L^{p}(\Omega,\mathbb{R}^{1}):=L^{p}(\Omega)$. Note that each of the above spaces $\mathcal{X}$ is a normed linear space. In particular $H^{1}(\Omega):=W^{1,2}(\Omega)$ is a Hilbert space with the inner product $(u,v)_{H^{1}(\Omega)}=(u,v)_{L^{2}(\Omega)}+(\nabla u,\nabla v)_{L^{2}(\Omega)}$. The subspace $H_{0}^{1}(\Omega)$ of $H^{1}(\Omega)$ consists of those $u\in H^{1}(\Omega)$ that vanish (i.e., have zero trace) on the boundary $\partial\Omega$. Note that the space $C_{0}^{1}(\Omega)$ consisting of continuously differentiable functions vanishing on $\partial\Omega$ is dense in $H^{1}_{0}(\Omega)$ (see [2]). Next, $|\cdot|$ denotes the absolute value of real quantities, $\|\cdot\|$ will denote the Euclidean norm of $\mathbb{R}^{D}$ and $\|\cdot\|_{\mathcal{X}}$ will denote the norm in any of the function spaces $\mathcal{X}$ above.

Following [35, 10], we will fix a sequence of radial functions $\rho_{n}:\mathbb{R}^{D}\rightarrow\mathbb{R},n=1,\ldots$ which satisfy the following properties:

Note that by radial we mean that $\rho_{n}(x)=\rho_{n}(Ox)$ for any orthogonal matrix $O:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ with unit determinant. Alternatively, this means that $\rho_{n}(x)=\hat{\rho}_{n}(\|x\|)$ for some function $\hat{\rho}_{n}:\mathbb{R}\rightarrow\mathbb{R}$. We can specify (but will not insist) that the $\rho_{n}$ have compact support or even be of class $C^{\infty}_{c}(\Omega)$. Note that for our analysis of the nonlocal stochastic gradient descent method, we shall interpret the $\rho_{n}$ as a sequence of probability densities on $\mathbb{R}^{D}$ that approach the Dirac density $\delta_{0}$ centered at the origin. Prototypical examples of $\rho_{n}$ are given by sequences of Gaussian distributions with increasing variance (the non-compact support case) and the so-called bump functions with increasing height (the compact support case). In both these cases, we can show that the sequence of $\rho_{n}$ converge as Schwartz distributions to the Dirac density $\delta_{0}$. Details can be found in [23]. Finally, the set of all radial probability densities on $\mathbb{R}^{D}$, i.e., those $\rho:\mathbb{R}^{D}\rightarrow\mathbb{R}$ that satisfy the first two conditions of equation 4 will be denoted by $\mathcal{P}$.

Let us now focus on the main objects of study in this paper, nonlocal differential operators. We note that there is no unique notion of “the" nonlocal gradient (resp. Hessian). Indeed, differing choices of kernels may give rise to differing notions of gradients (resp. Hessians).

After this introductory section, we consider $1^{\text{st}}$ order nonlocal optimization in section 2. section 3 deals with a stylized numerical example of the $1^{\text{st}}$ order theory. In section 4, we develop the second order nonlocal optimization theory and the nonlocal analog of Newton’s method. We conclude in section 5. The proofs of all theoretical results are provided in the supplementary appendix section 6.

Consider first the standard affine approximant

to $u(x)\in C^{1}(\Omega)$ at $x_{0}$ with the corresponding remainder term $r(x_{0},x)=u(x)-A(x_{0},x).$ The nonlocal analog of the affine approximant is

and the corresponding nonlocal remainder term is $r_{n}(x_{0},x)=u(x)-A_{n}(x_{0},x).$ Standard calculus dictates that the remainder term $r(x_{0},x)=o(\|x-x_{0}\|)$, and by uniqueness of the derivative, we can conclude that the nonlocal remainder $r_{n}(x_{0},x)$ will in general not be $o(\|x-x_{0}\|)$. Indeed, the standard (local) affine approximant is asymptotically the unique best fit linear polynomial and hence, the nonlocal affine approximant must necessarily be sub-optimal. It is therefore important to understand how $r_{n}(x_{0},x)$ relates to $r(x_{0},x)$ as $n\rightarrow\infty$, and this is the objective of the next result.

We define the nonlocal gradient descent as follows. We first fix the scale of nonlocality $n$, and initialize $x^{0,n}\in\Omega$. We then perform the update:

Although one can easily conceive of an adaptive procedure with $n$ varying with step $k$, we emphasize that for our analysis, the scale $n$ of nonlocality does not change between the steps of the descent. Henceforth, we denote by $\nabla u^{k}:=\nabla u(x^{k})$ and $\nabla_{n}u^{k,n}:=\nabla_{n}u(x^{k,n})$. The step size $\alpha_{n}^{k}$ and initial point $x^{0,n}\in\Omega$ are chosen small enough to ensure that the iterates $x^{k,n}$ all remain in $\Omega$. We emphasize that nonlocal gradient descent is applicable so long as $\nabla_{n}u$ exists, in particular, even if $\nabla u$ does not exist.

In this section, we consider the nonlocal version of stochastic gradient descent (SGD). A generic descent algorithm with a proxy-gradient $g^{k}\in\mathbb{R}^{D}$:

After $K$ iterations, we output the averaged vector $\bar{x}=\frac{1}{K}\sum_{k=1}^{K}x^{k}$. If the expected value of $g^{k}$ equals the true gradient $\nabla u^{k}$, then, roughly, the difference between the expected value of $u(\bar{x})$ and the true minimum $u(x^{*})$ approaches $0$ as the number of iterations approaches $\infty$. This assumes that we sample $g^{k}$ from a distribution such that $\mathbb{E}(g^{k}|x^{k})=\nabla u(x^{k})$, or, more generally, $\mathbb{E}(g^{k}|x^{k})\in\partial u(x^{k})$, where $\partial u(x^{k})$ is the subdifferential set of $u$ at $x^{k}$ (see [41]).

As motivation for our nonlocal stochastic gradient descent analysis, let $X,Y$ be $\Omega-$valued random variables, with a conditional density $p_{Y|X}(y|x)=\rho(x-y).$ Here, $\rho(\cdot)$ is a radial density in the class $\mathcal{P}$ defined in section 1.3. Note that since we interpret $\rho(x-y)$ to be $p_{Y|X}(y|x)$, we can choose a density $p_{X}(x)$ and arrive at a joint distribution $p_{X,Y}(x,y)=\rho(x-y)p_{X}(x)$, and upon integration, arrive at the marginal distribution $p_{Y}(y)$ of $Y$. Let us denote by

If we let $\mathbb{E}_{\rho}[\cdot]$ denote the expectation with respect to the density $\rho$, then it is evident that

or if we interpret $\nabla_{\rho}u(X)$ as a random variable, we can write $\nabla_{\rho}u(X)=\mathbb{E}_{\rho}[k_{u}(X,Y)|X].$ We will set $g^{k}=k_{u}(x^{k},y)$, where $y$ is sampled from $p_{Y}(y)$ and perform the standard SGD with the knowledge that $\mathbb{E}_{\rho}[g^{k}|x^{k}]=\nabla_{\rho}u(x^{k})$. We can therefore consider a sequence of nonlocal SGD with kernels $\rho_{n}(\cdot)$, and as $n\rightarrow\infty$, we prove that the nonlocal optima converge to the usual SGD optima. In the following, we assume $u$ is a convex, Lipschitz function with Lipschitz constant $M$.

Lemma 2.7 is a motivation of the following (well-known, see [36, 26, 48, 7]) generalization of the notion of the subderivative set of a function.

It is obvious from the definition of $\partial_{\epsilon}u(x)$ and Lemma 2.7 that given an $\epsilon>0$, there is an $N>0$ such that for all $n>N$, we have that $\nabla_{n}u(x)\in\partial_{\epsilon}u(x),$ and in particular, $\nabla u(x)\in\partial_{\epsilon}u(x)$ for any $\epsilon>0.$ If we choose to run SGD with an $\epsilon-$subgradient $g^{k}$ in place of a usual subgradient, we call this an $\epsilon$-stochastic subgradient descent ($\epsilon$-SGD) algorithm.

We come now to the main result of this section. The result and proof are modeled after Theorem 14.8 of [41].

Note that Theorem 2.9 shows that the estimate of $\mathbb{E}[u(\bar{x})]-u(x^{*})$ is composed of two parts: $\frac{BM}{\sqrt{K}}$ corresponds to the usual SGD estimate as in [41] while the $\epsilon$ term corresponds to the effects of nonlocal interactions. We can now specialize the general $\epsilon-SGD$ theory above to the case of the nonlocal gradient operator. Recall our observation that given an $\epsilon>0$, we have $\nabla_{n}u(x)\in\partial_{\epsilon}u(x)$ for $n$ large enough. Thus, we can set $g^{k}=k_{u}(x^{k},y)$ to be the $\epsilon-$ subgradients used in $\epsilon-SGD$ and choose $p_{Y}(y)=\rho_{n}(x-y)p_{X}(x)$ to be the distribution from which the samples $y$ are chosen. We then have a nonlocal $\epsilon-SGD$ procedure.