# Learning Hard Distributions with Quantum-enhanced Variational Autoencoders

Anantha S Rao  Department of Physics
Indian Institute of Science Education and Research Pune
Pune, India
<EMAIL_ADDRESS>Dhiraj Madan  IBM Research
India
<EMAIL_ADDRESS>Anupama Ray  IBM Research
India
<EMAIL_ADDRESS>Dhinakaran Vinayagamurthy  IBM Research
India
<EMAIL_ADDRESS>M.S.Santhanam  Department of Physics
Indian Institute of Science Education and Research Pune
Pune, India
<EMAIL_ADDRESS>
###### Abstract

An important task in quantum generative machine learning is to model the
probability distribution of measurements of quantum mechanical systems.
Classical generative models, such as generative adversarial networks (GANs)
and variational autoencoders (VAEs), can learn the distributions of product
states with high fidelity, but fail or require an exponential number of
parameters to model entangled states. In this paper, we introduce a quantum-
enhanced VAE (QeVAE), a generative quantum-classical hybrid model that uses
quantum correlations to improve the fidelity over classical VAEs, while
requiring only a linear number of parameters. We provide a closed form
expression for the output distributions of the QeVAE. We also empirically show
that the QeVAE outperforms classical models on several classes of quantum
states, such as 4-qubit and 8-qubit quantum circuit states, haar random
states, and quantum kicked rotor states, with a more than 2x increase in
fidelity for some states. Finally, we find that the trained model outperforms
the classical model when executed on the IBMq Manila quantum computer. Our
work paves the way for new applications of quantum generative learning
algorithms and characterizing measurement distributions of high-dimensional
quantum states.

##  I Introduction

Research in quantum information science promises to enable the development of
a fault-tolerant quantum computer that can perform certain tasks faster and
more efficiently than any classical computer. To this end, several scientists
have demonstrated that quantum algorithms can, in theory, outperform the best-
known conventional algorithms when tackling specific problems and, in some
situations, deliver a ‘quantum speedup’[1, 2]. For instance, certain quantum
algorithms can take exponentially fewer resources for tasks such as
factorization and eigenvalue decomposition, and quadratically fewer resources
to search through unsorted databases [3, 4, 5]. This pursuit of ‘quantum
speedup’ has motivated generations of physicists and engineers to discover
novel algorithms that leverage the properties of superposition, entanglement
and interference.

Through advances in processing power and algorithmic ability of computing
devices, machine learning techniques have evolved into fundamental tools for
detecting patterns in data. Theoretically, models like deep neural networks
have the potential to learn some of the most complex patterns that exist in
nature or human-made systems and have been demonstrated to be highly competent
at complex tasks like playing Go, identifying protein structures, and self-
driving cars [6, 7, 8]. However, many tasks are still intractable or very
expensive for these methods. Some learning tasks, for example, include
sampling from complex distributions or estimating the average values of
numerous parameters under a complicated distribution, both of which are
typically intractable (requiring exponential time or space resources).
Moreover, certain distributions derived from quantum-mechanical systems are
fundamentally intractable to conventional approaches [9, 10, 11]. Enhancing
and augmenting classical machine-learning methods using quantum correlations
has been the focus of quantum-enhanced machine learning and our work.

One of the important tasks for machine learning is that of generative
learning, that involves modeling or learning a distribution given independent
samples from the same. Previously, models such as Generative Adversarial
networks (GANs [12]) and Variational Autoencoders (VAEs [13]) have been used
to learn distributions of classical data such as texts and images.

In this paper, we study the problem of modeling the measurement distributions
obtained from unknown quantum states. This problem is fundamental in quantum
information science, as it can reveal useful information about the properties
and dynamics of quantum systems. Moreover, it can enable applications such as
quantum state reconstruction, and entanglement quantification [14, 15, 16,
17]. These applications can help us understand and manipulate quantum systems
in various fields, such as chemistry, materials science, and cybersecurity.

For our problem setting, Variational Autoencoders have also been shown to
learn and compress the measurement distribution obtained from product states,
termed as ‘easy states’ but require an exponential number of parameters to
learn the distribution of Haar random states (‘hard states’) [18, 19]. In this
work, we propose a quantum-enhanced VAE (QeVAE) which incorporates quantum
correlations through quantum circuits to enhance the performance of classical
VAEs. We show that a QeVAE constructed by substituting the decoder (generator)
of a conventional VAE with a parameterized quantum circuit can learn these
‘hard states’ with only a linear number of parameters in the number of qubits
and outperform conventional methods when the dataset contains quantum-
correlations.

In this paper, we make three main contributions to the field of quantum
generative learning. Firstly, we propose the Quantum-enhanced VAE (QeVAE),
which can enhance the expressive power of classical VAEs and produce
distributions that are classically intractable. Secondly, we provide a
mathematical closed-form expression to theoretically analyze the class of
models that QeVAE can model better than classical models. Finally, we
demonstrate experimentally that our QeVAE outperforms classical VAEs on
modeling measurement samples of different classes of quantum states, such as
haar random states, quantum circuit states, and quantum kicked rotor states,
and report an increase in fidelity by more than 2x for an 8-qubit quantum
circuit state.

Our work is organized as follows: In Section II, we discuss the background
work on generative learning, variational autoencoders, and the problem of
learning measurement distributions of quantum states. In Section III, we
propose the QeVAE and mathematically characterize the output distribution of
the model. We describe our experiments and results in section IV and V
respectively. Finally, in Section VI, we discuss our conclusions and future
outlook in the greater context of developing quantum algorithms for generative
learning.

##  II Background and Related work

###  II-A Generative learning and Variational Auto Encoders

The task of generative modeling involves modeling a parameterized distribution
$p_{\theta}(\mathbf{x})$, given independent and identically distributed (iid)
samples from the distribution $\\{\mathbf{x_{i}}\\}_{i=1}^{m}$ where each
$\mathbf{x_{i}}\sim p(\mathbf{x})$. Here the goal is to maximize the log
likelihood of the data given by
$L(\theta)=\sum_{i}\log{p_{\theta}(\mathbf{x_{i}})}$. One of the common
approaches for generative modeling is that of adversarial learning (for
example, GANs, [12]), which involves learning a parameterized generative
network and a parameterized discriminative network. While the generative
network maps samples from a fixed distribution to the data distribution, the
discriminative network distinguishes between data samples and the generated
samples.

On the other hand, variational learning involves modeling the distribution
through a latent vector $\mathbf{z}$. One can define a prior on the latent
variables $\mathbf{z}\sim p(\mathbf{z})$ and a parameterized likelihood
function $p_{\theta}(\mathbf{x}|\mathbf{z})$. The joint distribution over the
variables $\mathbf{x}$ and $\mathbf{z}$ can then be written as
$p(\mathbf{x},\mathbf{z})=p(\mathbf{z})p_{\theta}(\mathbf{x}|\mathbf{z})$. The
log likelihood can then be expressed as
$L(\theta)=\log(p_{\theta}(\mathbf{x}))=\log(\sum_{\mathbf{z}}p(\mathbf{x},\mathbf{z}))$.
Variational learning considers a lower bound for the above as:-

$\log(p_{\theta}(\mathbf{x}))\geq\operatorname*{\mathbb{E}}_{\mathbf{z}\sim
q(\mathbf{z})}\left(\frac{\log(p(\mathbf{x},\mathbf{z}))}{q(\mathbf{z})}\right)$
(1)

which holds for any distribution $q(\mathbf{z})$. The above lower bound called
Evidence Lower Bound (ELBO) holds tight when
$q(\mathbf{z})=p(\mathbf{z}|\mathbf{x})$. Variational Autoencoders [13] seek
to maximize the ELBO where $q_{\phi}(\mathbf{z}|\mathbf{x})$ can be defined
through a parameterized neural network also called the encoder. Here we can
split the ELBO as

$\text{ELBO}=\operatorname*{\mathbb{E}}_{q(\mathbf{z}|\mathbf{x})}\left(\log
p(\mathbf{x}|\mathbf{z})\right)-KL(q(\mathbf{z}|\mathbf{x})||p(\mathbf{z})))$
(2)

Figure 1: A schematic for a classical VAE. A classical VAE consists of an
encoder, decoder and a continuous latent dimension that can be modeled as a
multivariate gaussian with diagonal covariance. The dataset x is mapped to a
distribution p(z) and after training random samples x can be generated by
sampling from the decoder.

The algorithm involves using the posterior network (encoder) to create the
distribution $q_{\phi}(\mathbf{z}|\mathbf{x})$ and provides samples
$\mathbf{z}$, which are then passed through the likelihood network (decoder)
to produce a distribution over $\mathbf{x}$ conditioned on sampled vector
$\mathbf{z}$. One then seeks to maximize the difference of log likelihood of
generated distribution $p(\mathbf{x}|\mathbf{z})$ and the KL divergence
between posterior ($q(\mathbf{z}|\mathbf{x})$) and prior distributions
$p(\mathbf{z})$. The first term (reconstruction term) tries to maximize the
likelihood of recovering back $\bf{x}$ from the latent variable. The second
term is a regularization term that tries to ensure that posterior distribution
is close to prior. In order to be able to backpropagate through the sampling
step, one uses a trick called the reparametrization trick. For example, when
$q_{\phi}(\mathbf{z}|\mathbf{x})$ is a Gaussian with mean
$\mu(\mathbf{x};\theta)$ and a covariance matrix (typically diagonal)
$\Sigma(\mathbf{x};\theta)$, to sample
$\mathbf{z}\sim\mathcal{N}(\mathbf{\mu,\Sigma})$, one redefines
$\mathbf{z}=\mu+\Sigma^{1/2}\epsilon$ where $\epsilon\sim\mathcal{N}(0,I)$.
Separating out the noise $\epsilon$, enables us to backpropagate through the
$\mu$ and $\Sigma$ which depend on the parameters of the network.

Giving a higher weight to the second term leads to posterior becoming
independent of $\bf{x}$ (known as the posterior collapse problem[20]). In
practice one often weighs the KL divergence term by a hyperparameter $\beta$
to control the effect of regularization term [21], producing the overall cost
function as:

$\operatorname*{\mathbb{E}}_{q_{\phi}(\mathbf{z}|\mathbf{x})}\left(\log
p_{\theta}(\mathbf{x}|\mathbf{z})\right)-\beta
KL(q_{\phi}(\mathbf{z}|\mathbf{x})||p(\mathbf{z})))$ (3)

where $\beta$ is a hyperparameter that indicates the relative importance of
the regularization term with respect to reconstruction term.

###  II-B Quantum computation and quantum machine learning

Quantum computing is a model of computation with a potential to provide a
speedup over its classical counterpart. The fundamental units in quantum
computing framework consist of quantum bits or qubits. A qubit can take value
as a unit vector in two dimensional complex Hilbert Space $\mathbb{C}^{2}$.
The basis states {$\ket{0}$ and $\ket{1}$} correspond to the classical bits
{0,1}. An arbitrary state $\ket{\psi}$ can be considered as a superposition of
the basis states $\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$, where
$|\alpha|^{2}+|\beta|^{2}=1$.

