A deep learning theory
for neural networks grounded in physics

Artificial intelligence (AI) emerged as a research discipline in the 1950s from the idea that every aspect of human intelligence can in principle be discovered, understood, and built into machines. The idea of AI started after Alan Turing formalized the notion of computation and began to study how computers can solve problems on their own (Turing, 1950), but the term artificial intelligence was first coined by John McCarthy in 1956. Today, AI is usually used in a broader sense, to refer more broadly to the field of study concerned with designing computational systems to solve practical tasks that we want to automate, whether or not such tasks require some form of human-like intelligence, and whether or not such computational systems are inspired by the brain. However, human intelligence is a natural benchmark for us to build ‘intelligent’ machines. This is an arbitrary choice, and implies by no means that we humans should be regarded as the ‘perfect’ or ‘ultimate’ form of intelligence. But, until we are able to build machines that can do all the incredible things that we humans can do, as easily and as effortlessly as we do them, it seems rather natural to take human intelligence as a benchmark. In the quest of building intelligent machines, Turing proposed that the goal would be reached when our machines can exhibit intelligent behaviours indistinguishable from that of humans.

Consider the task of classifying images of cats and dogs mentioned earlier. Say that each image is made of $1000$ by $1000$ pixels, each pixel being described by three numbers (in the RGB color representation). Thus, each image can be represented by a vector of three million numbers. The goal is to come up with a program which, given such a vector $x$ as input, produces ‘dog’ or ‘cat’ as output, accordingly. Because there exist many very different vectors $x$ associated to the concept of ‘dog’, there is no obvious, simple and reliable rule to recognize a dog. To solve this task, the program must combine a very large number of ‘weak rules’.

Early forms of AI consisted of explicit, manually-crafted rules, e.g. depending on formal logic. However, using this methodology to figure out all the weak rules that are necessary to correctly classify images is a really daunting task, given the complexity of real-world images. One of the key features of the brain, which these traditional programs did not have, is its ability to learn from experience and to adapt to the environment. Arthur Samuel introduced a new approach to AI, called machine learning (ML) (Samuel, 1959), that takes inspiration from how we learn. In the ML approach, instead of operating with predetermined (i.e. immutable) instructions, the program is made of flexible rules that depend on adjustable parameters. As we modify the parameters, the program changes. The goal is then to tune these parameters so that the resulting program solves the task we want.

To solve the task of image classification, the ML approach requires to collect lots of examples that specify the correct output (the label) for a given input (the image). Such a collection of examples is called a dataset. Then we use these examples to adjust the parameters of the ML program so that, for each input image, the program produces as output the label associated to that image. Such a procedure to adjust the parameters is called a learning algorithm: the ML program learns from examples to solve the problem. Once trained, the program obtained can then be used to predict outputs for new unseen inputs. The performance of the program is assessed on a separate set of examples called the test dataset.

In the setting of image classification, the data is such that the labels are provided together with the corresponding images. This setting, where the expected result is known in advance for the available data, is called supervised learning. This type of learning is currently the most widely used and successful approach to ML. Depending on the task that we want to solve, and the type and the amount of data that is available, there are two other main machine learning paradigms for training an ML program: unsupervised and reinforcement learning. Unsupervised learning refers to data for which no explicitly identified labels exist. Reinforcement learning refers to the case where no exact labels exist, but a scalar value is available (usually called ‘reward’) that provides some knowledge on whether a proposed output is good or bad.

The foundations of modern neuroscience were laid by Santiago Ramon y Cajal, several decades before AI research started. Cajal was the first to observe the brain’s micro-organisation with a microscope. He observed that the brain consists of disjoint nerve cells (the neurons), not of a continuous network as the proponents of the reticular theory thought before him. Neurons have a very particular shape. Each neuron is composed of three main parts (Figure 2): a large ‘tree’ composed of thousands of branches (the dendrites¹¹1In Greek, the word dendron means tree), a cell body (also called the soma), and a long fiber which extends out of the cell body towards other neurons (the axon). A neuron collects information from other neurons through its dendritic tree. The messages collected in the dendrites converge to the cell body, where they are compiled. After compilation, the neuron sends a unique message, called action potential (or spike), which is carried along its axon away from the cell body. In turn, this message is delivered to other neurons.

While neurons are distinct cells, they come into contact at certain points called synapses (Figure 3). Synapses are junction zones through which neurons communicate. Specifically, each synapse is the point of contact of the axon of a neuron (called pre-synaptic neuron) and the dendrite of another neuron (called post-synaptic neuron). The message traveling through the axon of the pre-synaptic neuron is electrical, but the synapse turns it into a chemical message. The axon terminal of the pre-synaptic neuron contains some sorts of pockets (the vesicles) filled with molecules (the neurotransmitters). When the electrical signal reaches the axon terminal, the vesicles open and the neurotransmitters flow in the small synaptic gap between the two neurons. The neurotransmitters then bind with the membrane of the post-synaptic neuron at specific points (the receptors). A neurotransmitter acts on a receptor as a key in a lock: they open ‘gates’ (called channels) in the post-synaptic membrane. As a result, ions flow from the extra-cellular fluid through these channels and generate a current in the post-synaptic neuron. To sum up, the message coming from the pre-synaptic neuron went from electrical to chemical, back to electrical, and in the process, the message was transmitted to the post-synaptic neuron.

