A Simultaneous Perturbation Weak Derivative Estimator for Stochastic Neural Networks

TL;DR

This paper introduces a novel gradient estimation method for stochastic neural networks, specifically the Little model, using measure-valued differentiation and simultaneous perturbation, enabling training without explicit stationary distribution solutions.

Contribution

It presents a new derivative estimation technique applicable to general stochastic neural networks, extending previous methods limited to specific network structures.

Findings

Effective gradient estimation for the Little model.

Applicable to networks with complex, nonlinear stochastic units.

Extends previous algorithms to more general network architectures.

Abstract

In this paper we study gradient estimation for a network of nonlinear stochastic units known as the Little model. Many machine learning systems can be described as networks of homogeneous units, and the Little model is of a particularly general form, which includes as special cases several popular machine learning architectures. However, since a closed form solution for the stationary distribution is not known, gradient methods which work for similar models such as the Boltzmann machine or sigmoid belief network cannot be used. To address this we introduce a method to calculate derivatives for this system based on measure-valued differentiation and simultaneous perturbation. This extends previous works in which gradient estimation algorithms were presented for networks with restrictive features like symmetry or acyclic connectivity.

Equations

$$F({\xi }_i(t)<x)=\sigma (x):=\frac{1}{1+\exp (-x)}.$$

$$f_i(x,(w,b),\xi )=\{\begin{array}{l}1\text{ if }\sum _{j=1}^n{w}_{i,j}{x}_j+b_i>{\xi }_i,\\ 0\text{ otherwise.}\end{array}$$

$$x(t+1)=f(x(t),\theta ,\xi (t+1))$$

$$u_i(x,(w,b))=\sum _{j=1}^n{w}_{i,j}{x}_j+b_i.$$

$$P_{\theta }(x^0,x^1)=\prod _{i=1}^n\sigma (u_i(x^0,\theta ){)}^{x_i^1}(1-\sigma (u_i(x^0,\theta )){)}^{1-x_i^1}.$$

$$x^{†}=2x-1.$$

$$P_{\theta }(x^0,x^1)=\prod _{i=1}^n\sigma ((x_i^1{)}^{†}{u}_i(x^0,\theta )).$$

$$\underset{\theta \in \Theta }{\min }\left[J(\theta ):={\int }_{X}e(x)\mathrm{d}{\pi }_{\theta }(x)\right]$$

$$d_{TV}({\mu }_1,{\mu }_2)=\underset{B\subseteq X}{\sup }|{\mu }_1(B)-{\mu }_2(B)|$$

$$d_{TV}({\mu }_{1}P,{\mu }_{2}P)\le (1-ϵ)d_{TV}({\mu }_1,{\mu }_2).$$

$$d_{TV}(\mu P^t,\pi )\le (1-ϵ{)}^t{d}_{TV}(\mu ,\pi ).$$

$${\nabla }_{\theta }{\mu }_{\theta }(e)=c_{\theta }[{\mu }_{\theta }^{+}(e)-{\mu }_{\theta }^{-}(e)].$$

$$\underset{M^0\to \infty ,M^1\to \infty }{\lim }\mathbb{E}\left[{\Delta }_{M^0,M^1}^{MVD}\right]={\nabla }_{\theta }{\pi }_{\theta }(e).$$

$$P(v)=\prod _{i=1}^n[\frac{1}{2}{\delta }_{-1}(v_i)+\frac{1}{2}{\delta }_1(v_i)]$$

$$\sigma \left((x_i^1{)}^{†}{u}_i(x^0,\theta )\right)$$

$$\ge \sigma \left(-\sum _{j=1}^n|w_{i,j}|-|b_i|\right)$$

$$\ge \sigma \left(-\Vert w{\Vert }_{\infty }-\Vert b{\Vert }_{\infty }\right).$$

$$d_{TV}({\mu }_{1}P,{\mu }_{2}P)\le (1-\sigma (-\Vert w{\Vert }_{\infty }-\Vert b{\Vert }_{\infty }{)}^n)d_{TV}({\mu }_1,{\mu }_2)$$

$${\nabla }_{\theta }{\mu }_{\theta }(e)v=c_{\theta ,v}[{\mu }_{\theta ,v}^{+}(e)-{\mu }_{\theta ,v}^{-}(e)].$$

$$\nabla_{\theta}\mu_{\theta}(e)\>v=\nabla_{\lambda}\big{|}_{\lambda=0}\,\mu_{\theta+\lambda v}(e).$$

$${\mu }_{\theta }=\sum _{i=1}^m\frac{e^{-{\theta }_i}}{Z_{\theta }}{\mu }_i$$

$${\mu }_{\theta +\lambda v}(e)=\sum _{i=1}^m\frac{e^{-({\theta }_i+\lambda v_i)}}{Z_{\theta +\lambda v_i}}{\mu }_i(e).$$

$$\nabla_{\theta}\mu_{\theta}(e)\>v=\nabla_{\lambda}\big{|}_{\lambda=0}\,\mu_{\theta+\lambda v}(e)=c_{\theta,v}\left[\sum\limits_{i=1}^{m}\alpha_{i,v}\mu_{i}(e)-\sum\limits_{i=1}^{m}\beta_{i,v}\mu_{i}(e)\right]$$

$$c_{\theta ,v}=\frac{K_{\theta ,v}}{Z_{\theta }},$$

$${\alpha }_{i,v}=\frac{1}{Z_{\theta }{K}_{\theta ,v}}\left[(v_i{)}_{-}{Z}_{\theta }+\sum _{j=1}^m{e}^{-{\theta }_i}(v_j{)}_{+}\right],$$

$${\beta }_{i,v}=\frac{1}{Z_{\theta }{K}_{\theta ,v}}\left[(v_i{)}_{+}{Z}_{\theta }+\sum _{j=1}^m{e}^{-{\theta }_i}(v_j{)}_{-}\right].$$

$$p(v)=\prod _{i=1}^m[\frac{1}{2}{\delta }_{-1}(v_i)+\frac{1}{2}{\delta }_1(v_i)]$$

$${\nabla }_{\theta }\left[\sum _{x^1}e(x^1)P_{\theta }(x^0,x^1)\right]v=c_{\theta ,v}(x^0)\left[\sum _{x^1}e(x^1)P_{\theta ,v}^{+}(x^0,x^1)-\sum _{x^1}e(x^1)P_{\theta ,v}^{-}(x^0,x^1)\right],$$

$$c_{\theta ,v}(x)=\sum _{i=1}^n|v_i|\sigma (u_i(x^0,\theta ))+\sum _{i=1}^n\sum _{j=1}^n|v_{i,j}|x_j^0\sigma (u_i(x^0,\theta ))),$$