More generally, an $n$ qubit state lies in a Hilbert space spanned by basis
states corresponding to $2^{n}$ classical bit strings, $\ket{0...0}$ through
$\ket{1...1}$. An arbitrary state can be a unit vector spanned by the $2^{n}$
basis strings as $\ket{\psi}=\sum_{x\in\\{0,1\\}^{n}}\alpha_{x}\ket{x}$, where
$\sum_{x}|\alpha_{x}|^{2}=1$.

If there are $n$ qubits each in state $\\{\ket{\psi_{i}}\\}_{i=1}^{m}$, then
the overall state of the n qubit system (product state) is
$\otimes_{i=1}^{n}\ket{\psi_{i}}$.

However there are $n$ qubits states that cannot be expressed in this form and
are refered to as entangled states. The measurement of these states yields
correlated bit strings.

The basic operations on qubits include quantum gates, which are unitary
operators acting on one or more qubits. Common examples of single qubits gates
include the Hadamard, Pauli Gates(X,Y,Z), and S gate.

The exponentiated Pauli Gates also provide us a set of parameterized gates,
$R_{x}(\theta)=exp(-i\frac{\theta}{2}X)$ ,
$R_{y}(\theta)=exp(-i\frac{\theta}{2}Y)$ and
$R_{z}(\theta)=exp(-i\frac{\theta}{2}Z)$. A 2 qubit gate can be described by a
4 dimensional unitaries. A common 2 qubit gate is $CNOT$ which acts as
$CNOT\ket{x,y}=\ket{x,x\oplus y}$ on basis states. A quantum circuit consists
of a sequence of single qubit and two-qubit gates acting on an initial state
(typically $\ket{0...0}$) followed by measurement of the final state. A
measurement of a state $\ket{\psi}=\sum_{x}\alpha_{x}\ket{x}$ yields one of
the bit strings $x$ with probabilities $|\alpha_{x}|^{2}$.

A parameterized quantum circuit (PQC) is a quantum circuit that has some gates
that depend on parameters $\theta$. The unitary operator of the PQC is denoted
by $\hat{U}(\theta)$, which can produce a state
$\ket{\psi(\theta)}=\hat{U}(\theta)\ket{0}^{\otimes n}$. The PQC can also be
conditioned on an input $\mathbf{x}$ as $\hat{U}(\mathbf{x},\theta)$. This can
be decomposed into a feature map $\hat{U}_{\phi}(\mathbf{x})$ and a trainable
ansatz $\hat{V}(\theta)$. The feature map is a data encoding circuit that
transforms input data $\bar{x}\in\mathbb{R}^{n}$ into a quantum state using
single qubit and two qubit gates, potentially parameterized by the input
variables. The ansatz is a parameterized circuit that consists of alternating
rotation layers and entanglement layers. The rotation layers are single-qubit
gates applied on all qubits. The entanglement layer uses two-qubit gates to
entangle the qubits according to a predefined scheme. The parameters of the
PQC can be optimized to achieve a certain goal, such as minimizing the energy
of a quantum system or maximizing the accuracy of a machine learning task. The
parameters can be updated iteratively using classical optimization methods
such as gradient-based (ADAM [22], SPSA [23]) or gradient-free (COBYLA)
algorithms. The measurement distribution of
$\ket{\psi(\mathbf{x},\theta)}=\hat{U}(\mathbf{x},\theta)\ket{0}$ in Z-basis
defines the conditional distribution
$p(\mathbf{y}|\mathbf{x})=|\braket{y}{\psi(\mathbf{x},\theta)}|^{2}$ $\forall
y\in\\{0,1\\}^{n}$.

A variational quantum algorithm (VQA) is then commonly defined as a hybrid
quantum-classical algorithm that utilizes a PQC to optimize a cost function
[24]. The cost function can be based on the sampled measurement distribution
of the state (QML, or our work), or computing the expectation value of a
Hamiltonian (for quantum chemistry problems). In the context of Quantum ML,
this model is also called a quantum neural network. Moreover, a quantum
circuit born machine (QCBM) is a generative quantum neural network that uses a
PQC to represent the probability distributions over bit strings as the
measurement distribution of the quantum state
$\ket{\psi(\theta)}=\hat{U}(\theta)\ket{0}^{\otimes n}$. The parameters of the
model can be trained by minimizing the cross entropy loss with the generated
samples. After training, the model generates samples from the desired
distribution.

###  II-C Prior work in Learning distributions of quantum states

In this section, we review prior work on learning the measurement distribution
of quantum states, i.e., given samples from measuring an $n$-qubit quantum
state

$\displaystyle\ket{\psi}=\sum_{x\in\\{0,1\\}^{n}}\alpha_{x}\ket{x},$ (4)

with probabilities $p(x)=|\alpha_{x}|^{2}$ for each bitstring $x$, we want to
find a model that approximates the distribution with parameters $\theta$ and
yields a distribution $p_{\theta}(x)$. This problem can have potential
applications in quantum state compression, quantum state transfer, and quantum
state discrimination. Previous works have shown that classical generative
models such as restricted Boltzmann machines (RBMs), variational autoencoders
(VAEs), and autoregressive models can model such distributions but require an
exponential number of parameters to learn and generate samples [25, 18, 19].
Moreover, a recent work based on the probably approximately correct (PAC)
framework has shown that these distributions can be efficiently learned with
quantum resources, but not with purely classical approaches [26]. However, the
quantum learner proposed in that work assumes the availability of a fault-
tolerant quantum computer.

There have been some previous works on employing quantum algorithms and the
variational autoencoder framework for different problems using classical data.
The first proposed algorithm within the VAE framework utilized the annealing-
based framework [27]. A recent work focused on improving the latent space
representation of classical VAEs through parameterized quantum circuits (PQCs)
[28]. Only recently, PQCs within the VAE framework were proposed for the
problem of drug discovery [29]. To the best of our knowledge, this is the
first work to consider quantum VAEs for data generation from quantum states.

We propose a quantum-enhanced variational autoencoder (QeVAE) that can learn
the measurement distribution of an unknown quantum state using noisy quantum
devices. Our QeVAE can reconstruct the measurement distribution through an
iterative learning process that involves a parameterized quantum circuit (PQC)
as the generative model and a classical neural network as the inference model.
We also show that our QeVAE reduces to a quantum circuit born machine (QCBM)
in the zero latentsize limit. We hypothesize that our QeVAE can leverage the
quantum properties of superposition and entanglement to learn the complex and
high-dimensional measurement distributions of quantum states.

##  III Quantum-enhanced Variational Autoencoders

In this section, we define a quantum-enhanced variational autoencoder (QeVAE)
to model the measurement distribution of an unknown $n$-qubit quantum state.
As a baseline, classical VAEs have already been used to model such a
distribution [18]. The hybrid model that we propose consists of a feedforward
classical encoder, a continuous latent space, and a parametrized quantum
circuit as a decoder. We model the approximate posterior (encoder network)
$Q_{\phi}(\mathbf{z}|\mathbf{x})$ through a classical feedforward neural
network and the latent variable as $\mathbf{z}\sim\mathcal{N}(0,I)$. The
likelihood (generator) distribution $p_{\theta}(\mathbf{x}|\mathbf{z})$ is
defined via a quantum circuit i.e.
$p_{\theta}(\mathbf{x}|\mathbf{z})=|\braket{x}{\hat{U}(\theta,\mathbf{z})}{0\rangle^{\otimes
n}}|^{2}$. The Evidence Lower bound loss (ELBO) is optimized through a
classical optimizer such as ADAM. The model is trained to mimic the given
measurement distribution of states. During training, the parameters of the
encoder and the rotation gates in the decoder (with a pre-selected
entanglement type) are iteratively varied and learned. Such a model has
multiple applications: It will enable scientists to generate certain quantum
states in different physical quantum computers just by knowing the set of
rotation and entangling gates to perform. The algorithm also has applications
in state compression and transferring a state from one system to another upto
a phase (entire phase information can be learnt with additional models[18]).
We now present a theorem that allows us to mathematically characterize the
class of distributions that can be obtained via the above model.

######  Theorem 1.

Consider a latent variable model with $\mathbf{z}\sim p(\mathbf{z})$ and
$p(\mathbf{x}|\mathbf{z})=|\bra{\mathbf{x}}V_{\theta}U_{\phi}(\mathbf{z})\ket{0^{n}}|^{2}$,
where $U_{\phi}(\mathbf{z})$ is a feature map and $V_{\theta}$ is a
parameterized ansatz. Then

  1. 1.

$\exists$ a density matrix $\rho$, such that
$p\mathbf{(}x)=\bra{x}V_{\theta}\rho V_{\theta}^{\dagger}\ket{x}$. In other
words the distribution of $\mathbf{x}$ can be obtained by evolving a density
matrix $\rho$ under unitary ansatz $V_{\theta}$ followed by measurement in
standard basis.

  2. 2.

Conversely, for each density matrix $\rho$, there exists a prior
$p(\mathbf{z})$ and a feature map $U_{\phi}(\mathbf{z})$, such that
$p\mathbf{(}x)=\bra{x}V_{\theta}\rho V_{\theta}^{\dagger}\ket{x}$.

###### Proof.

  1. 1.

$\displaystyle p(\mathbf{x})$ $\displaystyle=\int
p(\mathbf{x},\mathbf{z})d\mathbf{z}$ $\displaystyle=\int
p(\mathbf{z})|\bra{\mathbf{x}}V_{\theta}U_{\phi}(\mathbf{z})\ket{0^{n}}|^{2}d\mathbf{z}$
$\displaystyle=\int
p(\mathbf{z})\bra{\mathbf{x}}V_{\theta}U_{\phi}(\mathbf{z})\ket{0^{n}}\bra{0^{n}}U_{\phi}(\mathbf{z})^{\dagger}V_{\theta}^{\dagger}\ket{\mathbf{x}}d\mathbf{z}$
$\displaystyle=\bra{\mathbf{x}}V_{\theta}(\int
p(\mathbf{z})U_{\phi}(\mathbf{z})\ket{0^{n}}\bra{0^{n}}U_{\phi}(\mathbf{z})^{\dagger}d\mathbf{z})V_{\theta}^{\dagger}\ket{\mathbf{x}}$
( By linearity of inner product)

We define $\rho\coloneqq\int
p(\mathbf{z})U_{\phi}(\mathbf{z})\ket{0^{n}}\bra{0^{n}}U_{\phi}(\mathbf{z})^{\dagger}d\mathbf{z}$
Note

    1. (a)

$\rho$ is a valid density matrix i.e. $\rho^{\dagger}=\rho$, $\rho\geq 0$ and
$Tr(\rho)=I$

    2. (b)

$\rho$ is independent of both $\mathbf{x}$ and $\mathbf{z}$.

Then,

$\displaystyle p(\mathbf{x})$ $\displaystyle=\bra{\mathbf{x}}V_{\theta}\rho
V_{\theta}^{\dagger}\ket{\mathbf{x}}$

Thus, such the $p(x)$ is equivalent to preparing a density state $\rho$,
evolving under unitary $V_{\theta}$, and performing a measurement in standard
basis.

  2. 2.