Each synapse is a chemical factory in which numerous elements can be modified: the number of vesicles and their size, the number of receptors and their efficacy, as well as the size and the shape of the synapse itself. All these elements affect the strength with which the pre-synaptic electrical message is transmitted to the post-synaptic neuron. Synapses are constantly modified and these modifications reflect what we learn.

The human brain is composed of around 100 billion ($10^{11}$) neurons, interconnected by a total of around a quadrillion ($10^{15}$) synapses. The brain is a huge parallel computer: in this incredibly complex machine, all the synapses work in parallel – like independent nanoprocessors – to process the messages sent between the neurons. Besides, synapses are modified in response to experience, and in turn these modifications alter our behaviours. Thus, synapses are both the computing units and the memory units of the brain. For every task we do, all our thoughts, memories and all our behaviours emerge from the neural activity generated by this machinery.

One of the fundamental questions of neuroscience is that of figuring out the learning algorithms of the brain: what is the set of rules which translate the experiences we have into synaptic changes, and how do these synaptic changes modify our behaviour? Understanding the brain’s learning algorithms is not only key to understanding the biological basis of intelligence, but would also unlock the development of truly intelligent machines.

Artificial neural networks are ML models that draw inspiration from real brains. Artificial neurons imitate the functionality of biological neurons. These models are highly simplified: they keep some essential ideas from real neurons and synapses but they discard many details of their working mechanisms. The first neuron model was introduced by McCulloch and Pitts (1943), but the idea to use such artificial neurons in machine learning was proposed by Rosenblatt (1958) and Widrow and Hoff (1960). Artificial neurons used today in deep learning are essentially unchanged and rely on the same basic math algebra.

Each neuron $i$ is described by a single number $y_{i}$. This number can be thought of as the firing rate of neuron $i$, that is the rate of spikes sent along its axon. Each synapse is also described by a single number, representing its strength. The strength of the synapse connecting pre-synaptic neuron $j$ to post-synaptic neuron $i$ is denoted $W_{ij}$. These artificial synapses can transmit signals with different efficacies depending on their strength. The neuron calculates a nonlinear function of the weighted sum of its inputs:

The pre-activation $x_{i}=\sum_{j}W_{ij}y_{j}$ is a weighted sum of the messages received from other neurons, weighted by the corresponding synaptic strengths. $x_{i}$ can be though of as the membrane voltage of neuron $i$. $\sigma$ is a function called an activation function, which maps $x_{i}$ onto the firing-rate $y_{i}$.

Such artificial neurons can be combined to form an artificial neural network (ANN). Each neuron in the network receives messages from other neurons (the $y_{j}$’s), compiles them ($x_{i}$), and sends in turn a message to other neurons ($y_{i}$). Thus, a network of interconnected neurons exploits the composition of many elementary operations to form more complex computations. The synaptic strengths (the $W_{ij}$’s), also called weights, play the role of adjustable parameters that parameterize this computation.

Deep learning refers to ANNs composed of multiple layers of neurons (LeCun et al., 2015; Goodfellow et al., 2016). These deep neural networks were inspired by the structure of the visual cortex in the brain, each layer corresponding to a different brain region. One of the core ideas of neural networks is that of distributed representations, the idea that the vector of neuron’s states can represent abstract concepts, by opposition to other approaches to AI that use discrete symbols to represent concepts. In a deep network, each layer of neurons applies specialized operations and transformations on its inputs, with the intuition that each layer builds up more abstract concepts than the previous (Bengio, 2009).

Several families of neural networks emerged in the 1980s. One of these families is that of energy-based models, which includes the Hopfield network (Hopfield, 1982) and the Boltzmann machine (Ackley et al., 1985). In these models, under the assumption that the synaptic weights are symmetric (i.e. $W_{ij}=W_{ji}$ for every pair of neurons $i$ and $j$), the dynamics of the network converges to an equilibrium state, after iterating Eq. 1.1 a large number of times for every neuron $i$. Because of the large number of iterations required, these models tend to be slow. This is one of the reasons why these neural networks have been mostly discontinued today. However, by reinterpreting the equilibrium equation of energy-based models as a variational principle of physics, I believe that these models could be the basis of a new generation of fast, efficient and scalable neural networks grounded in physics. We will come back to this point later in the discussion (Section 1.3).

