Hopfield Neural Network Flow: A Geometric Viewpoint

TL;DR

This paper offers a geometric perspective on Hopfield neural networks by interpreting their dynamics as gradient flows on Riemannian manifolds, linking activation functions to underlying geometry and enabling efficient computation methods.

Contribution

It establishes a connection between deterministic HNNs and natural gradient descent, and characterizes stochastic HNNs as gradient flows in probability measure spaces, revealing new geometric insights.

Findings

Deterministic HNNs can be viewed as natural gradient descent.

Stochastic HNNs induce gradient flows in Wasserstein space.

Gradient flow interpretation enables fast proximal algorithms.

Abstract

We provide gradient flow interpretations for the continuous-time continuous-state Hopfield neural network (HNN). The ordinary and stochastic differential equations associated with the HNN were introduced in the literature as analog optimizers, and were reported to exhibit good performance in numerical experiments. In this work, we point out that the deterministic HNN can be transcribed into Amari's natural gradient descent, and thereby uncover the explicit relation between the underlying Riemannian metric and the activation functions. By exploiting an equivalence between the natural gradient descent and the mirror descent, we show how the choice of activation function governs the geometry of the HNN dynamics. For the stochastic HNN, we show that the so-called "diffusion machine", while not a gradient flow itself, induces a gradient flow when lifted in the space of probability measures. We characterize this infinite dimensional flow as the gradient descent of certain free energy with respect to a Wasserstein metric that depends on the geodesic distance on the ground manifold. Furthermore, we demonstrate how this gradient flow interpretation can be used for fast computation via recently developed proximal algorithms.

Tables

Table 1. TABLE I : Comparison between the natural (finite dimensional) gradient descent for the deterministic HNN and the measure-valued (infinite dimensional) gradient descent for the stochastic HNN. The graphical illustrations in the first row show that the finite dimensional gradient flow for the deterministic HNN evolves on ℳ ℳ \mathcal{M} , and is a gradient descent of f 𝑓 f w.r.t. the distance d 𝑮 subscript 𝑑 𝑮 d_{\bm{G}} . On the other hand, the infinite dimensional gradient flow for the stochastic HNN evolves on 𝒫 2 ​ ( ℳ ) subscript 𝒫 2 ℳ \mathcal{P}_{2}(\mathcal{M}) , and is a gradient descent of F 𝐹 F w.r.t. the distance W 𝑮 subscript 𝑊 𝑮 W_{\bm{G}} . Notice that for stochastic HNN, the sample paths of the SDE ( 45 ) in ℳ ℳ \mathcal{M} , induce a flow of PDFs in 𝒫 2 ​ ( ℳ ) subscript 𝒫 2 ℳ \mathcal{P}_{2}(\mathcal{M}) .

Attributes	Gradient Flow in Deterministic HNN	Gradient Flow in Stochastic HNN
Graphical illustration
Variational problem	$\underset{𝒙\in ℳ}{inf}f\left(𝒙\right)$	$\underset{\rho \in {𝒫}_2\left(ℳ\right)}{inf}F\left(\rho \right)$
Objective	$f\left(𝒙\right)$	$F\left(\rho \right)$ given by (51)
Distance	$d_{𝑮}$ given by (13)	$W_{𝑮}$ given by (53)
Differential equation	$\frac{d𝒙}{dt}=-{\left(𝑮\left(𝒙\right)\right)}^{-1}\nabla f$	$\frac{\partial \rho }{\partial t}=\nabla \cdot \left(\rho {\left(𝑮\left(𝒙\right)\right)}^{-1}\frac{\delta F}{\delta \rho }\right)$
	$:=-{\nabla }^{d_{𝑮}}f\left(𝒙\right)$	$:=-{\nabla }^{W_{𝑮}}F\left(\rho \right)$
Algebraic recursion	${𝒙}_k={𝒙}_{k-1}-{h{\nabla }^{d_{𝑮}}f\left(𝒙\right)|}_{𝒙={𝒙}_{k-1}}$	${\rho }_k={\rho }_{k-1}-{h{\nabla }^{W_{𝑮}}F\left(\rho \right)|}_{\rho ={\rho }_{k-1}}$
Proximal recursion	${𝒙}_k=\underset{𝒙\in ℳ}{\arg inf}{d}_{𝑮}^2(𝒙,{𝒙}_{k-1})+hf\left(𝒙\right)$	${\rho }_k=\underset{\rho \in {𝒫}_2\left(ℳ\right)}{\arg inf}{W}_{𝑮}^2(\rho ,{\rho }_{k-1})+hF\left(\rho \right)$
	$:={\mathrm{prox}}_{hf}^{d_{𝑮}}\left({𝒙}_{k-1}\right)$	$:={\mathrm{prox}}_{hF}^{W_{𝑮}}\left({\rho }_{k-1}\right)$

Equations

$${\dot{\mathbf{x}}}_{\mathrm{H}}$$

$$\mathbf{x}$$

$${\left({\sigma }_1({x_{\mathrm{H}}}_1),{\sigma }_2({x_{\mathrm{H}}}_2),\dots ,{\sigma }_n({x_{\mathrm{H}}}_n)\right)}^{\top },$$

$$\displaystyle\dot{\bm{x}}=-\left(\nabla_{\bm{x}_{\rm{H}}}\bm{\sigma}\right)\big{|}_{\bm{x}_{\rm{H}}=\bm{\sigma}^{-1}(\bm{x})}\nabla f(\bm{x}).$$