For the second part, given a density matrix $\rho$, consider the sepctral
decomposition of $\rho$ as
$\rho=\sum_{z}\lambda_{z}\ket{\psi_{z}}\bra{\psi_{z}}$. Since $\rho$ is a
valid density matrix, $\lambda_{z}$’s define a distribution
$p(\mathbf{z})=\lambda_{z}$. Now, for each $\mathbf{z}$, choose a unitary such
that $U_{\phi(\mathbf{z})}\ket{0}=\ket{\psi_{z}}$. Now, one can verify that
such a feature map satisfies the required equation.

∎

Moreover, such a distribution defined using PQCs can encode information about
the entanglement structure of the states via the entangling ansatz
$V_{\theta}$ and density matrix $\rho$ which is inaccessible to purely
classical models. Note that this fully includes the set of QCBM distributions
since setting $\rho=\ket{0}\bra{0}$, will enable us to get QCBMs.

##  IV Methods

Figure 2: Ansatz and Feature maps for learning quantum distributions (a) A
two-local ansatz on three qubits with two repeating layers of $Rx$ and $Ry$
gates along with linear entanglement (b) A Pauli-Z feature map that embeds a
three-dimensional vector. ‘H’ represents the Hadamard gate. (c) A Pauli-ZZ
feature Map (d) A Pauli-(X, XY) feature map on a three qubit system.

In this section, we provide the methods followed to setup a hybrid quantum-
classical neural network that follows the variational autoencoder framework.
We provide details on the architecture of the model, the training algorithm,
metrics followed, datasets utilized and end the section with the
hyperparameters used for various models.

###  IV-A Architecture Details

The hybrid quantum-classical neural network consists of three components: An
encoder, a continuous latent-space and a decoder. In section we discuss the
structure of each component.

The _encoder_ is modeled through a classical network which defines the mean
$\mu(\mathbf{x},\theta)$, and diagonal covariance $\Sigma(\mathbf{x},\theta)$
for the posterior distribution $Q(\mathbf{z}|\mathbf{x})$, modeled as a
multivariate Gaussian. In our experiments, we use feed-forward neural networks
where the input bitstring is transformed through successive operations of
linear layers each followed by a non-linearity. Here, one starts with
$h^{0}=x$, which is the initial input. The $l^{th}$ layer transforms this as
$h^{l}=f^{l}(h^{l-1})=\phi(W^{l}h^{l-1}+b^{l})$, where $\Phi$ is the Leaky
ReLu activation function and $\\{W^{l},b^{l}\\}$ are the parameters of the
$l^{th}$ layer of the neural network (weights, biases). The final layer gives
us the vectors corresponding to mean $\mu$ and the diagonal entries of log
covariance matrix. Since we benchmark the performance of classical VAEs with
QeVAEs, the architecture of the encoder is constant over all simulations. The
encoder consists of 2 hidden layers each containing 8 and 7 neurons with the
leak of the activation function set to 0.01. The output of the encoder is
propagated to the latent space that is modelled as a continuous gaussian with
diagonal covariance. The encoded sample is obtained from the latent space
through the reparametrization trick.

The _decoder_ network in the QeVAE is modeled via a variational quantum
circuit with trainable rotation gates. We use a Pauli feature map (such as the
Z or ZZ featuremap given in figure 2) to encode the sampled latent variable
$\mathbf{z}$ onto the circuit. The Pauli featuremap is a data encoding circuit
that transforms input data $\bar{x}\in\mathbb{R}^{n}$
$U_{\Phi(\vec{x})}\ket{0}=\exp\left(i\sum_{S\subseteq[n]}\phi_{S}(\vec{x})\prod_{i\in
S}P_{i}\right)\ket{0}$ [30]. The variable $P_{i}\in\\{X,Y,Z,I\\}$ denotes the
Pauli matrices. The index $S$ describes connectivities between different
qubits in a given circuit. If $S=\\{Z\\}$, we obtain a $Z$ feature map. For
$S=\\{Z,ZZ\\}$ we have ZZ feature map.

We then use a two-local of ansatz with linear entanglement due to their
compatibility with current-era hardware. The two-local circuit is a
parameterized circuit consisting of alternating rotation layers and
entanglement layers. The rotation layers are single qubit gates applied on all
qubits while the entanglement layer uses two-qubit gates to entangle the
qubits according to a predefined scheme. In addition, while training, we find
it useful to add an additional linear feedforward layer before the quantum
circuit, after sampling a latent vector. This provides two benefits: (a)
Flexibility in choosing a latent size. The layer can linearly transforming the
latent vector of a different size to fit the input requirements of the quantum
circuit; (b) This also adds power to the network by introducing additional
parameters.

To benchmark the performance of the QeVAE, we train classical VAEs on
different 4 qubit and 8 qubit states, and different number of layers. The
number of trainable parameters is kept fixed close to $2^{\text{No. of
qubits}}$. Table I shows the sizes of classical decoder layers for circuits
corresponding to 8 qubit measurements.

TABLE I: No of neurons and number of hidden layers for different classical VAE benchmarks for 8 qubit states #Hidden layers | # parameters | #Neurons in each layer
---|---|---
1 | 256 | 8,15,8
2 | 269 | 8,10,9,8
3 | 255 | 8,6,7,9,8
4 | 266 | 8, 7, 6, 6, 7, 8
5 | 258 | 8, 7, 5, 4, 5, 7, 8

###  IV-B Training details

We train both the classical and hybrid models using ADAM optimizer. The
training is iterative where each iteration involves forward propagating a
single bitstring (‘x’) (obtained from the distribution) through the encoder,
obtaining a sample from the latent space (‘z’), preprocessing the sample
through a linear layer, embedding onto the PQC, propagating through the
quantum circuit, and computing the measurement distribution with multiple
shots. The loss function utilizes the quasi-probability of $p(x|z)$ and the KL
divergence of the latent space. To avoid overfitting on the training data, we
perform early stopping on the validation set loss with a patience factor
$\delta$. We use python packages Pytorch and the Torchconnector module in
qiskit to build these hybrid models[31, 32]. In table II, we summarize the
various hyperparameters used in the QeVAE model.

TABLE II:  No | Hyperparameter | Value
---|---|---
1 | Feature Map | Z - ZZ Pauli map
2 | Ansatz | Two-local
3 | Entanglement | Linear
4 | Measurement basis | Z
5 | Encoder learning rate | 0.001 - 0.01
6 | Decoder learning rate | 0.001 - 0.009
7 | No of latent variables | 0 - number of qubits
8 | Batchsize | 16-64
9 | KL-term weight ($\beta$) | 0.5-2
10 | Training type | Annealing, Fixed, Stepfunction
11 | Patience factor ($\delta$) | 5-7

###  IV-C Metrics

The metric we use to quantify the generated distribution is the fidelity
between two discrete distributions. If $\rho$ and $\sigma$ are $n$-qubit
states, we say that $\sigma$ is a good representation of $\rho$ if the
fidelity F = Tr($\sqrt{\rho^{1/2}\sigma\rho^{1/2}}$) $>1-\epsilon$ for an
$\epsilon>$ 0\. Through the result of [33], the fidelity can be expressed in
terms of the probability distributions over a measurement that maximally
distinguishes the two states. Thus given two random variables X, Y with
probabilities $p=(p_{1},p_{2},..p_{n})$ and $q=(q_{1},q_{2},..q_{n})$, the
fidelity of X and Y is defined to be the quantity:

$F(X,Y)=(\sum_{i}\sqrt{p_{i}q_{i}})^{2}$ (5)

where the measure $\sum_{i}\sqrt{p_{i}q_{i}}$ is the Bhattacharyya coefficient
between the two distributions. To compute the fidelity between the original
and learnt distributions, we propagate 5000 random samples from the latent
space through the decoder and construct the output distribution by computing
the average distribution.

###  IV-D Datasets

Figure 3: Different types of measurement datasets: (a) Product states obtained
by a combination of arbitrary single qubit gates. The orange box represents
repeatable layers of gates. (b) Quantum circuit states obtained from circuits
with local (nearest-neighbor) entanglement. (c) Haar states obtained from a
pure-quantum state by normalizing a $2^{n}$-dimensional complex vector. (d)
Quantum-kicked rotor states obtained by the time-evolution of an initial state
$|0\rangle$. Figure 4: QeVAE for learning quantum state distributions: (a)
Multiple copies of a quantum state $\rho$ are obtained naturally from a
quantum system and are measured through different measurement operators. The
measurement dataset is stored on a classical computer. (b) The QeVAE with a
parameterized quantum circuit as the generative network and a classical feed-
forward neural network as the inference network can be used to recrete the
distribution. After training, the circuit can be used to generate the original
distribution through any quantum computer and generate states amenable for
downstream processing. (c) A classical VAE with continuous Gaussian latent
variables that perform the same task.

In order to evaluate our proposed approach, we try to model measurements from
multiple families of states as described below. For each family we generate
multiple 4 and 8 qubit states with different random seeds. Multiple copies of
a state are measured in standard basis to yield a dataset consisting of
samples in $\\{0,1\\}^{n}$ for each n qubit state. Here, each measurement
dataset contains 1024 samples. We use 70% samples for training and the 30% for
validation. We benchmark the performance of QeVAE on several datasets as shown
in figure 3. We train both classical and quantum models to reproduce these
distributions. We are interested in discovering if a classical learner can
learn the same distribution and how the number of parameters required scales
with the size of the system. In the following section, we describe the
families of states that have been considered for creating our datasets.

_Random product states_ : Product states (i.e. Tensor Product of single qubit
states) are classically easy to simulate and are empirically found to be
_classically easy_ to learn. We generate random product states by simulating
quantum circuits with only single qubit gates with arbitrary angles of
rotation, generated according to a random seed (figure 3(a)). The state
prepared is on the form
$|\psi\rangle=\otimes_{i=0}^{n-1}\left\\{\alpha_{i}|0\rangle+\beta_{i}|1\rangle\right\\}$,
where $n$ is the number of qubits and projective Z-basis measurements generate
the samples, $x\in\\{0,1\\}^{n}$.

_Haar random states_ are quantum states that are uniformly distributed over
the Hilbert space according to the Haar measure and represent _classically
hard states_ i.e, they require exponential number of parameters in the number
of qubits to learn. These states can either be generated by first creating a
Haar unitary $U$ and then applying it on an initial state of dimension $2^{n}$
or by normalizing a complex-valued vector of dimension $2^{n}$. We use the
later method, where a complex-valued vector of
$|\psi\rangle=\sum_{l=0}^{2^{n}-1}(c_{1l}+ic_{2l})|l\rangle$ is initialized
with $|l\rangle$ corresponding to the orthonormal basis vector in the
$2^{n}$-dimensional Hilbert space, $\mathbb{C}^{2^{n}}$, and $c_{1l},c_{2l}$
are real numbers chosen independently from a standard Gaussian distribution.
This vector is normalized to yield a quantum state by using the constraint:
$\langle\psi|\psi\rangle=1$. After normalization, the states are uniformly
distributed on a unit hyper-sphere.

_Random quantum circuit states_ are obtained from random quantum circuits with
a pre-defined entanglement structure and circuit depth, as shown in figure
4(b). These states are useful for circuit compression and circuit compilation.

_Quantum kicked rotor states_ are obtained from the quantum kicked rotor
(QKR), a well-known model in quantum chaos and quantum information that
exhibits rich dynamics under time evolution. The QKR is defined by the
Hamiltonian:

$\displaystyle\hat{H}=\frac{\hat{p}^{2}}{2I}+k\cos{\hat{x}}\sum_{n}\delta(t-nT)$
(6)

where $\hat{x}$ and $\hat{p}$ are the position and momentum operators, $I$ is
the moment of inertia, $k$ is the kick strength, and $T$ is the kick period.
We set $I=1$ and introduce dimensionless parameters $\hbar_{s}=\hbar T/I$ and
$\kappa=k/\hbar$, where $\hbar$ is the reduced Planck’s constant and
$[\hat{x},\hat{p}]=i\hbar_{s}$. The time-evolution operator for one period is
then given by:

$\displaystyle\hat{U}=\hat{U}_{kick}\hat{U}_{free}=\exp{\left(-i\kappa\cos{\hat{x}}\right)}\exp{\left(-\frac{i}{2\hbar_{s}}\hat{p}^{2}\right)}$
(7)

To simulate the QKR, we apply $\hat{U}$ repeatedly to an initial state
$|\psi_{p}(0)\rangle=0$ and perform forward and inverse Fourier transforms
between each kick. The resulting wavefunction shows different behaviors
depending on the value of $\kappa$. For weak kicking ($\kappa\lessapprox
5.95$), the system exhibits quantum diffusion until a break time, where the
variance of momentum $\langle p^{2}\rangle$ grows linearly with time, and
after the break time, $\langle p^{2}\rangle$ saturates. For strong kicking
($\kappa\gtrapprox 5.95$), the system exhibits dynamical localization, where
$\langle p^{2}\rangle$ reaches a finite value. In contrast, the classical
kicked rotor shows chaotic behavior, where the future trajectory is highly
sensitive to initial conditions. For a comprehensive review of the QKR, we
refer the reader to [34]. In this work, we investigate whether a generative
model can learn the probability distribution of the wavefunction
$|\psi_{p}|^{2}$ after 1000 kicks for different values of
$\kappa\in\\{0.5,6\\}$ and $\hbar_{s}=1$.

##  V Results

In the tables presented below, we summarize the best fidelity obtained across
each type of measurement dataset. We compare the final fidelity between the
target distribution and that produced by a random uniform guess, a classical
variational autoencoder (CVAE), and a Quantum-enhanced variational autoencoder
(QeVAE). For each type of state, we consider five different random seeds.
QeVAE results include the best fidelity observed across different hyper-
parameters like latent size, feature-map, prepossessing-layer, and relative
KL-divergence term $\beta$.

TABLE III: Fidelity for Product states No qubits | 4 | 8
---|---|---
Seed | 12 | 16 | 27 | 44 | 102 | Mean $\downarrow$ | 12 | 16 | 27 | 44 | 102 | Mean $\downarrow$
Uniform | .471 | .356 | .302 | .156 | .436 | .344 | .164 | .306 | .081 | .306 | .106 | .193
CVAE | .998 | .995 | .998 | .995 | .995 | .996 | .983 | .979 | .982 | .983 | .987 | .983
QeVAE | .995 | .827 | .882 | .973 | .892 | .914 | .875 | .947 | .870 | .823 | .957 | .894
TABLE IV: Fidelity for Quantum circuit states No qubits | 4 | 8
---|---|---
Seed | 12 | 16 | 27 | 44 | 102 | Mean $\downarrow$ | 12 | 16 | 27 | 44 | 102 | Mean $\downarrow$
Uniform | .477 | .366 | .315 | .154 | .431 | .348 | .158 | .279 | .080 | .309 | .107 | .187
CVAE | .501 | .667 | .597 | .925 | .758 | .690 | .229 | .423 | .079 | .304 | .181 | .243
QeVAE | .981 | .976 | .950 | .873 | .912 | .938 | .665 | .654 | .388 | .548 | .591 | .569
TABLE V: Fidelity for Haar random states No qubits | 4 | 8
---|---|---
Seed | 42 | 96 | 27 | 101 | 102 | Mean $\downarrow$ | 12 | 43 | 16 | 27 | 2 | Mean $\downarrow$
Uniform | .772 | .776 | .777 | .771 | .768 | .773 | .766 | .770 | .773 | .781 | .772 | .772
CVAE | .795 | .798 | .800 | .788 | .788 | .794 | .754 | .755 | .757 | .766 | .763 | .759
QeVAE | .839 | .983 | .913 | .988 | .932 | .931 | .876 | .878 | .887 | .887 | .887 | .883
TABLE VI: Fidelity of Quantum-kicked rotor states No qubits | 4 | 8
---|---|---
Type | Localized (k=6) | Diffusive (k=0.5) | Mean $\downarrow$ | Localized (k=6) | Diffusive (k=0.5) | Mean $\downarrow$
Uniform | .175 | .838 | .506 | .053 | .418 | .236
CVAE | .723 | .908 | .815 | .061 | .406 | .233
QeVAE | .991 | .992 | .991 | .912 | .616 | .764
TABLE VII: Hardware results for a 4 qubit quantum circuit state State | Fidelity | Simulator | Hardware | Suppression | Mitigation
---|---|---|---|---|---
Uniform | 0.477 | $\checkmark$ |  |  |
CVAE | 0.501 | $\checkmark$ |  |  |
QeVAE | 0.981 | $\checkmark$ |  |  |
QeVAE | 0.658 |  | $\checkmark$ |  |
QeVAE | 0.642 |  | $\checkmark$ | $\checkmark$ | $\checkmark$

From tables(I-V), we observe that across all entangled quantum states we
considered (i.e. all clases of states other than product states) the final
fidelity obtained from a QeVAE outperforms the classical VAE and a random
guess. In addition, the number of learnable parameters in the classical VAE is
typically of the order $2^{(n)}$ (Table I) while those in a QeVAE is
$4n+\epsilon$ where $n$ is the number of qubits and $\epsilon$ is a constant
($\epsilon<$ 4). To further validate our findings, we run the best QeVAE
models on real quantum devices and see that the obtained fidelity is higher
than those achieved by classical methods (Table VII). We know from literature
that the measurement distribution obtained from product states are classically
_easy_ and that measurements obtained from states with quantum correlations is
classically _hard_. Our results not only corroborates with the above
observation but also shows that quantum-enhanced classical models can overcome
the drawbacks of purely classical models. Our main results and observations
from data is as follows:

  1. 1.

Our proposed algorithm achieves the highest fidelity across all types of
datasets, other than product states. The inherent ability of our model to
learn quantum correlations i.e we are able to produce entangled multi-qubit
states through variational quantum circuits and tailor the rotation gates to
reproduce a desired distribution allows them to outperform the classical model
for the quantum circuit, haar random and kicked rotor states.

  2. 2.

_Fidelity for Product states_ : Since product states do not employ quantum-
correlations between different qubits in a many-body system, a distribution on
$n$ bits obtained by measuring a product of $n$ qubits, can be characterized
by $n$ independent Bernoulli variables, thus requiring only $n$ parameters.
This polynomial dependence and presence of only classical correlations (each
bitstring in the distribution is independent) in the output distributions
enables classical VAEs to efficiently learn and reproduce the distribution
with very high fidelity. In addition, the ansatz in our model is not always
tailored for product states, i.e the ansatz has an entanglement structure and
the quantum state produced from the QeVAE will exhibit entanglement. Such a
class of distribution produced cannot in general be similar to that produced
by product states. This might explain the reduced fidelities for QVAEs.

  3. 3.

_Fidelity for Haar, Quantum circuit and Quantum kicked rotor states_ : We
notice that for measurement distributions obtained from more generic quantum
states, our quantum models outperform the classical models in terms of final
fidelity. In some cases (like seed 44 for quantum circuit states), the score
of the classical model is higher because the resulting measurement
distribution is highly concentrated around a single output ($>80\%$). In such
cases, it becomes easier to just predict the output statistically, rather than
learning the intrinsic structure of the quantum circuit producing the state.
Nonetheless, we find that QeVAE models require only $O(n)$ rotation gates to
learn the output distribution and achieve a fidelity score that is
unattainable to other methods.

  4. 4.

Since QeVAE is a hybrid model, it operates through the synergy of quantum and
classical resources. Classical models are well-versed at producing non-linear
transformations whereas quantum models are restricted to unitary or linear
transformations. Conversely, quantum models can encode quantum correlations
like entanglement that inaccessible to classical methods. An important feature
of our model is the ability to leverage the merits of both the classical and
quantum models, and outperform both models.

  5. 5.

The QeVAE model contains many hyper-parameters that can be tuned to achieve
optimal performance. In the above tables, we have presented the result for the
best hyper-parameter setup. When the latent-size is set to zero, there is no
contribution from the classical encoder and the output distribution is
produced by the circuit alone. In such cases, the resulting model is the
Quantum circuit Born Machine (QCBM), wherein the circuit parameters are
iteratively updated to minimize the difference between the output and target
distributions. Thus in the latent-size=0 limit, our model results in the QCBM
generative model. The KL divergence term in the ELBO loss function can be
neglected, and minimizing the negative expected log-likelihood becomes
equivalent to minimizing the KL-divergence between the output and target
distributions.

Hardware run: We execute the best model trained on the simulator on IBM
Hardware by transpiling the decoder circuit. To generate the output
distribution, random samples from $N(0,1)$ are propagated through the
preprocessor linear layer and then through the circuit (executed on hardware).
An average over all initial random points yields the desired output
distribution. Our results are depicted in figure 5 and in table VII. We find
that final fidelity is lesser than that on a simulator. Using error mitigation
and suppression techniques, we are able to perform better than the classical
VAE.

Figure 5: Inference on IBM hardware: The true distribution produced by an
unknown 4-qubit system (green) is used to train a QeVAE on IBM’s qasm-
simulator (blue) that results in a final of 0.97. After training, the decoder
is executed on the IBM-Manila. The measurement distribution produced from
hardware has a total fidelity of 0.658 which changes to 0.642 with error-
mitigation (Zero noise extrapolation) and error-suppression (Dynamic
decoupling)

###  V-A Circuit compilation

Figure 6: Circuit compilation with QeVAEs: (a) Original circuit structure that
produces a measurement distribution. The orange box represents a single layer
of rotation and entangling layer that is repeated 20 times. (b) The compiled
circuits requires very few gates when compared to the original circuit that
produces the state. Particularly, the number of CNOT gates is reduced by
$\sim$14X (c) The variational quantum circuit containing the Pauli-Z feature
Map that embeds data x from the latent space and the ansatz with only 16
parameters. (d) Results of training the ansatz to produce the measurement
distribution. We achieve a final fidelity of 0.956 (e) Original and the
distribution produced by the QeVAE after training.

We now present a practical application of QeVAEs on the task of circuit
compilation. Circuit compilation or circuit compression is an important area
of focus in the NISQ era. Since the depth of circuits executable on hardware
is limited, there is a need to transform deep circuits into shallow ones by
altering the sequence of gates and reducing the overall size and complexity.
Furthermore, multi-qubit gates like the CNOT gate, T-gate, and SWAP gates are
expensive to implement on hardware. Efficient circuit compilation assists in
faster computations, and allows accurate simulation of quantum systems. This
has immense applications in many-body physics, condensed matter physics, and
quantum chemistry. Through simulation studies, we show show QeVAEs can help in
reducing the complexity of circuits by learning to reproduce measurement
distributions with fewer gates.