The family of neural networks that is at the heart of the on-going deep learning revolution is that of differentiable neural networks, which became popular thanks to the discovery of the backpropagation algorithm to train them (Rumelhart et al., 1988). In such neural networks, each operation in the process of computation is differentiable. The earliest models of this kind were feedforward neural networks (e.g. the multi-layer perceptron), wherein the connections between the neurons do not form loops. Recurrent neural networks can also be cast to this category of differentiable neural networks, by unfolding the graph of their computations in time. Since their inception, differentiable neural networks have come a long way. Many novel architectures have been introduced, in particular: convolutional neural networks (Fukushima, 1980; LeCun et al., 1989), Long Short-term Memory (Hochreiter and Schmidhuber, 1997; Graves, 2013), and attention mechanisms (Bahdanau et al., 2014; Vaswani et al., 2017).

The computations performed by a neural network are parameterized by its synaptic weights. The goal is to find weight values for which the computations of the neural network solve the task of interest. One essential idea of machine learning is to introduce a loss function, which provides a numerical measure of how good or bad the computations of the model are, with respect to the task that we want to solve. The goal is then to minimize the loss function with respect to the model weights. For example, in image classification, the computations of the model produce an output which represents a ‘guess’ for the class of that image, and the loss provides a graded measure of ‘wrongness’ between that guess and the actual image label. A smaller value of the loss function means that the model produces an output closer to the desired target. The loss function is minimal when the output is equal to the desired target.

One of the most important ideas of deep learning today is stochastic gradient descent (SGD). Provided that the loss function is differentiable with respect to the network weights, we can use the gradient of the loss function to indicate the direction of the minimum of this function. SGD consists in taking examples from the training set one at a time, and adjusting the network weights iteratively in proportion to the negative of the gradient of the loss function. At each iteration, the network performance (as measured per the loss value) slightly improves.

A key discovery that has greatly eased and accelerated deep learning research is the following. Given a computer program that computes a differentiable scalar function $f:\mathbb{R}^{n}\to\mathbb{R}$, it is possible to automatically transform the program into another program that computes the gradient operator $\nabla f:\mathbb{R}^{n}\to\mathbb{R}^{n}$. The gradient $\nabla f(\theta)$ can then be evaluated at any given point $\theta\in\mathbb{R}^{n}$, with a computational overhead that scales linearly with the complexity of the original program. This technique, known as reverse-mode automatic differentiation (Speelpenning, 1980), provides a general framework for ‘backpropagating loss gradients’ (Rumelhart et al., 1988) in any differentiable computational graph. In the last decade, dozens of deep learning frameworks and libraries have been developed, which exploit reverse-mode automatic differentiation to compute gradients in arbitrary differentiable computational graphs. This includes Theano (Bergstra et al., 2010) – a framework that was developed at Université de Montréal – and the more recent Tensorflow (Abadi et al., 2016) and PyTorch (Paszke et al., 2017) frameworks. The emergence of these deep learning frameworks has considerably accelerated deep learning research, by enabling researchers to quickly design neural network architectures and train them by SGD. Thanks to these frameworks, deep learning researchers can explore the space of differentiable computational graphs much more rapidly, as they seek novel and more effective neural architectures.

It was not trivial to discover that deep neural networks can be trained at all. Until the seminal work of Hinton et al. (2006), common belief was that neural networks with more than two layers were essentially impossible to train. In particular, because the landscape of the loss function associated to a deep network is typically highly non-convex, a common misconception was that gradient-descent methods would likely get stuck at bad local minima. In terms of generalization performance, the large over-parameterization of neural networks was also against general prescriptions from classical statistics and learning theory. One of the surprising discoveries of the deep learning revolution was that, in such highly non-convex and over-parameterized statistical models, provided that the loss landscape has appropriate shape, SGD can solve complex tasks by finding excellent parameter values that generalize well to unseen examples.

Several elements contributed to unlock the training of deep neural networks ; among others: the discovery (Glorot et al., 2011) that the ReLU (‘Rectified Linear Unbounded’) activation function usually outperforms the sigmoid activation function, the discovery of better weight initialization schemes (Glorot and Bengio, 2010; Saxe et al., 2013; He et al., 2015), and the batch-normalization technique to systematically normalize signal amplitudes at each layer of a network (Ioffe and Szegedy, 2015). Besides, fundamental advances have come from the introduction of new network architectures such as the ones mentioned earlier, and from novel machine learning paradigms such as that of generative adversarial networks (Goodfellow et al., 2014). All these techniques have as an effect to modify the landscape of the loss function, as well as the starting parameter (before training) in the parameter space. Understanding how the landscape of the loss function can be appropriately shaped to ease optimization by SGD is an active area of research (Poggio et al., 2017; Arora, 2018).