$$\begin{array}{rl}& P_{\theta ,v}^{+}(x^0,x^1)=\\ & \frac{P_{\theta }(x^0,x^1)}{c_{\theta }(x^0)}\sum _{i=1}^n\left[x_i^1\left((v_i{)}_{+}+\sum _{j=1}^n(v_{i,j}{)}_{+}{x}_j^0\right)+\sigma (u_i(x^0,\theta ))\left((v_i{)}_{-}+\sum _{j=1}^n(v_{i,j}{)}_{-}{x}_j^0\right)\right].\end{array}$$

$$\begin{array}{rl}& P_{\theta ,v}^{-}(x^0,x^1)=\\ & \frac{P_{\theta }(x^0,x^1)}{c_{\theta }(x^0)}\sum _{i=1}^n\left[x_i^1\left((v_i{)}_{-}+\sum _{j=1}^n(v_{i,j}{)}_{-}{x}_j^0\right)+\sigma (u_i(x^0,\theta ))\left((v_i{)}_{+}+\sum _{j=1}^n(v_{i,j}{)}_{+}{x}_j^0\right)\right].\end{array}$$

Full text

1 Introduction

Many computational models in machine learning can be described as networks of homogeneous units. Although each unit may be very simple, capable of only trivial logical or mathematical operations, the hope is that through a training procedure the interconnection of these units can be arranged in order that the overall network can perform useful tasks, such as classification or regression. In machine learning, “training” means the optimization of the parameters of the model in order to minimize an average loss over available data. If gradient-based methods are to be used for the optimization, then a fundamental step will be computing the derivative of an appropriate cost function with respect to some network parameters. The exact means to compute the derivatives will depend on the details of the network - including whether the units are deterministic or stochastic, whether they satisfy any smoothness parameters, and network properties like symmetry, or the presence of cycles. While certain combinations of the aforementioned features lend themselves to easy treatment from the perspective of gradient estimation, in other cases it is a difficult problem. In this work we show how measure-valued differentiation can be fruitfully applied to this problem, to yield gradient estimators for very general network structures. This extends previous works in which gradient estimation algorithm’s were presented for networks with restrictive features like symmetry or acyclic connectivity.

The networks we are interested in are probabilistic and operate on a finite state space. Each unit can be in one of two states, [math] or $1$, and if there are $n$ nodes in total, then the state space is $X=\{0,1\}^{n}$. In general for stochastic neural networks, at each time step one or more units may change their state. The particular update rule in the network studied here was first defined in [14], and hence we sometimes refer to it as the Little model. Let $\xi(1),\xi(2)\ldots$ be a sequence of noise vectors in $\mathbb{R}^{n}$, with the entire collection $\{\xi_{i}(t);i=1,\ldots,n,t=1,2,\ldots\}$ independent and distributed according to the logistic distribution. That is, the cumulative distribution function of $\xi_{i}(t)$ is

[TABLE]

The parameters of the model are denoted by $\theta\in\Theta\subset\mathbb{R}^{n\times n}\times\mathbb{R}^{n}$, where $\theta=(w,b)$ comprises the weights $w_{i,j}$ (in network terminology, these are the values of the corresponding links $(i,j)$) and the bias $b_{i}$ associated with node $i$.

Define $f:X\times\Theta\times\Xi\to X$ as

[TABLE]

This function $f$ and the noise $\xi(1),\xi(2)\ldots$ determines the operation of the network; from the initial point $x(0)$ the states follows the recursion

[TABLE]

to generate the next state. They can be interpreted as threshold networks, where the thresholds are random at each time step.

We let $P_{\theta}$ be the Markov kernel corresponding to this recursion, and $P_{\theta}(x^{0},x^{1})$ be the probability of going to state $x^{1}\in X$ from state $x^{0}\in X$. The function $u_{i}(x)$, that determines the input to each node at the state $x$ is defined as

[TABLE]

The above equation describes the input “flow” to node $i$ in terms of the network model, adding the bias to the total flow from incoming nodes. Then

[TABLE]

Alternatively, we can use the following notation of [17]: for $x\in\{0,1\}$, define

[TABLE]