We simulate an unknown quantum state by considering a deep quantum circuit
with twenty layers of rotation and entangling gates. In reality, the form of
the circuit is unknown and one only has access to the measurement data. A
projective measurement on such a state produces a measurement distribution as
shown in figure 4.5(e). Note that our goal here is to reproduce the
measurement distribution and not to learn the original state itself. Through
the QeVAE learning approach, we can learn the measurement distribution with
high fidelity. After training, we can discard the encoder part of the circuit,
and the decoder provides a sequence of gates that can be implemented on
hardware to generate the same distribution. We achieve a final fidelity of
0.956 (figure 6(d)) and a multi-fold reduction in the number of gates as seen
in figure 6(b). The initial circuit contains: 80 Rx gates, 80 Ry gates, 63
CNOT gates which are reduced to 28, 20, and 6 respectively, yielding an
overall 4x reduction in total number of gates and 3x reduction in the number
of parameters.

##  VI Conclusion and Outlook

The major outcome of our work is a working hybrid quantum-classical machine
learning model for generative learning and its bench-marking against purely
classical and quantum models. In addition, our models are particularly
suitable for quantum devices with small number of noisy qubits with limited
connectivity. With this algorithm, we have investigated the ability of quantum
generative models to learn the distributions obtained from the measurement of
quantum many-body systems. Such distributions are known to be demanding for
classical generative models and we verify the same in our experiments. We
further go on to show that our proposed hybrid-model can learn these
distributions with a much higher final fidelity. We find this trend to be
universal across a range of different types of quantum states, from generic
haar random states to dynamic kicked rotor states. We find that our model has
multiple hyper-parameters to tune and in one such case, when the latent size
is set to zero, we obtain at the Quantum circuit Born Machine, the most
popular pure quantum-mechanical generative model. In addition, we have shown
that QeVAEs can be useful for the practical task of circuit compilation.

Our work highlights the ability of quantum-enhanced models to perform better
than classical models within the variational autoencoder framework. There are
some drawbacks of the VAE approach and they include: (1) information loss in
the encoding and decoding process that often results in blurry outputs; (2)
posterior collapse: if the likelihood function is much more complex than the
posterior, then the encoder ignores the input data and outputs a trivial
latent space, leading the decoder to reconstruct the data from noise. Other
classical generative methods like GANs, Diffusion Models, and Normalizing
flows have gained much traction for producing better quality images. Although,
they also contain their fair share of disadvantages like mode collapse in GANs
and the large number of parameters in Diffusion models, they have been shown
to produce better quality images. In future, quantized versions of such models
can also be utilized to model measurement distributions.

The QeVAE algorithm contains many tunable parameters as mentioned in Table II.
In our calculations, we observed that obtaining an accurate final distribution
requires careful fine-tuning of these hyperparameters. For example, we found
out that initializing the preprocessing layer near $N(0,0.5)$ for our QeVAE is
very essential to outperform the classical model for learning measurement
distributions. In the absence of the preprocessing layer, the models quickly
over-fit the training data and occasionally performed better than classical
VAE. Moreover, we found that feature maps with entanglement (like the Pauli ZZ
map) produced distributions with the same accuracy but require fewer epochs.
The correct hyper-parameters have to be found out through an extensive grid-
search, and through semi-empirical means.

In the future, work can focus on verifying our results for learning the
classical distribution on diverse datasets. Particularly, text datasets can be
used for discrete distributions. In addition, one can also examine if
modifications in the ansatz like incorporating a non-linearity (mid-circuit
measurements) can enhance the performance of the QeVAE and provide a better
bound to the ELBO. Such hybrid quantum models can also be used to model the
distribution of observables used in particle physics and condensed matter
physics.

##  VII Acknowledgemnts

ASR would like to thank the KVPY fellowship and IISER Pune computing resources
for support. The corresponding code and datafiles are available on request.

## References

  * [1] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, “Quantum machine learning,” Nature, vol. 549, pp. 195–202, Sept. 2017\.
  * [2] S. Aaronson, “Read the fine print,” Nature Physics, vol. 11, pp. 291–293, Apr. 2015.
  * [3] P. W. Shor, “Algorithms for quantum computation: discrete logarithms and factoring,” in Proceedings 35th annual symposium on foundations of computer science, pp. 124–134, Ieee, 1994.
  * [4] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” Physical review letters, vol. 103, no. 15, p. 150502, 2009\.
  * [5] L. K. Grover, “A fast quantum mechanical algorithm for database search,” in Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, pp. 212–219, 1996.
  * [6] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, Y. Chen, T. Lillicrap, F. Hui, L. Sifre, G. van den Driessche, T. Graepel, and D. Hassabis, “Mastering the game of go without human knowledge,” Nature, vol. 550, pp. 354–359, Oct 2017.
  * [7] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera-Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, “Highly accurate protein structure prediction with AlphaFold,” Nature, vol. 596, pp. 583–589, July 2021.
  * [8] “Safe driving cars,” Nature Machine Intelligence, vol. 4, pp. 95–96, Feb. 2022.
  * [9] M. J. Bremner, R. Jozsa, and D. J. Shepherd, “Classical simulation of commuting quantum computations implies collapse of the polynomial hierarchy,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 467, pp. 459–472, Aug. 2010.
  * [10] B. Fefferman and C. Umans, “The power of quantum fourier sampling,” arXiv preprint arXiv:1507.05592, 2015.
  * [11] R. Sweke, J.-P. Seifert, D. Hangleiter, and J. Eisert, “On the Quantum versus Classical Learnability of Discrete Distributions,” Quantum, vol. 5, p. 417, Mar. 2021.
  * [12] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” Communications of the ACM, vol. 63, no. 11, pp. 139–144, 2020.
  * [13] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” stat, vol. 1050, p. 1, 2014.
  * [14] H.-Y. Huang, “Learning quantum states from their classical shadows,” Nature Reviews Physics, vol. 4, pp. 81–81, Jan. 2022.
  * [15] G. García-Pérez, M. A. Rossi, B. Sokolov, F. Tacchino, P. K. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco, “Learning to measure: Adaptive informationally complete generalized measurements for quantum algorithms,” PRX Quantum, vol. 2, p. 040342, Nov 2021.
  * [16] A. Rocchetto, S. Aaronson, S. Severini, G. Carvacho, D. Poderini, I. Agresti, M. Bentivegna, and F. Sciarrino, “Experimental learning of quantum states,” Science Advances, vol. 5, Mar. 2019.
  * [17] H.-Y. Huang, R. Kueng, and J. Preskill, “Predicting many properties of a quantum system from very few measurements,” Nature Physics, vol. 16, pp. 1050–1057, June 2020.
  * [18] A. Rocchetto, E. Grant, S. Strelchuk, G. Carleo, and S. Severini, “Learning hard quantum distributions with variational autoencoders,” npj Quantum Information, vol. 4, no. 1, pp. 1–7, 2018.
  * [19] M. Y. Niu, A. M. Dai, L. Li, A. Odena, Z. Zhao, V. Smelyanskyi, H. Neven, and S. Boixo, “Learnability and complexity of quantum samples,” arXiv preprint arXiv:2010.11983, 2020.
  * [20] Y. Wang, D. Blei, and J. P. Cunningham, “Posterior collapse and latent variable non-identifiability,” Advances in Neural Information Processing Systems, vol. 34, pp. 5443–5455, 2021.
  * [21] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, “beta-vae: Learning basic visual concepts with a constrained variational framework,” in International conference on learning representations, 2017.
  * [22] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [23] J. C. Spall, “An overview of the simultaneous perturbation method for efficient optimization,” Johns Hopkins apl technical digest, vol. 19, no. 4, p. 482–492, 1998.
  * [24] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, “Noisy intermediate-scale quantum algorithms,” Rev. Mod. Phys., vol. 94, p. 015004, Feb 2022.
  * [25] G. Carleo and M. Troyer, “Solving the quantum many-body problem with artificial neural networks,” Science, vol. 355, no. 6325, pp. 602–606, 2017.
  * [26] M. Hinsche, M. Ioannou, A. Nietner, J. Haferkamp, Y. Quek, D. Hangleiter, J.-P. Seifert, J. Eisert, and R. Sweke, “Learnability of the output distributions of local quantum circuits,” arXiv preprint arXiv:2110.05517, 2021.
  * [27] A. Khoshaman, W. Vinci, B. Denis, E. Andriyash, H. Sadeghi, and M. H. Amin, “Quantum variational autoencoder,” Quantum Science and Technology, vol. 4, no. 1, p. 014001, 2018.
  * [28] P. Rivas, L. Zhao, and J. Orduz, “Hybrid quantum variational autoencoders for representation learning,” in 2021 International Conference on Computational Science and Computational Intelligence (CSCI), pp. 52–57, 2021\.
  * [29] J. Li and S. Ghosh, “Scalable variational quantum circuits for autoencoder-based drug discovery,” in 2022 Design, Automation & Test in Europe Conference & Exhibition (DATE), pp. 340–345, IEEE, 2022.
  * [30] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, “Supervised learning with quantum-enhanced feature spaces,” Nature, vol. 567, no. 7747, pp. 209–212, 2019.
  * [31] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., “Pytorch: An imperative style, high-performance deep learning library,” Advances in neural information processing systems, vol. 32, 2019.
  * [32] H. Abraham, I. Adachi, G. Aguilar, A. Ahmed, J. Albert, G. Aleksandrowicz, T. Alexander, G. Alexandrowics, J. Aliaga, D. Allen, V. Altomare, C. Amico, A. Anastasiou, A. Annovi, R. Ansari, M. Antonelli, M. Antonello, D. Arbelaez, F. Arute, A. Asfaw, A. Azizkhanian, I. Babushkin, L. Baiardi, G. Balint-Kurti, S. Banerjee, P. Barkoutsos, G. Barron, S. Bartlett, L. Bassman, A. Behera, B. Behera, and P. . Behera, “Qiskit: An open-source framework for quantum computing,” 2019.
  * [33] C. A. Fuchs and C. M. Caves, “Ensemble-dependent bounds for accessible information in quantum mechanics,” Physical Review Letters, vol. 73, no. 23, p. 3047, 1994.
  * [34] M. Santhanam, S. Paul, and J. B. Kannan, “Quantum kicked rotor and its variants: Chaos, localization and beyond,” Physics Reports, vol. 956, pp. 1–87, 2022.  Quantum kicked rotor and its variants: Chaos, localization and beyond.


# Analysis of biologically plausible neuron models for regression with spiking
neural networks ††thanks: This work was supported by the Vannevar Bush Faculty
Fellowship award (GEK) from ONR (N00014-22-1-2795). The work of GEK is
supported by the U.S. Department of Energy, Advanced Scientific Computing
Research program, under the Scalable, Efficient and Accelerated Causal
Reasoning Operators, Graphs and Spikes for Earth and Embedded Systems (SEA-
CROGS) project, DE-SC0023191.