In recent years, deep neural networks have proved capable of solving very complex problems across a wide variety of domains. Today, they achieve state-of-the-art performance in image recognition (He et al., 2016), speech recognition (Hinton et al., 2012; Chan et al., 2016), machine translation (Vaswani et al., 2017), image-to-text (Sharma et al., 2018), text-to-speech (Oord et al., 2016), text generation (Brown et al., 2020), and synthesis of realistic-looking portraits (Karras et al., 2019), among many other applications. Neural networks have become better than humans at playing Atari games (Mnih et al., 2013), playing the game of Go (Silver et al., 2018) and playing Starcraft (Vinyals et al., 2019). Perhaps what is most exciting is that, although these neural networks are designed to solve very different tasks and deal with different types of data, they are all trained using the same handful of basic principles. As we scale these neural networks, more advanced aspects of intelligence seem to emerge from this handful of principles.

While the pace of progress in neural network research is breathtaking, we should also emphasize that, without any question, neural networks are nowhere close to ‘surpass’ humans. In their current form, they miss key elements of human intelligence. Whenever neural networks beat humans, they beat us at a very specific task and/or under very specific conditions. We may need new learning paradigms to move away from ‘task-specific’ neural networks towards ‘multi-functional’ and continually learning neural networks. Besides, as of today, neural networks cannot handle and combine abstract concepts nearly as flexibly as humans do. Making progress along these lines may require new breakthroughs, to give these neural networks the ability to develop a ‘thought language’ and a sense of causal reasoning, among others.

The on-going deep learning revolution owes to other technological developments too. In the past decades, the amount of data available has greatly increased, and powerful multi-core parallel processors for general purpose computing, such as Graphics Processing Units (GPUs), have emerged (Owens et al., 2007). Training large models on large datasets – here is part of the recipe to make a deep neural network solve a challenging task. In the 1980s, when deep neural networks were first conceived, the lack of data and computational power to train them made it practically unfeasible to demonstrate their effectiveness.

The last decade of neural network research seems to hint at a simple and straightforward strategy to further improve performance of our AI systems: scaling. Using more memory and more compute to train larger models on larger datasets – here is the current trend to build state-of-the-art deep learning systems. The largest model ever built thus far, GPT-3 (Brown et al., 2020), has a capacity of 175 billion parameters.

Training these neural networks requires very large amounts of computations. Standard practice today is to distribute the computations across more and more GPUs, to train larger and larger models. Still, even using thousands of GPUs working in parallel, training these neural networks can take months. For example, AlphaZero learnt to play the game of Go by playing 140 million games, which took 5000 processors and two weeks. Moreover, training such models can cost millions of dollars, just for electricity consumption, not to mention their ecological impact. Yet, even these large neural networks are only a tiny fraction of the size of the human brain. What causes such inefficiency, preventing us from building models of the size of the human brain?

If at the conceptual level the neural networks used today take their overall strategy from the brain, on the hardware implementation level however, they use little of the cleverness of nature. Our current processors, on which these neural networks are trained and run, operate in fundamentally different ways than brains. They rely on the von Neumann architecture. In this computer architecture, the memory unit where information is stored, is separated from the processing unit where calculations are done. A so-called bus moves information back and forth between these two units. Over the course of history of computing, this computer architecture has become the norm, and today, the von Neumann architecture is used in virtually all computer systems: laptops, smartphones, and all kinds of embedded systems. The GPUs, massively used for neural network training today, also rely on the von Neuman architecture, where memory is separated from computing.

The brain on the other hand deeply merges memory and computing by using the same functional unit: the synapse. The human brain is composed of a quadrillion ($10^{15}$) synapses. In other words, the human brain has $10^{15}$ nanoprocessors working in parallel.

Training neural networks on von Neumann hardware as we do it today is extraordinarily energy inefficient in comparison with the way brains operate. The necessity to move the data back and forth between the memory and processing units in the von Neumann architecture is energy intensive and creates considerable latency. This limitation is known as the von Neumann bottleneck. How inefficient is it compared to biological systems like the brain? The human brain is composed of $10^{11}$ neurons and consumes around 20W to conduct all of its activities (Attwell and Laughlin, 2001). In comparison, training a BERT model (a state-of-the-art natural language processing model) on a modern supercomputer requires 1500 kW.h (Strubell et al., 2019), which is the total amount of energy consumed by a brain in nine years. Besides, a GPU running real-time object detection with YOLO (Redmon et al., 2016), a network smaller than the brain by four orders of magnitude, consumes around 200W.

This striking mismatch holds more broadly with biological systems in general. For example, Kempes et al. (2017) study the energy efficiency of ‘cellular computation’ in the process of biological translation. What is the amount of energy required (in ATP equivalents) by ribosomes to assemble amino acids into proteins ? They point out that "the best supercomputers perform a bit operation at roughly $5.27\times 10^{-13}J$, […] which is about five orders of magnitude less efficient than biological translation."

To sum up, our current neural networks are orders of magnitude less energy efficient than biological systems at processing information, and the von Neumann bottleneck is largely responsible for this inefficiency. If using more and more GPUs may increase speed, and thereby speed up training and inference of neural networks, this strategy however can’t improve energy efficiency.