Using this, together with the identity $1-\sigma(x)=\sigma(-x)$, an equivalent expression to (4) when $x_{i}\in\{0,1\}$ is given by:

[TABLE]

In practice, a user can compute with this stochastic network by fixing an input and iterating the update rule for a large number of steps before observing the network state. To ensure that the long-run statistical behavior of the network is independent of the initial conditions, we should establish ergodicity of the Markov kernel. In Section 4 we show that $P_{\theta}$ is ergodic and find the convergence rate of the Markov chain in terms of $\|\theta\|$ and $n$. Denote by $\pi_{\theta}$ be the stationary measure of this Markov chain.

Given a cost function $e:X\to\mathbb{R}$, the optimization problem is then

[TABLE]

and $\pi_{\theta}$ is defined as the solution to $\pi_{\theta}P_{\theta}=\pi_{\theta}$, for the Markov kernel $P_{\theta}$ defined in (4). In this work the focus is on how to compute the gradient $\nabla_{\theta}J$.

2 Related work

Several stochastic neural networks on discrete state spaces have been studied, and their gradient estimation procedures are based on having closed form solutions for the resulting probability distributions. The works [11, 18, 3] assume specific constraints on network connectivity - for instance symmetry, or prohibiting cycles. The earliest neural network models to be studied from the computational view were the deterministic threshold networks [15, 22]. In this model, each unit senses the states of its neighbors, takes a weighted sum of the values, and applies a threshold to determine its next state (either on or off). For single layer versions of these networks, where the units are partitioned into input and output groups, with connections only from input to output nodes, the corresponding optimization problem can be solved by the perceptron algorithm [22]. Any iterative algorithm for optimizing threshold networks has to address the credit assignment problem [16]. This means that during optimization, the algorithm must identify which internal components of the network are not working correctly, and adjust those units to improve the output. The difficulty in solving the credit assignment problem for threshold networks with multiple layers prevents simple deterministic threshold models from being used in complex problems like image recognition. Although there is yet no universal method to train the deterministic threshold networks, recent research focuses on specific types of networks that can be trained using gradient-descent optimization. For instance, one can abandon the threshold units, and work with units that have a smooth, graded, response such as the sigmoid neural networks [24]. In this case methods of calculus are available to determine unit sensitivities. These new networks are still deterministic but now operate on a continuous state space.

Another approach is to keep the space discrete but make the network probabilistic, and use the smoothing effects of the noise to obtain a model one can apply methods of calculus to. One can interpret the sigmoid belief networks in this way. These networks were introduced in [17] and so named because they combine features of sigmoid neural networks and Bayesian networks. In these networks, when a unit receives a large positive input it is very likely to turn on, while a large negative input means the unit is likely to remain off. In fact, these networks can be interpreted as threshold networks with random thresholds. The use of the sigmoid function, which is the cumulative distribution function (CDF) of the logistic distribution, leads to an interpretation of a network with thresholds drawn from the logistic distribution. In [17], the author derived formulas for the gradient in these networks, and showed how Markov chain Monte Carlo (MCMC) techniques can be used to implement gradient estimators. The networks studied in [17] had a feed-forward architecture, but one could also define variants that allow cycles among the connections. If the connectivity graph in the Little model is acyclic, then one obtains a model resembling the sigmoid belief network. We can enforce this by requiring $w_{i,j}=0$ if $i<j$. In this way one is lead to the random threshold networks. In this case, one would be interested in the long-term average behavior of the network. Such a generalization would resemble the random threshold networks that are our focus. It would be interesting to obtain a gradient estimator for these new networks.

Another motivation to study general random threshold networks comes from the Boltzmann machine [11, 1, 25]. This is a network of stochastic units that are connected symmetrically. This means there is feed-back in the network, and the problem in these networks is to optimize the long-term behavior. The Boltzmann machine follows an update rule similar to (1), except that the weights are constrained to be symmetric ($w_{i,j}=w_{j,i}$) and the nodes are updated one at a time, instead of all at once. This model was an important ingredient in many machine learning systems [10, 27]. The symmetry in the network, and the use of the sigmoid function to calculate the probabilities, leads to a nice closed form solution for the stationary measure in this model. Based on formulas for the stationary distribution, expressions for the gradient of long-term costs can be obtained, leading to MCMC based gradient estimators. However, if one changes the model, by for instance using non-symmetric connections, or changing the type of nonlinearity, closed formulas are no longer available. Instead, one winds up with a model of random threshold networks. Specifically, a model like the Boltzmann machine is obtained if the weights are symmetric, meaning $w_{i,j}=w_{j,i}$. Technically, if one puts a symmetry requirement on our threshold networks, one does not exactly recover the Boltzmann machine, but a variant known as the synchronous or parallel Boltzmann machine [18]. The synchronous Boltzmann machine also has a known, simple, stationary distribution [18]. This provides another motivation for studying gradient estimation in the Little model.

In the case of networks where cycles are allowed, the work [2] considered gradient estimation in the finite horizon setting. Our interest is in the long-term average cost in networks that have general connectivity, where only knowledge of the transition probabilities is available. Methods such as forward sensitivity analysis cannot be used in this case, as they rely on the differential structure of the underlying state space. Instead, we propose an algorithm that computes descent directions based on simultaneous perturbation analysis and measure-valued differentiation (MVD).

