2023

[1]Jérémie Laydevant

[1]Julie Grollier

1]Unité Mixte de Physique CNRS, Thales, Université Paris-Saclay, 91767
Palaiseau, France

# Training an Ising Machine with Equilibrium Propagation

<EMAIL_ADDRESS>Danijela Marković  danijela.markovic@cnrs-
thales.fr<EMAIL_ADDRESS>[

###### Abstract

Ising machines, which are hardware implementations of the Ising model of
coupled spins, have been influential in the development of unsupervised
learning algorithms at the origins of Artificial Intelligence (AI). However,
their application to AI has been limited due to the complexities in matching
supervised training methods with Ising machine physics, even though these
methods are essential for achieving high accuracy. In this study, we
demonstrate a novel approach to train Ising machines in a supervised way
through the Equilibrium Propagation algorithm, achieving comparable results to
software-based implementations. We employ the quantum annealing procedure of
the D-Wave Ising machine to train a fully-connected neural network on the
MNIST dataset. Furthermore, we demonstrate that the machine’s connectivity
supports convolution operations, enabling the training of a compact
convolutional network with minimal spins per neuron. Our findings establish
Ising machines as a promising trainable hardware platform for AI, with the
potential to enhance machine learning applications.

###### keywords:

Machine learning, Ising machines, Equilibrium Propagation, Neuromorphic
computing

##  1 Introduction

Investigating physical systems that can execute cognitive tasks based on their
dynamic behaviors or statistical features has long been a topic of interest in
physics, motivated by the quest to unravel the brain’s learning capabilities
review-physics-neuro . The Ising system ising_beitrag_1925  of coupled spins,
described by the Ising energy function111We write here the Energy function
that the D-Wave Ising machine really minimizes harris_experimental_2010 ,
which differs from the classical Ising Energy function by a minus sign, in
order to be consistent with all the equations and learning rules in the rest
of the paper.:

$E_{Ising}=\sum_{i>j}J_{ij}\sigma_{i}\sigma_{j}+\sum_{i}h_{i}\sigma_{i},$ (1)

has played a significant role in these developments little_existence_1974 ;
Hopfield1982 ; Amit1985 ; Mezard1986 . It can indeed be likened to a neural
network, where the state of a spin $\sigma_{i}$ (up or down) corresponds to
the activity of a binary neuron, the value of the coupling between spins
$J_{ij}$ corresponds to the strength of the synaptic connection between the
neurons they emulate, and the bias fields $h_{i}$ applied to individual spins
correspond to the biases of the artificial neurons. The majority of learning
demonstrations on Ising machines adachi_application_nodate ; Benedetti2017 ;
dorband_boltzmann_2015 ; liu_adiabatic_2018 ; adachi ; dixit_training_2021 ;
niazi2023training  have focused on implementing Boltzmann machines methods
hinton_optimal_nodate ; hinton_training_2002 . This algorithm capitalizes on
the properties of physical systems, characterized by an energy function, to
evolve towards an equilibrium state governed by Boltzmann statistics. However,
Boltzmann machines, as generative models, underperform in difficult
classification tasks when compared to standard supervised learning algorithms
like backpropagation krizhevsky_convolutional_nodate .

The recent boom in AI, driven by the advent of these highly effective
supervised algorithms, has led to the development of various new hardware
platforms for AI applications. These platforms aim to address the growing
power consumption and computational demands associated with both training and
inference phases in AI systems. Emerging AI hardware solutions exploit local
memory and leverage the unique physical phenomena exhibited by novel
components wang-resistive-switching ; xia_memristive_2019 ;
markovic_physics_2020  . However, these new platforms face compatibility
challenges with the most efficient supervised training methods, which rely on
the minimization of a global cost function, such as error backpropagation.
This is due to the intrinsically non-local nature of these methods and the
fact that the calculation of associated gradients is based on mathematical
procedures that do not correspond to the physics of the emerging devices used.
As a result, a significant research effort is currently underway to develop
algorithms capable of training these novel, AI-dedicated hardware platforms
efficiently nokland_direct_2016 ; neftci_surrogate_2019 ; martin_eqspike_2020
; kendall_training_2020 ; Frenkel2021 ; ernoult_towards_2022 ;
wright_deep_2022 ; schuman2022 .

Equilibrium Propagation scellier_equilibrium_2017 , introduced in 2017, has
garnered significant attention in this field of research for its ability to
train physical systems kendall_training_2020 ; martin_eqspike_2020 ; dillavou-
PRA ; Yi2022  in a supervised way using a local learning rule that closely
approximates the gradients of state-of-the-art software methods, such as
backpropagation through time NEURIPS2019_67974233 . The algorithm requires
that the physical system evolves towards a stable equilibrium state through
the minimization of an energy function such as:

$E_{EP}=\frac{1}{2}\sum_{i}s_{i}^{2}-\sum_{i>j}W_{ij}\rho(s_{i})\rho(s_{j})-\sum_{i}b_{i}\rho(s_{i}),$
(2)

where $s_{i}$, $s_{j}$ are the real and continuous states of the neurons,
$W_{ij}$ are the symmetric synaptic weights connecting neurons $i$ and $j$,
$b_{i}$ are individual biases applied to the neurons, and $\rho$ is a non-
linear activation function such as tanh. The first term in Eq. (2) is a
damping term that allows the system to reach a stable equilibrium state. In
Equilibrium Propagation, the inference is performed by conditioning the steady
state of the system with input values (free phase), while learning is achieved
by dynamically perturbing the outputs to align them with the desired values
(nudge phase) and thus minimize the objective loss function $\mathcal{L}$. The
parameter changes required for learning are derived from local measurements of
the equilibrium states, as opposed to a complex non-local analytic
mathematical procedure like backpropagation. The corresponding learning rule
for synaptic weights writes:

$-\frac{\partial\mathcal{L}}{\partial W_{ij}}=\Delta
W_{ij}\propto\frac{1}{\beta}\left[(\rho(s_{i})\rho(s_{j}))^{*,nudge}-(\rho(s_{i})\rho(s_{j}))^{*,free}\right],$
(3)

with a term $\beta$ that characterizes the strength of the nudging force. It
favors the equilibrium state obtained after the nudge phase, whose outputs are
closer to the target value, and destabilizes the one obtained after the free
phase. This is accomplished by decreasing the energy of the nudge state and
increasing the energy of the free state. Numerical simulations of Equilibrium
Propagation have demonstrated state-of-the-art performance on benchmark tasks
such as MNIST scellier_equilibrium_2017 ; NEURIPS2019_67974233 , CIFAR-10
laborieux_scaling_2020 , and Image-net-32 laborieux-holo , using both fully-
connected networks and convolutional architectures.

Equilibrium Propagation is therefore an excellent candidate for training
physical systems described by an energy function dillavou-PRA ; Yi2022 . Ising
machines mohseni_ising_2022  are particularly interesting hardware systems for
this purpose, as they are designed to find the ground state of the Ising
model. Moreover, they offer thousands of spins, and their reconfigurable
coupling parameters facilitate training. However, their applications are
mostly limited to solving combinatorial problems with fixed parameters
lucas_ising_2014 . Training Ising machines using Equilibrium Propagation would
broaden their scope to supervised classification tasks and leverage their
adjustable parameters. Nonetheless, three fundamental differences exist
between the Ising model (Eq. (1)) and the EP model (Eq. (2)), which present
challenges for training.

First, the Ising energy function (Eq. (1)) lacks the damping term present in
the Equilibrium Propagation (EP) model (Eq. (2)), which allows the latter to
reach steady state equilibrium intrinsically. Ising machines can approach the
ground state using various annealing methods cai_power-efficient_2020 ;
Farhi2001 ; harris_experimental_2010  or minimum gain principle
yamamoto_coherent_2017 , but destabilizing it for the nudge phase in EP
remains challenging. Developing methods to gently manipulate the equilibrium
state is necessary.

Second, the difference between Ising spins’ two-state nature and EP’s
continuous-state neurons poses a challenge. Approaches must be devised to
create smooth modifications of the spin system by the outputs, emulating the
gradual changes in a neural network’s learning process. This would bridge the
gap between the two-state Ising machine and continuous-state EP neurons,
enabling efficient training.

Third, implementing EP on Ising machines involves addressing the physical
connectivity challenge. Neural networks have dense connectivity, while
physical spin systems have sparser connections. Strategies need to be
established to compensate for reduced connectivity or adapt the network
architecture to the Ising hardware’s connectivity.

Figure 1: Training the Ising Machine with EP. a. Illustration of the free
phase and nudge phase of the Equilibrium Propagation algorithm applied to an
Ising spin system. For both phases, the input is fed to the chip through bias
fields (see Section 2.2 and 2.3), with a strength that depends on the task.
The steady spins states obtained at equilibrium after the free and the nudge
phases can be directly measured on the chip to compute the parameters updates.
b. Annealing schedule used to drive the Ising machine during the two
sequential phases of EP. At the end of both phases, the probability of
transition between states ends at 0, in the steady state where we measure the
states of all the spins. The small bump in the probability during the nudge
phase, achieved through reverse annealing, allows the system to be sensitive
to the nudge signal applied to the output neurons c. Binary activations - such
as spins - in a dynamical neural network trained with EP can cause the
vanishing gradient issue if the input of the neuron is only weakly modified
between the free and the nudge phase d. Schematic of the D-Wave chip used with
the specific Chimera architecture where spins are arranged as small 4x4 fully
connected square lattices and laterally coupled to 2 neighbors. The spins
$\sigma_{i}$ are represented as horizontal and vertical lines whereas the
couplings $J_{ij}$ are represented as plain circles for intra-cluster coupling
and dotted bridges for cluster-to-cluster coupling.

In this study, we report a critical advancement towards utilizing Ising
machines for machine learning applications. Employing the commercial D-Wave
Ising machine harris_experimental_2010 , composed of thousands of two-state
components, and the Equilibrium Propagation algorithm
scellier_equilibrium_2017 , we successfully recognize handwritten digits from
the MNIST/100 database truncated-mnist , achieving recognition rates
comparable to a fully-connected network trained using software simulations on
standard digital hardware.

We train the Ising machine with Equilibrium Propagation, taking advantage of
its capacity to reach the ground state of its energy function through
annealing during the free phase and reverse-annealing during the nudge phase.
Moreover, we implement gradual changes to the spin system’s ground state, as
required for smooth steps in gradient-based optimization, by increasing the
number of outputs and adapting procedures developed for training dynamic
binary activation networks Laydevant_2021_CVPR .

Finally, we demonstrate that the connectivity between near-neighbor spins on
the DW-2000 chip, featuring the Chimera architecture harris_experimental_2010
, is inherently compatible with convolutional operations. We successfully
train a compact convolutional network entirely on the chip, achieving
recognition rates on par with software performance.

Our results indicate that Ising machines hold significant potential as machine
learning hardware, with their physics allowing for inference, error
backpropagation, and gradient computation. Furthermore, our findings highlight
the promise of physics-based learning algorithms, such as Equilibrium
Propagation, in training fully connected and convolutional networks on
emerging hardware.

##  2 Results

###  2.1 Training an Ising Machine with EP through annealing

In this study, we have chosen the commercial Ising machine D-Wave as the
demonstration platform for our algorithm, owing to its distinct advantages
compared to other publicly available IMs. Specifically, D-Wave offers a large
number of spins, ranging from 2000-5000 depending on the version (2000Q or
Advantage 4-5), high-precision coupling parameters (4-6 bits), and the ability
to control these parameters online through a Python interface. This interface
is fully compatible with the code developed for the training algorithm, which
is essential, as the parameters need iterative adjustments during the training
process.

We illustrate the implemented training procedure in Fig. 1a and provide a
comprehensive description of how we execute the free phase, nudge phase, and
learning rule in the subsequent sections. For both the free and nudge phases,
we introduce the input data to the spin system through bias fields. In the
following sections, we elaborate on how the bias fields are set according to
the input data and the architecture implemented on the chip.

In each phase, the input bias fields are maintained constant, allowing the
system to stabilize at an equilibrium state conditioned by the input data.
Ising machines employ extrinsic mechanisms, such as simulated annealing
kirkpatrick_optimization_1983 , noise annealing cai_power-efficient_2020 , or
quantum annealing harris_experimental_2010 , to guide the spin system towards
its ground state rather than local energy minima. These annealing procedures
regulate the system’s exploration of its energy landscape by adjusting a
probability parameter that dictates the system’s capacity to escape a given
configuration. This probability is controlled by the temperature of the system
in simulated annealing, by the noise in the system in noise annealing, and by
the tunneling rate between states in quantum annealing. For this
demonstration, we utilized the quantum annealing procedure of the D-Wave chip
to achieve the free phase of EP (Fig. 1a (free phase) and Fig. 1b
($0<time<20\mu s$)). The ground state is obtained by progressively reducing
the probability from a high value where the system explores all possible
configurations, towards zero probability (Fig. 1b: $0<time<20\mu s$).

The nudge phase of Equilibrium Propagation involves adding to the energy of
the system a term proportional to the cost function ($C$) that describes the
discrepancy between the output neurons and the target. Here, we use as a cost
function the Mean Square Error (see Methods 5.2) that represents the deviation
between the activation functions of the output neurons and their target
states:

$C(\rho(y),\rho(\hat{y}))=\sum_{i\in Y}(\rho(y_{i})-\rho(\hat{y})_{i})^{2},$
(4)

where $Y$ denotes the set of output neurons, the subscript $i$ refers to the
index of a specific output neuron and $\rho(\hat{y}_{i})$ represents the
target value sought for their activation function. For an Ising system, this
corresponds to minimizing the spin energy

$E_{Ising}+\beta
C(\sigma^{y},\hat{\sigma}^{y})=\sum_{i>j}J_{ij}\sigma_{i}\sigma_{j}+\sum_{i}h_{i}\sigma_{i}+\beta\sum_{i\in
Y}(\sigma_{i}-\hat{\sigma}_{i})^{2},$ (5)

where the output of the activation function is the binary value of the state
$\sigma$. Since the output spin $\sigma^{y}$ and their corresponding target
state $\hat{\sigma}^{y}$ always take a value equal to $\pm 1$ , this formula
can be rewritten as (see Supplementary Note LABEL:sec:suppl-nudge-ising):

$E_{Ising}+\beta
C(\sigma^{y},\hat{\sigma}^{y})=\sum_{i>j}J_{ij}\sigma_{i}\sigma_{j}+\sum_{i\not\in
Y}h_{i}\sigma_{i}+\sum_{i\in Y}(h_{i}-\beta\hat{\sigma}_{i})\sigma_{i}.$ (6)

This encodes the nudge as an additional bias term $\beta\hat{\sigma}_{i}$ that
is only applied to the output spins.

At the end of the free phase, the system is frozen in its ground state 1b
($time=20~{}\mu$s). We found that the application of biases at the output
alone is not sufficient to destabilize it, which prevents the error from being
backpropagated. To overcome this problem, we employed the Reverse Quantum
Annealing procedure PerdomoOrtiz2010 . As depicted in Fig. 1b ($20~{}\mu$s
$<time<40~{}\mu$s), we slightly increase the interstate tunneling probability,
allowing the system to evolve to a state that is close to the equilibrium
state of the free phase but closer to the desired output state. Finally, we
decrease the tunneling probability to 0 again, ensuring that the system
reaches the new nudged steady state, as illustrated in Fig. 1b ($40~{}\mu$s
$<time<60~{}\mu$s).

While reverse annealing allows for the nudge phase, EP requires additional
adaptation to effectively train systems that have abrupt ON/OFF flip-flop
activation functions, such as binary neurons or Ising spins. For example, for
the same range of input values, a flip-flop neuron is statistically much less
likely to change state than a neuron with a continuous activation function
when the same nudge bias is applied to them. As illustrated in Fig. 1c, for
small nudge biases applied to output spins, the network does not change its
state in practice and does not learn. Applying too high nudge biases also
poses a problem, as the switching of the output neurons between their two
extreme values leads to very strong changes in the states of the neurons of
the whole network, hindering the learning process which needs to be
progressive. To solve this problem, Laydevant et al.Laydevant_2021_CVPR
proposed using several neurons to represent each output class instead of using
only one as is commonly done. We have adopted this method, which allows us to
induce more flip-flop events in the output and back-propagate them more easily
to the rest of the network.

The steady states of the spins at the end of the free and nudge phase
(respectively $\sigma^{*,0}$ and $\sigma^{*,\beta}$) are measured, recorded,
and used to calculate the gradient of the loss function with respect to the
couplings according to the following learning rule:

$-\frac{\partial C}{\partial J_{ij}}=\Delta
J_{ij}\propto-\frac{1}{\beta}\left[(\sigma_{i}\sigma_{j})^{*,nudge}-(\sigma_{i}\sigma_{j})^{*,free}\right]$
(7)

for a fully-connected architecture. The updates are then applied to the
weights using the standard stochastic gradient descent algorithm.

The number of layers, of neurons per layer, and the connectivity of a neural
network depend on the task to be solved. Additionally, in our case, this
architecture is constrained by the specific connectivity of the hardware
system. The D-Wave machine’s spins are organized in locally connected sub-
networks, as illustrated in n Fig. 1d. As a result, mapping a fully-connected
architecture onto the chip is not straightforward. We employ the embedding
procedure provided by D-Wave to map the neural network architecture to the
chip’s architecture. This procedure relies heavily on the ”chaining” process
of spins, allowing a spin to couple with more spins than its direct six
neighbors. The chaining process involves strongly coupling a chain of physical
spins on the chip so that they maintain the same value at each time step of
the annealing procedure, effectively implementing a single spin. Through the
embedding procedure, we have successfully mapped our fully-connected
architecture at the expense of using chains of approximately six physical
spins per neuron.

###  2.2 Training a Fully-connected neural network

We now apply the previously described methods to train a neural network
featuring a fully-connected architecture for recognizing handwritten digits
from the MNIST database mnist , a widely used benchmark for evaluating the
performance of hardware-based neural networks.

Fully-connected neural networks for solving the MNIST problem typically
consist of several layers of neurons, including an input layer with 784
neurons (one neuron per image pixel), one or more hidden layers, and an output
layer. Here we chose to embed the largest possible neural network with one
hidden layer as described below. Since Equilibrium Propagation-based training
relies on the dynamics of the neurons within the network and input neurons
have fixed values, we do not implement them on the chip (see Eq. (8)).
Instead, we calculate the product of the input data $X$ (an input image, for
instance) and a trainable weight matrix $W_{input}$ using a digital computer.
The resulting vector of fixed bias fields $h_{input}$ encodes the input on the
chip:

$h_{input}=X*W_{input}.$ (8)

This bias is then augmented by a second bias that is the equivalent of the
standard bias of artificial neural networks: $h_{hidden}=h_{input}+h_{bias}$.

We find that using 4 output neurons per class of digits to be recognized (from
0 to 9, resulting in 40 output neurons) allows the nudge phase to function
effectively. When mapping this fully-connected network architecture onto the
5000 locally connected spins of D-Wave, we determined that the maximum number
of neurons we can implement in the hidden layer is 120 (Fig. 2a). Due to the
limited access time to the D-wave machine, we train this network with only a
part of the MNIST data - using 1000 images for training and 100 images for
testing (see Methods). We refer to this task as MNIST/100, following the
notation of truncated-mnist . As shown in Fig. 2b, we obtained a recognition
rate of 98.8% ($\pm 0.8$) on the training data and 88.8% ($\pm 1.5$) on the
test data.

Figure 2: Training a fully-connected neural network on MNIST. a. Embedding the
fully-connected architecture on the Ising machine. We first compute in
software the product between the input vector (an MNIST image) and the first
weight matrix. The result is a vector of small constant bias fields that are
directly applied to the hidden spins on the chip. The hidden and the output
layer are embedded, coupled and eventually nudged on the actual chip b.
Training on D-Wave: training and testing accuracy as a function of the number
of epochs. c. Training and testing accuracy as a function of the number of
hidden neurons. We plot the accuracy obtained with the D-Wave Ising machine
vs. those obtained with Simulated Annealing (SA) and with the deterministic
Artificial Neural Network based on binary activations and real-value weights
(ANN-EP). d. Training and testing accuracy as a function of the size of the
training dataset. We plot the accuracy obtained on the D-Wave Ising machine
vs. that obtained with Simulated annealing (SA).

We benchmark the physical system in Fig. 2c-d by comparing the obtained
accuracies to numerical simulations carried out on a digital processor for two
identically designed networks. The first network (dashed lines in Fig. 2c-d)
is a spin network identical to the one on the chip, and trained in the same
way by replacing the quantum annealing with Simulated Annealing (SA). The
second network (solid lines in Fig. 2c) is a software Artificial Neural
Network with binary activations and real-value weights (ANN-EP) evolving
according to a Hopfield energy and trained by Equilibrium Propagation. As
shown in Fig. 2c-d, the accuracy that the hardware spin network reaches on
MNIST/100 (denoted by the dots with error bars on 3 repetitions) is always
higher or equal to that of the simulated ideal networks for both the training
and test databases, demonstrating the quality of the training performed on the
Ising machine.

The MNIST/100 task is the most complex task that the D-Wave machine has been
trained to solve to date. Previous results, based on unsupervised contrastive
divergence learning, have been limited to smaller subsets of the MNIST
dataset, such as MNIST/20 dorband_boltzmann_2015 , or images reduced to (6x6)
pixels instead of the 28x28 pixels of the original database
adachi_application_nodate ; adachi . The simulations in Fig. 2d indicate that
the recognition rate on the test data of the implemented 784-120-40 network
can be further improved to 97% ($\pm 1\%$) when trained with more images, in
this case, 10000.

The simulations with simulated annealing in Fig.2c also reveal that the
recognition rate on the reduced MNIST/100 database can be further improved to
94% by increasing the number of hidden layer neurons to 1920. While the total
number of neurons required for this network in hardware (1960 in total for the
two layers, hidden and output) is lower than the total number of spins
available on D-Wave, the sparse connectivity of the chip and the need for an
embedding step make it impossible to implement in practice. This is because it
requires the use of several spins per neuron, $\sim 6$ per neuron for the
implemented network (Fig. 2a).

Therefore, it is essential to consider neural network architectures that are
congruent with the local connectivity of most Ising machines to make the best
use of their resources and, in particular, to use fewer spins per neuron in
order to embed larger, and thus more powerful, architectures. In the following
section, we demonstrate that it is possible to map and train a complete
convolutional neural network on the Chimera graph structure proposed by D-Wave
(Fig. 1d), with less than 1.6 spins per neuron on average.

###  2.3 Training a Convolutional neural network

Figure 3: Training a convolutional neural network. a. Detailed view of the
convolutional neural network trained on the D-Wave Ising machine. The inputs
are 3x3 pixel images. The convolution layer applies four different sets of 2x2
weight filters to the input images, generating four 2x2 feature maps. Each
feature map is then condensed to a single value through an average pooling
operation. After flattening, a fully-connected classifier provides the output
of the network, a vector of dimension 4. b. Schematic of the convolutional
neural network’s implementation on the Chimera architecture of a D-Wave Ising
machine. The four peripheral crossbar arrays perform the convolution operation
in parallel. Each array receives a distinct patch of pixels from the input
data, yet all share identical couplings, ensuring uniform filter application
across different patches. The blue spins on the four crossbar arrays
individually represent the values $x_{1},x_{2},x_{3},x_{4}$ of the feature map
highlighted in (a), obtained by convolving the input with the filter encoded
in the couplings depicted in light blue circles. The output of the
convolutional operation is then down-sampled via the averaged pooling coupling
($J=\frac{1}{4}$) and linked through identity couplings ($J=1$). The dotted
blue chain represents the output of the averaged pooling operation applied to
the feature map depicted with the blue spins. Finally, the results of the
averaged pooling operation are fed into the fully-connected classifier, which
predicts the input’s class. Here we have four output neurons as we use two
output neurons to encode a class. The convolutional neural network is mapped
as-is on the chip, eliminating the need for an embedding step. c. The training
dataset consisting in the 2 patterns used for training the CNN implemented on
the D-Wave Ising machine. d. Training curve (mean squared error) related to
training the CNN on the D-Wave Ising machine. e. Training curve (accuracy (%)
related to training the CNN on the D-Wave Ising machine.

In this section we show that we can directly map a convolutional neural
network (CNN) to the Chimera connectivity graph of the D-Wave Ising machine
employed in our study.

Convolutional neural networks, which are currently one of the state-of-the-art
architectures for image classification, function by sliding convolutional
filters over input data (Conv2D in Fig. 3a) in order to extract different
learnable features. The resulting feature maps are then down-sampled (Pooling
operation in Fig. 3a), combined and fed into a final fully-connected
classifier (Flattening and Fully-connected in Fig. 3a) to assign a class to
the input data.

Typically, the convolutional operation is sequential since the filters must
move over the input data to compute multiple local dot products (Fig. 3a).
However, with the D-Wave Ising machine’s Chimera connectivity graph, we can
perform this operation in a completely parallel manner by utilizing the local
clusters consisting of 4 spins connected to 4 other spins by $4\times 4$
adjustable couplings represented in Fig. 1d).

The convolutional layer applies a filter, a set of four weights, to extract
2x2 pixel patches of the image (e.g. $p_{1},p_{2},p_{4},p_{5}$) and apply a
non-linearity (the binary value of the spin), generating a single output value
for each patch (e.g., $x_{1}$). This operation creates four output values for
a given filter, constituting a feature map.

We employ four different crossbars, labelled $W_{c}onv$ in Fig. 3b to process
the four different 2x2 pixel patches of the input. In the schematic crossbar
view, the spins are physically represented by horizontal and vertical lines.
The pixel values $p_{i}$ are applied to the input spins through strong biases
$h_{i}$ that set their direction. The output spins in the crossbar then align
according to the sum of the inputs weighted by the coupling values. The set of
couplings highlighted in blue across the crossbars correspond to the same
filter, generating the values $x_{1},x_{2},x_{3},x_{4}$ in the blue output
spin of each crossbar. This method allows us to apply four different filters
to the full input image simultaneously, thereby executing the entire
convolutional operation in parallel.

After the convolution operation, we apply a pooling operation to reduce the
dimensionality of the output. Standard CNNs typically employ max pooling,
which would, in our example, correspond to keeping the maximum of the four
outputs in blue, and coupling only this one to the next layers in the network.
We found more convenient and relevant here, given the layout of the chip and
the potential degeneracy in the maximum values of a set of binary spins, to
implement an average pooling, which averages out the four outputs before
connecting them to the next layers. To achieve this, we connect the outputs of
the different convolution operations (e.g. the blue ”spins” in Fig. 3c) to
another chain of spins (e.g. the blue horizontal dotted line in Fig. 3c)
through couplings with the value $J=1/4$, highlighted in yellow in Fig. 3c.
This arrangement enables the computation of the weighted average of the
convolution outputs and concurrently performs the ”flatten” operation of CNNs.

The final stage of the neural network, subsequent to the pooling and
flattening operations, consists of a fully-connected layer with four outputs.
his stage is realized using the central crossbar array depicted in Fig. 3c,
where the array’s couplings implement the layer’s weights.

Here, we use two spins/per class to classify the $3\times 3$ pixel thumbnails
depicted in Fig. 3d. As shown in Fig. 3e-f, the network implemented and
trained on the D-Wave’s Ising machine achieves a 100% success rate on the
training patterns. This result demonstrates the feasibility of training a
convolutional network by performing the two phases as well as gradient
computations directly on a quantum computer.,This method utilizes connectivity
far more efficiently than a fully-connected network, requiring on average only
$1.6$ spins per neuron as opposed to approximately $6$. These results also
demonstrate the power and flexibility of Equilibrium Propagation to train
hardware systems with constrained connectivity.

##  3 Discussion

In the past, several attempts have been made to implement neural networks on
spin systems by using their physics and taking into account hardware
constraints. Boltzmann machines, in particular, can take advantage of the
annealing procedures available in Ising machines. However, the size of the
networks that can be embedded on the chip is limited in practice because in
Boltzmann machines, all the input neurons (784 for MNIST) must be physically
present on the chip. Most implementations up to date were made with D-Wave,
which offers a few thousand spins. References adachi_application_nodate ;
Benedetti2017 ; dorband_boltzmann_2015 ; liu_adiabatic_2018 ; adachi  used
D-Wave to train a Restricted Boltzmann Machine (RBM) on a coarse-grained
version of MNIST downsampled to $6\times 6$ pixel images in order to fit the
chip that was available at the time of publication. They then fine-tuned the
network with backpropagation outside of the chip. Reference
dixit_training_2021  only trains some couplings between the hidden nodes (80
nodes) of the Boltzmann Machine. The couplings between the visible and the
hidden layer are computed on a side computer. Only dorband_boltzmann_2015
trains the Ising Machine on the standard version of MNIST but restricted to
200 training images. However, it is still trained layer-wise and shows limited
performances on MNIST/20 (maximum 67% test accuracy with 479 sparsely
connected hidden nodes).

In a recent study niazi2023training , researchers trained a 2-hidden layer
Deep Boltzmann machine on an FPGA, featuring connectivity similar to the one
used for the D-Wave Ising machine. The Boltzmann machine was kept sparse,
avoiding the use of embedding for a fully-connected architecture, which is
reminiscent of our approach to leverage the Ising machine’s connectivity for
convolution operations. As we do with Equilibrium Propagation, the authors
trained the network as a whole instead of layer-by-layer, which is typically
done for Boltzmann machines. They also employed multiple neurons per encoded
class (five spins per class, compared to our four spins per class), and showed
that their approach converges for systems composed of two-state flip-flop
systems such as Ising machines. Nonetheless, the test accuracy achieved on the
full MNIST dataset (90% for both train and test accuracy) underscores the
advantage of using Equilibrium Propagation and reverse annealing to train
Ising systems. As we showed in the simulations of Fig. 2d, a test accuracy of
97% is theoretically attainable by training the D-Wave Ising machine with
Equilibrium propagation on MNIST/1000, a performance that is likely to improve
when employing the full database.

The D-Wave Ising machine has also been employed to train specific components
of other neural network types. For instance, nguyen_generating_2017 ;
nguyen_image_2018 ; sleeman-dwave  leverage the probabilistic nature of this
hardware to generate a sparse latent representation of an auto-encoder.
However, again, the authors do not train the entire auto-encoder on the D-Wave
Ising machine.

In this work, we demonstrate the feasibility of performing inference,
backpropagation of errors, and computation of gradients of a global cost
function solely through the dynamics of a spin system. Our approach paves the
way for training Ising machines using modern supervised learning algorithms,
employing standard gradient-based methods.

Our experiments were carried out using the D-Wave machine, whose quantum
properties are a topic of ongoing discussion in the scientific community
experimental-quantum-annealing ; Rnnow2014 . We found that the accuracy
obtained on the hardware through quantum annealing was slightly better than
software simulations, particularly simulated annealing. However, it is
challenging to conclusively determine a quantum advantage as the exact
training conditions of the hardware cannot be easily replicated in simulation
(see Methods). To establish more definitive conclusions on this matter, it
would be valuable to investigate whether the D-Wave machine maintains higher
accuracy compared to software simulations when trained on other tasks
(especially those necessitating an optimal ground state for the Ising system),
and ideally, to compare the hardware accuracy with that of classical Ising
machines with the same connectivity.

Future generations of D-Wave Ising machines will be capable of modeling larger
spin systems, allowing to embed more complex neural networks. In particular,
they will be equipped with larger locally-connected crossbar arrays zephyr-
archi , enabling the direct application of $3\times 3$ filters and the
performance of convolutions on benchmark images, such as those from MNIST or
CIFAR-10. One constraint of our current convolution implementation for
progressing in this direction is the requirement for binary inputs. However,
there are at least two possible ways to tackle this issue in the future when
implementing convolutional neural networks on a D-Wave Ising machine. The
first approach is similar to the method we employed in this article for
training the fully-connected architecture, where we compute the initial
vector-matrix product on a digital computer and send the results as inputs to
the chip. We could use the same technique with a first convolutional layer
computed on a digital computer with real-value inputs, as is commonly done
with binary neural networks courbariaux_binarized_2016 ; rastegari_xnor-
net_2016 . A second way to manage binary inputs directly on the chip would be
to take inspiration from hirtzlin_stochastic_2019 , where the authors
stochastically binarize the inputs to accommodate binary inputs, even though
the original inputs are real-valued (CIFAR-10).

Our results and the algorithm employed to obtain them can be applied to any
type of Ising machine. Those with an ultra-low power consumption are
particularly appealing for reducing the overall electrical consumption of AI
and deploying it in embedded systems. Memristors or spintronic nano-components
are currently being extensively researched as building blocks for such
systems, as they enable the co-integration of memory, novel physical
functionality and computing. This greatly enhances the efficiency and
scalability of the system, making it more suitable for real-world
applications. In line with this trend, a recent study Yi2022  successfully
trained a crossbar array of memristors to emulate the couplings of an Ising-
like model using a learning law similar to Equilibrium Propagation. Contrary
to our study, the system in this research is not intrinsically dynamic.
Instead, the dynamics are emulated by iteratively and recursively using the
memristive system. Combining our approach, in which an Ising spin system
computes gradients through its intrinsic dynamics, with the use of memristive
couplings presents significant potential for future low-power, embedded Ising
machines that can be trained on the edge.

Figure 4: Training the stochastic Ising machine to behave in a deterministic
way. Distribution of the steady-state energy of the spin system modeled by the
D-Wave Ising machine when the exact same inputs are applied 1000 times to the
same neural network architecture. The light pink distribution corresponds to
the energy of the system before training it on MNIST/100. The purple
distribution corresponds to the energy of the system after training it with
EP. The samples before training have a flat distribution with a large standard
deviation - which is the signature of a highly stochastic system. The samples
before training exhibit a flat distribution with a substantial standard
deviation, indicative of a highly stochastic system. In contrast, the samples
after training are concentrated within a narrower region with a considerably
smaller standard deviation. This demonstrates that the training procedure
causes the Ising machine to exhibit more deterministic behavior.

Finally, Ising machines, designed to reach the ground state of an Ising
system, are inherently stochastic in nature hamerly_experimental_2019 . D-Wave
for instance operates at finite non-zero temperature, resulting in thermal
excitation competing with quantum annealing. When the machine is used for
solving combinatorial problems, the state needs to be sampled multiple times
to obtain an accurate solution. We also had to sample the state of the machine
10 times at the end of each phase in order to achieve successful trainings.
However, we found that after training with our method, the initially
stochastic D-Wave Ising machine is brought into a much more deterministic
regime (as seen in Fig. 4). This means that after training, the solution is
given in a single call to the Ising machine, greatly reducing the inference
time. These results also open up the possibility of learning combinatorial
problems through a data-driven approach that can provide faster and more
accurate solutions than the traditional approach where Ising system parameters
are defined by the problem.

##  4 Conclusion

In conclusion, this study presents a significant advancement in the field of
utilizing Ising machines as hardware platforms for Artificial Intelligence.
Leveraging Equilibrium Propagation together with annealing methods, we have
successfully demonstrated that Ising machines can be trained using modern
supervised learning algorithms, overcoming the limitations of previous
attempts to implement neural networks on spin systems. Our experiments,
conducted on the D-Wave Ising machine, show that the accuracy obtained through
quantum annealing is on par with that of software simulations. Additionally,
our results indicate that the local connectivity of Ising machines can be
effectively harnessed for performing convolutions. The potential for future
developments by combining our approach, where Ising spin systems compute
gradients through their intrinsic dynamics, with low-power hardware that
employs nanotechnologies such as memristors for implementing local couplings,
presents exciting opportunities for the future of embedded AI.

Acknowledgments This work was supported by the European Research Council
advanced grant GrenaDyn (reference: 101020684). The text of the article was
partially edited by a large language model (OpenAI ChatGPT). The authors would
like to thank D. Querlioz for discussion and invaluable feedback.

Author contributions statement J.G., J.L. and D.M. devised the study. J.L.
performed all the simulations and experiments. J.G and J.L. wrote the initial
version of the manuscript. All authors discussed the results and reviewed the
manuscript.

Competing interests The authors declare no competing interests.

Data availability The datasets analyzed and all data measured in this study
are available from the corresponding author upon reasonable request.

Code availability The code to reproduce the results is available on github at
the following link: https://github.com/jlaydevant/Ising-Machine-EqProp.

Supplementary information

We provide additional information in the Supplementary Materials file.

## References

  *   * (1) Marković, D., Mizrahi, A., Querlioz, D., Grollier, J.: Physics for neuromorphic computing. Nature Reviews Physics 2(9), 499–510 (2020). https://doi.org/10.1038/s42254-020-0208-2
  * (2) Ising, E.: Beitrag zur Theorie des Ferromagnetismus. Zeitschrift für Physik 31(1), 253–258 (1925). https://doi.org/10.1007/BF02980577. Accessed 2022-08-24
  * (3) Harris, R., Johnson, M.W., Lanting, T., Berkley, A.J., Johansson, J., Bunyk, P., Tolkacheva, E., Ladizinsky, E., Ladizinsky, N., Oh, T., Cioata, F., Perminov, I., Spear, P., Enderud, C., Rich, C., Uchaikin, S., Thom, M.C., Chapple, E.M., Wang, J., Wilson, B., Amin, M.H.S., Dickson, N., Karimi, K., Macready, B., Truncik, C.J.S., Rose, G.: Experimental Investigation of an Eight Qubit Unit Cell in a Superconducting Optimization Processor. Physical Review B 82(2), 024511 (2010). https://doi.org/10.1103/PhysRevB.82.024511. arXiv: 1004.1628. Accessed 2022-02-04
  * (4) Little, W.A.: The existence of persistent states in the brain. Mathematical Biosciences 19(1), 101–120 (1974). https://doi.org/10.1016/0025-5564(74)90031-5. Accessed 2022-07-15
  * (5) Hopfield, J.J.: Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences 79(8), 2554–2558 (1982). https://doi.org/10.1073/pnas.79.8.2554
  * (6) Amit, D.J., Gutfreund, H., Sompolinsky, H.: Spin-glass models of neural networks. Physical Review A 32(2), 1007–1018 (1985). https://doi.org/10.1103/physreva.32.1007
  * (7) Mezard, M., Parisi, G., Virasoro, M.: Spin Glass Theory and Beyond. WORLD SCIENTIFIC, ??? (1986). https://doi.org/10.1142/0271. https://doi.org/10.1142/0271
  * (8) Adachi, S.H., Henderson, M.P.: Application of Quantum Annealing to Training of Deep Neural Networks (2015)
  * (9) Benedetti, M., Realpe-Gómez, J., Biswas, R., Perdomo-Ortiz, A.: Quantum-assisted learning of hardware-embedded probabilistic graphical models. Physical Review X 7(4) (2017). https://doi.org/10.1103/physrevx.7.041052
  * (10) Dorband, J.E.: A Boltzmann Machine Implementation for the D-Wave. In: 2015 12th International Conference on Information Technology - New Generations, pp. 703–707. IEEE, Las Vegas, NV, USA (2015). https://doi.org/10.1109/ITNG.2015.118. http://ieeexplore.ieee.org/document/7113557/ Accessed 2022-02-03
  * (11) Liu, J., Spedalieri, F.M., Yao, K.-T., Potok, T.E., Schuman, C., Young, S., Patton, R., Rose, G.S., Chamka, G.: Adiabatic Quantum Computation Applied to Deep Learning Networks. Entropy 20(5), 380 (2018). https://doi.org/10.3390/e20050380. Number: 5 Publisher: Multidisciplinary Digital Publishing Institute. Accessed 2022-02-03
  * (12) Job, J., Adachi, S.: Systematic comparison of deep belief network training using quantum annealing vs. classical techniques (2020)
  * (13) Dixit, V., Selvarajan, R., Alam, M.A., Humble, T.S., Kais, S.: Training Restricted Boltzmann Machines With a D-Wave Quantum Annealer. Frontiers in Physics 9 (2021). Accessed 2022-02-03
  * (14) Niazi, S., Aadit, N.A., Mohseni, M., Chowdhury, S., Qin, Y., Camsari, K.Y.: Training Deep Boltzmann Networks with Sparse Ising Machines (2023)
  * (15) Hinton, G.E., Sejnowski, T.J.: Optimal Perceptual Inference, 6
  * (16) Hinton, G.E.: Training Products of Experts by Minimizing Contrastive Divergence. Neural Computation 14(8), 1771–1800 (2002). https://doi.org/10.1162/089976602760128018. Accessed 2022-09-01
  * (17) Krizhevsky, A.: Convolutional Deep Belief Networks on CIFAR-10, 9
  * (18) Wang, Z., Wu, H., Burr, G.W., Hwang, C.S., Wang, K.L., Xia, Q., Yang, J.J.: Resistive switching materials for information processing. Nature Reviews Materials 5(3), 173–195 (2020). https://doi.org/10.1038/s41578-019-0159-3
  * (19) Xia, Q., Yang, J.J.: Memristive crossbar arrays for brain-inspired computing. Nature Materials 18(4), 309–323 (2019). https://doi.org/10.1038/s41563-019-0291-x. Accessed 2021-04-22
  * (20) Markovifá, D., Mizrahi, A., Querlioz, D., Grollier, J.: Physics for neuromorphic computing. Nature Reviews Physics, 1–12 (2020). https://doi.org/10.1038/s42254-020-0208-2. Accessed 2020-08-19
  * (21) Nøkland, A.: Direct Feedback Alignment Provides Learning in Deep Neural Networks. arXiv. arXiv:1609.01596 [cs, stat] (2016). https://doi.org/%****␣main.bbl␣Line␣375␣****10.48550/arXiv.1609.01596. http://arxiv.org/abs/1609.01596 Accessed 2022-08-23
  * (22) Neftci, E.O., Mostafa, H., Zenke, F.: Surrogate Gradient Learning in Spiking Neural Networks: Bringing the Power of Gradient-Based Optimization to Spiking Neural Networks. IEEE Signal Processing Magazine 36(6), 51–63 (2019). https://doi.org/10.1109/MSP.2019.2931595. Accessed 2022-08-24
  * (23) Martin, E., Ernoult, M., Laydevant, J., Li, S., Querlioz, D., Petrisor, T., Grollier, J.: EqSpike: Spike-driven Equilibrium Propagation for Neuromorphic Implementations. arXiv:2010.07859 [cs] (2020). arXiv: 2010.07859. Accessed 2020-11-25
  * (24) Kendall, J., Pantone, R., Manickavasagam, K., Bengio, Y., Scellier, B.: Training End-to-End Analog Neural Networks with Equilibrium Propagation. arXiv:2006.01981 [cs] (2020). arXiv: 2006.01981. Accessed 2021-04-09
  * (25) Frenkel, C., Lefebvre, M., Bol, D.: Learning without feedback: Fixed random learning signals allow for feedforward training of deep neural networks. Frontiers in Neuroscience 15 (2021). https://doi.org/10.3389/fnins.2021.629892
  * (26) Ernoult, M., Normandin, F., Moudgil, A., Spinney, S., Belilovsky, E., Rish, I., Richards, B., Bengio, Y.: Towards Scaling Difference Target Propagation by Learning Backprop Targets. arXiv:2201.13415 [cs] (2022). arXiv: 2201.13415. Accessed 2022-02-07
  * (27) Wright, L.G., Onodera, T., Stein, M.M., Wang, T., Schachter, D.T., Hu, Z., McMahon, P.L.: Deep physical neural networks trained with backpropagation. Nature 601(7894), 549–555 (2022). https://doi.org/10.1038/s41586-021-04223-6. Number: 7894 Publisher: Nature Publishing Group. Accessed 2022-08-30
  * (28) Schuman, C.D., Kulkarni, S.R., Parsa, M., Mitchell, J.P., Date, P., Kay, B.: Opportunities for neuromorphic computing algorithms and applications. Nature Computational Science 2(1), 10–19 (2022). https://doi.org/10.1038/s43588-021-00184-y
  * (29) Scellier, B., Bengio, Y.: Equilibrium Propagation: Bridging the Gap between Energy-Based Models and Backpropagation. Frontiers in Computational Neuroscience 11 (2017). Accessed 2022-06-21
  * (30) Dillavou, S., Stern, M., Liu, A.J., Durian, D.J.: Demonstration of decentralized physics-driven learning. Phys. Rev. Appl. 18, 014040 (2022). https://doi.org/10.1103/PhysRevApplied.18.014040
  * (31) Yi, S.-i., Kendall, J.D., Williams, R.S., Kumar, S.: Activity-difference training of deep neural networks using memristor crossbars. Nature Electronics (2022). https://doi.org/%****␣main.bbl␣Line␣550␣****10.1038/s41928-022-00869-w
  * (32) Ernoult, M., Grollier, J., Querlioz, D., Bengio, Y., Scellier, B.: Updates of equilibrium prop match gradients of backprop through time in an rnn with static input. In: Wallach, H., Larochelle, H., Beygelzimer, A., d' Alché-Buc, F., Fox, E., Garnett, R. (eds.) Advances in Neural Information Processing Systems, vol. 32. Curran Associates, Inc., ??? (2019). https://proceedings.neurips.cc/paper/2019/file/67974233917cea0e42a49a2fb7eb4cf4-Paper.pdf
  * (33) Laborieux, A., Ernoult, M., Scellier, B., Bengio, Y., Grollier, J., Querlioz, D.: Scaling Equilibrium Propagation to Deep ConvNets by Drastically Reducing its Gradient Estimator Bias. arXiv:2006.03824 [cs] (2020). arXiv: 2006.03824. Accessed 2020-11-24
  * (34) Laborieux, A., Zenke, F.: Holomorphic Equilibrium Propagation Computes Exact Gradients Through Finite Size Oscillations. arXiv (2022). https://doi.org/10.48550/ARXIV.2209.00530. https://arxiv.org/abs/2209.00530
  * (35) Mohseni, N., McMahon, P.L., Byrnes, T.: Ising machines as hardware solvers of combinatorial optimization problems. Nature Reviews Physics 4(6), 363–379 (2022)
  * (36) Lucas, A.: Ising formulations of many NP problems. Frontiers in Physics 2 (2014). Accessed 2022-06-24
  * (37) Cai, F., Kumar, S., Van Vaerenbergh, T., Sheng, X., Liu, R., Li, C., Liu, Z., Foltin, M., Yu, S., Xia, Q., Yang, J.J., Beausoleil, R., Lu, W.D., Strachan, J.P.: Power-efficient combinatorial optimization using intrinsic noise in memristor Hopfield neural networks. Nature Electronics 3(7), 409–418 (2020). https://doi.org/10.1038/s41928-020-0436-6. Number: 7 Publisher: Nature Publishing Group. Accessed 2022-07-11
  * (38) Farhi, E., Goldstone, J., Gutmann, S., Lapan, J., Lundgren, A., Preda, D.: A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem. Science 292(5516), 472–475 (2001). https://doi.org/10.1126/science.1057726
  * (39) Yamamoto, Y., Aihara, K., Leleu, T., Kawarabayashi, K.-i., Kako, S., Fejer, M., Inoue, K., Takesue, H.: Coherent Ising machines—optical neural networks operating at the quantum limit. npj Quantum Information 3(1), 49 (2017). https://doi.org/10.1038/s41534-017-0048-9. Accessed 2022-03-14
  * (40) Nielsen, M.: Reduced MNIST: how well can machines learn from small data? https://cognitivemedium.com/rmnist. Accessed: 2022-09-03
  * (41) Laydevant, J., Ernoult, M., Querlioz, D., Grollier, J.: Training dynamical binary neural networks with equilibrium propagation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, pp. 4640–4649 (2021)
  * (42) Kirkpatrick, S., Gelatt, C.D., Vecchi, M.P.: Optimization by Simulated Annealing. Science 220(4598), 671–680 (1983). https://doi.org/10.1126/science.220.4598.671. Publisher: American Association for the Advancement of Science. Accessed 2022-08-24
  * (43) Perdomo-Ortiz, A., Venegas-Andraca, S.E., Aspuru-Guzik, A.: A study of heuristic guesses for adiabatic quantum computation. Quantum Information Processing 10(1), 33–52 (2010). https://doi.org/10.1007/s11128-010-0168-z
  * (44) LeCun, Y., Cortes, C.: MNIST handwritten digit database (2010)
  * (45) Nguyen, N.T.T., Larson, A.E., Kenyon, G.T.: Generating Sparse Representations Using Quantum Annealing: Comparison to Classical Algorithms. In: 2017 IEEE International Conference on Rebooting Computing (ICRC), pp. 1–6. IEEE, Washington, DC (2017). https://doi.org/10.1109/ICRC.2017.8123653. http://ieeexplore.ieee.org/document/8123653/ Accessed 2022-03-14
  * (46) Nguyen, N.T.T., Kenyon, G.T.: Image classification using quantum inference on the D-Wave 2X. 2018 IEEE International Conference on Rebooting Computing (ICRC), 1–7 (2018). https://doi.org/10.1109/ICRC.2018.8638596. arXiv: 1905.13215. Accessed 2022-02-03
  * (47) Sleeman, J., Dorband, J., Halem, M.: A hybrid quantum enabled RBM advantage: convolutional autoencoders for quantum image compression and generative learning. In: Donkor, E., Hayduk, M. (eds.) Quantum Information Science, Sensing, and Computation XII, vol. 11391, p. 113910. SPIE, ??? (2020). https://doi.org/%****␣main.bbl␣Line␣825␣****10.1117/12.2558832. International Society for Optics and Photonics. https://doi.org/10.1117/12.2558832
  * (48) Boixo, S., Albash, T., Spedalieri, F.M., Chancellor, N., Lidar, D.A.: Experimental signature of programmable quantum annealing. Nature Communications 4(1), 2067 (2013). https://doi.org/10.1038/ncomms3067
  * (49) Rønnow, T.F., Wang, Z., Job, J., Boixo, S., Isakov, S.V., Wecker, D., Martinis, J.M., Lidar, D.A., Troyer, M.: Defining and detecting quantum speedup. Science 345(6195), 420–424 (2014). https://doi.org/10.1126/science.1252319
  * (50) Zephyr Topology of D-Wave Quantum Processors. https://www.dwavesys.com/media/2uznec4s/14-1056a-a_zephyr_topology_of_d-wave_quantum_processors.pdf. Accessed: 2023-02-02
  * (51) Courbariaux, M., Hubara, I., Soudry, D., El-Yaniv, R., Bengio, Y.: Binarized Neural Networks: Training Deep Neural Networks with Weights and Activations Constrained to +1 or -1. arXiv:1602.02830 [cs] (2016). arXiv: 1602.02830. Accessed 2021-03-23
  * (52) Rastegari, M., Ordonez, V., Redmon, J., Farhadi, A.: XNOR-Net: ImageNet Classification Using Binary Convolutional Neural Networks. In: Leibe, B., Matas, J., Sebe, N., Welling, M. (eds.) Computer Vision – ECCV 2016. Lecture Notes in Computer Science, pp. 525–542. Springer, Cham (2016). https://doi.org/10.1007/978-3-319-46493-0_32
  * (53) Hirtzlin, T., Penkovsky, B., Bocquet, M., Klein, J.-O., Portal, J.-M., Querlioz, D.: Stochastic Computing for Hardware Implementation of Binarized Neural Networks. IEEE Access 7, 76394–76403 (2019). https://doi.org/10.1109/ACCESS.2019.2921104. Conference Name: IEEE Access
  * (54) Hamerly, R., Inagaki, T., McMahon, P.L., Venturelli, D., Marandi, A., Onodera, T., Ng, E., Langrock, C., Inaba, K., Honjo, T., Enbutsu, K., Umeki, T., Kasahara, R., Utsunomiya, S., Kako, S., Kawarabayashi, K.-i., Byer, R.L., Fejer, M.M., Mabuchi, H., Englund, D., Rieffel, E., Takesue, H., Yamamoto, Y.: Experimental investigation of performance differences between coherent Ising machines and a quantum annealer. Science Advances 5(5), 0823 (2019). https://doi.org/10.1126/sciadv.aau0823. Accessed 2022-05-19
  * (55) Ocean API - D-Wave. https://docs.ocean.dwavesys.com/en/stable/. Accessed: 2023-02-02
  * (56) Park, J., Lee, J., Jeon, D.: A 65-nm neuromorphic image classification processor with energy-efficient training through direct spike-only feedback. IEEE Journal of Solid-State Circuits 55(1), 108–119 (2020). https://doi.org/10.1109/JSSC.2019.2942367

##  5 Methods

In this section, we outline the essential steps of the methods employed to
generate the results presented in the paper. Further details can be found in
the Supplementary Materials file. We begin by discussing how to interact with
the D-Wave Ising machine and identifying the most crucial features to
facilitate training with Equilibrium Propagation (EP) on this specific Ising
machine. Following this, we provide a brief overview of the requirements and
learning rules of EP. Lastly, we delve into the training methods for both the
fully-connected and convolutional architectures.

Additionally, we have made our code available at at the following link:
https://github.com/jlaydevant/Ising-Machine-EqProp to facilitate the
reproduction of our results. While access to the D-Wave Ising machine is
required for some of the results, similar outcomes can be achieved using
Simulated Annealing, which is more easily reproducible. In particular, we
provide our modified Simulated Annealing sampler, based on D-Wave’s Simulated
Annealing code, which we employed to perform both the free and nudge phases
with the same annealing schedule baseline, akin to our approach with the
D-Wave Ising machine.

###  5.1 Methods for handling the D-Wave Ising machine

The training of neural networks on the D-Wave Ising machine were conducted as
follows.

We accessed the Ising machine virtually through a Python API ocean-api-dwave
provided by D-Wave. Initially, we needed to define a sampler, which refers to
the type of D-Wave machine we intended to use (e.g., Chimera or Pegasus
architecture) and the embedding procedure we employed to map the neural
network architecture onto the actual architecture of the Ising machine.

Following the initial implementation of EP scellier_equilibrium_2017 , we
performed the two sequential phases of EP on the Ising machine. The input
image remained static during all the annealing process, as required for the
system to converge given a fixed input. We had to adjust the duration of the
annealing for the free phase. We chose the native duration which is 20$\mu s$
(see Fig. 1b). For the reverse annealing we specified the duration and the
schedule. We set the duration to 40 $\mu s$, the initial annealed fraction to
0, the annealed fraction at time 20$\mu s$ to $0.25$ and the final annealed
fraction to $0$ (see Fig. 1b). This enabled the system to change its state
according to the nudging signal. It seems that the annealed fraction needed at
the middle of the reverse annealing depends on the size of the architecture to
train (see Supplementary Note LABEL:sec:supp-sa with Simulated annealing).

The D-Wave Ising machine is not perfect and thus does not return the ground
state of the problem at each call to the solver (see Fig. 4 \- distribution
before training). This alters the training procedure and we had to sample the
problem multiple times per input image to get a good estimate of the ground
state. For that purpose, we sampled the same problem 10 times and selected the
sample with the lowest energy.

We used the Ising formalism to sample our problems - i.e. the spins are $\pm
1$ and not $0/1$ as with a QUBO (Quadratic Unconstrained Binary Optimization)
- as the D-Wave Ising machine really optimizes the Ising problem. We can
submit problems in the QUBO formulation to the D-Wave Ising machine but the
parameters (couplings and bias) are scaled non-trivially, which we found to
affect the training procedure.

An essential feature to disable when executing a two-phase training algorithm,
such as EP, on the D-Wave Ising machine, is the auto-scale feature:

$scaling=\max(\\\ \max(\frac{\max(\mathbf{h},0)}{\max(h_{range},0)}),\\\
\max(\frac{\min(\mathbf{h},0)}{\min(h_{range},0)}),\\\
\max(\frac{\max(\mathbf{J},0)}{\max(J_{range},0)}),\\\
\max(\frac{\min(\mathbf{J},0)}{\min(J_{range},0)}))\\\ $ (9)

This feature enables the scaling (by a factor $scaling$) of problem parameters
to optimally fit within the accessible value range for the parameters on the
chip, in order to gain in performance. However, when training the Ising
machine with EP, we introduce nudging biases on the output neurons (which
always have a larger magnitude than the network biases). We immediately see
from Eq. (LABEL:eq:auto-scale) that the scaling for all the parameters will
change for the second phase. As a result, the optimization problem will be
significantly different, making the resulting gradient irrelevant. We found
that disabling this feature and manually scaling the parameters (see
Supplementary Table LABEL:tab:parameters-dwave-hyper) allows for the
successful training of a neural network on the D-Wave Ising machine. One
drawback of this approach is the need to clip the parameters within the
available range of values accessible on the chip, which we do after the
stochastic gradient descent (SGD) update step.

###  5.2 Objective function for training a neural network on an Ising machine
with Equilibrium Propagation

The particularities of Equilibrium Propagation on Ising spin models have been
discussed in the main text. Here, we revisit the choice of the objective
function to minimize. In standard software-based simulations of EP, we use the
Mean Squared Error (MSE) between the internal state of the output neurons and
their target states:

$C(y,\hat{y})=\sum_{i\in Y}(y_{i}-\hat{y}_{i})^{2}$ (10)

However, when using the Ising machine, we no longer have access to the
internal states $y_{i}$. The only measurement available to us is the
activation - the spin’s state $\sigma(y_{i})$ \- of the output neurons, as
$y_{i}$ is implicitly computed through the couplings. . This is the reason
behind our choice of MSE for training the system (Eq.(4)).

Nevertheless, the choice of using the MSE as the cost function naturally fits
the Ising machine frameworks as the minimization of the MSE during the nudge
phase simply translates in a nudge bias (see Supplementary Note
LABEL:sec:suppl-nudge-ising).

###  5.3 Methods for training a fully-connected architecture

We now discuss the two primary challenges associated with training a fully-
connected neural network on the D-Wave Ising machine. The first challenge
involves embedding a dense graph (the fully-connected architecture) onto a
sparsely connected graph, i.e., the architecture of the D-Wave Ising machine.
The second challenge concerns efficiently feeding inputs to the Ising machine.

####  5.3.1 Embedding the fully-connected architecture onto the chip layout

To map our problem, which is the underlying graph of the neural network we
want to train, onto the actual architecture of the chip, we used the
LazyFixedEmbeddingComposite function. This function operates as follows:
during the first annealing, it calls a heuristic function minorminer, which
aims to find a way to embed the neural network onto the chip layout. However,
this procedure is itself an NP-hard problem, resulting in a time-consuming
process. LazyFixedEmbeddingComposite finds the embedding only for the first
annealing and reuses this embedding for subsequent training examples, which
significantly reduces the annealing time. In addition to this advantage, it
allows the training process to account for the local imperfections of the chip
(faulty spins, noise, etc.) so that the parameters are updated based on the
actual dynamics of the chip.

The embedding process relies heavily on the chaining procedure - using
multiple hardware spins to represent a single spin in order to couple it to
more neighbors than the architecture allows for - and the chaining strength
can be adjusted. In practice, we set the chaining strength (i.e., the value of
the couplings between the hardware spins that represent the same spin) to $-1$
which has proven effective.

For training a fully-connected architecture, we used the Pegasus architecture
on the Advantage 4 chip available at the time of the simulations.

####  5.3.2 Feeding the input data to the chip

As mentioned earlier, the input image remains static throughout the entire
annealing process, as required for EP. This enables us to perform the vector-
matrix product between the input image and the first weight matrix in silico
only once, and apply the result as an input bias to the hidden neurons. This
approach allows us to train the Ising machine with a large input (MNIST: 28x28
pixels), which was previously impossible with existing training methods that
required embedding the input data on the chipadachi_application_nodate ;
niazi2023training . This weight matrix is also trained, given only the steady-
state spin states in the Ising machine and the input data.

####  5.3.3 Training data

Training a neural network reduces to an iterative process where the training
data is shown multiple times to the network. So training a neural network on
the D-Wave Ising machine requires multiple calls to the solver, each being
billed according to the time the solver is used.

For instance, a free phase is billed the following way: time to initialize the
parameters on the chip (couplings and field biases): $\approx 7\;ms$, time to
thermalize the chip after the initialization of the parameters: $\approx
1\;ms$, annealing time: $\approx 20\;\mu s$ being done 10 times per example in
our case, readout and thermalization time before next annealing: $\approx
200\;\mu s$ which result in $\approx 10\;ms$ per free phase. The nudge phase
is even more costly in time as we need to re-initialize the spins for each
reverse annealing which adds $10*7\;ms$ per example and gives $\approx
80\;ms$. In total this results in $\approx 9\;ms$ per training example (free
and nudge phases). An epoch on full MNIST (60k training images and 10k testing
images) would be $\approx 5525\;s$ on the solver. And 50 epochs would take
$\approx 5k\;\min$ which is $\approx 77\;h$. This would result in an insane
bill, that is why we chose to use a subset of MNIST: MNIST/100. An epoch with
1000 training examples would be $\approx 91\;s$ on the solver, 50 epochs would
take $\approx 76\;min$ which is $\approx 1.27\;h$. So to train one neural
network on the D-Wave Ising machine for 50 epochs on MNIST/100 we have to pay
for $\approx 1.27\;h$ access time, which is much more accessible than the
$76h$ for full MNIST.

Following the notation of truncated-mnist , our dataset is called MNIST/100,
which contains 1000 training images with 100 training images per class, and
100 testing images with 10 testing images per class. To create our training
(resp. testing) dataset, we select the first 100 (resp. 10) images of each
class from the MNIST dataset downloaded with Pytorch. This equidistribution
helps avoid training (resp. testing) bias in such a small training (resp.
testing) dataset. We did not use data augmentation techniques.

Furthermore, to save access time to the Ising machine, we adopted a simple
scheme similar to refs martin_eqspike_2020 ; park-skip-nudge . At the end of
the free phase, we compared the state of the output layer of the network to
the target. If the output layer was equal to the target vector, we skipped the
nudge phase for that particular input data, meaning no gradient had to be
computed for this image. This approach significantly saved access time and
accelerated the training.

####  5.3.4 Learning rule for the fully-connected architecture

The D-Wave Ising machine minimizes the following Ising energy function:

$E_{Ising}=\sum_{i>j}J_{ij}\sigma_{i}\sigma_{j}+\sum_{i}h_{i}\sigma_{i}$ (11)

Thus, the learning rule for the couplings is directly derived from Eq. (11):

$\frac{\partial\mathcal{L}}{\partial
J_{ij}}=\frac{1}{\beta}\left(\frac{\partial E_{Ising}}{\partial
J_{ij}}(x,\theta,\sigma^{*,\beta},\beta,\hat{\sigma})-\frac{\partial
E_{Ising}}{\partial J_{ij}}(x,\theta,\sigma^{*,0},0,0)\right)$ (12)

where $\sigma^{*,0}$ and $\sigma^{*,\beta}$ stands for the two sequential free
and nudge equilibrium states.

Finally, Eq. (12) directly read as:

$-\frac{\partial\mathcal{L}}{\partial J_{ij}}=\Delta
J_{ij}=-\frac{1}{\beta}\left((\sigma_{i}\sigma_{j})^{*,\beta}-(\sigma_{i}\sigma_{j})^{*,0}\right)$
(13)

where $\mathcal{L}$ stands for the loss function to be minimized during the
training procedure - here the Mean Squared Error function.

This learning rule differs from a negative sign to the conventional learning
rule of fully-connected layers with EP as conversely to the standard energy
minimized in EP, the D-Wave Ising machine minimizes this specific Ising energy
function (Eq. (11)) where there is a positive sign in front of the coupling
sum.

Similarly, we can derive the learning rule for the bias fields:

$-\frac{\partial\mathcal{L}}{\partial h_{i}}=\Delta
h=-\frac{1}{\beta}\left(\sigma_{i}^{*,\beta}-\sigma_{i}^{*,0}\right)$ (14)

We feed these gradients to a SGD optimizer without momentum and no mini-batch:

$J\leftarrow J+\eta\cdot\Delta J$ (15)

where $\eta$ stands for the learning rate which is a tunable parameter.

###  5.4 Methods for training a convolutional architecture

In this section, the main challenge, besides training the system with EP, is
to find the correct embedding that implements the convolutional neural network
we want to train.

####  5.4.1 Handcrafting the embedding dictionary

Contrary to the embedding procedure for the fully-connected architecture, we
handcrafted the embedding for this case, as we leveraged the actual
architecture to do specific computations. To achieve this, we manually created
a dictionary that maps the index of the spins in our architecture to a
specific site on the chip. For the spins requiring chaining (e.g., for
implementing the results of the average pooling operation), the dictionary
maps the spin to a list of hardware sites.

For the convolutional architecture, we used the Chimera architecture on the
chip DW-2000.

A typical embedding dictionary is:

        embedding={0: [560], 1: [561], ..., 39: [579]}

where the key of the dictionary is the index of a specific neuron in the
architecture and the linked list is the corresponding spin(s) on the chip. The
embedding is shown in Supplementary Fig. LABEL:fig:figsm11.

####  5.4.2 Feeding the input data to the chip

The concept behind the convolutional architecture is centered around the
”local” product of input data and the coupling. To demonstrate that this works
in practice for realizing convolutional operations, we had to embed the inputs
directly on the chip, unlike the fully-connected architecture. For that
purpose, we set strong biases (i.e., $\pm 4$, the largest possible value for
biases on DW-2000) on the spins corresponding to the inputs. This ensures that
they remain constant throughout the entire annealing procedure and have the
binary value corresponding to the sign of the bias applied to them.

####  5.4.3 Learning rule for the convolutional architecture

The learning process for the convolutional architecture differs from that of
the fully-connected model. The operations involved in the energy function have
changed, resulting in a different learning rule for the convolutional weights.
However, the learning for the last weights involved in the classifier remains
the same as for the standard fully-connected architecture.

The learning rule for the convolutional weights $\mathbf{J}$ between the input
$x$ and the ”hidden” layer $h$ \- after activation but before the pooling
operation - is:

$\Delta\mathbf{J}=\frac{1}{\beta}\left[\left(x\star
h\right)^{*,\beta}-\left(x\star h\right)^{*,0}\right]$ (16)

where $\star$ denotes the convolution operation between $x$ and $h$. Here the
gradient $\Delta\mathbf{J}$ has the same dimension as the convolutional
weights tensor $\mathbf{J}$, so we can directly apply the update through
stochastic gradient descent:

$\mathbf{J}\leftarrow\mathbf{J}+\eta\cdot\Delta\mathbf{J}$ (17)

The learning rule Eq. (16) can be easily extended to the layers that are not
directly linked to the fixed inputs:

$\Delta\mathbf{J}=\frac{1}{\beta}\left[\left(h^{\ell}\star
h^{\ell+1}\right)^{*,\beta}-\left(h^{\ell}\star
h^{\ell+1}\right)^{*,0}\right]$ (18)

where $h^{\ell}$ and $h^{\ell+1}$ stands for the states of the two consecutive
layers $\ell$ and $\ell+1$.

###  5.5 Methods for the Simulated Annealing simulations

We performed digital simulations using Simulated Annealing (SA) to benchmark
the results obtained on the D-Wave Ising machine. SA shares the same goal as
quantum annealing algorithms, but relies on temperature to control the
probability of a system to escape from a given configuration. The temperature
is initially set to a high value, allowing the system to explore various
configurations. It is then gradually decreased, so that y the end of the
annealing process, the system reaches a steady state (see Alg.
LABEL:alg:simulated-annealing).

We used a code provided by D-Wave to perform those simulations.

However, similarly to the ”auto-scale” feature that had to be disabled on the
D-Wave Ising machine, we manually set the temperature range, which is natively
auto-scaled to match the range of parameters for the problem to be solved
through simulated annealing. The temperature range is chosen based on the
initial parameter values and remains fixed throughout the training. Although
the temperature range may not be ideal, it works well in practice.

We also modified the original code to implement a similar kind of reverse
annealing to the one of D-wave, but with temperature. In this case, we
reversed the temperature schedule during the first annealing for a certain
duration up to a tunable value (as with simulations on D-Wave) and then
decreased it back to zero. The temperature schedule - i.e., the schedule of
$p_{1\rightarrow 2}$ as here $p_{1\rightarrow 2}\propto e^{-\frac{\Delta
E_{1\rightarrow 2}}{T}}$\- follows the same curve as in Fig. 1.b.

To run the simulations, we integrated the Simulated Annealing sampler into the
code developed for the training on the D-Wave Ising machine. The user simply
needs to specify which sampler to use for a specific training.

Similar to the trainings on the D-Wave Ising machine, the states reached with
SA are stochastic. As a result, we also had to sample the states multiple
times per image to obtain a reliable estimate of the ground state.
Additionally, we skipped the nudge phase when the free phase ad already
produced the correct output state.

###  5.6 Methods for the deterministic simulations

We also conducted digital simulations using a similar type of neural network:
one with binary activations and real-valued weights but with deterministic
dynamics, i.e., a gradient dynamics on the energy function.

For these simulations, we adapted the code from Laydevant_2021_CVPR , in which
the neurons are described by the standard EP energy function:

$E_{EP}=\sum_{i}s_{i}^{2}-\sum_{i>j}W_{ij}\rho(s_{i})\rho(s_{j})-\sum_{i}b_{i}\rho(s_{i})$
(19)

where the weights $W_{ij}$ are real-valued and $\rho$ is the binary Heaviside
step function (see Fig. 1.c):

$\rho(s)=\begin{cases}0\;if\;s<0\\\ 1\;else\end{cases}$ (20)

A dynamical equation for the neurons that minimizes this energy function is
given by the following gradient dynamics:

$\frac{ds}{dt}=-\frac{\partial E_{EP}}{\partial s}$ (21)

where $s$ stands for the internal state of a specific neuron in the network.
For a particular neuron $s_{i}$, we can simply rewrite this equation as:

$\frac{ds_{i}}{dt}=-s_{i}+\rho^{\prime}(s_{i})\left[\sum_{j}W_{ij}\rho(s_{j})+b_{i}\right]$
(22)

This equation governs the internal dynamics of the binary neurons - with the
particular choice of $\rho^{\prime}(s_{i})=1_{0<s_{i}<1}$ which is arbitrary
as the derivative of the Heaviside step function is almost zero everywhere.

We use an Euler scheme to solve this dynamics:

$s_{t+dt}=s_{t}+dt*\frac{ds}{dt}$ (23)

where the time step $dt$, the number of time steps for the free and nudge
phases resp. $T$ and $K$ are hyperparameters to tune (see Supplementary Table
LABEL:tab:parameters-dwave-hyper-deterministic). No annealing is used here so
given an initial state (always 1 here) and a set of parameters, the dynamics
always converges toward the same steady state so we do not need to repeat the
simulation for each image.

For the nudge phase, the energy function is augmented by a cost term as
follows:

$\frac{ds}{dt}=-\frac{\partial E_{EP}}{\partial s}-\beta\frac{\partial
C}{\partial s}$ (24)

where $C$ is the cost function to be minimized and $\beta$ the nudging
parameter.

The output neurons $y$ now have the following dynamics:

$\frac{dy_{i}}{dt}=-y_{i}+\rho^{\prime}(y_{i})\left[\sum_{j}W_{ij}\rho(s_{j})+b_{i}\right]+\beta\cdot(y_{i}-\hat{y}_{i})$
(25)

where $\hat{y}_{i}$ is the target state for the corresponding neuron $y_{i}$.

We compute the gradients given the free and nudge equilibrium states:

$-\frac{\partial\mathcal{L}}{\partial W_{ij}}=\Delta
W=\frac{1}{\beta}\left((\rho(s_{i})\rho(s_{j}))^{*,\beta}-(\rho(s_{i})\rho(s_{j}))^{*,0}\right)$
(26)

and feed them to a SGD optimizer (Eq. (15)).


# The Diophantine equation $x^{4}\pm y^{4}=iz^{2}$ in Gaussian integers

Filip Najman

###### Abstract

In this note we find all the solutions of the Diophantine equation $x^{4}\pm
y^{4}=iz^{2}$ using elliptic curves over $\mathbb{Q}(i)$. Also, using the same
method we give a new proof of Hilbert’s result that the equation $x^{4}\pm
y^{4}=z^{2}$ has only trivial solutions in Gaussian integers.

##  1 Introduction.

The Diophantine equation $x^{4}\pm y^{4}=z^{2}$, where $x,y,$ and $z$ are
integers, was studied by Fermat, who proved that there exist no nontrivial
solutions. Fermat proved this using the _infinite descent_ method, proving
that if a solution can be found, then there exists a smaller solution (see for
example [1, Proposition 6.5.3]). This was the first particular case of
Fermat’s Last Theorem proven (the theorem was completely proven by Wiles in
[11]).

The same Diophantine equation, but now with $x,y,$ and $z$ being Gaussian
integers, i.e., elements of $\mathbb{Z}[i]$, was later examined by Hilbert
(see [3, Theorem 169]). Once again, it was proven that there exist no
nontrivial solutions. Other authors also examined similar problems. In [9]
equations of the form $ax^{4}+by^{4}=cz^{2}$ in Gaussian integers with only
trivial solutions were studied. In [2] a different proof than Hilbert’s is
given, using descent, that $x^{4}+y^{4}=z^{4}$ has only trivial solutions in
Gaussian integers. The equations $x^{4}\pm y^{4}=z^{2}$ over some quadratic
fields were also considered in [7] and [8], again proving that there exist no
nontrivial solutions. Some applications of Diophantine equations of this type
can be found in [12] and [13].

In this short note, we will examine the Diophantine equation

$x^{4}\pm y^{4}=iz^{2}$

in Gaussian integers and find all the solutions of this equation. Also, we
will give a new proof of Hilbert’s results. Our strategy will differ from the
one used by Hilbert and will be based on elliptic curves. Elliptic curves have
also been used in [4] to prove that the Diophantine equation
$x^{3}+y^{3}=z^{3}$ has only trivial solutions in Gaussian integers, but in a
somewhat different way than in this note. We will use elliptic curves over
quadratic fields that have nontrivial torsion, while in [4], an elliptic curve
with trivial torsion over the rationals was examined.

For an elliptic curve $E$ over a number field $K$, it is well known, by the
Mordell-Weil theorem, that the set $E(K)$ of points on the curve $E$ with
coordinates in the field $K$ is a finitely generated abelian group. The group
$E(K)$ is isomorphic to $T\oplus\mathbb{Z}^{r}$, where $r$, which is called
the _rank_ , is a nonnegative integer and $T$, which is called the _torsion
subgroup_ , is the group of all elements of finite order. Thus, there are
finitely many points on an elliptic curve over a field if and only if it has
rank 0.

We will be interested in the case when $K=\mathbb{Q}(i)$. We will work only
with elliptic curves with rational coefficients, and by a recent result of the
author (see [5]), if an elliptic curve has rational coefficients, then the
torsion of the elliptic curve over $\mathbb{Q}(i)$ is either cyclic of order
$m$, where $1\leq m\leq 10$ or $m=12$, of the form
$\mathbb{Z}_{2}\oplus\mathbb{Z}_{2m}$, where $1\leq m\leq 4$, or
$\mathbb{Z}_{4}\oplus\mathbb{Z}_{4}$.

Throughout this note, the following extension of the Lutz-Nagell Theorem is
used to compute torsion groups of elliptic curves.

######  Theorem (Extended Lutz-Nagell Theorem).

Let $E$ be an elliptic curve in the form $E:y^{2}=x^{3}+Ax+B$ with
$A,B\in\mathbb{Z}[i]$. If a point $(x,y)\in E(\mathbb{Q}(i))$ has finite
order, then:

  1. 1.

$x,y\in\mathbb{Z}[i].$

  2. 2.

Either $y=0$ or $y^{2}|4A^{3}+27B^{2}.$

The main step in the proof of the Lutz-Nagell Theorem (for curves over
$\mathbb{Q}$) is to show that all the torsion points have integer coordinates.
This is done by showing that no prime can divide the denominators of the
coordinates of the torsion points. The proof of the Lutz-Nagell Theorem can
easily be extended to elliptic curves over $\mathbb{Q}(i)$. Details of the
proof can be found in [10, Chapter 3]. An implementation in Maple can be found
in [10, Appendix A].

Note that every elliptic curve can be put in the form $E:y^{2}=x^{3}+Ax+B$
(this is the _short Weierstrass form_) over any field of characteristic zero,
and thus in particular over $\mathbb{Q}(i)$. The Extended Lutz-Nagell Theorem
enables us to easily get a finite list of possible candidates for the torsion
points, and then check which ones are actually torsion points. Although when
all the torsion points are found, one could easily compute the group structure
of the torsion subgroup using addition laws, all we need is the list of all
the torsion points.

##  2 The Diophantine equation $x^{4}\pm y^{4}=iz^{2}$.

######  Definition 1.

We call a solution $(x,y,z)$ of the Diophantine equation

$x^{4}\pm y^{4}=cz^{2},$

for some given $c\in\mathbb{C}$, _trivial_ if $xyz=0$.

We are now ready to prove our main result.

######  Theorem 1.

  * (i)

The equation $x^{4}-y^{4}=iz^{2}$ has only trivial solutions in Gaussian
integers.

  * (ii)

The only nontrivial solutions satisfying $\gcd(x,y,z)=1$ in Gaussian integers
of the equation $x^{4}+y^{4}=iz^{2}$ are $(x,y,z)$, where $x,y\in\\{\pm i,\pm
1\\},\ z=\pm i(1+i)$.

_Proof:_
(i) Suppose $(x,y,z)$ is a nontrivial solution. Dividing the equation by
$y^{4}$ and making the variable change $s=x/y,\ t=z/y^{2}$, we obtain the
equation $s^{4}-1=it^{2}$, where $s,t\in\mathbb{Q}(i)$. We can rewrite this
equation as

$r=s^{2},$ (1) $r^{2}-1=it^{2}.$ (2)

Multiplying these equations we obtain $i(st)^{2}=r^{3}-r$. Again, making the
variable change $a=st,\ b=-ir$ and dividing by $i$, we obtain an equation
defining an elliptic curve

$E:a^{2}=b^{3}+b.$

Using the program [6], written in PARI, we compute that the rank of this curve
is 0. It is easy to compute, using the Extended Lutz-Nagell Theorem, that
$E(\mathbb{Q}(i))_{tors}=\mathbb{Z}_{2}\oplus\mathbb{Z}_{2}$ and that
$b\in\\{0,\pm i\\}$. It is obvious that all the possibilities lead to trivial
solutions.
(ii) Suppose $(x,y,z)$ is a nontrivial solution satisfying $\gcd(x,y,z)=1$.
Dividing the equation by $y^{4}$ and making the variable change $s=x/y,\
t=z/y^{2}$, we obtain the equation $s^{4}+1=it^{2}$, where
$s,t\in\mathbb{Q}(i)$. We can rewrite this equation as

$r=s^{2},$ (3) $r^{2}+1=it^{2}.$ (4)

Multiplying these equations we obtain $i(st)^{2}=r^{3}+r$. Again, making the
variable change $a=st,\ b=-ir$ and dividing by $i$, we obtain an equation
defining an elliptic curve

$E:a^{2}=b^{3}-b.$

Using the program [6], we compute that the rank of this curve is 0. Using the
Extended Lutz-Nagell Theorem we compute that
$E(\mathbb{Q}(i))_{tors}=\mathbb{Z}_{2}\oplus\mathbb{Z}_{4}$ and that
$b\in\\{0,\pm i,\pm 1\\}$. Obviously $b=0$ leads to a trivial solution. It is
easy to see that $b=\pm 1$ leads to $r=\pm i$ and this is impossible, since
$r$ has to be a square by (3). This leaves us the possibility $b=\pm i$. Since
we can suppose that $x$ and $y$ are coprime, this case leads us to the
solutions stated in the theorem. ∎

##  3 A new proof of Hilbert’s results.

We now give a new proof of Hilbert’s result, which is very similar to Theorem
1.

######  Theorem 2.

The equation $x^{4}\pm y^{4}=z^{2}$ has only trivial solutions in Gaussian
integers.

_Proof:_
Suppose $(x,y,z)$ is a nontrivial solution. Dividing the equation by $y^{4}$
and making the variable change $s=x/y,\ t=z/y^{2}$, we obtain the equation
$s^{4}\pm 1=t^{2}$, where $s,t\in\mathbb{Q}(i)$. We can rewrite this equation
as

$r=s^{2},$ (5) $r^{2}\pm 1=t^{2},$ (6)

and by multiplying these two equations, and making the variable change $a=st$,
we get the two elliptic curves

$a^{2}=r^{3}\pm r.$

As in the proof of Theorem 1, both elliptic curves have rank 0 and it is easy
to check that all the torsion points on both curves lead to trivial
solutions.∎

######  Remark.

Note that from the proofs of Theorems 1 and 2 it follows that the mentioned
solutions are actually the only solutions over $\mathbb{Q}(i)$, not just
$\mathbb{Z}[i]$.

### Acknowledgments.

The author would like to thank the referees for many helpful suggestions. The
author was supported by the Ministry of Science, Education and Sports,
Republic of Croatia, Grant 037-0372781-2821.

## References

  * [1] H. Cohen, _Number Theory, Volume I: Tools and Diophantine Equations_ , Graduate Texts in Mathematics, Springer, New York, 2007.
  * [2] J. T. Cross, In the Gaussian integers $\alpha^{4}+\beta^{4}\neq\gamma^{4}$, _Math. Mag._ 66 (1993) 105–108.
  * [3] D. Hilbert, Die Theorie der algebraischen Zahlkörper, _Jahresber. Deutsch. Math.-Verein._ 4 (1894/1895) 517–525.
  * [4] E. Lampakis, In Gaussian integers $x^{3}+y^{3}=z^{3}$ has only trivial solutions - a new approach, _Integers_ 8 (2008) A32.
  * [5] F. Najman, Torsion of elliptic curves over quadratic cyclotomic fields (to appear).
  * [6] D. Simon, Le fichier gp (2007), available at http://www.math.unicaen.fr/~simon/ell.gp.
  * [7] S. Szabó, The Diophantine equation $x^{4}+y^{4}=z^{2}$ in $\mathbb{Q}(\sqrt{-2})$, _Indian J. Pure Appl. Math._ 30 (1999) 857–861.
  * [8] ———, The Diophantine equation $x^{4}-y^{4}=z^{2}$ in three quadratic fields, _Acta Math. Acad. Paedagog. Nyházi. (N.S.)_ 20 (2004) 1–10.
  * [9] ———, Some fourth degree diophantine equations in Gaussian integers, _Integers_ 4 (2004) A16.
  * [10] T. Thongjunthug, _Elliptic Curves Over $\mathbb{Q}(i)$,_ Honours thesis, University of New South Wales, Sydney, 2006.
  * [11] A. Wiles, Modular elliptic curves and Fermat’s last theorem, _Ann. of Math. (2)_ 141 (1995) 443–541.
  * [12] K. Xu and H. Qin, Some Diophantine equations over $\mathbb{Z}[i]$ and $\mathbb{Z}[\sqrt{-2}]$ with applications to $K_{2}$ of a field, _Comm. Algebra_ 30 (2002) 353–367.
  * [13] K. Xu and Y. Wang, Several Diophantine equations in some ring of integers of quadratic imaginary fields, _Algebra Colloq._ 14 (2007) 661–668.

Department of Mathematics, University of Zagreb, Bijenička cesta 30, 10000
Zagreb, Croatia
<EMAIL_ADDRESS>