In order to build massively parallel neural networks that are energy efficient and can scale to the size of the human brain, we need to fundamentally rethink the underlying computing hardware. We need to design neural networks so that computations are performed at the physical location of the synapses, where the strength of the connections (the weights of the neural network) are stored and adjusted, just like in the brain. The concept of hardware that merges memory and computing is called in-memory computing (or in-memory processing), and the field tackling this problem is called neuromorphic computing. This field of research, started by Carver Mead in the 1980s (Mead, 1989) aims at mimicking brains at a hardware level, by building physical neurons and synapses onto a chip.

The most common approach to in-memory computing today is to use programmable resistors as synapses. Programmable resistors, such as memristors (Chua, 1971), are resistors whose conductance can be changed (or ‘programmed’). The weights of a neural network can be encoded in the conductance of such devices. In the last decade, important advances in nanotechnology were made, and a number of new technologies have emerged and have been studied as potentially promising programmable resistors (Burr et al., 2017; Xia and Yang, 2019).

Neuromorphic computing thus explores analog computations that fundamentally depart from the standard digital computing paradigm.

Analog processing differs from digital processing in important ways. Whereas digital circuits manipulate binary signals with reliably distinguishable on and off-states, analog circuits on the other hand manipulate real-valued currents and voltages that are subject to analog noise that bounds the precision with which computation may be performed. More importantly, analog devices suffer from mismatches, i.e. small random variations in the physical characteristics of devices, which occurs during their manufacturing. No two devices are exactly alike in their characteristics, and it is impossible to make a perfect clone of one. These variations result in behavioral differences between identically designed devices. Due to the accumulation of the mismatch errors from individual devices, it is very difficult to analytically predict the behavior of a large analog circuit.

A growing field of research in the neuromorphic literature attempts to perform in analog the operations that we normally do in software, so as to implement feedforward neural networks and the backpropagation algorithm efficiently. In this approach, the starting point is an equation of the kind of Eq. 1.1. Many of these operations are then performed in analog and combined to form the computations of a feedforward network. However, because of device mismatches, it is hard to perform such idealized operations, and as we combine many of these operations, the resulting computation may be different from the desired one. Either these idealized operations are performed with low precision, or we may spend a lot of energy trying to improve precision, e.g. by using analog-to-digital conversion.

Not coincidentally, the constraint of device nonidealities and device variability is shared with biology too. No two neurons are exactly the same. This realization demonstrates, in principle, that it is possible to train (biological) neural networks even in the presence of noise and imperfect ‘devices’. It invites us to rethink the learning algorithm for neural networks (and the notion of computation altogether).

The general idea of the manuscript is the following. We consider a physical system composed of multiple parts whose characteristics and working mechanisms may be only partially known. The system has some ‘adjustable parameters’, some of them playing the role of ‘inputs’, and we may read or measure a ‘response’ on some other ‘output’ part of the system. We can think of this black box system as performing computations and implementing a nonlinear input-to-output mapping function (which may be analytically unknown). We wish to tune the adjustable parameters of this system by stochastic gradient descent (SGD), as we normally do in deep learning. The question is: how can we compute or estimate the gradients in such a physical system, by relying on the physics of the system?

The main theoretical result of the thesis is that, for a large class of physical systems (those whose state or dynamics derive from a variational principle), there is a simple procedure to estimate the parameter gradients, which in many practical situations requires only locally available information for each parameter. This procedure, called equilibrium propagation (EqProp), preserves the key benefit of being compatible with SGD, while offering the possibility to directly exploit physics to implement and train neural networks.

Rather than Eq. 1.1, our starting point is a variational equation of the form

where $E$ is a scalar function. If $s$ is the state of the system, then Eq. 1.2 is an equilibrium condition. In this case, we say that the system is an energy-based model (EBM) and we call $E$ the energy function.

In this thesis, we also introduce the concept of Lagrangian-based model (LBM). Variational equations exist not just to characterize equilibrium states, but also entire trajectories. Many physical systems are such that their trajectory derives from a principle of stationary action (e.g. a principle of least action). Denoting $\mathrm{s_{t}}$ the state of the system at time $t$, this means that the (continuous-time) trajectory $\mathrm{s}=\{\mathrm{s_{t}}\}_{0\leq t\leq T}$ over a time interval $[0,T]$ minimizes a functional of the form

where $\mathrm{\dot{s}_{t}}$ is the time derivative of $\mathrm{s_{t}}$, $L$ is a function called the Lagrangian function of the system, and $\mathcal{S}$ is a scalar functional called the action. The stationarity of the action tells us that $\frac{\delta\mathcal{S}}{\delta\mathrm{s}}=0$, which is another variational equation of the kind of Eq. 1.2. These systems, which we call Lagrangian-based models (LBMs), are suitable in particular in the setting with time-varying inputs and can thus play the role of ‘recurrent neural networks’.