In machine learning, most of the optimization problems are solved either using a closed form of the stationary averages and using a deterministic gradient-based method, or using a stochastic gradient when the cost function is smooth and one can interchange derivative and expectation. However, this limits the class of models that can be used in practice. On the other hand, networks with less constraints may have a more powerful modeling capability but estimating their gradients is not as straight forward and so far, little has been done to apply more sophisticated gradient estimation techniques in part because the typical practical applications require very high dimensional gradients and the code must be very efficient to be of practical use.

The main contribution of this paper to the area of machine learning is the application of MVD to a feedback neural network in order to train the machine using gradient descent stochastic approximation. The main contribution of the paper to the area of gradient estimation is the description of a novel methodology to efficiently calculate high dimensional gradients. Our proposed method combines the main idea of directional derivatives that has enjoyed unprecedented success with SPSA using finite differences, and MVD gradient estimation, which yields unbiased estimation and thus, bounded variance.

3 Preliminaries

Since we are interested in performance of Markov chains, we should introduce our ergodicity framework.

3.1 Ergodicity Framework

We will work with the total variation metric: The total variation distance between two probability measures $\mu_{1},\mu_{2}$ on a discrete space $X$ is

[TABLE]

where the supremum is over all subsets of $X$. The abstract result we use is the following:

Proposition 3.1**.**

Let a Markov kernel $P$ have the property that $\min_{x,y}P(x,y)\geq\epsilon>0$. Then $P$ is a contraction for $d_{TV}$: For any two measures $\mu_{1},\mu_{2}$,

[TABLE]

In particular, if $\pi$ is the stationary measure of $P$ then for any initial measure $\mu$,

[TABLE]

For instance, see Lemma 2.20 in [21]. That lemma says that the convergence rate is at least $1-\min_{i,j}\sum_{k}\min\{p_{i,k},p_{j,k}\}$, and it’s clear from our assumption that this is greater than $1-\epsilon$.

Next let us recall some of the basic ideas from measure-valued differentiation and simultaneous perturbation.

3.2 Measure-Valued Differentiation

The idea of measure-valued differentiation is to express the derivative of an expectation as the difference of two expectations. Each of these expectations involves the same cost function, but the underlying measures are different. This enables simple, unbiased derivative estimators if these measures are easy to sample from. For simplicity, in the following definitions we will only consider the (relevant to our purposes) case of a finite state space $X$. The method was pioneered by the “weak derivative” estimation method of [19, 20], and extended to cover discontinuous, and unbounded cost functions in the introductory papers [8, 7].

Definition 3.2**.**

Consider a measure $\mu_{\theta}$ on a finite state space $X$ that depends on a real parameter $\theta$. The measure $\mu_{\theta}$ is said to be differentiable at $\theta\in\mathbb{R}$ if there is a triple $(c_{\theta},\mu^{+}_{\theta},\mu^{-}_{\theta})$, consisting of a real number $c_{\theta}$, and two probability measures $\mu^{+}_{\theta},\mu^{-}_{\theta}$ on $X$ such that, for any function $e:X\to\mathbb{R}$,

[TABLE]

An MVD gradient estimator consists of two parts: First, sample a random variable $Y^{+}$ distributed according to $\mu^{+}_{\theta}$, then sample a random variable $Y^{-}$ according to $\mu^{-}_{\theta}$, and finally form the estimate $\Delta^{MVD}=c_{\theta}[e(Y^{+})-e(Y^{-})]$. For more background see [9]. The concept of MVD can be extended from measures to Markov kernels, and then applied to derivatives of stationary costs [19, 20, 7].

Definition 3.3**.**

A Markov kernel $P_{\theta}$ that is defined on a finite space $X$ and which depends on a real parameter $\theta$ is said to be differentiable at $\theta$ if for each $x\in X$ there is a triple $(c_{\theta}(x),P^{+}_{\theta}(x,\cdot),P^{-}_{\theta}(x,\cdot))$ which is the measure-valued derivative of the measure $P_{\theta}(x,\cdot)$ at the parameter $\theta$ in the sense of Definition 3.2.

If the Markov kernel $P_{\theta}$ is ergodic with stationary measure $\pi_{\theta}$, then in certain cases it is possible to use the $(c_{\theta},P^{+}_{\theta},P^{-}_{\theta})$ to compute the stationary derivatives $\nabla_{\theta}\pi_{\theta}(e)$ [20]. The procedure is shown in Algorithm 1. The correctness is expressed in Theorem 3.4, a result proved in [20]. It gives a condition on a Markov chain $P_{\theta}$ that guarantees the corresponding stationary costs $\pi_{\theta}(e)$ are differentiable, and establishes that the procedure presented in Algorithm 1 (see below) can be used to compute the derivative of the stationary cost. Note that we are recalling a simplified version of the Theorem in the case of a finite state space. Similar results on differentiability for discrete chains have been developed by other authors, sometimes leading to different computational methods; see related results in [4, 12].

Observe that Algorithm 1 has two settings $M^{0}$ and $M^{1}$. To emphasize the dependence of the random variable $\Delta^{MVD}$ on these parameters, we use the notation $\Delta^{MVD}_{M^{0},M^{1}}$.

Theorem 3.4** ([20]).**