$$g_{ij}=\{\begin{array}{ll}1/{\sigma }_i^{\prime }({\sigma }_i^{-1}(x_i))& \text{for}i=j,\\ 0& \text{otherwise,}\end{array}$$

$$\displaystyle\bm{x}(k+1)=\bm{x}(k)-h\left(\bm{G}(\bm{x})\right)^{-1}\nabla f(\bm{x})\big{|}_{\bm{x}=\bm{x}(k)},$$

$$\dot{\mathbf{x}}=-{\left(\mathbf{G}(\mathbf{x})\right)}^{-1}\nabla f(\mathbf{x}).$$

$$\mathrm{d}{s}^2=\sum _{i=1}^n\frac{\mathrm{d}{x}_i^2}{{\sigma }_i^{\prime }({\sigma }_i^{-1}(x_i))},$$

$${\ddot{\gamma }}_k(t)+\sum _{i,j=1}^n{\Gamma }_{ij}^k(\mathbf{\gamma }(t)){\dot{\gamma }}_i(t){\dot{\gamma }}_j(t)=0,k=1,\dots ,n,$$

$${\Gamma }_{ij}^k:=\frac{1}{2}\sum _{l=1}^n{g}^{kl}\left(\frac{\partial g_{jl}}{\partial x_i}+\frac{\partial g_{il}}{\partial x_j}-\frac{\partial g_{ij}}{\partial x_l}\right),$$

$${\Gamma }_{ij}^k=0,$$

$${\Gamma }_{ii}^k=-\frac{1}{2}{g}^{kk}\frac{\partial g_{ii}}{\partial x_k},$$

$${\Gamma }_{ik}^k=\frac{\partial }{\partial x_i}\left(\log \sqrt{|g_{kk}|}\right),$$

$${\Gamma }_{kk}^k=\frac{\partial }{\partial x_k}\left(\log \sqrt{|g_{kk}|}\right).$$

$${\ddot{\gamma }}_i(t)+{\Gamma }_{ii}^i(\mathbf{\gamma }(t)){\left({\dot{\gamma }}_i(t)\right)}^2=0,i=1,\dots ,n.$$

$${\ddot{\gamma }}_i(t)+\frac{{\sigma }_i^{\prime \prime }\left({\sigma }_i^{-1}\left({\gamma }_i\right)\right){\sigma }_i^{\prime }\left({\gamma }_i\right)}{2{\sigma }_i^{\prime }\left({\sigma }_i^{-1}\left({\gamma }_i\right)\right){\left({\sigma }_i\left({\gamma }_i\right)\right)}^2}{\left({\dot{\gamma }}_i(t)\right)}^2=0,$$

$$d_{\mathbf{G}}(\mathbf{x},\mathbf{y})$$

$$={\int }_0^1{\left(\sum _{i=1}^n\frac{{\left({\dot{\gamma }}_i(t)\right)}^2}{{\sigma }_i^{\prime }({\sigma }_i^{-1}({\gamma }_i(t)))}\right)}^{1/2}\mathrm{d}t.$$

$$(\nabla \psi {)}^{-1}:\mathcal{M}\mapsto \mathcal{N}.$$

$$D_{\psi }\left(\mathbf{z},\tilde{\mathbf{z}}\right):=\psi (\mathbf{z})-\psi (\tilde{\mathbf{z}})-(\mathbf{z}-\tilde{\mathbf{z}}{)}^{\top }\nabla \psi (\tilde{\mathbf{z}}).$$

$$\nabla \psi (\mathbf{y}(k+1))$$

$$\mathbf{z}(k+1)$$

$${\mathrm{proj}}_{\mathcal{Z}}^{D_{\psi }}\left(\mathbf{\eta }\right):=\underset{\mathbf{\xi }\in \mathcal{Z}}{\mathrm{argmin}}{D}_{\psi }\left(\mathbf{\xi },\mathbf{\eta }\right).$$

$${\psi }^{\ast }\left(\mathbf{x}\right):=\underset{\mathbf{z}}{\sup }\left({\mathbf{x}}^{\top }\mathbf{z}-\psi (\mathbf{z})\right),$$

$$\nabla {\psi }^{\ast }={\left(\nabla \psi \right)}^{-1},$$

$$D_{\psi }\left(\mathbf{z},\tilde{\mathbf{z}}\right)=D_{{\psi }^{\ast }}\left(\nabla \psi (\tilde{\mathbf{z}}),\nabla \psi (\mathbf{z})\right).$$

$$z_i$$

$$=\frac{\partial {\psi }^{\ast }}{\partial x_i}=\underset{=:\varphi (x_i)}{\underbrace{\int \frac{\mathrm{d}{x}_i}{{\sigma }_i^{\prime }({\sigma }_i^{-1}(x_i))}}}+\underset{\begin{array}{c}\text{arbitrary function that}\\ \text{does not depend on xi}\end{array}}{\underbrace{k_i\left({\mathbf{x}}_{-i}\right)}},$$

$${\psi }^{\ast }(\mathbf{x})=\int \varphi (x_i)\mathrm{d}{x}_i+\kappa \prod _{i=1}^n{x}_i+c,$$

$${\dot{\mathbf{x}}}_{\mathrm{H}}=-\nabla f(\mathbf{x}),\mathbf{x}=\nabla {\psi }^{\ast }\left({\mathbf{x}}_{\mathrm{H}}\right),$$

$${\sigma }_i(x_{{\mathrm{H}}_i}):=\frac{1}{2}\tanh \left({\beta }_i\left(x_{{\mathrm{H}}_i}-\frac{1}{2}\right)\right)+\frac{1}{2},{\beta }_i>0;$$

$${\sigma }_i^{\prime }({\sigma }_i^{-1}(x_i))=\frac{1}{2}{\beta }_i{\mathrm{sech}}^2\left({\tanh }^{-1}(2x_i-1)\right),$$

$$({\nabla }^2{\psi }^{\ast }{)}_{ii}=g_{ii}=1/{\sigma }_i^{\prime }({\sigma }_i^{-1}(x_i))$$

$$=\frac{2}{{\beta }_i}{\cosh }^2\left({\tanh }^{-1}(2x_i-1)\right)=\frac{2}{{\beta }_i}\frac{1}{1-(1-2x_i{)}^2},$$

$${\ddot{\gamma }}_i+\frac{2{\gamma }_i-1}{2{\gamma }_i\left(1-{\gamma }_i\right)}{\left({\dot{\gamma }}_i\right)}^2=0,$$

Full text

For the stochastic HNN, we show that the so-called “diffusion machine”, while not a gradient flow itself, induces a gradient flow when lifted in the space of probability measures. We characterize this infinite dimensional flow as the gradient descent of certain free energy with respect to a Wasserstein metric that depends on the geodesic distance on the ground manifold. Furthermore, we demonstrate how this gradient flow interpretation can be used for fast computation via recently developed proximal algorithms.

Keywords: Natural gradient descent, mirror descent, proximal operator, stochastic neural network, optimal mass transport, Wasserstein metric.

I Introduction

Computational models such as the Hopfield neural networks (HNN) [1, 2] have long been motivated as analog machines to solve optimization problems. Both deterministic and stochastic differential equations for the HNN have appeared in the literature, and the numerical experiments using these models have reported appealing performance [3, 4, 5, 6, 7]. In [8], introducing the stochastic HNN as the “diffusion machine”, Wong compared the same with other stochastic optimization algorithms, and commented: “As such it is closely related to both the Boltzmann machine and the Langevin algorithm, but superior to both in terms of possible integrated circuit realization”, and that “A substantial speed advantage for diffusion machines is likely”. Despite these promising features often backed up via numerical simulations, studies seeking basic understanding of the evolution equations remain sparse [9, 10]. The purpose of this paper is to revisit the differential equations for the HNN from a geometric standpoint.

Our main contribution is to clarify the interplay between the Riemannian geometry induced by the choice of HNN activation functions, and the associated variational interpretations for the HNN flow. Specifically, in the deterministic case, we argue that the HNN flow is natural gradient descent (Section II) of a function over a manifold $\mathcal{M}$ with respect to (w.r.t.) a distance induced by certain metric tensor $\bm{G}$ that depends on the activation functions of the HNN. This identification leads to an equivalent mirror descent interpretation (Section III). We work out an example and a case study (Section IV) to highlight these ideas.

In the stochastic HNN, the random fluctuations resulting from the noise do not allow the sample paths of the stochastic differential equation to be interpreted as gradient flow on $\mathcal{M}$. However, we argue (Section V) that the stochastic sample path evolution engenders a (deterministic) flow of probability measures supported on $\mathcal{M}$ which does have an infinite dimensional gradient descent interpretation. In other words, the gradient flow interpretation for the stochastic HNN holds in the macroscopic or ensemble sense. This viewpoint seems novel in the context of stochastic HNN, and indeed leads to proximal algorithms enabling fast computation (Section V-C). A comparative summary of the gradient flow interpretations for the deterministic and the stochastic HNN is provided in Table I (after Section V-C).

Notations

The notation $\nabla$ stands for the Euclidean gradient operator. We sometimes put a subscript as in $\nabla_{\bm{q}}$ to denote the gradient w.r.t. vector $\bm{q}$, and omit the subscript when its meaning is obvious from the context. As usual, $\nabla$ acting on vector-valued function returns the Jacobian. The symbol $\nabla^{2}$ is used for the Hessian. We use the superscript ⊤ to denote matrix transposition, and $\mathbb{E}_{\rho}\left[\cdot\right]$ to denote the expectation operator w.r.t. the probability density function (PDF) $\rho$, i.e., $\mathbb{E}_{\rho}\left[\cdot\right]:=\int_{\mathcal{M}}(\cdot)\rho(\bm{x})\>{\rm{d}}\bm{x}$. In the development that follows, we assume that all probability measures are absolutely continuous, i.e., the corresponding PDFs exist. The notation $\mathcal{P}_{2}\left(\mathcal{M}\right)$ stands for the set of all joint PDFs supported on $\mathcal{M}$ with finite second raw moments, i.e., $\mathcal{P}_{2}\left(\mathcal{M}\right):=\{\rho:\mathcal{M}\mapsto\mathbb{R}_{\geq 0}\mid\int_{\mathcal{M}}\rho\>{\rm{d}}\bm{x}=1,\mathbb{E}_{\rho}\left[\bm{x}^{\top}\bm{x}\right]<\infty\}$. The inequality $\bm{A}\succ\bm{0}$ means that the matrix $\bm{A}$ is symmetric positive definite. The symbol $\bm{I}$ denotes an identity matrix of appropriate dimension; $\bm{1}$ denotes a column vector of ones. For $\bm{x}\in\mathcal{M}$, the symbol $\mathcal{T}_{\bm{x}}\mathcal{M}$ denotes the tangent space of $\mathcal{M}$ at $\bm{x}$; we use $\mathcal{T}\mathcal{M}$ to denote the tangent space at a generic point in $\mathcal{M}$. The set of natural numbers is denoted as $\mathbb{N}$. We use the symbols $\odot$ and $\oslash$ for denoting element-wise (i.e., Hadamard) product and division, respectively.

II HNN Flow as Natural Gradient Descent

We consider the deterministic HNN dynamics given by

[TABLE]

with initial condition $\bm{x}(t=0)\in(0,1)^{n}$. Here, the function $f:\mathbb{R}^{n}\mapsto\mathbb{R}$ is continuously differentiable, and therefore admits a minimum on $[0,1]^{n}$. Without loss of generality, we assume $f$ to be non-negative on $[0,1]^{n}$. The vectors $\bm{x},\bm{x}_{\rm{H}}\in\mathbb{R}^{n}$ are referred to as the state and hidden state, respectively. The so called “activation function” $\bm{\sigma}:\mathbb{R}^{n}\mapsto[0,1]^{n}$ is assumed to satisfy the following properties:

$\bullet$

homeomorphism,

$\bullet$

differentiable almost everywhere,

$\bullet$

has structure

[TABLE]

i.e., its $i$-th component depends only on the $i$-th component of the argument,

$\bullet$

$\sigma_{i}(\cdot)$ for all $i=1,\ldots,n$ are strictly increasing.

Examples of activation function include the logistic function, hyperbolic tangent function, among others.

To view the continuous-time dynamics (1) from a geometric standpoint, we start by rewriting it as

[TABLE]

In words, the flow of $\bm{x}(t)$ is generated by a vector field that is negative of the Jacobian of $\bm{\sigma}$ evaluated at $\bm{x}_{\rm{H}}=\bm{\sigma}^{-1}(\bm{x})$, times the gradient of $f$ w.r.t. $\bm{x}$. Due to (2), the Jacobian in (3) is a diagonal matrix. Furthermore, since each $\sigma_{i}$ is strictly increasing and differentiable almost everywhere, the diagonal terms are positive. Therefore, the right-hand-side (RHS) of (3) has the following form: negative of a positive definite matrix times the gradient of $f$. This leads to the result below.

Theorem 1**.**

The flow $\bm{x}(t)$ governed by (3) (equivalently by (1)) is a natural gradient descent for function $f(\bm{x})$ on the manifold $\mathcal{M}=(0,1)^{n}$ with metric tensor $\bm{G}=[g_{ij}]_{i,j=1}^{n}$, given by

[TABLE]

where $i,j=1,\ldots,n$, and ′ denotes derivative.

Proof.

Amari’s natural gradient descent [11] describes the steepest descent of $f(\bm{x})$ on a Riemannian manifold $\mathcal{M}$ with metric tensor $\bm{G}\succ\bm{0}$, and is given by the recursion

[TABLE]

where the discrete iteration index $k\in\mathbb{N}$, and the step-size $h>0$ is small. For $h\downarrow 0$, we get the natural gradient dynamics

[TABLE]

Notice that due to the assumptions listed on the activation function $\bm{\sigma}$, its Jacobian appearing in (3) is a diagonal matrix with positive entries. Therefore, (3) is of the form (6) with $\bm{G}(\bm{x})={\rm{diag}}(1/\sigma^{\prime}_{i}(\sigma_{i}^{-1}(x_{i})))$, $i=1,\ldots,n$. In particular, $\bm{G}(\bm{x})\succ\bm{0}$ for all $\bm{x}\in(0,1)^{n}$. This completes the proof. $\blacksquare$

Theorem 1 allows us to interpret the flow generated by (3) as the steepest descent of $f$ in a Riemannian manifold $\mathcal{M}$ with length element ${\rm{d}}s$ given by

[TABLE]

i.e., an orthogonal coordinate system can be associated with $\mathcal{M}$. Further insights can be obtained by interpreting $\mathcal{M}$ as a Hessian manifold [12]. To this end, we would like to view the positive definite metric tensor $\bm{G}={\rm{diag}}(1/\sigma^{\prime}_{i}(\sigma_{i}^{-1}(x_{i})))$ as the Hessian of a (twice differentiable) strictly convex function. Such an identification will allow us to associate a mirror descent [13, 14] in the dual manifold [15] corresponding to the natural gradient descent in manifold $\mathcal{M}$. Before delving into mirror descent, we next point out that (7) allows us to compute geodesics on the manifold $\mathcal{M}$, which completes the natural gradient descent interpretation.

II-A Geodesics

The (minimal) geodesic curve $\bm{\gamma}(t)$ parameterized via $t\in[0,1]$ that connects $\bm{x},\bm{y}\in(0,1)^{n}$, solves the well-known Euler-Lagrange equation (see e.g., [16, Ch. 3.2])

[TABLE]

subject to the boundary conditions $\bm{\gamma}(0)=\bm{x}$, $\bm{\gamma}(1)=\bm{y}$. Here $\Gamma_{ij}^{k}$ denote the Cristoffel symbols of the second kind for $i,j,k=1,\ldots,n$, and are given by

[TABLE]

wherein $g^{kl}$ denotes the $(k,l)$-th element of the inverse metric tensor $\bm{G}^{-1}$. For $i\neq j\neq k$, and diagonal $\bm{G}$, we get

[TABLE]

For our diagonal metric tensor (4), since the $i$-th diagonal element of $\bm{G}$ only depends on $x_{i}$, the terms (10a)–(10c) are all zero, i.e., the only non-zero Cristoffel symbols in (10) are (10d). This simplifies (8) as a set of $n$ decoupled ODEs:

[TABLE]

Using (4) and (10d), for $i=1,\ldots,n$, we can rewrite (11) as

[TABLE]

which solved together with the endpoint conditions $\bm{\gamma}(0)=\bm{x}$, $\bm{\gamma}(1)=\bm{y}$, yields the geodesic curve $\bm{\gamma}(t)$. In Section IV, we will see an explicitly solvable example for the system of nonlinear ODEs (12).

Once the geodesic $\bm{\gamma}(t)$ is obtained from (12), the geodesic distance $d_{\bm{G}}(\bm{x},\bm{y})$ induced by the metric tensor $\bm{G}$, can be computed as

[TABLE]

This allows us to formally conclude that the flow generated by (3), or equivalently by (6), can be seen as the gradient descent of $f$ on the manifold $\mathcal{M}$, measured w.r.t. the distance $d_{\bm{G}}$ given by (13b).

Remark 1**.**

In view of the notation $d_{\bm{G}}$, we can define the Riemannian gradient $\nabla^{d_{\bm{G}}}\left(\cdot\right):=(\bm{G}(\bm{x}))^{-1}\nabla\left(\cdot\right)$, and rewrite the dynamics (6) as $\dot{\bm{x}}=-\nabla^{d_{\bm{G}}}f(\bm{x})$.

III HNN Flow as Mirror Descent

We now provide an alternative way to interpret the flow generated by (3) as mirror descent on a Hessian manifold $\mathcal{N}$ that is dual to $\mathcal{M}$. For $\bm{x}\in\mathcal{M}$ and $\bm{z}\in\mathcal{N}$, consider the mapping $\mathcal{N}\mapsto\mathcal{M}$ given by the map $\bm{x}=\nabla_{\bm{z}}\psi(\bm{z})$ for some strictly convex and differentiable map $\psi(\cdot)$, i.e.,

[TABLE]

The map $\psi(\cdot)$ is referred to as the “mirror map”, formally defined below. In the following, the notation ${\rm{dom}}(\psi)$ stands for the domain of $\psi(\cdot)$.

Definition 1**.**

(Mirror map) A differentiable, strictly convex function $\psi:{\rm{dom}}(\psi)\mapsto\mathbb{R}$, on an open convex set ${\rm{dom}}(\psi)\subseteq\mathbb{R}^{n}$, satisfying $\parallel\nabla\psi\parallel_{2}\rightarrow\infty$ as the argument of $\psi$ approaches the boundary of the closure of ${\rm{dom}}(\psi)$, is called a mirror map.

Given a mirror map, one can introduce a Bregman divergence [17] associated with it as follows.

Definition 2**.**

(Bregman divergence associated with a mirror map) For a mirror map $\psi$, the associated Bregman divergence $D_{\psi}:{\rm{dom}}(\psi)\times{\rm{dom}}(\psi)\mapsto\mathbb{R}_{\geq 0}$, is given by

[TABLE]

Geometrically, $D_{\psi}\left(\bm{z},\widetilde{\bm{z}}\right)$ is the error at $\bm{z}$ due to first order Taylor approximation of $\psi$ about $\widetilde{\bm{z}}$. It is easy to verify that $D_{\psi}\left(\bm{z},\widetilde{\bm{z}}\right)=0$ iff $\bm{z}=\widetilde{\bm{z}}$. However, the Bregman divergence is not symmetric in general, i.e., $D_{\psi}\left(\bm{z},\widetilde{\bm{z}}\right)\neq D_{\psi}\left(\widetilde{\bm{z}},\bm{z}\right)$.

In the sequel, we assume that the mirror map $\psi$ is twice differentiable. Then, $\psi$ being strictly convex, its Hessian $\nabla^{2}\psi$ is positive definite. This allows one to think of $\mathcal{N}\equiv{\rm{dom}}(\psi)$ as a Hessian manifold with metric tensor $\bm{H}\equiv\nabla^{2}\psi$. The mirror descent for function $f$ on a convex set $\mathcal{Z}\subset\mathcal{N}$, is a recursion of the form

[TABLE]

where $h>0$ is the step-size, and the Bregman projection

[TABLE]

For background on mirror descent, we refer the readers to [13, 14]. For a recent reference, see [18, Ch. 4].

Raskutti and Mukherjee [15] pointed out that if $\psi$ is twice differentiable, then the map $\mathcal{N}\mapsto\mathcal{M}$ given by $\bm{x}=\nabla_{\bm{z}}\psi(\bm{z})$, establishes a one-to-one correspondence between the natural gradient descent (5) and the mirror descent (15). Specifically, the mirror descent (15) on $\mathcal{Z}\subset\mathcal{N}$ associated with the mirror map $\psi$, is equivalent to the natural gradient descent (5) on $\mathcal{M}$ with metric tensor $\nabla^{2}\psi^{*}$, where

[TABLE]

is the Legendre-Fenchel conjugate [19, Section 12] of $\psi(\bm{z})$. In our context, the Hessian $\nabla^{2}\psi^{*}\equiv\bm{G}$ is given by (4). To proceed further, the following Lemma will be useful.

Lemma 1**.**

(Duality) Consider a mirror map $\psi$ and its associated Bregman divergence $D_{\psi}$ as in Definitions 1 and 2. It holds that

[TABLE]

Proof.

See for example, [22, Section 2.2]. $\blacksquare$

Combining (4) with (18a), we get

[TABLE]

for $i=1,\ldots,n$, where the notation $k_{i}\left(\bm{x}_{-i}\right)$ means that the function $k_{i}(\cdot)$ does not depend on the $i$-th component of the vector $\bm{x}$. Therefore, the map $\mathcal{M}\mapsto\mathcal{N}$ associated with the HNN flow is of the form (19).

To derive the mirror descent (15) associated with the HNN flow, it remains to determine $D_{\psi}$. Thanks to (18b), we can do so without explicitly computing $\psi$. Specifically, integrating (19) yields $\psi^{*}$, i.e.,

[TABLE]

where $\phi(\cdot)$ is defined in (19), and $\kappa,c$ are constants. The constants $\kappa,c$ define a parametrized family of maps $\psi^{*}$. Combining (14) and (20) yields $D_{\psi^{*}}$, and hence $D_{\psi}$ (due to (18b)). In Section IV, we illustrate these ideas on a concrete example.

Remark 2**.**

We note that one may rewrite (1) in the form of a mirror descent ODE [20, Section 2.1], given by

[TABLE]

for some to-be-determined mirror map $\psi$, thereby interpreting (1) as mirror descent in $\mathcal{M}$ itself, with $\bm{x}_{\rm{H}}\in\mathbb{R}^{n}$ as the dual variable, and $\bm{x}\in[0,1]^{n}$ as the primal variable. This interpretation, however, is contingent on the assumption that the monotone vector function $\bm{\sigma}$ is expressible as the gradient of a convex function (here, $\psi^{*}$), for which the necessary and sufficient condition is that the map $\bm{\sigma}$ be maximally cyclically monotone [21, Section 2, Corollary 1 and 2]. This indeed holds under the structural assumption (2).

IV Examples

In this Section, we first give a simple example (Section IV-A) to elucidate the ideas presented in Sections II and III. Then, we provide a case study (Section IV-B) to exemplify their application in the context of an engineering problem.

IV-A An Illustrative Example

Consider the activation function $\bm{\sigma}:\mathbb{R}^{n}\mapsto[0,1]^{n}$ given by the soft-projection operator

[TABLE]

see Fig. 1. We next illustrate the natural gradient descent and the mirror descent interpretations of (3) with activation function (22).

IV-A1 Natural Gradient Descent Interpretation

Direct calculations yield

[TABLE]

for all $i=1,\ldots,n$. Notice that $x_{i}\in[0,1]$ implies $(2x_{i}-1)\in[-1,1]$, and hence the RHS of (23) is positive. We have

[TABLE]

which allows us to deduce the following: the HNN flow with activation function (22) is the natural gradient descent (5) with diagonal metric tensor $\bm{G}$ given by (24), i.e., $g_{ii}=1/\left(2\beta_{i}x_{i}(1-x_{i})\right)$, $i=1,\ldots,n$.

In this case, (12) reduces to the nonlinear ODE

[TABLE]

which solved with the boundary conditions $\gamma_{i}(0)=x_{i}$, $\gamma_{i}(1)=y_{i}$, $i=1,\ldots,n$, determines the geodesic curve $\bm{\gamma}(t)$ connecting $\bm{x},\bm{y}\in(0,1)^{n}$. Introducing the change-of-variable $v_{i}:=\dot{\gamma}_{i}$, noting that $\ddot{\gamma}_{i}=v_{i}\dfrac{{\rm{d}}v_{i}}{{\rm{d}}\gamma_{i}}$, and enforcing $v_{i}\not\equiv 0$, we can transcribe (25) into the separable ODE

[TABLE]

which gives

[TABLE]

and consequently,

[TABLE]

where $\{\left(a_{i},b_{i}\right)\}_{i=1}^{n}$ are constants of integration. Using the boundary conditions $\gamma_{i}(0)=x_{i}$, $\gamma_{i}(1)=y_{i}$ in (28) yields the geodesic curve $\bm{\gamma}(t)$ as

[TABLE]

where $i=1,\ldots,n$, and $t\in[0,1]$. From (29), it is easy to verify that $\bm{\gamma}(t)$ is component-wise in $[0,1]$, and satisfies the boundary conditions $\bm{\gamma}(0)=\bm{x}$, $\bm{\gamma}(1)=\bm{y}$. From (13b) and (29), the geodesic distance associated with the activation function (22) is

[TABLE]

where all vector operands such as square-root, arcsin, division (denoted by $\oslash$), are element-wise. The contour plots in Fig. 2 show the geodesic balls of $d_{\bm{G}}$ given by (30), centered at $(0.5,0.5)$. We summarize: the HNN flow generated by (3) with activation function (22) is natural gradient descent of $f$ measured w.r.t. the geodesic distance (30).

IV-A2 Mirror Descent Interpretation

To derive a mirror descent interpretation, integrating (24) we obtain

[TABLE]

where $i=1,\ldots,n$, and as before, $k_{i}\left(\bm{x}_{-i}\right)$ means that the function $k_{i}(\cdot)$ does not depend on the $i$-th component of the vector $\bm{x}$. Equation (31) implies that

[TABLE]

where $\kappa,c$ are constants. Therefore, (31) can be re-written as

[TABLE]

Comparing the RHS of (33) with (19), for this specific example, we have that

[TABLE]

Let us consider a special case of (32), namely $\kappa=c=0$, that is particularly insightful. In this case,

[TABLE]

i.e., a separable (weighted) sum of the bit entropy (see [22, Table 1] and [23, Table 1]). The associated Bregman divergence

[TABLE]

can be interpreted as the weighted sum of logistic loss functions associated with the Bernoulli random variables. Furthermore, $z_{i}=\log(x_{i}/(1-x_{i}))/2\beta_{i}$ yields

[TABLE]

Consequently, using (18b) and (36), we get

[TABLE]

which is the dual Logistic loss [22, Table 1], and the corresponding mirror map can be identified as

[TABLE]

Thus, an alternative interpretation of the HNN dynamics with activation function (22) is as follows: it can be seen as mirror descent of the form (15) on $\mathcal{Z}\equiv\mathbb{R}^{n}$ in variables $z_{i}$ given by (33); a candidate mirror map is given by (38) with associated Bregman divergence (37).

IV-B Case Study: Economic Load Dispatch Problem

We now apply the HNN gradient descent to solve an economic load dispatch problem in power system – a context in which HNN methods have been used before [7, 24, 25, 26]. Our objective, however, is to illustrate the geometric ideas presented herein.

Specifically, we consider $n_{\text{G}}$ generators which at time $t=0$ are either ON or OFF, denoted via state vector $\bm{x}_{0}\in\{0,1\}^{n_{\text{G}}}$. For $i=1,\ldots,n_{\text{G}}$, if the $i$th generator is ON (resp. OFF) at $t=0$, then the $i$th component of $\bm{x}_{0}$ equals 1 (resp. 0). Let $\bm{y}_{0}\in[0,1]^{n_{\text{G}}}$ be the vector of the generators’ power output magnitudes at $t=0$. A dispatcher simultaneously decides for the next time step $t=1$, which of the generators should be ON or OFF ($\bm{x}\in\{0,1\}^{n_{\text{G}}}$) and the corresponding power output ($\bm{y}\in[0,1]^{n_{\text{G}}}$) in order to satisfy a known power demand $\pi_{d}>0$ while minimizing a convex cost of the form

[TABLE]

Here, $\bm{p}$ denotes the (known) marginal production cost of the generators; the second summand in (39) models the symmetric cost for switching the generator states from ON to OFF or from OFF to ON; the last summand in (39) models the generator ramp-up or ramp-down cost. The weights $c_{1},c_{2}>0$ are known to the dispatcher. The optimization problem the dispatcher needs to solve (over single time step) reads

[TABLE]

The constraint (40c) enforces that the supply equals the demand; the constraint (40b) ensures that only those generators which are ON, can produce electricity.

To solve the Mixed Integer Quadratically Constrained Quadratic Programming (MIQCQP) problem (40), we construct the augmented Lagrangian

[TABLE]

where the parameter $r>0$, and the dual variables $\lambda_{1},\lambda_{2}\in\mathbb{R}$. Then, the augmented Lagrange dual function

[TABLE]

For some initial guess of the pair $(\lambda_{1},\lambda_{2})$, we apply the HNN natural gradient descent (5) to perform the minimization (42), and then do the dual ascent updates

[TABLE]

where $h_{\text{dual}}>0$ is the step size for the dual updates. Specifically, for the iteration index $k=0,1,2,...$, we perform the recursions

[TABLE]

where (44a) is solved via Hopfield sub-iterations with step size $h_{\text{Hopfield}}>0$. In [7], this algorithm was referred to as the dual Hopfield method, and was shown to have superior numerical performance compared to other standard approaches such as semidefinite relaxation.

For numerical simulation, we implement (44) with $n_{\text{G}}=40$ generators, uniform random vector $\bm{p}\in[0,1]^{n_{\text{G}}}$, the pair $(c_{1},c_{2})$ uniform random in $(0,1]^{2}$, $r=1$, $h_{\text{Hopfield}}=10^{-2}$, and $h_{\text{dual}}=10^{-1}$. The vectors $\bm{x}_{0}$ and $\bm{y}_{0}$ in (39) are chosen as follows: $\bm{x}_{0}$ is chosen to be uniform random in $[0,1]^{n_{\text{G}}}$ and then element-wise rounded to the nearest integers; $\bm{y}_{0}$ is chosen uniform random in $[0,1]^{n_{\text{G}}}$ and then element-wise multiplied with $\bm{x}_{0}$. For $u$ chosen uniform random in $[0,1]$, we set $\pi_{d}=\left(1+0.1(u-0.5)\right)\bm{1}^{\top}\bm{y}_{0}$. The activation functions are as in (22) with $\bm{\beta}=\bm{1}$. The Hopfield sub-iterations were run until error tolerance $10^{-7}$ was achieved with maximum number of sub-iterations being $10^{4}$.

We fix the problem data generated as above, and perform 100 Monte Carlo simulations, each solving an instance of the dual Hopfield method for the same problem data with randomly chosen $(\lambda_{1}(0),\lambda_{2}(0))$. Fig. 3 shows the convergence trends for the joint vector $\bm{z}:=\left(\bm{x};\bm{y}\right)\in\mathbb{R}^{2n_{\text{G}}}$ associated with the Monte Carlo simulations by plotting the distance between the current iterate $\bm{z}(k)$ and the converged vector $\bm{z}^{*}$ over iteration index $k\in\mathbb{N}$. In particular, Fig. 3 highlights how the geometry induced by the HNN governs the convergence trend: the geodesic distance to minimum $d_{\bm{G}}\left(\bm{z}(k),\bm{z}^{*}\right)$ (blue curves in Fig. 3), where $d_{\bm{G}}$ is given by (30), decays monotonically, as expected in gradient descent. However, $\|\bm{z}(k)-\bm{z}^{*}\|_{2}$ (red curves in Fig. 3) does not decay monotonically since the Euclidean distance does not encode the Riemannian geometry of the HNN recursion (5).

V Stochastic HNN Flow As Gradient Descent In the Space of Probability Measures

We next consider the stochastic HNN flow, also known as the “diffusion machine” [8, Section 4], given by the Itô stochastic differential equation (SDE)

[TABLE]

where the state***We make the standard assumption that $g_{ii}^{-1}=\sigma_{i}^{\prime}\left(\sigma_{i}(x_{i})\right)=0$ at $x_{i}=0,1$. If the activation functions $\sigma_{i}(\cdot)$ do not satisfy this, then a reflecting boundary is needed at each $x_{i}=0,1$, to keep the sample paths within the unit hypercube; see e.g., [27]. $\bm{x}\in(0,1)^{n}$, the standard Wiener process $\bm{w}\in\mathbb{R}^{n}$, the notation $\bm{1}$ stands for the $n\times 1$ vector of ones, and the parameter $T>0$ denotes the thermodynamic/annealing temperature. Recall that the metric tensor $\bm{G}(\bm{x})$ given by (4) is diagonal, and hence its inverse and square roots are diagonal too with respective diagonal elements being the inverse and square roots of the original diagonal elements.

The diffusion machine (45) was proposed in [8, 28] to solve global optimization problems on the unit hypercube, and is a stochastic version of the deterministic HNN flow discussed in Sections II–IV. This can be seen as a generalization of the Langevin model [29, 30] for the Boltzmann machine [31]. Alternatively, one can think of the diffusion machine as the continuous-state continuous-time version of the HNN in which Gaussian white noise is injected at each node.

Let $\mathcal{L}_{\rm{FPK}}$ be the Fokker-Planck-Kolmogorov (FPK) operator [32] associated with (45) governing the evolution of the joint PDFs $\rho(\bm{x},t)$, i.e.,

[TABLE]

where the given initial joint PDF $\rho_{0}\in\mathcal{P}_{2}\left(\mathcal{M}\right)$. The drift and diffusion coefficients in (45) are tailored so that $\mathcal{L}_{\rm{FPK}}$ admits the stationary joint PDF

[TABLE]

i.e., a Gibbs PDF where $Z(T)$ is a normalizing constant (known as the “partition function”) to ensure $\int\rho_{\infty}{\rm{d}}\bm{x}=1$. This follows from noting that for (45), we have†††We remind the readers that the raised indices denote the elements of the inverse of the metric tensor, i.e., $\bm{G}^{-1}=[g^{ij}]_{i,j=1}^{n}$.

[TABLE]

and that $g^{ii}(x_{i})>0$ for each $x_{i}\in(0,1)$. From (47), the critical points of $f$ coincides with that of $\rho_{\infty}$, which is what allows to interpret (45) as an analog machine for globally optimizing the function $f$. In fact, for the FPK operator (48), the solution $\rho(\bm{x},t)$ of (46) enjoys an exponential rate of convergence [33, p. 1358-1359] to (47). We remark here that the idea of using SDEs for solving global optimization problems have also appeared in [34, 35].

In the context of diffusion machine, the key observation that part of the drift term can be canceled by part of the diffusion term in the FPK operator (48a), thereby guaranteeing that the stationary PDF associated with (48) is (47), was first pointed out by Wong [8]. While this renders the stationary solution of the FPK operator meaningful for globally minimizing $f$, it is not known if the transient PDFs generated by (46) admit a natural interpretation. We next argue that (46) is indeed a gradient flow in the space of probability measures supported on the manifold $\mathcal{M}$.

V-A Wasserstein Gradient Flow

Using (48b), we rewrite the FPK PDE (46) as

[TABLE]

where

[TABLE]

For $\rho\in\mathcal{P}_{2}\left(\mathcal{M}\right)$, consider the free energy functional $F(\rho)$ given by

[TABLE]

which is a sum of the potential energy $\int_{\mathcal{M}}f\rho\>{\rm{d}}\bm{x}$, and the internal energy $T\int_{\mathcal{M}}\rho\log\rho\>{\rm{d}}\bm{x}$. The free energy (51) serves as the Lyapunov functional for (49) since direct calculation yields (Appendix -A)

[TABLE]

along any solution $\rho(\bm{x},t)$ generated by (49), thanks to $\bm{G}(\bm{x})\succ\bm{0}$ for all $\bm{x}\in\mathcal{M}=(0,1)^{n}$. In particular, the RHS of (52) equals zero iff $\nabla\zeta=0$, i.e., at the stationary PDF (47). For $t<\infty$, the RHS of (52) is $<0$ for any transient PDF $\rho(\bm{x},t)$ solving (49).

To view (49) as gradient flow over $\mathcal{P}_{2}(\mathcal{M})$, we will need the notion of (quadratic) Wasserstein metric $W_{\bm{G}}$ between two probability measures $\mu$ and $\nu$ on $\mathcal{M}$, defined as

[TABLE]

where $\pi(\bm{x},\bm{y})$ denotes a joint probability measure supported on $\mathcal{M}\times\mathcal{M}$, and the geodesic distance $d_{\bm{G}}$ is given by (13). The infimum in (53) is taken over the set $\Pi_{2}\left(\mu,\nu\right)$, which we define as the set of all joint probability measures having finite second moment that are supported on $\mathcal{M}\times\mathcal{M}$, with prescribed $\bm{x}$-marginal $\mu$, and prescribed $\bm{y}$-marginal $\nu$. If $\mu$ and $\nu$ have respective PDFs $\rho_{\bm{x}}$ and $\rho_{\bm{y}}$, then we can use the notation $W_{\bm{G}}(\rho_{\bm{x}},\rho_{\bm{y}})$ in lieu of $W_{\bm{G}}(\mu,\nu)$. The subscript in (53) is indicative of its dependence on the ground Riemannian metric $\bm{G}(\cdot)$, and generalizes the Euclidean notion of Wasserstein metric [36], i.e., $\bm{G}\equiv\bm{I}$ case. See also [37, Section 3] and [38]. Following [36, Ch. 7] and using the positive definiteness of $\bm{G}$, it is easy to verify that (53) indeed defines a metric on $\mathcal{P}_{2}\left(\mathcal{M}\right)$, i.e., $W_{\bm{G}}$ is non-negative, zero iff $\mu=\nu$, symmetric in its arguments, and satisfies the triangle inequality.

The square of (53) is referred to as the “optimal transport cost” (see [39] for the Euclidean case) that quantifies the minimum amount of work needed to reshape the PDF $\rho_{\bm{x}}$ to $\rho_{\bm{y}}$ (or equivalently, $\rho_{\bm{y}}$ to $\rho_{\bm{x}}$). This can be formalized via the following dynamic variational formula [40, Corollary 2.5] for (53):

[TABLE]

The infimum above is taken over the PDF-vector field pairs $(\rho,\bm{u})\in\mathcal{P}_{2}(\mathcal{M})\times\mathcal{U}$, where $\mathcal{U}:=\{\bm{u}:\mathcal{M}\times[0,\infty)\mapsto\mathcal{TM}\}$. Recognizing that (54a) is the continuity/Liouville equation for the single integrator dynamics $\dot{\bm{x}}=\bm{u}$, one can interpret (54) as a fixed horizon “minimum weighted energy” stochastic control problem subject to endpoint PDF constraints (54b)-(54c).

Following [41] and [40, Section 1.2], we can endow $\mathcal{P}_{2}(\mathcal{M})$ with the metric $W_{\bm{G}}$, resulting in an infinite dimensional Riemannian structure. Specifically, the tangent space $\mathcal{T}_{\rho}\mathcal{P}_{2}(\mathcal{M})$ can be equipped with the inner product $L^{2}_{\bm{G}}$ given by

[TABLE]

for tangent vectors $\bm{u},\bm{v}$. Given a smooth curve $t\mapsto\rho(\bm{x},t)$ in $\mathcal{P}_{2}(\mathcal{M})$, the associated tangent vector $\bm{u}(\bm{x},t)$ is characterized as the vector field that solves (54). We refer the readers to [42, Ch. 13] for further details. In particular, the “Wasserstein gradient” [43, Ch. 8] of a functional $\Phi(\rho)$ at $\rho\in\mathcal{P}_{2}\left(\mathcal{M}\right)$, can be defined as

[TABLE]

where $\delta\Phi/\delta\rho$ denotes the functional derivative of $\Phi(\rho)$. Consequently, the FPK PDE (49) can be transcribed into the familiar gradient flow form:

[TABLE]

Furthermore, the Euler (i.e., first order in time) discretization of the LHS of (57) yields the familiar gradient descent form:

[TABLE]

where $\rho_{k}:=\rho(\bm{x},t=kh)$ for $k\in\mathbb{N}$. The functional recursion (58) evolves on the metric space $\left(\mathcal{P}_{2}\left(\mathcal{M}\right),W_{\bm{G}}\right)$. We will see next that (58) can be recast as an equivalent proximal recursion, which will prove helpful for computation.

Remark 3**.**

We clarify here that the FPK operator appearing in the RHS of (49), or equivalently in the RHS of (57), is not quite same as the Laplace-Beltrami operator [44]

[TABLE]

This is because the Itô SDE (45) was hand-crafted in [8] such that the associated stationary PDF $\rho_{\infty}$ becomes (47) w.r.t. the volume measure ${\rm{d}}\bm{x}$, which is desired from an optimization perspective since then, the local minima of $f$ would coincide with the location of the modes of $\rho_{\infty}$. In contrast, (59) admits stationary PDF $\exp(-f/T)$ w.r.t. the volume measure $\sqrt{{\mathrm{det}}(\bm{G}(\bm{x}))}\>{\rm{d}}\bm{x}$ (up to a normalization constant), and in general, would not correspond to the local minima of $f$.

V-B Proximal Recursion on $\mathcal{P}_{2}\left(\mathcal{M}\right)$

A consequence of identifying (49) as gradient flow in $\mathcal{P}_{2}(\mathcal{M})$ w.r.t. the metric $W_{\bm{G}}$ is that the solution $\rho(\bm{x},t)$ of (57) can be approximated [40] via variational recursion

[TABLE]

with the initial PDF $\varrho_{0}\equiv\rho_{0}(\bm{x})$, and small step-size $h>0$. In particular, $\varrho_{k}\rightarrow\rho(\bm{x},t=kh)$ in strong $L^{1}(\mathcal{M})$ sense as $h\downarrow 0$. That a discrete time-stepping procedure like (60) can approximate the solution of FPK PDE, was first proved in [45] for $\bm{G}\equiv\bm{I}$, and has since become a topic of burgeoning research; see e.g., [43, 46, 47].

The RHS in (60) is an infinite dimensional version of the Moreau-Yosida proximal operator [48, 49, 50], denoted as ${\mathrm{prox}}_{hF}^{W_{\bm{G}}}\left(\cdot\right)$, i.e., (60) can be written succinctly as

[TABLE]

and (47) is the fixed point of this recursion. Just like the proximal viewpoint of finite dimensional gradient descent, (60) can be taken as an alternative definition of gradient descent of $F(\cdot)$ on the manifold $\mathcal{P}_{2}\left(\mathcal{M}\right)$ w.r.t. the metric $W_{\bm{G}}$; see e.g., [51, 52, 53, 54].

The utilitarian value of (60) over (58) is as follows. While the algebraic recursion (58) involves first order differential operator (w.r.t. $\bm{x}$), the variational recursion (60) is zero-th order. In the optimization community, the proximal operator is known [55] to be amenable for large scale computation. In our context too, we will see in Section V-C that (58) allows scattered weighted point cloud-based computation avoiding function approximation or spatial discretization, which will otherwise be impossible had we pursued a direct numerical method for the algebraic recursion (58).

Notice from (51) that $F(\rho)$ is strictly convex in $\rho\in\mathcal{P}_{2}(\mathcal{M})$. Furthermore, the map $\rho\mapsto W^{2}_{\bm{G}}(\widetilde{\rho},\rho)$ is convex in $\rho$ for any $\widetilde{\rho}\in\mathcal{P}_{2}(\mathcal{M})$. Therefore, (61) is a strictly convex functional recursion wherein each step admits unique minimizer.

V-B1 Interpretation for small $T$

In the global optimization context, the parameter $T$ in the free energy (51) plays an important role. For $T\downarrow 0$, the SDE (45) reduces to the ODE (6), and (51) simplifies to $\mathbb{E}_{\rho}[f]$ with $\rho$ being the time-varying Dirac PDF along the solution of (6). In that case, $W_{\bm{G}}\equiv d_{\bm{G}}$, and hence (60) reduces to the proximal recursion for finite dimensional natural gradient descent (Section II):

[TABLE]

Furthermore, the FPK PDE (49) reduces to the Liouville PDE

[TABLE]

whose stationary PDF becomes an weighted sum of Diracs located at the stationary points of $f$. The stationary solution of (62) converges to this stationary PDF. Thus, (60) with small $T$ approximates the deterministic natural gradient descent.

V-B2 Interpretation for large $T$

For $T$ large, the second summand in (51) dominates, i.e., (60) reduces to maximum entropy descent in $\mathcal{P}_{2}(\mathcal{M})$ w.r.t. the metric $W_{\bm{G}}$. The resulting proximal recursion admits uniform PDF as the fixed point, which is indeed the $T\rightarrow\infty$ limit of the stationary PDF (47). Thus, (60) with large $T$ approximates pure diffusion.

V-C Proximal Algorithm

The recursion (60) is not only appealing from theoretical perspective in that it reveals the metric geometry of gradient flow, but also from an algorithmic perspective in that it opens up the possibility to compute the transient solutions of the PDE (49) in a scalable manner via convex optimization. Similar to the finite dimensional proximal algorithms [55], the recently introduced infinite dimensional proximal algorithms [53, 54] were shown to be amenable for large scale implementation. Next, we adapt the proximal algorithm in [54, Section III.B] to solve (60) for an example problem.

V-C1 Example for stochastic HNN

We consider the activation functions $\sigma_{i}(\cdot)$ as in Section IV, resulting in

[TABLE]

which in turn gives the following instance of (45):

[TABLE]

To compute the transient joint PDF flow associated with (64), we solve (60) with $W_{\bm{G}}$ given by (53), wherein $d_{\bm{G}}$ is given by (30).

For numerical simulation, we set $n=2$, $T=25$, $\beta_{i}=1/4$ for $i=1,2$, and

[TABLE]

where $(x_{1},x_{2})\in[0,1]^{2}$. The function $f$ in (65) is a re-scaled version of the so-called Himmelblau’s function [56], often used as a benchmark for non-convex optimization. As seen in Fig. 4, the function $f$ has four local minima, and the corresponding stationary PDF (47) has four modes (i.e., local maxima) at the arg mins of $f$.

We generate $N=500$ samples from the uniform initial joint PDF $\rho_{0}$ supported on $[0,1]^{2}$ in the form of the point cloud $\{\bm{x}_{0}^{i},\varrho_{0}^{i}\}_{i=1}^{N}$. Here, $\{\bm{x}_{0}^{i}\}_{i=1}^{N}$ denotes the location of the samples in $[0,1]^{2}$ while $\{\varrho_{0}^{i}\}_{i=1}^{N}$ denotes the values of the initial joint PDF evaluated at these locations. We apply the gradient flow algorithm in [54, Section III.B] to recursively compute the updated point clouds $\{\bm{x}_{k}^{i},\varrho_{k}^{i}\}_{i=1}^{N}$ for $k\in\mathbb{N}$. Specifically, we update the sample locations $\{\bm{x}_{k}^{i}\}_{i=1}^{N}$ via the Euler-Maruyama‡‡‡Instead of the Euler-Maruyama, it is possible to use any other stochastic integrator for this part of computation; see [54, Remark 1 in Section III.B]. scheme [57, Ch. 10] applied to (64) with time-step $h=10^{-4}$, i.e., for $k\in\mathbb{N}$,

[TABLE]

where $\{\Delta\bm{w}_{k}^{i}\}_{i=1}^{N}$ are independent and identically distributed samples from a Gaussian joint PDF with zero mean and covariance $h\bm{I}$. To compute the joint PDF values $\{\varrho_{k}^{i}\}_{i=1}^{N}$ at the updated sample locations, we apply the gradient flow algorithm in [54, Section III.B] to solve the proximal recursion (60) with entropic-regularization parameter $\epsilon=0.1$. This results in a scattered weighted point cloud evolution as depicted in Fig. 5, and thanks to the exponential rate-of-convergence [33] to the stationary PDF, the local minima of $f$ can be quickly ascertained from the evolution of the modes of the transient joint PDFs. Fig. 6 shows that the resulting computation has extremely fast runtime due to certain nonlinear contraction mapping property established in [54].

An appealing feature of this proximal gradient descent computation is that it does not involve spatial discretization or function approximation, and is therefore scalable to high dimensions. In particular, we discretize time but not the state space, and evolve weighted point clouds over time. We eschew the algorithmic details and refer the readers to [54].

This example demonstrates that the HNN SDE (45) and its associated proximal recursion (60) for gradient descent in $\mathcal{P}_{2}(\mathcal{M})$ can be a practical computational tool for global optimization.

VI Conclusions

In this paper, we highlight the metric geometry of the dynamics associated with the deterministic and the stochastic Hopfield Neural Networks (HNN). We show that the deterministic HNN flow is a natural gradient descent with respect to a metric tensor that depends on the choice of activation functions. Alternatively, the same deterministic HNN flow can be interpreted as mirror descent on a suitable manifold that again is governed by the activation functions. The stochastic HNN, also known as the “diffusion machine”, defines a Wasserstein gradient flow on the space of probability measures induced by the underlying stochastic differential equation on the ground Riemannian manifold. This particular viewpoint leads to infinite dimensional proximal recursion on the space of probability density functions that provably approximates the same Wasserstein gradient flow in the small time step limit. We provide a numerical example to illustrate how this proximal recursion can be implemented over probability weighted scattered point cloud, obviating spatial discretization or function approximation, and hence is scalable. Our numerical experiments reveal that the resulting implementation is computationally fast, and can serve as a practical stochastic optimization algorithm. A conceptual summary of the gradient descent interpretations for the deterministic and the stochastic HNN is outlined in Table I.

We clarify here that our intent in this paper has been understanding the basic geometry underlying the differential equations for the HNN flow. We hope that the results of this paper will be insightful to the practitioners in the optimization and machine learning community, and motivate further theoretical and algorithmic development of global optimization via neural networks.

-A Derivation of (52)

We start by noting that

[TABLE]

where the last step follows from integration-by-parts w.r.t. $\bm{x}$, and zero flux of probability mass across the boundary of $\mathcal{M}$.

For $F(\rho)$ given by (51), $\delta F/\delta\rho=f+T\left(1+\log\rho\right)$, and therefore,

[TABLE]

Combining (67) and (68), we arrive at (52).