Mario De Florio  Division of Applied Mathematics
Brown University
Providence, RI
<EMAIL_ADDRESS>Adar Kahana  Division of Applied Mathematics
Brown University
Providence, RI
<EMAIL_ADDRESS>George Em Karniadakis  Division of Applied Mathematics
Brown University
Providence, RI
<EMAIL_ADDRESS>
###### Abstract

This paper explores the impact of biologically plausible neuron models on the
performance of Spiking Neural Networks (SNNs) for regression tasks. While SNNs
are widely recognized for classification tasks, their application to
Scientific Machine Learning and regression remains underexplored. We focus on
the membrane component of SNNs, comparing four neuron models: Leaky Integrate-
and-Fire, FitzHugh–Nagumo, Izhikevich, and Hodgkin-Huxley. We investigate
their effect on SNN accuracy and efficiency for function regression tasks, by
using Euler and Runge-Kutta 4th-order approximation schemes. We show how more
biologically plausible neuron models improve the accuracy of SNNs while
reducing the number of spikes in the system. The latter represents an
energetic gain on actual neuromorphic chips since it directly reflects the
amount of energy required for the computations.

###### Index Terms:

Spiking Neural Networks, LIF model, FitzHugh-Nagumo model, Izhikevich model,
Hodgkin-Huxley model, Regression, Scientific Machine Learning

##  I Introduction

The neuromorphic research has been gaining a lot of traction recently. The
main reasons are the potential efficiency, that unlocks faster training and
inference when using neural networks [1]. Another reason is the scalability.
With more neuromorphic hardware [2] and software [3] made available every day,
researchers are inclined to explore applying their methods developed on
traditional hardware to neuromorphic chips. Finally, the biological
plausibility of neuromorphic research is a backwind that drives this research
[4]. The core of neuromorphic research is currently Spiking Neural Networks
(SNNs). Several works [5, 7, 6, 8, 9] have shown the superiority of SNNs in
terms of efficiency, while also maintaining sufficient accuracy compared to
their Artificial Neural Network (ANN) counterparts. Most works so far that
utilize SNNs for Machine Learning (ML) focus on classification tasks. This is
due to the natural usage of neuromorphic hardware for classification, where
the inference is usually binary - $0$ if the sample belongs to the specific
class and $1$ otherwise (also referred to as one-hot classification). This
aligns with the neuronal model well - $0$ for no spike and $1$ for spike.

A less explored area of neuromorphic research is the use of SNNs for
Scientific Machine Learning (SciML). An example of a SciML task would be
approximating a continuous function using ML. These problems are often using
regression instead of classification, which is a challenge for SNNs since they
do not naturally fit the paradigm like in classification. However, to some
extent, regression is still possible using SNNs as presented in [10].
Following regression, we have shown that SNNs can be used for operator
regression, numerical schemes for dynamical systems [11], physics-informed
machine learning [12], and more [13, 14, 15]. There is still a lot to study,
revolving around neuromorphic algorithms for SciML, and this work presents a
fundamental step.

SNNs have three main components. The first is the input data, which is in the
form of spikes (binary zeros and ones). Data is encoded into spikes before
being input into the SNN, and outputs are often decoded back to at the end
(depending on the task). The second element is the synapse, which can be
thought of as a fully connected layer in the classical deep learning sense.
The synapses are trained using learning rules [16]. Last are the membranes,
which are analogous to the activation functions. An illustration of the
structure of an SNN is presented in Fig. 1 The membranes are the focus of this
work.

From the definition of the biological neuron [17], the membrane is a voltage
gate at the entrance of the neuronal cell. The membrane “decides” how much
voltage enters, and when a certain voltage threshold is crossed the neuron
spikes for a short period of time, and the membrane voltage resets to its
resting voltage. The first discovery of the membrane model was in 1907 by
Louis Lapicque [18], where he injected a frog’s leg with voltage until it
twitched, inventing the Leaky-Integrate and Fire (LIF) neuron model that is
commonly used in many studies. This model is famous for its simplicity, and
also biological accuracy. Many neuroscience studies [19, 20, 21, 22] use this
model, and all SNN software and hardware simulators we are aware of have
implementations of a LIF model. The LIF neuron is modeled using a Resistor-
Capacitor circuit (hence the term neuron model), as explained in detail in the
following section.

Figure 1: Sketch of an SNN architecture. The spiking data is used as inputs to
Hodgkin-Huxley neurons (membranes) and then fed into a synapse. In this
architecture, we have a membrane potential that creates output spikes,
followed by a synapse, and another membrane that generates output spikes,
followed by a second synapse. Finally, a classification or regression
mechanism should be applied. In this study, we employ only one membrane
followed by a single synapse to perform regression tasks, as elaborated in the
results section.

Several neuroscience-related studies have been conducted to develop more
accurate membrane formulations. The FitzHugh-Nagumo model [23, 24] introduces
a set of two coupled Ordinary Differential Equations (ODEs) to better capture
the oscillatory behavior of the neuronal cell. The Izhikevich model [25] is
also a set of ODEs, which captures better the spiking patterns. The Hodgkin-
Huxley model [26] is a set of four non-linear coupled ODEs that goes deeper
into the transfer of potassium and sodium in and out of the neuronal cell.
While the aforementioned three models and LIF are explored in this work, we
note that several other models have been proposed [27, 28, 29, 30]. The trade-
off when estimating the eligibility of the models for SciML is between the
accuracy of the overall process as described in the Results section and the
complexity of the model, which directly impacts its efficiency. These were the
guidelines for choosing the models to compare in this work, where LIF is the
simplest and most efficient and Hodgkin-Huxley on the other accuracy end of
the spectrum. The intermediate ones are, in addition to their qualities, very
popular amongst this scientific community.

The main question we explore in this work is: how do the different membrane
models contribute to the success of an SNN? As mentioned, SNNs are measured by
their overall accuracy, as well as their efficiency. We first develop a
generic implementation of a SNN with a synapse and a membrane, and simply
alternate between the different membrane models. To do this, we implemented
all neuron models and fine-tuned the parameters to achieve the desired
results. The findings we report in this manuscript bring valuable insights for
those experimenting with SNNs. In addition, we find that more complex models,
although taking longer to compute, reduce the total number of spikes
registered, which means lower overall energy (higher efficiency).

##  II Biological Neuron Models

###  II-A Leaky Integrate-and-Fire model

The Leaky Integrate-and-Fire (LIF) model is a simplified neuronal model widely
used for its computational efficiency and ease of analysis. It provides a
basic representation of neuronal dynamics, capturing the essential features of
membrane potential integration and action potential generation. The earliest
form of the Integrate-and-Fire model was originally introduced by Lapicque in
1907 [18], and subsequently, the LIF model was developed after adding an
attenuation term [31].
The dynamics of the LIF model is governed by a single ODE that describes the
evolution of the membrane potential $V(t)$ over time:

$\dfrac{dV(t)}{dt}=\frac{V(t)-V_{rest}}{RC}+I(t),$ (1)

where $C$ is the membrane capacitance, $R$ is the membrane resistance, $I(t)$
is the input current, and $V_{rest}$ is the membrane resting voltage. In our
simulation, $RC=0.2$ and $I_{c}(t)$ are input currents of magnitude $1$.
The LIF model is commonly used to study the response of a neuron to input
currents. A spike is generated when the membrane potential reaches a threshold
value $V_{\text{th}}=1$. Upon reaching this threshold, the membrane potential
is reset to a resting value $V_{rest}$. The thresholding mechanism is
expressed as follows:

$\text{if }V\geq V_{th},\text{ then spike and }V\rightarrow V_{rest}$ (2)

###  II-B FitzHugh–Nagumo model

The FitzHugh–Nagumo (FHN) model is a simplified two-variable system of ODEs
that was introduced by Richard FitzHugh in 1961 [23] and later refined by J.
Nagumo et al. in 1962 [24]. The FHN model serves as a mathematical abstraction
to capture the characteristics of neuronal excitability and oscillatory
behavior, and it has the following form

$\begin{cases}\dfrac{dV}{dt}=V-\dfrac{V^{3}}{3}-W+I\\\
\dfrac{dW}{dt}=\dfrac{1}{\gamma}(V+\alpha-\beta
W)\end{cases}\;\text{s.t.}\qquad\begin{cases}V(0)=0\\\ W(0)=0\end{cases}$ (3)

where $V$ is the membrane potential, $W$ is a recovery variable, and $I$ is
the input current. The parameters governing the system’s behavior are the
bifurcation parameter $\alpha=0.7$ determining the onset of excitability,
$\beta=0.8$ which controls the recovery variable’s sensitivity to the membrane
potential, and $\gamma=12.5$ governs the timescale of the recovery variable. A
spike is generated and the membrane resets to the resting value 0 once the
threshold value of 1 is reached.

###  II-C Izhikevich model

The Izhikevich model, introduced by Eugene M. Izhikevich in 2003 [25], is a
simple yet powerful two-variable neuronal model that captures a wide range of
spiking behaviors observed in real neurons. This model is particularly
valuable for its computational efficiency and ability to replicate diverse
neural firing patterns.
The dynamics of the Izhikevich model are governed by the following set of
equations:

$\begin{cases}\dfrac{dv}{dt}=0.04v^{2}+5v+140-u+I\\\
\dfrac{du}{dt}=a(bv-u)\end{cases}\;\text{s.t.}\quad\begin{cases}v(0)=-70mV\\\
u(0)=-14\end{cases}$ (4)

with the auxiliary condition

$\text{if }v\geq v_{th},\text{ then spike and }\begin{cases}v\rightarrow c\\\
u\rightarrow u+d\end{cases}$ (5)

where $v$ represents the membrane potential, $u$ is a recovery variable, and
$I=10mA$ is the input current. The parameters shaping the spiking dynamics are
$a=0.02\frac{1}{ms}$, which controls the time scale of the recovery variable,
$b=0.2$ that influences the sensitivity of the recovery variable to the
membrane potential, $c=-50mV$, which is the after-spike reset value of the
membrane potential caused by the fast high-threshold $K^{+}$ conductances, and
$d=2mV$, which is the after-spike reset of the recovery variable caused by
slow high-threshold $Na^{+}$ and $K^{+}$ conductances.

###  II-D Hodgkin-Huxley model

The Hodgkin-Huxley model, proposed in 1952 by British physiologists Alan Lloyd
Hodgkin and Andrew Huxley [26], revolutionized our understanding of the
mechanisms underlying the generation and propagation of action potentials in
neurons. This mathematical model provides a framework for describing the
complex dynamics of ion channels and membrane potentials in excitable cells,
particularly neurons. The activation of the membrane potential is described by
the following system of four non-linear Ordinary Differential Equations (ODEs)
subject to certain initial conditions [32]:

$\begin{cases}C_{m}\dfrac{dV_{m}}{dt}=-g_{l}(V_{m}-e_{l})-g_{k}n^{4}(V_{m}-e_{k})\\\
\qquad\qquad-g_{Na}m^{3}h(V_{m}-e_{Na})+I\\\
\dfrac{dn}{dt}=\alpha_{n}(V_{m})(1-n)-\beta_{n}(V_{m})n\\\
\dfrac{dm}{dt}=\alpha_{m}(V_{m})(1-m)-\beta_{m}(V_{m})m\\\
\dfrac{dh}{dt}=\alpha_{h}(V_{m})(1-h)-\beta_{h}(V_{m})h\\\ \end{cases}$ (6)