Let $\theta\in\mathbb{R}$ be a scalar parameter of a family of Markov kernels $P_{\theta}$, and let $\delta_{x}P_{\theta}$ denote the (conditional) probability measure $P_{\theta}(x,\cdot)$, given state $x\in X$. Assume that $(\delta_{x}P_{\theta})(e)$ is differentiable in $\theta$ for each cost function $e$. Suppose further that $P_{\theta}$ is a contraction on the space of probability measures in the sense of Inequality 8. Then the stationary cost $\pi_{\theta}(e)$ is differentiable, and Algorithm 1 can be used to estimate the derivatives. Specifically, if we let $\Delta_{M^{0},M^{1}}^{MVD}$ be the output of the algorithm, then

[TABLE]

More general results on MVD for stationary measures can be found in [6, 7], including the extension to discontinuous, and unbounded cost functions.

Observe that we call for the use of common random numbers in the above algorithm. This ensures that the variance of the estimate $\Delta^{MVD}$ can be bounded independently of the parameters $M^{1}$ and $M^{0}$ [21, Theorem 4.36]. This can be achieved through any coupling with a contraction property - in the case of a finite state space, using common random numbers is one way - but this is not the only way; see [21] for details.

For measure-valued differentiation, like finite differences, it seems that in order to estimate the full derivative for a system with $m$ parameters, $m$ parallel simulations are required although the variance characteristics are much more favorable compared to finite differences (see [21], Section 4.3). In finite differences, one must trade off bias for variance, but for MVD the variance can be shown to be bounded independently of the parameters $M^{1},M^{0}$ which determine the bias.

3.3 Simultaneous Perturbation

One interesting solution to the prohibitive computational complexity of estimating derivatives with finite differences in high-dimensions is known as simultaneous perturbation [26]. In this scheme, one picks a random direction $v$, and then approximates the directional derivative using stochastic finite differences. In this case only two simulations are needed for a system with $m$ parameters. Using random directions in stochastic approximation was also studied in [13, Section 2.3.5] and [5].

The variance issues of finite differences remain with this approach; in order to decrease the bias of the estimator, one has to deal with a larger variance. For generating the directions, one possibility is to let $v$ be a random point on the hypercube $\{-1,1\}^{m}$, as suggested in [26]. For the theoretical analysis, it is important that the directions have zero mean, and that the random variable $\frac{1}{\|v\|}$ be integrable. The procedure is shown in Algorithm 2.

In the above algorithm, the random direction is generated according to Spall’s original SPSA method, namely each component of the random vector is either $-1$ or $1$ with equal probability, and components are independent [26, Section V-A]. Other ways to generate random directions have been studied further, and our method can easily be generalized in this way as well.

This algorithm will form the basis for the gradient estimator we derive for discrete attractor networks below.

4 Ergodicity of Little model

In this section we consider the ergodicity of the Little model [14]. We will show that $P_{\theta}$ has a contraction property with respect to $d_{TV}$.

The contraction coefficient $\alpha$ will depend on the weights of the network; larger weights will lead to a worse bound on the convergence. Let us calculate the lower bound $\epsilon$ required by Proposition 3.1 for the transition probability $P_{\theta}(x_{0},x_{1})$. According to the product form of $P_{\theta}$, (6), it suffices to find an $\epsilon_{1}$ such that $\sigma((x_{i}^{1})^{{\dagger}}u_{i}(x^{0},\theta))\geq\epsilon_{1}$ and then we can set $\epsilon=\epsilon_{1}^{n}$. Note that for any value of $x_{i}^{1}$, we have

[TABLE]

where $\|w\|_{\infty}$ is the $\infty$-norm of the matrix $w$.111This norm for matrices is defined as as $\|w\|_{\infty}=\sup_{\|u\|_{\infty}=1}\|wu\|_{\infty}$, where for the vector $u$ and $wu$, the norm $\|\cdot\|_{\infty}$ is defined in the usual way. Letting $\epsilon=\sigma\left(-\|w\|_{\infty}-\|b\|_{\infty}\right)^{n}$, results in the formula

[TABLE]

It follows that the predicted speed of convergence decreases as the weights increase.

Given the contraction property, we now consider the differentiability of the stationary distribution. Firstly, we are dealing with a finite state space and a Markov kernel with smooth transition probabilities (due to the smoothness of $\sigma$). Hence we are in a simple setting and the differentiability is relatively easy to establish. For instance Lemma 4 in [20] can be applied. The conclusion is that for any cost function $e$, the stationary expectation $\pi_{\theta}(e)$ is a differentiable function of $\theta$.

Remark: When using gradient-descent or other updating methods to find optimal $\theta$, one must ensure that either $\theta$ stays always inside a compact set, or that the sequence of consecutive updates is tight. This will follow from the convexity of the cost function, from Lyapunov stability, or from truncations performed in practice when applying such algorithms. This in turn establishes uniform ergodicity and differentiability of the stationary measure at each step of the updates.

5 Gradient Estimation

Gradient estimation has been studied for closely related models, such as the Boltzmann machine and sigmoid belief networks, as we discuss below. The case of the Little model is somewhat more challenging since there is not a known closed-form solution for the stationary distribution. One has to focus on gradient estimation methods that only use the Markov kernel associated to the process.

We propose an estimator which combines features of SP and MVD. The algorithm generates a random direction, as in SP, and then uses measure-valued differentiation to approximate this directional derivative. In this way one deals with a small number of simulations, as in SP, while avoiding the variance issues with finite differences. We call the method simultaneous perturbation measure-valued differentiation. The only requirement is that one can compute the measure-valued derivative along arbitrary directions. After developing the method in an abstract setting, we will consider the method in the context of the Little model.

5.1 SPMVD