Equilibrium propagation (EqProp) allows to compute gradients with respect to arbitrary loss functions in these EBMs and LBMs. Furthermore, if the energy function (resp. Lagrangian function) of the system has a property called sum-separability, meaning that it is the sum of the energies (resp. Lagrangians) of its parts, then computing the loss gradients with EqProp requires only information that is locally available for each parameter (i.e. the learning rule is local).

In the 1650s, Pierre de Fermat proposed the principle of least time which states that, between two given points, the light travels the path which takes the least time. He showed that both the laws of reflection and refraction can be derived from this principle. Fermat’s least time principle is an instance of what we call more generally a variational principle. Today, the variational approach pervades much of modern physics and engineering (Lanczos, 1949), with applications not only in optics, but also in mechanics, electromagnetism, thermodynamics, etc. Even at a fundamental level of description, our universe seems to behave according to variational principles: for example, Einstein’s equations of general relativity can be derived from the Einstein-Hilbert action, and in a sense, Feynman’s path integral formulation of quantum mechanics can be seen as a generalized principle of least action (Feynman, 1942).

Thus, many physical systems qualify as energy-based models or Lagrangian-based models. This offers in principle a lot of options for implementing our proposed method on physical substrates. In this manuscript we will present one such option in details: nonlinear resistive networks.

Interestingly, our theoretical framework invites us to rethink not only the von Neumann architecture on which our current processors rely, but also the notion of computation altogether.

Much of computer science deals with computation in the abstract, without worrying about physical implementation (Lee, 2017). Our computers today rely on the computing paradigm introduced by Turing, where computers operate on digital data and carry out computations algorithmically, via step-by-step (discrete-time) processes. The von Neumann architecture was invented to implement these computations (i.e. to bridge the gap between physics and these abstract computations), but suffers from the speed and energy efficiency problems mentioned earlier (Section 1.2.8).

The theoretical framework presented in this manuscript can be seen as an alternative approach to computing that takes advantage of the ways by which Nature operates. We suggest a novel computing paradigm, which uses the variational formulations of the laws of physics as first principles. Together with the equilibrium propagation (EqProp) training procedure, our approach suggests a way to implement the core principles of deep learning by exploiting physics directly. As will become apparent in Section 4.3.2, the process of ‘computations’ in EqProp is very different from the step-by-step processes of Turing’s conventional computing paradigm. Although we may call EqProp a learning ‘algorithm’, it is not an algorithm in the conventional sense (one that performs step-by-step computations).

More technically, the main idea of the manuscript can be summarized by the following mathematical identity – a novel differentiation method compatible with variational principles. Consider two functions $E(\theta,s)$ and $C(s)$ of the two variables $\theta$ and $s$. We wish to compute the gradient $\frac{d}{d\theta}C(s_{\theta})$, where $s_{\theta}$ is such that $\frac{\partial E}{\partial s}(\theta,s_{\theta})=0$. To do this, we introduce the function $F(\theta,\beta,s)=E(\theta,s)+\beta\;C(s)$, where $\beta$ is a scalar. For fixed $\theta$, we further define $s_{\star}^{\beta}$ by the relationship $\frac{\partial F}{\partial s}(\theta,\beta,s_{\star}^{\beta})=0$, for any $\beta$ in the neighborhood of $\beta=0$. In particular $s_{\star}^{0}=s_{\theta}$. Then we have the identity

where $\frac{\partial E}{\partial\theta}$ denotes the partial derivative of the function $E$ with respect to its first argument. This result is developed and proved in Chapter 2.

We have seen in Chapter 2 that EqProp is an algorithm to train systems that possess equilibrium states, i.e. states characterized by a variational equation of the form $\frac{\partial E}{\partial s}\left(\theta,x,s_{\star}\right)=0$, where $E\left(\theta,x,s\right)$ is a scalar function called energy function. Recall that $\theta$ is the set of adjustable parameters of the system, and $x$ is an input. We have called such systems energy-based models. The class of energy-based models that we study here is that of systems whose dynamics spontaneously minimizes the energy function $E$ by following its gradient. In such a system, called gradient system, the state follows the dynamics

Here $s_{t}$ denotes the state of the system at time $t$. The energy of the system decreases until $\frac{ds_{t}}{dt}=0$, and the equilibrium state $s_{\star}$ reached after convergence of the dynamics is an energy minimum (either local or global). The function $E$ is also sometimes called a Lyapunov function for the dynamics of $s_{t}$.

In this setting, equilibrium states are stable: if the state is slightly perturbed around equilibrium, the dynamics will tend to bring the system back to equilibrium. For this reason, such equilibrium states are also called ‘attractors’ or ‘retrieval states’, because the system’s dynamics can ‘retrieve’ them if they are only partially known. Thus, gradient systems can recover incomplete data, by storing ‘memories’ in their point attractors.