subject to $V_{m}(0)=-65mV,n(0)=0.3177,m(0)=0.0529$, and $h(0)=0.5960$. Each
ODE describes the time evolution of the four variables of the HH model. In
particular, $V_{m}$ is the membrane potential, $n$ and $m$ are the activation
gating variables for potassium and sodium channels, respectively, and $h$ is
the inactivation gating variable for the sodium channels. $C_{m}$ represents
the membrane capacitance, $g_{l}$, $g_{K}$ and $g_{Na}$ are the leak,
potassium and sodium conductances, respectively, $e_{l}$, $e_{k}$ and $e_{Na}$
are their reversal potentials, and $I$ is the input current per unit area. The
$\alpha$ and $\beta$ functions represent the rate constants for the ion
channels, and they can be expressed as

$\displaystyle\alpha_{n}(V_{m})=\dfrac{0.1-0.01U_{m0}}{e^{1-0.1U_{m0}}-1}\qquad$
$\displaystyle\beta_{n}(V_{m})=0.125e^{-\dfrac{U_{m0}}{80}}$
$\displaystyle\alpha_{m}(V_{m})=\dfrac{2.5-0.1U_{m0}}{e^{2.5-0.1U_{m0}}-1}\qquad$
$\displaystyle\beta_{m}(V_{m})=4e^{-\dfrac{U_{m0}}{18}}$
$\displaystyle\alpha_{h}(V_{m})=0.07e^{-\dfrac{U_{m0}}{20}}\qquad$
$\displaystyle\beta_{h}(V_{m})=\dfrac{1}{1+e^{3-\dfrac{U_{m0}}{10}}}$

where $U_{m0}=V_{m}-V_{m}(0)$, with $V_{m}(0)=-65mV$, and the values of the
parameters used in our simulations, assuming $1cm^{2}$ of membrane, are
$C_{m}=1\mu F$, $g_{l}=0.3mS$, $g_{k}=36mS$, $g_{Na}=120mS$, $e_{k}=-77mV$,
$e_{l}=-54mV$, $e_{Na}=50mV$. In our simulations, the input current $I$ is
represented by spikes of $5mA$ of intensity. The membrane has a threshold
parameter set at $30mV$. When the membrane potential reaches the threshold, a
spike is registered in the neuron, and the membrane resets to the value of
$V_{m}(0)$. An example of this membrane’s behavior is represented in Figure 2.

Figure 2: Example of the membrane potential over time for given input spikes.
The reset mechanism is set to a resting position without refraction. The top
plot presents the input spikes over time. The middle plot presents the
membrane activity modeled by the Hodgkin-Huxley model. Notice that when the
membrane potential reaches the threshold, an output spike is registered, and
the membrane resets. The bottom plot presents the registered output spikes
over time.

##  III Numerical methods

###  III-A Euler method

To simplify the calculations, we can employ the forward Euler finite
difference method to directly approximate the solution of the Hodgkin-Huxley
model [33] from equations (6), as follows:

$\displaystyle V_{m}(t+\Delta t)$ $\displaystyle=V_{m}(t)+\Delta
t\bigl{(}I(t)-g_{K}n^{4}(V_{m}(t)-e_{K})$ (7) $\displaystyle-g_{Na}m^{3}$
$\displaystyle h(V_{m}(t)-e_{Na})-g_{L}(V_{m}(t)-e_{L})\bigr{)}/C_{m}$
$\displaystyle n(t+\Delta t)$ $\displaystyle=n(t)+\Delta
t(\alpha_{n}(1-n(t))-\beta_{n}n(t))$ $\displaystyle m(t+\Delta t)$
$\displaystyle=m(t)+\Delta t(\alpha_{m}(1-m(t))-\beta_{m}m(t))$ $\displaystyle
h(t+\Delta t)$ $\displaystyle=h(t)+\Delta t(\alpha_{h}(1-h(t))-\beta_{h}h(t))$

where $\Delta t$ is set equal to 0.1. This serves as a simplified
approximation to the solution of the Hodgkin-Huxley model, and it’s
anticipated that the accuracy in computing the membrane potential will be
reduced. Nevertheless, this approach offers a more lightweight implementation,
enabling more efficient computations.

###  III-B Runge-Kutta 4th order method

In contrast to the Euler method, the Runge-Kutta 4th Order (RK4) numerical
integration scheme provides a more accurate approximation of the solution to
the Hodgkin-Huxley model [34]. The RK4 method involves multiple steps to
iteratively update the variables, providing enhanced accuracy in capturing the
intricate dynamics of the membrane potential and gating variables.
The slopes $k_{1},k_{2},k_{3}$, and $k_{4}$ for all the HH variables are

$\displaystyle k_{1V}$ $\displaystyle=\Delta
t\bigl{(}I-g_{K}n^{4}(V_{m}-e_{K})-g_{Na}m^{3}h(V_{m}-e_{Na})$ (8)
$\displaystyle-g_{L}(V_{m}-e_{L})\bigr{)}/C_{m}$ $\displaystyle k_{1n}$
$\displaystyle=\Delta t(\alpha_{n}(1-n)-\beta_{n}n)$ $\displaystyle k_{1m}$
$\displaystyle=\Delta t(\alpha_{m}(1-m)-\beta_{m}m)$ $\displaystyle k_{1h}$
$\displaystyle=\Delta t(\alpha_{h}(1-h)-\beta_{h}h)$

$\displaystyle k_{2V}=\Delta
t\bigl{(}I-g_{K}(n+0.5k_{1n})^{4}(V_{m}+0.5k_{1V}-e_{K})$ (9) $\displaystyle-
g_{Na}(m+0.5k_{1m})^{3}(h+0.5k_{1h})(V_{m}+0.5k_{1V}-e_{Na})$
$\displaystyle\qquad-g_{L}(V_{m}+0.5k_{1V}-e_{L})\bigr{)}/C_{m}$
$\displaystyle k_{2n}=\Delta
t(\alpha_{n}(1-(n+0.5k_{1n}))-\beta_{n}(n+0.5k_{1n}))$ $\displaystyle
k_{2m}=\Delta t(\alpha_{m}(1-(m+0.5k_{1m}))-\beta_{m}(m+0.5k_{1m}))$
$\displaystyle k_{2h}=\Delta
t(\alpha_{h}(1-(h+0.5k_{1h}))-\beta_{h}(h+0.5k_{1h}))$

$\displaystyle k_{3V}=\Delta
t\bigl{(}I-g_{K}(n+0.5k_{2n})^{4}(V_{m}+0.5k_{2V}-e_{K})$ (10) $\displaystyle-
g_{Na}(m+0.5k_{2m})^{3}(h+0.5k_{2h})(V_{m}+0.5k_{2V}-e_{Na})$
$\displaystyle\qquad-g_{L}(V_{m}+0.5k_{2V}-e_{L})\bigr{)}/C_{m}$
$\displaystyle k_{3n}=\Delta
t(\alpha_{n}(1-(n+0.5k_{2n}))-\beta_{n}(n+0.5k_{2n}))$ $\displaystyle
k_{3m}=\Delta t(\alpha_{m}(1-(m+0.5k_{2m}))-\beta_{m}(m+0.5k_{2m}))$
$\displaystyle k_{3h}=\Delta
t(\alpha_{h}(1-(h+0.5k_{2h}))-\beta_{h}(h+0.5k_{2h}))$

$\displaystyle k_{4V}=\Delta
t\bigl{(}I-g_{K}(n+k_{3n})^{4}(V_{m}+k_{3V}-e_{K})$ (11) $\displaystyle\qquad-
g_{Na}(m+k_{3m})^{3}(h+k_{3h})(V_{m}+k_{3V}-e_{Na})$ $\displaystyle\qquad-
g_{L}(V_{m}+k_{3V}-e_{L})\bigr{)}/C_{m}$ $\displaystyle k_{4n}=\Delta
t(\alpha_{n}(1-(n+k_{3n}))-\beta_{n}(n+k_{3n}))$ $\displaystyle k_{4m}=\Delta
t(\alpha_{m}(1-(m+k_{3m}))-\beta_{m}(m+k_{3m}))$ $\displaystyle k_{4h}=\Delta
t(\alpha_{h}(1-(h+k_{3h}))-\beta_{h}(h+k_{3h}))$

Thus, the updated equations for the RK4 method applied to the Hodgkin-Huxley
model are as follows:

$\displaystyle V_{m}(t+\Delta
t)=V_{m}+\frac{1}{6}(k_{1V}+2k_{2V}+2k_{3V}+k_{4V})$ (12) $\displaystyle
n(t+\Delta t)=n+\frac{1}{6}(k_{1n}+2k_{2n}+2k_{3n}+k_{4n})$ $\displaystyle
m(t+\Delta t)=m+\frac{1}{6}(k_{1m}+2k_{2m}+2k_{3m}+k_{4m})$ $\displaystyle
h(t+\Delta t)=h+\frac{1}{6}(k_{1h}+2k_{2h}+2k_{3h}+k_{4h})$

Here, $\Delta t$ is the time step, and $V_{m},n,m,$ and $h$ represent the
membrane potential and gating variables. The RK4 method ensures a more
accurate numerical approximation of the Hodgkin-Huxley model’s dynamics
compared to the Euler method, making it suitable for simulations that require
higher precision.

##  IV Results and Discussions

To test the performance of the different neural models computed with the Euler
and RK4 methods we follow the experimental setup presented in [10] for
function regression. We use the same codes, changing only the activation from
LIF neurons into the various ones presented in this work. Since we do not have
access to a neuromorphic chip for testing the method with spikes, we are using
a GPU-based training with TensorFlow 2 [35]. The energetic gains are measured
by the number of spikes in the system, rather than execution time. This is a
valid measure, as on an actual chip the number of spikes directly reflects the
amount of energy required for the computations.