We introduce the directional derivatives, as usual. Following standard notation, the gradient of the expected value of the cost function $e$ under measure $\mu_{\theta}$ is denoted by $\nabla_{\theta}\mu_{\theta}(e)$ and it is a row vector with as many components as the dimension of $\theta$. Therefore, the projection of this vector along a direction $v$ (a column vector with same number of components) is the (scalar) inner product $\nabla_{\theta}\mu_{\theta}(e)\>v$.

Definition 5.1**.**

Let $\mu_{\theta}$ be a measure depending on an $m$-dimensional vector parameter $\theta$. Let $v\in\mathbb{R}^{m}$ be a direction. A triple $(c_{\theta,v},\mu^{+}_{\theta,v},\mu^{-}_{\theta,v})$ is called a measure-valued directional directional derivative at $\theta$ in the direction $v$ if for all $e:X\to\mathbb{R}$,

[TABLE]

Note that in this expression, the left hand side is the dot product between the $m$ dimensional vectors $\nabla_{\theta}\mu_{\theta}(e)\>$ and $v$. The three terms making up the expression on the right side are all real-numbers.

In practice, one can try to calculate the MVD in direction $v$ as follows. By basic calculus,

[TABLE]

Therefore, to find the MVD of $\mu_{\theta}$ in direction $v$ it suffices to find the usual scalar, MVD for $\mu_{\theta+\lambda v}$ at $\lambda=0$. This is the approach in the following example.

Example 5.2**.**

Let $\mu_{1},\mu_{2},\ldots,\mu_{m}$ be $m$ probability measures, and for a vector parameter $\theta\in\mathbb{R}^{m}$ define $\mu_{\theta}$ as

[TABLE]

where $Z_{\theta}=\sum\limits_{i=1}^{m}e^{-\theta_{i}}$. For any function $e:X\to\mathbb{R}$ and direction $v\in\mathbb{R}^{m}$, then, we have

[TABLE]

Introduce the notation $\gamma_{+}=\max\{0,\gamma\}$ and $\gamma_{-}=\max\{0,-\gamma\}$, for the positive and negative part of a scalar $\gamma$, respectively. Then we have the following identities: $\gamma=\gamma_{+}-\gamma_{-}$ and $|\gamma|=\gamma_{+}+\gamma_{-}$. Define $K_{\theta,v}=\sum\limits_{j=1}^{m}e^{-\theta_{j}}|v_{j}|$. Differentiating (10) at $\lambda=0$, and doing some algebra, one can get the following representation for the directional derivative:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Therefore the triple $\left(c_{\theta,v},\sum\limits_{i=1}^{m}\alpha_{i}\mu_{i},\sum\limits_{i=1}^{m}\beta_{i}\mu_{i}\right)$ is the measure-valued derivative of $\mu$ at $\theta$ in the direction $v$.

To generate samples from the mixture $\sum\limits_{i=1}^{m}\alpha_{i}\mu_{i}$ one needs to choose one of the components $i=1,..,m$ according to the mixture probabilities $\alpha_{i},i=1,\ldots,m$, and then generate a sample from $\mu_{i}$. Similarly for the negative measure. Thus, one needs to generate two random samples from the underlying collection of measures to get an unbiased derivative estimate. This concludes the example.

Defining the extension of directional derivatives to Markov chains is straightforward.

Definition 5.3**.**

A triple $(c_{\theta,v},P^{+}_{\theta,v},P^{-}_{\theta,v})$ is a measure-valued derivative for the Markov kernel $P$ at $\theta$ in the direction $v$ if for each $x$, $(c_{\theta,v}(x),P^{+}_{\theta,v}(x,\cdot),P^{-}_{\theta,v}(x,\cdot))$ is an MVD at $\theta$ in the direction $v$ for the measure $P_{\theta}(x,\cdot)$ in the sense of Definition 5.1.

The gradient estimator for stationary costs proceeds by choosing a random direction, and applying the stationary MVD procedure of Algorithm 1. The pseudocode is presented as Algorithm 3.

Note that $\Delta^{MVD}$ is a real-number, and the algorithm output, $\Delta^{MVD}v$, is a vector, of the same dimensionality as $\theta$. This estimator will be applied to our motivating example, the random threshold networks.

5.2 Application to the Little model

Let us now give the measure-valued directional derivatives for the Little model. Fix a parameter $\theta.$ A perturbation in the parameter space is represented as a vector $v\in\mathbb{R}^{n\times n}\times\mathbb{R}^{n}$, where $v_{i,j}$ is a perturbation along the weight from $j$ to $i$ and $v_{i}$ is a perturbation along the bias at node $i$.

Using definitions (3), (5), (6), and after some algebra, one can obtain the following expression for the directional MVD:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

See Section A.1 for a detailed derivation of Equation 12. This yields two Markov kernels $P^{+}_{\theta,v}$ and $P^{-}_{\theta,v}$, that depend not just on the parameter $\theta$ (as would be the case in scalar MVD) but also the direction of interest $v$.

In order for SPMVD to be useful, there must be a practical procedure for running the Markov kernels $P^{+}_{\theta,v}$ and $P^{-}_{\theta,v}$. The Little networks operate on a large state space, of size $2^{n}$ when $n$ nodes are used, so this is not necessarily trivial. The original Markov kernel $P_{\theta}$ has a relatively simple structure, allowing each node to be updated independently (see Eqn. 6), and one may hope for a similar situation with the MVD pair $P_{\theta,v}^{+},P_{\theta,v}^{-}$.