In this manuscript, we are more specifically interested in the supervised learning problem, where the loss to optimize is of the form

After training is complete, the model can be used to ‘retrieve’ a label $y$ associated to a given input $x$.

In a gradient system, EqProp takes the following form. In the first phase, or free phase, the state of the system follows the gradient of the energy (Eq. 3.1). At the end of the first phase the system is at equilibrium ($s_{\star}$). In the second phase, or nudged phase, starting from the equilibrium state $s_{\star}$, a term $-\beta\;\frac{\partial C}{\partial s}$ (where $\beta>0$ is a hyperparameter called nudging factor) is added to the dynamics of the state and acts as an external force nudging the system dynamics towards decreasing the cost $C$. Denoting $s_{t}^{\beta}$ the state of the system at time $t$ in the second phase (which depends on the value of the nudging factor $\beta$), the dynamics is defined as¹¹1As discussed in the case of Eq. 2.6, we can also define $\frac{ds_{t}^{\beta}}{dt}=-\frac{\partial E}{\partial s}\left(\theta,x,s_{t}^{\beta}\right)-\beta\;\frac{\partial C}{\partial s}\left(s_{\star},y\right)$ without changing the conclusions of the theoretical results.

The system eventually settles to a new equilibrium state $s_{\star}^{\beta}$. Recall from Theorem 2.1 that the gradient of the loss $\mathcal{L}$ can be estimated based on the two equilibrium states $s_{\star}$ and $s_{\star}^{\beta}$. Specifically, in the limit $\beta\to 0$,

Furthermore, if the energy function has the sum-separability property (as defined by Eq. 2.5), then the learning rule for each parameter $\theta_{k}$ is local:

Note that the learning rule of Eqs. 3.4-3.5 only depends on the equilibrium states $s_{\star}$ and $s_{\star}^{\beta}$, not on the specific trajectory that the system follows to reach them. Indeed, as we have seen in the previous chapter, EqProp applies to any energy based model, not just gradient systems. But under the assumption of a gradient dynamics, we can say more about the transient states, when the system gradually moves from the free state ($s_{\star}$) towards the nudged state ($s_{\star}^{\beta}$): we show that the transient states ($s_{t}^{\beta}$ for $t\geq 0$) perform gradient computation with respect to a function called the projected cost function.

Recall that in the free phase, the system follows the dynamics of Eq. 3.1. In particular, the state $s_{t}$ at time $t\geq 0$ depends not just on $\theta$ and $x$, but also on the initial state $s_{0}$ at time $t=0$. Let us define the projected cost function

where we omit $x$ and $y$ for brevity of notations. $L_{t}(\theta,s_{0})$ is the cost of the state projected a duration $t$ in the future, when the system starts from $s_{0}$ and follows the dynamics of the free phase. Note that $L_{t}$ depends on $\theta$ and $s_{0}$ (as well as $x$) implicitly through $s_{t}$. For fixed $s_{0}$, the process $\left(L_{t}(\theta,s_{0})\right)_{t\geq 0}$ represents the successive cost values taken by the state of the system along the free dynamics when it starts from the initial state $s_{0}$. In particular, for $t=0$, the projected cost is simply the cost of the initial state, i.e. $L_{0}(\theta,s_{0})=C\left(s_{0}\right)$. As $t\to\infty$, we have $s_{t}\to s_{\star}$ and therefore $L_{t}(\theta,s_{0})\to C(s_{\star})=\mathcal{L}(\theta)$, i.e. the projected cost converges to the loss at equilibrium.

The following result shows that the transient states of EqProp ($s_{t}^{\beta}$) can be expressed in terms of the projected cost function ($L_{t}$), when $s_{0}=s_{\star}$.

The left-hand-side of Eq. 3.7 represents the gradient provided by EqProp if we substitute $s_{t}^{\beta}$ to $s_{\star}^{\beta}$ in the gradient formula (Eq. 3.4). This corresponds to a truncated version of EqProp, where the second phase (nudged phase) is halted before convergence to the nudged equilibrium state. Eq. 3.7 provides an analytical formula for this truncated gradient in terms of the projected cost function, when $s_{0}=s_{\star}$

The left-hand side of Eq. 3.8 is the temporal derivative of $s_{t}^{\beta}$ rescaled by a factor $\frac{1}{\beta}$. In essence, Eq. 3.8 shows that, in the second phase of EqProp (nudged phase), the temporal derivative of the state codes for gradient information (namely the gradients of the projected cost function, when $s_{0}=s_{\star}$).

In this section, we prove Theorem 3.1. In doing so, we also establish a link between EqProp and the recurrent backpropagation algorithm of Almeida (1987) and Pineda (1987), which we briefly present below.

First, let us introduce the temporal processes $(\overline{S}_{t},\overline{\Theta}_{t})$ and $(\widetilde{S}_{t},\widetilde{\Theta}_{t})$ defined by