Figure 3: Comparison of function regressions using LIF model and Hodgkin-Huxley model with Euler method in the SNN implementation for three different functions, with and without noise. TABLE I: Relative $L_{2}$ error and number of output spikes for neuron membrane models and solvers, for both noiseless and noisy functions. The encoded neuron models are Leaky Integrate-and-Fire (LIF), FitzHugh–Nagumo (FHN), Izhikevich (IZH), and Hodgkin-Huxley (HH). |  | Relative $L_{2}$ error (# of output spikes)
---|---|---
Neuron model |  Numerical method | Noiseless | Noisy |  Time (s)
Discontinuity | Square | Sine | Discontinuity | Square | Sine
LIF | Euler | 4.50e-03 | 1.30e-02 | 1.43e-02 | 4.30e-03 | 1.18e-02 | 1.52e-02 | 12
(2271) | (2242) | (2312) | (2257) | (2212) | (2225)
RK | 1.864e-03 | 6.13e-03 | 5.82-03 | 1.66e-03 | 4.78e-03 | 5.53e-03 | 22
(904) | (886) | (896) | (908) | (902) | (897)
FHN | Euler | 6.63e-04 | 1.96e-03 | 1.87e-03 | 5.84e-04 | 2.03e-03 | 2.19e-03 | 16
(281) | (269) | (271) | (266) | (276) | (272)
RK | 6.16e-04 | 5.49e-04 | 6.10e-04 | 6.24e-04 | 5.02e-04 | 6.40e-04 | 48
(55) | (53) | (51) | (58) | (54) | (56)
IZH | Euler | 5.78e-04 | 2.20e-03 | 2.05e-03 | 6.58e-04 | 1.97e-03 | 1.91e-03 | 20
(219) | (222) | (217) | (216) | (218) | (221)
RK | 5.16e-04 | 5.67e-04 | 5.89e-04 | 5.52e-04 | 6.11e-04 | 7.32e-04 | 55
(58) | (62) | (64) | (65) | (67) | (66)
HH | Euler | 4.17e-04 | 1.26e-03 | 1.24e-03 | 4.61e-04 | 1.20e-03 | 1.37e-03 | 48
(118) | (114) | (122) | (113) | (109) | (121)
RK | 1.69e-03 | 1.02e-03 | 8.95e-04 | 9.20e-04 | 8.07e-04 | 9.04e-04 | 130
(71) | (71) | (62) | (74) | (74) | (86)

The SNN architecture we use is a membrane layer, followed by a synapse. Since
we are not using a neuromorphic chip, we use the back-propagation algorithm to
train the synapse and learn the weights of the network. The shallow network we
choose here is sufficient for demonstrating the advantages of using different
membrane models and analyzing the results. We point out that using a more
complex network architecture will likely yield better performance. In
addition, with this specific architecture, the back-propagation algorithm
operates only on the last layer, which is considered more “biologically
plausible”. The input data are encoded using $N_{t}=150$ time steps and
$N_{x}=100$ collocation points for spatial discretization. The outputs are
floating numbers (no decoding required). The three functions we attempt for
regression, with and without noise, are:

  * •

$Square$

$y(x)=x^{2}$ (13)

  * •

$Discontinuity$

$y(x)=\begin{cases}1\qquad\text{for}\qquad x\leq 0\\\ 2\qquad\text{for}\qquad
x>0\end{cases}$ (14)

  * •

$Sine$

$y(x)=\frac{\sin(kx)}{k^{2}}\qquad\text{with}\qquad k=1.2$ (15)

The square function is more of a sanity check. The discontinuous function is
challenging, since it has a sharp jump (not differentiable), and the sine
oscillates around 0 which often causes issues. We also consider the cases in
which the functions are perturbed with 0.1 additive normally distributed
noise. The performance of the regression is reported in terms of the
$L_{2}$-norm difference between the predicted and exact solution as

$L_{2}\text{-norm}=\sum_{j=1}^{N_{x}}\left(y_{pred}(x_{j})-y_{exact}(x_{j})\right)^{2},$
(16)

where $N_{x}$ is the number of collocation points of the spatial
discretization. In Figure 3, a qualitative visualization of the regressions
made by using the Euler method for the LIF and HH models is reported for the
three functions, with and without noise. One can see how with HH there is a
better overlap with the exact solutions. To quantitatively evaluate the
different performances in terms of the $L_{2}$-norm, efficiency (number of
spikes produced), and computational times we can refer to Table I. As
expected, the computational time increases from Euler to RK4, and from 1-ODE
model (LIF) to 4-ODE model (HH), with a decrease in the $L_{2}$-norm error.
One can see that in terms of accuracy, IZH and HH outperform LIF and FHN
neuron models for all the test cases, but with a significant cost in
computational time. It is interesting to notice how a greater accuracy
corresponds to a lower number of output spikes, which results in a more
efficient cost for the neuromorphic hardware in regression tasks.

## References

  * [1] D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier,  “Physics for neuromorphic computing,”  _Nature Reviews Physics_ , vol. 2, no. 9, pp. 499–510, 2020,  Nature Publishing Group UK London.
  * [2] J. Schemmel, D. Brüderle, A. Grübl, M. Hock, K. Meier, and S. Millner,  “A wafer-scale neuromorphic hardware system for large-scale neural modeling,”  in _2010 IEEE International Symposium on Circuits and Systems (ISCAS)_ , pp. 1947–1950, 2010,  IEEE.
  * [3] T. Titirsha, S. Song, A. Balaji, and A. Das,  “On the role of system software in energy management of neuromorphic computing,”  in _Proceedings of the 18th ACM International Conference on Computing Frontiers_ , pp. 124–132, 2021.
  * [4] J.-Q. Yang, R. Wang, Y. Ren, J.-Y. Mao, Z.-P. Wang, Y. Zhou, and S.-T. Han,  “Neuromorphic engineering: from biological to spike-based hardware nervous systems,”  _Advanced Materials_ , vol. 32, no. 52, p. 2003610, 2020,  Wiley Online Library.
  * [5] B. Han and K. Roy,  “Deep spiking neural network: Energy efficiency through time based coding,”  in _European Conference on Computer Vision_ , pp. 388–404, 2020,  Springer.
  * [6] S. Misra, C. Allemang, C. Arose, B. Taylor, A. Dubovskiy, A. Sidi El Valli, L. Rehm, A. Haas, A. Kent, L. Bland, S. Cardwell, D. Smith and J.B. Aimone,  _How to use stochastic devices in probabilistic calculations_ ,  Bulletin of the American Physical Society, 2024.
  * [7] E. Lemaire, L. Cordone, A. Castagnetti, P.-E. Novac, J. Courtois, and B. Miramond,  “An analytical estimation of spiking neural networks energy efficiency,”  in _International Conference on Neural Information Processing_ , pp. 574–587, 2022,  Springer.
  * [8] R. V. W. Putra and M. Shafique,  “Fspinn: An optimization framework for memory-efficient and energy-efficient spiking neural networks,”  _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 39, no. 11, pp. 3601–3613, 2020,  IEEE.
  * [9] S. Park, S. Kim, H. Choe, and S. Yoon,  “Fast and efficient information transmission with burst spikes in deep spiking neural networks,”  in _Proceedings of the 56th Annual Design Automation Conference 2019_ , pp. 1–6, 2019.
  * [10] A. Kahana, Q. Zhang, L. Gleyzer, and G. E. Karniadakis,  “Function regression using spiking deeponet,”  _arXiv preprint arXiv:2205.10130_ , 2022.
  * [11] S. O. Fatunla,  _Numerical methods for initial value problems in ordinary differential equations_ ,  Academic Press, 2014.
  * [12] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang,  “Physics-informed machine learning,”  _Nature Reviews Physics_ , vol. 3, no. 6, pp. 422–440, 2021,  Nature Publishing Group UK London.
  * [13] Q. Yu, H. Tang, K. C. Tan, and H. Yu,  “A brain-inspired spiking neural network model with temporal encoding and learning,”  _Neurocomputing_ , vol. 138, pp. 3–13, 2014,  Elsevier.
  * [14] H. Mostafa,  “Supervised learning based on temporal coding in spiking neural networks,”  _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 29, no. 7, pp. 3227–3235, 2017,  IEEE.
  * [15] I. M. Comsa, K. Potempa, L. Versari, T. Fischbacher, A. Gesmundo, and J. Alakuijala,  “Temporal coding in spiking neural networks with alpha synaptic function,”  in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pp. 8529–8533, 2020,  IEEE.
  * [16] Y. Bengio, S. Bengio, and J. Cloutier,  _Learning a synaptic learning rule_ ,  Citeseer, 1990.
  * [17] I. B. Levitan and L. K. Kaczmarek,  _The neuron: Cell and molecular biology_ ,  Oxford University Press, USA, 2015.
  * [18] L. M. Lapicque, “Recherches quantitatives sur l’excitation electrique des nerfs,” J Physiol Paris, vol. 9, pp. 620–635, 1907.
  * [19] C. Laing and G. J. Lord,  _Stochastic methods in neuroscience_ ,  Oxford University Press, 2009.
  * [20] S. Ostojic and N. Brunel,  “From spiking neuron models to linear-nonlinear models,”  _PLoS computational biology_ , vol. 7, no. 1, p. e1001056, 2011,  Public Library of Science San Francisco, USA.
  * [21] S. Krasovskaya and W. J. MacInnes,  “Salience models: A computational cognitive neuroscience review,”  _Vision_ , vol. 3, no. 4, p. 56, 2019,  MDPI.
  * [22] R. E. Kass, S.-I. Amari, K. Arai, E. N. Brown, C. O. Diekman, M. Diesmann, B. Doiron, U. T. Eden, A. L. Fairhall, G. M. Fiddyment, and others,  “Computational neuroscience: Mathematical and statistical perspectives,”  _Annual review of statistics and its application_ , vol. 5, pp. 183–214, 2018,  Annual Reviews.
  * [23] R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical journal, vol. 1, no. 6, pp. 445–466, 1961, Elsevier.
  * [24] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” Proceedings of the IRE, vol. 50, no. 10, pp. 2061–2070, 1962, IEEE.
  * [25] E. M. Izhikevich, “Simple model of spiking neurons,” IEEE Transactions on neural networks, vol. 14, no. 6, pp. 1569–1572, 2003, IEEE.
  * [26] A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” The Journal of physiology, vol. 117, no. 4, p. 500, 1952, Wiley-Blackwell.
  * [27] P. F. Pinsky and J. Rinzel,  “Intrinsic and network rhythmogenesis in a reduced Traub model for CA3 neurons,”  _Journal of computational neuroscience_ , vol. 1, pp. 39–60, 1994,  Springer.
  * [28] S. Hayman,  “The McCulloch-Pitts Model,”  in _IJCNN’99. International Joint Conference on Neural Networks. Proceedings (Cat. No. 99CH36339)_ , vol. 6, pp. 4438–4439, 1999,  IEEE.
  * [29] H. Lecar,  “Morris-Lecar model,”  _Scholarpedia_ , vol. 2, no. 10, p. 1333, 2007.
  * [30] M. Storace, D. Linaro, and E. de Lange,  “The Hindmarsh–Rose neuron model: Bifurcation analysis and piecewise-linear approximations,”  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 18, no. 3, 2008,  AIP Publishing.
  * [31] L. F. Abbott, “Lapicque’s introduction of the integrate-and-fire model neuron (1907),” Brain research bulletin, vol. 50, no. 5-6, pp. 303–304, 1999, Citeseer.
  * [32] W. Teka, D. Stockton, and F. Santamaria, “Power-law dynamics of membrane conductances increase spiking diversity in a Hodgkin-Huxley model,” PLoS computational biology, vol. 12, no. 3, p. e1004776, 2016, Public Library of Science San Francisco, CA USA.
  * [33] M. Hanslien, K. H. Karlsen, and A. Tveito, “A maximum principle for an explicit finite difference scheme approximating the Hodgkin-Huxley model,” BIT Numerical Mathematics, vol. 45, pp. 725–741, 2005, Springer.
  * [34] J. W. Moore and F. Ramon, “On numerical integration of the Hodgkin and Huxley equations for a membrane action potential,” Journal of theoretical Biology, vol. 45, no. 1, pp. 249–273, 1974, Elsevier.
  * [35] S. Raschka and V. Mirjalili,  _Python machine learning: Machine learning and deep learning with Python, scikit-learn, and TensorFlow 2_ ,  Packt Publishing Ltd, 2019.