To investigate this, let us consider the Markov kernel $P_{\theta,v}^{+}$. Fix an $x_{0}$ and a $\theta$. We want to see how to generate an $x^{1}$ sample from $P_{\theta,v}^{+}(x^{0},\cdot)$.

Define the following variables:

1. Pi.

$d=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sigma(u_{i}(x^{0},\theta))(v_{i,j})_{-}\,x_{j}^{0}+\sum\limits_{i=1}^{n}\sigma(u_{i}(x^{0},\theta))(v_{i})_{-}$ 2. Pii.

$a_{i}=(v_{i})_{+}+\sum\limits_{j=1}^{n}(v_{i,j})_{+}\,x_{j}^{0}$,   $i=1,\ldots,n$, 3. Piii.

$\beta_{i}=\sigma(u_{i}(x^{0},\theta))$,   $i=1,\ldots,n$, 4. Piv.

$c=c_{\theta,v}(x^{0})$.

Each of these depends on $x_{0},\theta$ and $v$. Then the probability $Q(x)=P_{\theta,v}^{+}(x^{0},x)$ has the representation

[TABLE]

We will sample from this distribution by sequentially generating the random variables $x_{1},\ldots,x_{n}$. First we will calculate and sample from the marginal distribution $Q(x_{1})$. Then we sample from the conditional distribution $Q(x_{2}\mid x_{1})$, followed by sampling from $Q(x_{3}\mid x_{1},x_{2})$ and so on until finally sampling from $Q(x_{n}\mid x_{1},\ldots,x_{n-1})$. This is a standard technique for generating random vectors (see Section 4.6 of [23]). As we shall see, it is feasible since the conditional probabilities are easy to compute.

To start, note that by the definition of conditional probability,

[TABLE]

Based on this equation, if we can compute the marginal probabilities quickly then we can compute the conditional probabilities quickly. Starting from the formula (15), one can show that for $x_{1}\in\{0,1\}$,

[TABLE]

and for any $(x_{k},\ldots,x_{1})\in\{0,1\}^{k}$,

[TABLE]

See Appendix A.2 for a derivation of (17) and (18). Based on these formulas, we construct the sampling algorithm (Algorithm 4). Using the identity (18), we can write (16) as

[TABLE]

The following proposition certifies the correctness of Algorithm 4.

Proposition 5.4**.**

For any data $c,d,a_{1},\ldots,a_{n},\beta_{1},\ldots,\beta_{n}$, the output of Algorithm 4 is distributed as follows:

[TABLE]

Proof.

This follows from the identity (16) and the computational formulas (17), (18), whose correctness is established in the appendix. ∎

Analogous to the variables (Pi) - (Piv), we have definitions for $P^{-}$:

1. Mi.

$d=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sigma(u_{i}(x^{0},\theta))(v_{i,j})_{+}\>x_{j}^{0}+\sum\limits_{i=1}^{n}\sigma(u_{i}(x^{0},\theta))(v_{i})_{+}$ 2. Mii.

$a_{i}=(v_{i})_{-}+\sum\limits_{j=1}^{n}(v_{i,j})_{-}\>x_{j}^{0}$,   $i=1,\ldots,n$, 3. Miii.

$\beta_{i}=\sigma(u_{i}(x^{0},\theta))$,   $i=1,\ldots,n$, 4. Miv.

$c=c_{\theta,v}(x^{0})$.

Corollary 5.5**.**

Let $c,d,a_{1},\ldots,a_{n},\beta_{1},\ldots,\beta_{n}$ be defined as in (Pi) - (Piv). Then the distribution of the output of Algorithm 4 is $P_{\theta,v}^{+}(x_{0},x_{1})$.

Alternatively, if the $c,d,a_{1},\ldots,a_{n},\beta_{1},\ldots,\beta_{n}$ are defined as in (Mi) - (Miv) then the distribution of the output of Algorithm 4 is $P_{\theta,v}^{-}(x_{0},x_{1})$.

One can use Algorithm 4 to simulate the Markov chains $P_{\theta,v}^{+}$ and $P_{\theta,v}^{-}$ that are needed in the directional MVD algorithm. Two copies would be used - one with the data for $P_{\theta,v}^{+}$ and one for $P_{\theta,v}^{-}$.

The gradient estimation procedure is summarized below. It is the MVD gradient estimation algorithm (Algorithm 1), customized to include random directions and the special sampling algorithm (Algorithm 4).

Observe that the algorithm uses a coupling by common random numbers - see the discussion immediately after Algorithm 1 for the discussion of why this is case, and alternatives.

The theoretical soundness of this procedure rests on Theorem 3.4. That theorem tells us that conditioned on the random direction $v$, the output $\Delta^{MVD}$ converges in expectation to the directional derivative of the stationary cost in direction $v$, as $M^{0}$ and $M^{1}$ tend to infinity. This is established more formally in the following Proposition.

Proposition 5.6**.**

Let $\Delta^{MVD}_{M^{0},M^{1}}$ and $v$ be as defined in Algorithm 5. Then

[TABLE]

Proof.

By Theorem 3.4, we know that for any sampled direction $v$,

[TABLE]

In what follows, we will omit the notation indicating the dependence on $M^{0},M^{1}$, and write $\Delta^{MVD}$ instead. By conditioning on $v$, then,

[TABLE]