and

The processes $\overline{S}_{t}$ and $\widetilde{S}_{t}$ take values in the state space (space of the state variable $s$). The processes $\overline{\Theta}_{t}$ and $\widetilde{\Theta}_{t}$ take values in the parameter space (space of the parameter variable $\theta$). Using these notations, Theorem 3.1 states that $\overline{S}_{t}=\widetilde{S}_{t}$ and $\overline{\Theta}_{t}=\widetilde{\Theta}_{t}$ for every $t\geq 0$. These identities are a direct consequence of the following lemma.

We refer to Scellier and Bengio (2019) for a proof of this result. Since the differential equation of Lemma 3.2 is linear with constant coefficients, we can express $S_{t}$ and $\Theta_{t}$ using closed-form formulas. Specifically $S_{t}=\exp\left(-t\frac{\partial^{2}E}{\partial s^{2}}\left(\theta,s_{\star}\right)\right)\cdot\frac{\partial C}{\partial s}\left(s_{\star}\right)$ and $\Theta_{t}=\frac{\partial^{2}E}{\partial\theta\partial s}\left(\theta,s_{\star}\right)\cdot\left(\frac{\partial^{2}E}{\partial s^{2}}\left(\theta,s_{\star}\right)\right)^{-1}\cdot\left[\textrm{Id}-\exp\left(-t\frac{\partial^{2}E}{\partial s^{2}}\left(\theta,s_{\star}\right)\right)\right]\cdot\frac{\partial C}{\partial s}\left(s_{\star}\right)$.

Lemma 3.2 also suggests an alternative procedure to compute the parameter gradients of the loss $\mathcal{L}$ numerically. This procedure, known as Recurrent Backpropagation (RBP), was introduced independently by Almeida (1987) and Pineda (1987). Specifically, RBP consists of the following two phases. The first phase is the same as the free phase of EqProp: $s_{t}$ follows the free dynamics (Eq. 3.1) and relaxes to the equilibrium state $s_{\star}$. The state $s_{\star}$ is necessary for evaluating $\frac{\partial^{2}E}{\partial s^{2}}\left(\theta,s_{\star}\right)$ and $\frac{\partial^{2}E}{\partial\theta\partial s}\left(\theta,s_{\star}\right)$, which the second phase requires. In the second phase, $S_{t}$ and $\Theta_{t}$ are computed iteratively for increasing values of $t$ using Eq. 3.12-3.15. Finally, $\Theta_{t}$ provides the desired loss gradient in the limit $t\to\infty$. To see this, we first note that Lemma 3.2 tells us that the vector $\Theta_{t}$ computed by this procedure is equal to $\widetilde{\Theta}_{t}$ for any $t\geq 0$. Then by definition, $\widetilde{\Theta}_{t}=\lim_{\beta\to 0}\frac{1}{\beta}\left(\frac{\partial E}{\partial\theta}\left(\theta,s_{t}^{\beta}\right)-\frac{\partial E}{\partial\theta}\left(\theta,s_{\star}\right)\right)$. It follows that, as $t\to\infty$, we have $\widetilde{\Theta}_{t}\to\lim_{\beta\to 0}\frac{1}{\beta}\left(\frac{\partial E}{\partial\theta}\left(\theta,s_{\star}^{\beta}\right)-\frac{\partial E}{\partial\theta}\left(\theta,s_{\star}\right)\right)=\frac{\partial\mathcal{L}}{\partial\theta}$.

An important benefit of EqProp over RBP it that EqProp requires only one kind of dynamics for both phases of training. RBP requires a special computational circuit in the second phase for computing the gradients.

The original RBP algorithm was described for a general state-to-state dynamics. Here, we have presented RBP in the particular case of gradient dynamics. We refer to LeCun et al. (1988) for a more general derivation of RBP based on the adjoint method.

A Hopfield network is a neural network with the following characteristics. The state of a neuron $i$ is described by a scalar $s_{i}$, loosely representing its membrane voltage. The state of a synapse connecting neuron $i$ to neuron $j$ is described by a real number $W_{ij}$ representing its efficacy (or ‘strength’). The notation $\sigma(s_{i})$ is further used to denote the firing rate of neuron $i$. The function $\sigma$ is called activation function ; it takes a real number as input, and returns a real number as output. Using the formalism of the previous section, the state of the system is the vector $s=\left(s_{1},s_{2},\ldots,s_{N}\right)$ where $N$ is the number of neurons in the network, and the set of parameters to be adjusted is $\theta=\{W_{ij}\}_{ij}$. As we will see shortly, one biologically unrealistic requirement of the Hopfield model is that synapses are assumed to be bidirectional and symmetric: the synapse connecting $i$ to $j$ shares the same weight value as the synapse connecting $j$ to $i$, i.e. $W_{ij}=W_{ji}$.