The interchange of the limit and expectation in the above equalities is allowed since $v$ is a discrete random variable that can only take values in the finite set $\{-1,1\}^{n}$. Continuing, then, we can combine Equation 20 with Equation 21 to see that

[TABLE]

In terms of the Euclidean basis vectors $e_{i}$, we have $v=\sum\limits_{i=1}^{n}v_{i}e_{i}$, and

[TABLE]

In the third equality we used the independence of the $v_{i}$, and that $\mathbb{E}[v_{i}]=0$ for all $i$. In the fourth inequality, we used that $\mathbb{E}[v_{i}^{2}]=1$. Combining 22 with 23 yields the result. ∎

This result can be compared with Lemma 1 of [26], which concerns the bias in the traditional SPSA gradient estimates from Algorithm 2.

A complexity analysis reveals that running the algorithm requires $\mathcal{O}\left((M^{0}+M^{1}+1)N^{2}\right)$ operations. Note that this assumes full connectivity of the underlying network, which results in $N^{2}$ parameters. This complexity follows from the fact that it takes $\mathcal{O}\left(N^{2}\right)$ operations to iterate the underlying chain for one step, and also that our intermediate step of preparing and running Algorithm 4 takes $\mathcal{O}(N^{2})$ operations.

In the next section we empirically investigate the behavior of Algorithm 5 when it used inside an optimization procedure.

6 Numerical Experiments

We have implemented Algorithm 5 to estimate gradients as part of a neural network training procedure. The network was trained to classify digits from the MNIST dataset222http://yann.lecun.com/exdb/mnist/. The network was trained using the first 5000 images digits of the dataset. We report how the empirical error evolved during training.

Formally, our objective function is the sum of $k=5000$ terms $J^{1},J^{2},\ldots,J^{k}$, one for each image. For simplicity, let us consider the term for a single image. In this case, each training datum consists of a $28\times 28$ pixel image $v\in\{0,1\}^{784}$ and a label $\ell\in\{0,1\}^{10}$. The vector $\ell$ is an indicator vector, with a single entry equal to one in the position corresponding to the class of the image and zeros elsewhere. The first $784$ nodes of the network receive the image pixels as an external input. The last ten nodes are output units. Hence the network has $784+10=794$ units (e.g. $n=794$). Among the units we allow arbitrary connectivity, resulting in $794\times 794$ connections. The error function $e$ is

[TABLE]

and the stationary cost is $J_{\theta}=\mathbb{E}_{x\sim\pi_{\theta}}[e(x)]$. Note that $\pi_{\theta}$ depends on the parameters $\theta$ and also the external input $v$.

In words, the goal of training is to find the parameters $\theta$ which guarantee that when the image $v$ is fixed as the input to the network, the network will compute the class of the image and indicate this class at the output nodes.

Each weight $w_{i,j}$ is initialized with a sample from a uniform distribution on $[-0.01,0.01]$. The biases $b_{i}$ are initialized with a sample from the same distribution. The parameters $\theta_{n+1}$ follow the recursion

[TABLE]

where at time $n$, the perturbation $\Delta_{n}$ is computed via Algorithm 5 as described above.

The parameters of Algorithm 5 were set up as follows. The value of $M^{0}$ was 10 and we experimented with two possible values for $M^{1}$, either $M^{1}=10$ or $M^{1}=50$. Common random numbers were used in the routines for updating $x^{+}$ and $x^{-}$ for variance reduction. The algorithm ran for 30000 parameter updates and the empirical error is reported every 500 updates. Results are shown in Figure 1. The trend is similar but using $M^{1}=50$ iterations inside Algorithm 5 seems to have less variance.

7 Discussion

In this paper we studied gradient estimation for the Little model. Since a closed form solution for the stationary distribution is not known, methods which work for similar models such as the Boltzmann machine or sigmoid belief network cannot be used. To address this we introduced a method to calculate derivatives in the Little model based on measure-valued differentiation. To get a gradient estimate this way, one has to run the Little model for some time to get an initial state, then generate a random direction, and then run two Markov kernels $P_{\theta,v}^{+}$ and $P_{\theta,v}^{-}$ and use the observed errors in both chains to form the gradient estimate. There are several parameters of the method - the $M^{0}$ and $M^{1}$ of Algorithm 5. Future work will study the dependence of the algorithm performance, such as bias and variance, on these parameters. A more detailed numerical study will also aid in tuning these parameters. The authors believe the general idea of pairing random directions with measure-valued-differentiation could enable optimization in other models as well.

Appendix A Derivations Related to the Little Model

A.1 Derivation of Equation 12

Fix an $x^{0}$ and a direction $v\in\mathbb{R}^{n\times n}\times\mathbb{R}^{n}$. Then

[TABLE]

Evaluating this at $\delta=0$ we find that

[TABLE]

Note also that

[TABLE]

which means

[TABLE]

Combining (24) and (25),

[TABLE]

Splitting each $v_{i,j}$ and $v_{i}$ into positive and negative parts,

[TABLE]

Note that

[TABLE]

Combining (27) with (28) and the definitions (13) and (14) we obtain (12).

A.2 Derivation of Equations 17 and 18

We have

[TABLE]

To simplify this equation, note that for $n>1$,

[TABLE]

and if $n=1$ then

[TABLE]

Combining Equation 30 and Equation 31, we see that for any $n\geq 1$,

[TABLE]

Combining Equation 29 with Equation 32,

[TABLE]

In general,

[TABLE]