# Regressions on quantum neural networks at maximal expressivity

I. Panadero1,2,3, Y. Ban1,4, H. Espinós1, R. Puebla1, J. Casanova3,5,6 E.
Torrontegui1 1 Departamento de Física, Universidad Carlos III de Madrid, Avda.
de la Universidad 30, 28911 Leganés, Spain  2 Arquimea Research Center, Camino
las Mantecas s/n, 38320 Santa Cruz de Tenerife, Spain  3 Department of
Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644,
48080 Bilbao, Spain  4 TECNALIA, Basque Research and Technology Alliance
(BRTA), 48160 Derio, Spain  5 EHU Quantum Center, University of the Basque
Country UPV/EHU, Leioa, Spain  6 IKERBASQUE, Basque Foundation for Science,
Plaza Euskadi 5, 48009 Bilbao, Spain<EMAIL_ADDRESS>/
<EMAIL_ADDRESS>
###### Abstract

We analyze the expressivity of a universal deep neural network that can be
organized as a series of nested qubit rotations, accomplished by adjustable
data re-uploads. While the maximal expressive power increases with the depth
of the network and the number of qubits, it is fundamentally bounded by the
data encoding mechanism. Focusing on regression problems, we systematically
investigate the expressivity limits for different measurements and
architectures. The presence of entanglement, either by entangling layers or
global measurements, saturate towards this bound. In these cases, entanglement
leads to an enhancement of the approximation capabilities of the network
compared to local readouts of the individual qubits in non-entangling
networks. We attribute this enhancement to a larger survival set of Fourier
harmonics when decomposing the output signal.

Keywords: quantum neural networks, quantum machine learning.

##  1 Introduction

The convergence of machine learning (ML) and quantum computing leads to the
birth of quantum machine learning [1, 2, 3], which exploits both information
processing paradigms to provide better algorithms. This is achieved either by
using ML to improve the efficiency of quantum protocols [4, 5, 6, 7, 8] or by
taking advantage of quantum resources to create new paradigms with better
performance [9, 10, 11, 12, 13], such as quantum neural networks (QNN). While
there are several studies comparing the expresssivity of classical and quantum
artificial neural networks [14, 15], the goal of beating classical ML is
undergoing critical revision [16]. Recent works, in the quest of finding
better quantum models for ML purposes, characterize and compare the
expressivity among different networks typologies [17, 18, 19, 20, 21]. This is
made possible by the emergence of unitary-based QNNs [22, 23, 24, 25, 26, 27,
28], which offer experimental feasibility across multiple quantum platforms,
such as trapped ions [29, 30] and superconducting circuits [25, 31, 32, 33].

On the one hand, a perceptron-type QNN, that can be implemented as a quasi-
adiabatic passage in an Ising model [24, 27, 28], expresses any continuous
function as a weighted linear combination of different adjustable sigmoids
based on the universal approximation theorem [34]. Correlations resulting from
the entanglement among various perceptrons provide these artificial networks
with at least as much power as their classical counterparts. However, the
entanglement is intrinsically generated along hidden layers, and understanding
its role in improving the network approximation capability in classification
and regression problems remains a challenge. On the other hand, QNNs based on
parametric quantum circuits, by virtue of the Solovay–Kitaev theorem [35], are
universal classifiers, and are implemented in repeat-until-success [22],
random circuits [23, 36], or by data re-uploading [25]. A single qubit has
already been employed to approximate continuous functions, specifically in the
context of regression problems [37]. Quantum resources may be introduced in an
easily controllable manner by adding entangling gates into the topology of the
network.

While the expressivity of all these unitary-based QNN models is fundamentally
bounded by the data encoding strategy [38], we show that the maximal
performance in regression tasks is conditioned by the final QNN state readout.
To this end, we assess the expressivity of a universal deep QNN that follows a
data re-uploading strategy, employing various measurements and architectures
that vary in terms of network depth and number of qubits. This construction
mitigates the discrepancies attributed to the use of different models,
enabling a fair comparison of networks with the same topology, but
distinguished by the presence of quantum resources, such as entanglement. We
make use of a partial Fourier series representation of the QNN outputs to
characterize and bound the maximal expression power of the network. In
particular, we find that the maximal performance is conditioned by the final
QNN state readout. With the aid of the teacher-student scheme [19], we
systematically benchmark the network accuracy in regression problems, showing
that global measurements of the multiqubit system or the use of a final
entangling layer leads to the saturation of the maximal expressive bound,
optimizing the performance of the network.

The article is structured as follows. We begin by introducing a universal
quantum regressor based on re-uploading continuous input data, which results
in a trainable QNN. Then, we quantify the QNN’s expressivity through partial
Fourier decomposition and numerically benchmark several QNN distinguished by
the presence of entanglement. The manuscript ends with the main conclusions
and prospective research lines.

##  2 Deep data re-uploading QNN

We propose an encoding scheme of arbitrary normalized real functions
$f(\vec{x})$ $\in$ $[0,1]$ through the degrees of freedom of a set of qubits.
These qubits are assumed to constitute a circuit that performs a unitary
transformation of the initial state $\ket{\phi_{0}}$ into the final output
state

$\ket{\phi(\vec{x})}=\hat{\mathcal{U}}(\vec{x})\ket{\phi_{0}},$ (1)

by the unitary transformer $\hat{\mathcal{U}}(\vec{x})$, where
$\vec{x}=(x_{1},x_{2},\dots,x_{n})\in\mathbb{R}^{n}$ are independent input
variables that set the topology of the circuit. The circuit is approximated by
a data re-uploading QNN with an adjustable set of parameters
$\vec{\theta}=\\{\vec{\theta}_{1},\vec{\theta}_{2},\cdots\vec{\theta}_{L}\\}$,

$\hat{\mathcal{U}}(\vec{x})\simeq\hat{U}(\vec{x},\vec{\theta})=\prod_{l=1}^{L}\hat{\mathcal{L}}_{l}(\vec{x},\vec{\theta}_{l}),$
(2)

where we adopt the convention $\prod_{n=1}^{N}A_{n}=A_{N}\cdot A_{N-1}...\cdot
A_{1}$, with the dot being the standard matrix multiplication. The $l$-th
layer is represented by the unitary transformation
$\hat{\mathcal{L}}_{l}(\vec{x},\vec{\theta}_{l})$, with input data $\vec{x}$
and undetermined parameters $\vec{\theta}_{l}$. Similar to perceptron-based
neural networks [24], the complexity (or interconectivity) of the network
increases with the number of layers $L$. We choose the unitary transformation
for each layer to be constituted by a single qubit SU(2) rotation with the
data encoded as the rotation angle,

$\hat{\mathcal{L}}_{l}(\vec{x},\vec{\theta}_{l})=\hat{R}_{x}(2(\vec{\omega}_{l}\cdot\vec{x}+\beta_{l}))\hat{R}_{y}(2\alpha_{l}),$
(3)

where $\hat{R}_{\gamma}(2\tau)=e^{i\tau\hat{\sigma}_{\gamma}}$ is the rotation
around the $\gamma-$axis, with $\gamma\equiv\\{x,y,z\\}$ and
$\hat{\sigma}_{\gamma}$ representing the standard Pauli matrices. For every
layer, there are $n+2$ adjustable parameters
$\vec{\theta}_{l}=\\{\vec{\omega}_{l},\beta_{l},\alpha_{l}\\}\in\\{\mathbb{R}^{n},\mathbb{R},\mathbb{R}\\}$,
and the scalar product is defined as
$\vec{\omega}_{l}\cdot\vec{x}=\omega_{l,1}x_{1}+\dots+\omega_{l,n}x_{n}$.
Thus, the complexity of the circuit increases linearly with the dimension of
the input data.

The output of the QNN, $F(\vec{x},\vec{\theta})$, is used to approximate a
target function $f(\vec{x})$, and it is obtained by calculating the
expectation value of an observable $\hat{\mathcal{M}}$ in the output state
$\ket{\phi(x)}$,

$F(\vec{x},\vec{\theta})=|\bra{\phi(\vec{x})}\hat{\mathcal{M}}\ket{\phi(\vec{x})}|\rightarrow{f(\vec{x})}.$
(4)

Initially, we choose $\hat{\mathcal{M}}=\hat{\sigma}_{z}$, such that the
mapping is done in the excitation amplitude of the qubit. However, more
sophisticated mappings are possible, e.g.,
$\hat{\mathcal{M}}=\hat{O}=u_{o}\mathbb{I}+u_{x}\hat{\sigma}_{x}+u_{y}\hat{\sigma}_{y}+u_{z}\hat{\sigma}_{z}$,
with $u_{o},u_{x},u_{y},u_{z}\in\mathbb{R}$ also trainable parameters. The
optimal network topology is constituted by the set of parameters
$\vec{\theta}$ that minimizes a particular cost function $\mathcal{C}$, i.e.,
$\vec{\theta}_{T}=\mbox{argmin}_{\vec{\theta}}\ \mathcal{C}(f;\vec{\theta})$.
For the considered regression problems, we use the squared error loss
function, which corresponds to

$\mathcal{C}(f;\vec{\theta})=\frac{1}{K}\sum_{k=1}^{K}\bigg{[}F(\vec{x}^{k},\vec{\theta})-f(\vec{x}^{k})\bigg{]}^{2},$
(5)

where $K$ is the total input samples of $\vec{x}$. The optimal parameters
$\vec{\theta}_{T}$ for a particular target function are obtained minimizing
Eq. (5) using gradient descent optimization, with $\vec{\theta}$ initialized
randomly. To explore a larger region of the parameter landscape, the
minimization is repeated $R$ times.

Note that the topology given by Eq. (3) is able to connect any two states on
the Bloch sphere. However, in order to capture non-trivial patterns in the
original data $\vec{x}$ and to approximate an arbitrary real function
$f(\vec{x})$, this data $\vec{x}$ needs to be reintroduced in the subsequent
layers [25]. Meanwhile, the complexity and approximation power of the network
can be increased by extending the model from one qubit to an arbitrary number
of qubits, $Q$. We first consider that the qubits forming the QNN do not
interact. Thus, the total transformation is the outer product of the unitary
transformations of all qubits after each layer, such that

$\hat{U}(\vec{x},\vec{\theta})=\prod_{l=1}^{L}\bigotimes_{q=1}^{Q}\hat{\mathcal{L}}_{l,q}(\vec{x},\vec{\theta}_{l,q}),$
(6)

where each qubit has a set of adjustable parameters
$\vec{\theta}_{l,q}=\\{\vec{\omega}_{l,q},\beta_{l,q},\alpha_{l,q}\\}$ and the
approximation of $f(\vec{x})$ is done via
$\hat{\mathcal{M}}=\hat{\sigma}_{z}^{tot}$, with
$\hat{\sigma}_{z}^{tot}=\sum_{q=1}^{Q}\hat{\sigma}_{z}^{q}$. The QNN has a
total number $N_{p}=LQ(n+2)$ of adjustable parameters $\vec{\theta}$. We now
introduce entangling gates at each layer among the different qubits in the QNN
model described by Eq. (6),

$\hat{U}(\vec{x},\vec{\theta})=\prod_{l=1}^{L}\hat{E}^{l}\cdot\bigg{(}\bigotimes_{q=1}^{Q}\hat{\mathcal{L}}_{l,q}(\vec{x},\vec{\theta}_{l,q})\bigg{)}.$
(7)

Although there are different ways of introducing entanglement to enrich the
network topology (by modifying either the entangling gates at each layer or
the entanglement allocation), for simplicity, we consider here a fixed gate
applied at each layer. In particular, to emphasize the role of the
entanglement, we define $\hat{E}^{l}$ as a maximally entangling gate, which
produces a generalization of Bell states when applied to the computational
basis [39],

$\hat{E}^{l}=\prod_{i=1}^{Q-1}\bigg{(}(\bigotimes_{i-1}\mathbb{I})\otimes\hat{E}_{2}\otimes(\bigotimes_{Q-i-1}\mathbb{I})\bigg{)},$
(8)

with

$\hat{E}_{2}=\frac{1}{\sqrt{2}}(\mathbb{I}\otimes\mathbb{I}+i\sigma_{x}\otimes\sigma_{y}).$
(9)

Note that the amount of entanglement is introduced as a fixed quantity, i.e.,
the gates in Eq. (8) do not introduce new adjustable parameters. To
distinguish different architectures, we label the topology of the QNN as
$\mathcal{A}_{0}(L,Q)$ ($\mathcal{A}_{1}(L,Q)$) with the subscripts $0$ ($1$)
indicating the absence (presence) of entangling gates. Two architectures,
$\mathcal{A}_{0}(L,Q)$ and $\mathcal{A}_{1}(L,Q)$, contain the same number
$N_{p}$ of variational parameters $\vec{\theta}$ and only differ by the
incorporation of entanglement. Consequently, we are able to evaluate
exclusively the role of this quantum resource in the performance of the QNN.

##  3 Approximation power

In order to characterize the expressivity of the QNN, we use the Fourier
decomposition of the unitary transformation produced by the QNN and compute
the number of harmonics that do not vanish after encoding the output state
through the observable $\hat{\mathcal{M}}$. We introduce a novel strategy for
bounding the expressivity of the QNN [38, 40] that enables to compute the
Fourier amplitudes and to identify zero coefficients.

For each qubit, the single layer unitary transformation given by Eq. (3) is
decomposed into the fundamental frequencies $\pm\vec{\omega}_{l,q}$ of the
qubit,

$\displaystyle\hat{\mathcal{L}}_{l,q}(\vec{x},\vec{\theta}_{l,q})=\hat{R}_{x}(2(\vec{\omega}_{l,q}\cdot\vec{x}+\beta_{l,q}))\hat{R}_{y}(2\alpha_{l,q})=\hat{T}_{l,q}^{+}e^{i\vec{\omega}_{l,q}\cdot\vec{x}}+\hat{T}_{l,q}^{-}e^{-i\vec{\omega}_{l,q}\cdot\vec{x}},$
(10)

where

$\hat{T}^{+}_{l,q}=\frac{1}{2}\left(\matrix{t_{l,q}^{-}&-t_{l,q}^{+}\cr-
t_{l,q}^{-}&t_{l,q}^{+}}\right)e^{i\beta_{l,q}},\quad\hat{T}^{-}_{l,q}=\frac{1}{2}\left(\matrix{t_{l,q}^{+}&t_{l,q}^{-}\cr
t_{l,q}^{+}&t_{l,q}^{-}}\right)e^{-i\beta_{l,q}},$ (11)

with $t_{l,q}^{\pm}=\cos(\alpha_{l,q})\pm\sin(\alpha_{l,q})$. The unitary
transformation of the network Eq. (6), is then given in the form of a partial
Fourier series,

$\hat{U}(\vec{x},\vec{\theta})=\prod_{l=1}^{L}\hat{E}^{l}\bigotimes_{q=1}^{Q}(\hat{T}_{l,q}^{+}e^{i\vec{\omega}_{l,q}\cdot\vec{x}}+\hat{T}_{l,q}^{-}e^{-i\vec{\omega}_{l,q}\cdot\vec{x}})=\sum_{j=1}^{\mathcal{N}_{h}}\hat{C}_{j}(\alpha,\beta)\exp{(i\vec{M}_{j}(\vec{\omega})\cdot\vec{x})},$
(12)

where $\alpha=\\{\alpha_{l,q}\\}$, $\beta=\\{\beta_{l,q}\\}$,
$\vec{\omega}=\\{\vec{\omega}_{l,q}\\}$, $\hat{C}_{j}$ are $2^{Q}\times 2^{Q}$
real matrices, and $\vec{M}_{j}(\vec{\omega})$ $\in$ $\mathbb{R}^{n}$ results
from the linear combination of the distinct frequencies $\vec{\omega}_{l,q}$.
The product of the unitaries in each layer has been rewritten as a single
summation that accounts for the total number of terms $\mathcal{N}_{h}$ in the
partial Fourier decomposition of the network. The spectrum of the QNN has a
size of $\mathcal{N}_{h}=2^{LQ}$, and depends solely on the data encoding
strategy, regardless of the presence of entanglement, which only affects the
values of the matrix coefficients $\hat{C}_{j}$ of the Fourier series [38].
Additional information regarding this derivation can be found in A.

Similarly, the output generated by the QNN can be decomposed into a partial
Fourier series,

$F(\vec{x},\vec{\theta})=\bra{0}\hat{U}^{\dagger}(\vec{x},\vec{\theta})\hat{\mathcal{M}}\hat{U}(\vec{x},\vec{\theta})\ket{0}=\sum_{j=1}^{N_{h}}c_{j}(\alpha,\beta)\exp\left(i\vec{\Omega}_{j}(\vec{\omega})\cdot\vec{x}\right).$
(13)

The details on how to perform this decomposition are given in B. The number of
terms in this expansion is $N_{h}\leq 3^{LQ}$. Equation (13) describes the
family of functions that a given QNN can learn through two interrelated
characteristics: the spectrum of Fourier frequencies $\vec{\Omega}_{j}$ $\in$
$\mathbb{R}^{n}$ accessible by the quantum model, and the expressivity, that
refers to the number of expansion coefficients $c_{j}$ $\in$ $\mathbb{R}$ that
a class of QNN can control independently, and defines the set of attainable
functions. It is important to notice that these two characteristics are fully
determined by the $N_{p}=LQ(n+2)$ topological hyper-parameters
$\vec{\theta}_{l,q}$, and each set plays separate roles; $\vec{\omega}_{l,q}$
controls the frequency spectrum, while $\alpha_{l,q}$ and $\beta_{l,q}$, the
amplitude of each Fourier term $c_{j}$. Crucially, the number of non-vanishing
expansion coefficients $N_{h}$ depends on the observable $\hat{\mathcal{M}}$
used to map the function $F(\vec{x},\vec{\theta})\rightarrow f(\vec{x})$. Note
that, in general, $N_{h}\gg N_{p}$. As a result, only a subset of both $c_{j}$
and $\Omega_{j}$ may be totally tunable, while the rest are constrained. The
variational formulation of a specific machine learning problem is equivalent
to finding the best fit for the Fourier decomposition. While the optimal
combination depends on the specific problem to be solved –i.e., $f(\vec{x})$–,
the expressivity of the network is determined by the topological parameters
$N_{p}$ and $N_{h}$. As we will show, the presence of entanglement, either by
a global readout or a final entangling layer, produces a larger set $N_{h}$
compared to a local readout in a non-entangling QNN. This enriches the
accessible functions that the QNN can access, leading to an enhancement of the
network accuracy when approximating $f(\vec{x})$.

In the upcoming sections, we determine the count of non-zero expansion
coefficients, denoted as $N_{h}$. This calculation is conducted both with and
without considering entanglement between layers, and depends on the choice of
observable utilized for the readout of the output state.

###  3.1 One qubit network

For the particular case of a single qubit architecture $\mathcal{A}_{0}(L,1)$,
we recursively count the number of vanishing $c_{j}(\alpha,\beta)$. According
to Eq. (10), $L$ layers acting on a single qubit produce the transformation

$\hat{v}=\prod_{l=1}^{L}(\hat{T}_{l}^{+}e^{i\vec{\omega}_{l}\cdot\vec{x}}+\hat{T}_{l}^{-}e^{-i\vec{\omega}_{l}\cdot\vec{x}}).$
(14)

For the particular layer architecture given by Eq. (3), the number of
expansion terms $N_{h}$ depends on the considered observable
$\hat{\mathcal{M}}$, being the maximum number of achivable terms
$\max(N_{h})=\chi=3^{L}$. As shown in C, the cases
$\hat{\mathcal{M}}=\hat{\sigma}_{z}$ and $\hat{\mathcal{M}}=\hat{\sigma}_{y}$
lead to $N_{h}\leq 2\cdot\chi/3$, while measuring the observable
$\hat{\mathcal{M}}=\hat{\sigma}_{x}$ gives $N_{h}\leq\chi/3$. In general, such
a layer topology is capable of reaching the maximal expressivity bound
$N_{h}=\chi$ by measuring the operator
$\hat{\mathcal{M}}=\hat{O}=u_{o}\mathbb{I}+u_{x}\hat{\sigma}_{x}+u_{y}\hat{\sigma}_{y}+u_{z}\hat{\sigma}_{z}$
along various directions defined by the coefficients
$u_{o},u_{x},u_{y},u_{z}\in\mathbb{R}$. Notice that using another layer
topology when constructing a data re-uploaded QNN might lead to a different
number of non-vanishing expansion coefficients. For instance, the architecture
given by

$\hat{\mathcal{L}}^{*}_{l}(\vec{x},\vec{\theta}_{l})=\hat{R}_{x}(2\alpha_{l})\hat{R}_{y}\left(2\left(\vec{\omega}_{l}\cdot\vec{x}+\beta_{l}\right)\right),$
(15)

results in $N_{h}\leq 2\cdot\chi/3+1$, regardless of the observable used for
the readout. Table 1 presents the different achievable expresivities in terms
of the observable $\hat{\mathcal{M}}$ and encoding layer.

$Q=1$ | $\hat{\mathcal{M}}$
---|---
| $\hat{\sigma}_{x}$ | $\hat{\sigma}_{y}$ | $\hat{\sigma}_{z}$ | $\hat{O}$
$N_{h}$ | $\hat{\mathcal{L}}$ | $\chi/3$ | $2\chi/3$ | $2\chi/3$ | $\chi$
$\hat{\mathcal{L}}^{*}$ | $2\chi/3+1$ | $2\chi/3+1$ | $2\chi/3+1$ | $2\chi/3+1$
Table 1: Maximum number $N_{h}$ of non-vanishing expansion coefficients when
expanding the output of the network using one qubit ($Q=1$) into a partial
Fourier series. Different observables $\hat{\mathcal{{M}}}$ and the encoding
architectures given by Eqs. (3) and (15) are considered. Here $\chi=3^{L}$ is
the spectrum bound.

###  3.2 Multiqubit architectures

The previous results can easily be generalized for multiqubit
$\mathcal{A}_{0}(L,Q)$ architectures without entanglement by taking the outer
product of the unitary acting on each qubit, i.e.,

$\hat{U}=\bigotimes_{q=1}^{Q}\hat{v}_{q},$ (16)

with
$\hat{v}_{q}=\prod_{l=1}^{L}(\hat{T}_{l,q}^{+}e^{i\vec{\omega}_{l,q}\cdot\vec{x}}+\hat{T}_{l,q}^{-}e^{-i\vec{\omega}_{l,q}\cdot\vec{x}})$.
The absence of the entangling gates between layers allows an exchange of
products over layers and qubits in Eq. (12) to recover Eq. (16).

Now, two different cases are considered depending on the election of the
observable $\hat{\mathcal{M}}$. At first, we consider that the mapping
$\hat{\mathcal{M}}$ of the whole system is composed of local measurements
$\hat{\mathcal{M}}_{q}$ of the individual qubits
$\hat{\mathcal{M}}=\sum_{q=1}^{Q}\hat{\mathcal{M}}_{q}$. As a result, the
generated output is decomposed into

$F(\vec{x},\vec{\theta})=\bra{0}(\bigotimes_{q=1}^{Q}\hat{v}_{q})^{\dagger}(\hat{\mathcal{M}}_{1}+...+\hat{\mathcal{M}}_{Q})(\bigotimes_{q=1}^{Q}\hat{v}_{q})\ket{0}=\sum_{q=1}^{Q}\bra{0}\hat{v}_{q}^{\dagger}\hat{\mathcal{M}}_{q}\hat{v}_{q}\ket{0}.$
(17)

Each term $\bra{0}v_{q}^{\dagger}\hat{\mathcal{M}}_{q}\hat{v}_{q}\ket{0}$ in
the expansion produces $N_{h}(L,1)$ coefficients, leading to a total of
$N_{h}(L,Q)=Q\cdot N_{h}(L,1)$ non-vanishing expansion coefficients of the
partial Fourier series for the $\mathcal{A}_{0}(L,Q)$ architecture.

Alternatively, we can perform a global measurement of the network by defining
the observable $\hat{\mathcal{M}}=\otimes_{q=1}^{Q}\hat{\mathcal{M}}_{q}$. In
this case,

$F(\vec{x},\vec{\theta})=\bra{0}(\bigotimes_{q=1}^{Q}\hat{v}_{q})^{\dagger}(\hat{\mathcal{M}}_{1}\otimes\hat{\mathcal{M}}_{2}\otimes...\otimes\hat{\mathcal{M}}_{Q})(\bigotimes_{q=1}^{Q}v_{q})\ket{0}=\prod_{q=1}^{Q}\bra{0}v_{q}^{\dagger}\hat{\mathcal{M}}_{q}\hat{v}_{q}\ket{0},$
(18)

leading to $N_{h}(L,Q)=N_{h}(L,1)^{Q}$ non-vanishing expansion coefficients.
The same architecture can generate more harmonics by just modifying the
measured observable.

$Q=2,\hat{\mathcal{L}}$ | $\hat{\mathcal{M}}$
---|---
| $\hat{\sigma}_{z}^{1}+\hat{\sigma}_{z}^{2}$ | $\hat{\sigma}_{z}^{1}\otimes\hat{\sigma}_{z}^{2}$ | $\hat{O}_{1}\otimes\hat{O}_{2}$
$N_{h}$ | without $\hat{E}_{2}$ | $4\chi/3$ | $4\chi^{2}/9$ | $\chi^{2}$
$\hat{E}_{2}$ all layers | $5\chi^{2}/9$ | $4\chi^{2}/9$ | $\chi^{2}$
|  $\hat{E}_{2}$ last layer
---
$5\chi^{2}/9$ | $4\chi^{2}/9$ | $\chi^{2}$
Table 2: Maximum number $N_{h}$ of non-vanishing expansion coefficients when
expanding the network’s output $F(\vec{x},\vec{\theta})$ into a partial
Fourier series using $Q=2$ qubits. Different observables $\hat{\mathcal{{M}}}$
are considered when using the encoding architecture. In particular, the number
$N_{h}$ is compared with and without entanglement. Here $N_{h}$ is compared
with the maximum number of harmonics of a single qubit network $\chi=3^{L}$.

###  3.3 Last layer entanglement

As commented at the beginning of Sec. 3 and further elaborated by Schuld et
al. [38], entanglement does not extend the frequency spectrum of the unitary
tranformation given by Eq. (12), but only the values of the matrix
coefficients $\hat{C}_{j}$. A systematic approach to analyze the full
tunability of $C_{j}$ is still needed. However, entanglement does affect the
size of the spectrum of $F(\vec{x},\vec{\theta})$ that involves the mapping of
$f(\vec{x})$ trough an observable $\hat{\mathcal{M}}$. To gain a physical
understanding of the effects of introducing entanglement exclusively in the
final layer, we examine a specific network architecture with $Q=2$ and the
observable
$\hat{\mathcal{M}}=\hat{\sigma}_{z}^{tot}=\hat{\sigma}_{z}^{1}+\hat{\sigma}_{z}^{2}$.
When the last output layer incorporates an entangling gate in the form of Eq.
(9), the resulting neural network is equivalent to a non-entangling
$\mathcal{A}_{0}(L,2)$ topology being mapped by the effective observable

$\hat{\mathcal{M}}^{\prime}=\hat{E}_{2}^{\dagger}\hat{\mathcal{M}}\hat{E}_{2}=\hat{E}_{2}^{\dagger}\hat{\sigma}_{z}^{tot}\hat{E}_{2}=\hat{\sigma}_{x}^{1}\otimes\hat{\sigma}_{x}^{2}-\hat{\sigma}_{y}^{1}\otimes\hat{\sigma}_{y}^{2}.$
(19)

In this case, $\hat{\mathcal{M}}^{\prime}$ is a global observable and
increases the number of non-vanishing expansion coefficients compared to the
local $\hat{\mathcal{M}}=\hat{\sigma}_{z}^{tot}$ measurements. More generally,
the exact number of generated expansion coefficients depends on the newly
transformed observable $\hat{\mathcal{M}}^{\prime}$ through $\hat{E}_{2}$. In
Table 2, we provide the values of the spectrum size of
$F(\vec{x},\vec{\theta})$ generated by various mapping strategies and
topologies for two-qubit QNNs. The situation becomes more complex when dealing
with architectures featuring entanglement across all layers, as the
factorization performed in Eq. (17) is not permitted, but it results in
similar scaling behaviors, as discussed in D.

Figure 1: Fourier transform of the output generated by a randomly initialized
QNN showing different expansion coefficients $c_{j}$, by using (a)
$\mathcal{A}_{0}^{T}(2,2)$ and (b) $\mathcal{A}_{1}^{T}(2,2)$. Sharing the
same $\vec{x}$ and $\vec{\theta}$ values, both QNNs have the same number of
variable parameters, however, the $\mathcal{A}_{1}^{T}(2,2)$ model generates
more harmonics. The measured operator is
$\hat{\mathcal{M}}=\hat{\sigma}_{z}^{tot}=\hat{\sigma}_{z}^{1}+\hat{\sigma}_{z}^{2}$.

##  4 Numerical results

For simple topologies, we were able to analytically characterize the
expressivity of a network by counting the number of non-vanishing expansion
coefficients in the partial Fourier decomposition. However, for complex deep
QNNs, this procedure becomes more intrincate. An alternative approach involves
directly Fourier transforming the generated output
$F(\vec{x},\vec{\theta})\rightarrow\mathcal{F}(\vec{\Omega},\vec{\theta})$,
with the number of coefficients $c_{j}$ determined by the count of peaks in
the Fourier transform. This numerical procedure allows us to extrapolate the
harmonics count to deep entangling $\mathcal{A}_{1}(L,Q)$ models
straightforwardly. In Fig. 1, we plot the Fourier transform of the outputs
generated by (a) $\mathcal{A}_{0}(2,2)$ and (b) $\mathcal{A}_{1}(2,2)$ using
the same encoding observable
$\hat{\mathcal{M}}=\sum_{q=1}^{Q}\hat{\sigma}_{z}^{q}$ and the random but
identically initialized set of parameters $\vec{\theta}$. The number of peaks
in both the architecture without and with entanglement ($12$ and
$(36)$,respectively) matches the analytical expressions in Table 2. Although
both $\mathcal{A}_{0}$ and $\mathcal{A}_{1}$ have the same number of tunable
parameters, the inclusion of entanglement in the QNN triples the spectrum.

To quantify the ability of the network to generate non-vanishing expansion
coefficients in the partial Fourier decomposition, we define the ratio
$\Gamma=N_{h}/N_{p}$. This index quantifies the number of harmonics per
adjustable parameter in the QNN, and facilitates the visualization of the QNN
expressivity. Considering the encoding observable described above
$\hat{\mathcal{M}}=\sum_{q=1}^{Q}\hat{\sigma}_{z}^{q}$, for
$\mathcal{A}_{0}(L,Q)$ models, $N_{h}=Q\cdot 2\chi/3$. Consequently, the ratio
$\Gamma$ grows as $3^{L}/L$, regardless of the number of qubits, as depicted
in Fig. 2. In this figure, we compare this result with the $\Gamma$ parameter
of $\mathcal{A}_{1}(L,Q)$ models, which incorporate entangling gates in all
layers. While both $\mathcal{A}_{0}$ and $\mathcal{A}_{1}$ have the same
number of adjustable parameters $N_{p}$, remarkably, the mere presence of
entanglement further amplifies the exponential relationship between $\Gamma$
and $L$ as the number of qubits increases.

Figure 2: Ratio $\Gamma$ as a function of the number of layers $L$ for
different number of qubits $Q$, computed numerically through the procedure
depicted in Fig. 1 (symbols). All the $\mathcal{A}_{0}$ models collapse to the
same ratio independently of the number of qubits (red-solid line). entangling
models $\mathcal{A}_{1}$ show a different exponential slope for $Q=2$ (cyan
triangles), $Q=3$ (blue squares), and $Q=4$ (black circles). The analytical
results $\Gamma=2(3^{L-1})/[(2+n)L]$ (red-solid line) and
$\Gamma=5(3^{L-1})^{2}/[(2+n)2L]$ (cyan-solid line) fit these trends. The
dashed lines facilitate the visualization of the linear dependency. In all
cases the measured operator is
$\hat{\mathcal{M}}=\sum_{q=1}^{Q}\hat{\sigma}_{z}^{q}$.

###  4.1 Teacher-student benchmark

In this section, we systematically evaluate the neural network approximation
capabilities for real-continuous functions when considering various
architectures, i.e., network depth, number of qubits, and presence or absence
of entanglement in all layers. To benchmark the QNNs and minimize any
potential dataset bias, we employ the teacher-student scheme [19, 41]. The
teacher model (T) generates datasets by mapping random inputs to outputs, and
these datasets are subsequently learned by the student models (S). To assess
the performance of the student models, we employ prediction maps to
qualitatively compare and visualize the similarity between T and S, see E.1.
For a rigorous quantitative comparison, we benchmark the average performance
[42] of both models through the loss function in Eq. (5).

####  4.1.1 The teacher output.

The role of the teacher is to map the input data $\vec{x}$ through a QNN into
a random output function $F(\vec{x},\vec{\theta})$ by the observable
$\hat{\mathcal{M}}=\sum_{q}^{Q}\hat{\sigma}_{z}^{q}$, which then becomes the
target $f(\vec{x})$ of S models. All the adjustable parameters
$\vec{\theta}\in[-\frac{\pi}{2},\frac{\pi}{2}]$ and $K$ input values
$\vec{x}\in[0,1]^{n}$ are randomly initialized to avoid any possible bias on
the datasets. The complexity of the target function $f(\vec{x})$ can be
controlled by the depth of the teacher architectures.

In Fig. 3, we show the output functions generated by several T-models
containing $4$ qubits, without entanglement (upper panel) and with it (lower
panel), when the input data is two-dimensional, $n=2$. Deeper architectures
(left to right) provide in general a ‘richer’ topology in
$F(\vec{x},\vec{\theta})$.

Figure 3: Prediction maps generated by different T models of a QNN containing
$Q=4$ qubits in different number of layers. Upper (lower) panels correspond to
QNNs of $\mathcal{A}_{0}^{T}(\mathcal{A}_{1}^{T})$ architectures. Any pair
$\\{\mathcal{A}_{0}^{T},\mathcal{A}_{1}^{T}\\}$ with same $L$ is equally
initialized in the network parameters $\vec{\theta}$ and input $\vec{x}$, but
differs in the presence (absence) of entanglement.

####  4.1.2 The student performance.

First, we analyze the performance of the student as a function of the depth of
the network, $L$. To this end, we focus on unidimensional ($n=1$) regression
problems. We create a $\mathcal{A}^{T}_{0}(10,10)$ T-model, which randomly
generates the output corresponding to the target function $f(x^{k})$ that is
depicted by the red-dashed line of Fig. 4a. Then, different
$\mathcal{A}_{0}^{S}(L,1)$ S-models are trained to approximate this function.
As inferred from Fig. 4b, the S-model’s performance improves as the number of
layers is increased, until it reaches a point of saturation, at which the
minimum value of the loss function $\mathcal{C}$ becomes independent of $L$.
For a fixed network topology, the particular value of $\min(\mathcal{C})$ and
layer at which it saturates depend on the specific target function and input-
data size. Note that a $\mathcal{A}_{0}^{S}(4,1)$ QNN programmed in ibmq_belem
(black circles of Fig. 4a) already reaches $\mathcal{C}\sim 10^{-5}$,
evidencing the capabilities of a single qubit as a universal approximator
[37].

Figure 4: (a) S-model outputs composed by a single qubit with different number
of layers $L$ (solid-lines) are listed to learn the targets $f(x^{k})$
function randomly generated by the T-model $\mathcal{A}^{T}_{0}(10,10)$ (red-
dashed). A single qubit QNN $\mathcal{A}^{S}_{0}(4,1)$ run in ibmq_belem shows
its learnability (black circles). (b) Corresponding loss error $\mathcal{C}$
(5) in terms of the data re-uploads $L$. (c) Output of the teacher
$\mathcal{A}_{0}^{T}(5,8)$ (continuous contour plot) for a two dimensional
problem $n=2$ and the fitting result for a smaller $\mathcal{A}_{0}^{S}(2,2)$
student (dashed contour plot). (d) Corresponding loss error $\mathcal{C}$ for
different student models $\mathcal{A}_{0}^{S}(2,Q)$ fitting the bi-dimensional
data (c) generated by $\mathcal{A}_{0}^{T}(5,8)$ (red-dashed line). Parameters
$K=400$ and $R=10$.

Secondly, we study the performance of QNNs in terms of the number of qubits,
$Q$, by fixing a $\mathcal{A}_{0}^{T}(5,8)$ T-model architecture for a two-
dimensional input data ($n=2$) and considering different
$\mathcal{A}_{0}^{S}(2,Q)$ students. In Fig. 4c, we observe that a student
with two qubits and two layers can produce an output function that closely
resembles that of the teacher. In Fig. 4d, we represent the cost function as a
function of the number qubits of the student. A similar trend is obtained, the
approximation power of the S-models increases with the number of qubits until
the loss function saturates.

Analogously to classical neural networks, where the approximation power
increases with the number of hidden layers and neurons [43], the accuracy of
the QNN increases with the amount of data re-uploadings (or layers) and number
of qubits. However, in contrast to its classical counterpart, the QNN
regression capabilities can be significantly enhanced with the incorporation
of quantum resources, such as global readouts or entangling gates.

Therefore, our objective is to examine how entanglement influences the network
approximation capabilities. For two-dimensional input data ($n=2$), we
generate different T-models, $\mathcal{A}^{T}_{0}(L_{Max},Q)$ or
$\mathcal{A}^{T}_{1}(L_{Max},Q)$. Subsequently, we evaluate the performance of
S-models, $\mathcal{A}^{S}_{0}(L,Q)$ and $\mathcal{A}^{S}_{1}(L,Q)$, for both
of these teacher models. The depth of the students ranges from 1 to $L_{Max}$,
where $L_{Max}$ represents the number of layers in the T-models. To enlarge
the benchmark, each T-model is randomly initialized $M$ times, creating
multiple target functions for the students. As a result, for fixed teacher and
student architectures, the QNN is evaluated $M\times R$ times, where $R$ is
the number of different initializations for S-models.

In Fig. 5, we plot the student average performance
$\bar{\mathcal{C}}=\sum_{i}^{M}\mathcal{C}_{i}/M$ over different realizations
as a function of the number of layers for both the $\mathcal{A}^{S}_{0}$ and
$\mathcal{A}^{S}_{1}$ models. We focus on the $\mathcal{A}_{1}^{T}$ target
function, since it generates a more complex data structure due to the effect
of an entangling readout. The $\mathcal{A}_{0}^{T}$ case is analyzed in E.2.
Nevertheless, both scenarios lead to similar conclusions. When comparing the
various plots in Fig. 5, it becomes apparent that, as the number of layers in
the QNN increase, the differentiation in the performance of both
$\mathcal{A}_{0}^{S}$ and $\mathcal{A}_{1}^{S}$ becomes increasingly
noticeable. This observation is readily apparent in Fig. 2; for two QNN pairs
with the same number of qubits, $\mathcal{A}_{0}^{S}$ and
$\mathcal{A}_{1}^{S}$, the difference in their $\Gamma$ values grows with $L$,
primarily owing to the greater number of generated harmonics in the entangling
network. A similar behavior is also observed when the number of qubits is
increased in the $\mathcal{A}_{1}^{S}$ topology, as it was already anticipated
in Fig. 2. For a fixed number of layers in the entangling architecture, the
steepness of $\Gamma$ increases with $Q$, leading to a larger set of Fourier
harmonics. However, this is not the case for $\mathcal{A}_{0}^{S}$ networks,
since they all lead to the same $\Gamma$ value, regardless of the number of
qubits. As a consequence, the efficiency of different $\mathcal{A}_{0}^{S}$
models saturate to a common $\bar{\mathcal{C}}$ value. This behavior is more
evident in the figures presented in E.2. As a result, a deep QNN involving a
global mapping of $f(\vec{x})$ leads to improvements of 1-2 orders of
magnitude in the approximation accuracy compared to local readouts in a non-
entangling network, see Fig. 5. Finally, note that, in general, deep
entangling QNNs also present higher robustness over different realizations in
the predicted output, as inferred from the smaller vertical lines of Figs. 5
and 7 showing the standard deviation
$\Delta\mathcal{C}=\sqrt{\sum_{i}^{M}(\mathcal{C}_{i}-\bar{\mathcal{C}})^{2}/M}$.

Figure 5: Mean loss error $\bar{\mathcal{C}}$ as a function of the number of
layers $L$ for different student models $\mathcal{A}_{0}^{S}(L,Q)$ (red
dashed-lines) and $\mathcal{A}_{1}^{S}(L,Q)$ (blue solid-lines). The target
function is prepared by teacher models $\mathcal{A}_{1}^{T}(L_{Max},Q)$
containing entanglement and (a) $Q=2$, (b) $Q=3$, and (c) $Q=4$ qubits.
Parameters: $n=2,K=400,R=5$, and $M=5$.

The expressive power of the QNN primarily depends on its topological
parameters. However, when considering regression tasks, other factors such the
QNN readout $\hat{\mathcal{M}}$ come into play. Generally, a larger $\Gamma$
parameter tends to result in a smaller average error $\bar{\mathcal{C}}$,
providing a higher accuracy when interpolating $f(\vec{x})$, except in two
distinct scenarios. Firstly, for small QNNs (as depicted in Fig. 5a), or when
dealing with simple $f(\vec{x})$ structures (as shown in Fig. 7), there can be
an overestimation of $f(\vec{x})$. The additional harmonics introduced by
$\mathcal{A}_{1}^{S}(L,Q)$ can lead to a less accurate estimation compared to
the non-entangling model $\mathcal{A}_{0}^{S}(L,Q)$, resulting in a decrease
of the network accuracy. Secondly, even though the expressive power grows
exponentially with the number of layers, we observe a saturation in the
accuracy of deep QNNs. Both phenomena can be attributed to the fact that the
$N_{h}$ Fourier coefficients have limited tunability due to the finite number
$N_{p}=LQ(n+2)$ of adjustable parameters.

## Outlook

We have analyzed the expressivity of a universal data re-uploading QNN and
systematically benchmarked its performance in regression tasks, showing the
relevance of the observable that maps the network generated output and the
target function. We show that the presence of entanglement, produced by a
global readout of the QNN or by the incorporation of a final entangling layer,
leads to the maximal expressivity of the network. As a result, the
approximation capabilities improve with respect to those based on local
readouts of the output qubits.

This work paves the way to several prospective research lines. For example,
(i) searching for alternative mapping and encoding strategies to expand the
partial Fourier decomposition of the generated output, and thus enlarge the
approximation efficiency of the network, such as using different building
blocks for the QNN like multi-level qudits [44] or unbounded quantum systems;
(ii) looking for applications where the designed QNN acts as a new element in
a larger quantum protocol. In this context, there is the potential to design
novel neural network-based quantum gates [29], which expand the quantum
computing toolkit and streamline the complexity of existing algorithms. These
gates can be engineered to influence the response of specific observables to
input signals [6], opening the door to new sensing schemes, as well as
providing efficient representation and classification of large quantum many-
body systems [45, 46].

The authors thank J. J. García-Ripoll for fruitful discussions. We also
acknowledge financial support from the Spanish Government PID2021-126694NA-C22
and PID2021-126694NB-C21, and by Comunidad de Madrid-EPUC3M14 and CAM/FEDER
Project No. S2018/TCS-4342 (QUITEMAD-CM). This work is also supported by
Arquimea Research Center and by Horizon Europe, Teaming for Excellence, under
grant agreement No 101059999, project QCircle. H. E. acknowledges the FPU
program (FPU20/03409). J. C. and E. T. acknowledge the Ramón y Cajal program
(RYC2018-025197-I and RYC2020-030060-I). Y. B. acknowledges the CDTI within
the Misiones 2021 program and the Ministry of Science and Innovation under the
Recovery, Transformation and Resilience Plan-Next Generation EU under the
project “CUCO: Quantum Computing and its Application to Strategic Industries”.

##  Appendix A Derivation of the circuit unitary transformation
$\hat{U}_{L,Q}(\vec{x},\vec{\theta})$

The unitary transformation in each layer can be decomposed as the sum of two
harmonics $\omega_{l,q}$ and $-\omega_{l,q}$

$\displaystyle\hat{\mathcal{L}}_{l,q}(\vec{x},\vec{\theta}_{l})=R_{x}(2(\vec{\omega}_{l,q}\cdot\vec{x}+\beta_{l,q}))R_{y}(2\alpha_{l,q})=$
$\displaystyle=\left(\matrix{\cos(\vec{\omega}_{l,q}\cdot\vec{x}+\beta_{l,q})&-i\sin(\vec{\omega}_{l,q}\cdot\vec{x}+\beta_{l,q})\cr-i\sin(\vec{\omega}_{l,q}\cdot\vec{x}+\beta_{l,q})&\cos(\vec{\omega}_{l,q}\cdot\vec{x}+\beta_{l,q})}\right)\left(\matrix{\cos(\alpha_{l,q})&-\sin(\alpha_{l,q})\cr\sin(\alpha_{l,q})&\cos(\alpha_{l,q})}\right)=$
$\displaystyle=\hat{T}_{l,q}^{+}e^{i\vec{\omega}_{l,q}\cdot\vec{x}}+\hat{T}_{l,q}^{-}e^{-i\vec{\omega}_{l,q}\cdot\vec{x}},$
(20)

where

$\hat{T}^{+}_{l,q}=\frac{1}{2}\left(\matrix{t_{l,q}^{-}&-t_{l,q}^{+}\cr-
t_{l,q}^{-}&t_{l,q}^{+}}\right)e^{i\beta_{l,q}},\quad\hat{T}^{-}_{l,q}=\frac{1}{2}\left(\matrix{t_{l,q}^{+}&t_{l,q}^{-}\cr
t_{l,q}^{+}&t_{l,q}^{-}}\right)e^{-i\beta_{l,q}},$ (21)

with $t_{l,q}^{\pm}=\cos(\alpha_{l,q})\pm\sin(\alpha_{l,q})$. We use this
expression to decompose a network into products of this sum. Then, we write
these products as a factorization of single layers gates:

$\hat{U}(\vec{x},\vec{\theta})=\prod_{l=1}^{L}\bigotimes_{q=1}^{Q}(\hat{T}_{l,q}^{+}e^{i\vec{\omega}_{l,q}\cdot\vec{x}}+\hat{T}_{l,q}^{-}e^{-i\vec{\omega}_{l,q}\cdot\vec{x}})=\prod_{l=1}^{L}\hat{u}_{l},$
(22)

with
$\hat{u}_{l}=\bigotimes_{q=1}^{Q}(\hat{T}_{l,q}^{+}e^{i\vec{\omega}_{l,q}\cdot\vec{x}}+\hat{T}_{l,q}^{-}e^{-i\vec{\omega}_{l,q}\cdot\vec{x}})$.
The harmonics generated by this transformation are combinations of
$\omega_{l,q}$ with positive or negative sign. To compute the Kronecker
product, we define $P=P(Q)$ as the set of all possible sign permutations of
$Q$ terms. For example, for $Q=3$ we have

$P(3)=\\{\\{+,+,+\\},\\{+,+,-\\},\\{+,-,-\\},\cdots,\\{-,-,-\\}\\}.$ (23)

We denote $P^{k}$ the $k$-th element in $P$, and $P^{k}_{n},$ the $n$-th
element in $P^{k}$ with $P^{k}_{n}=+$ or $-$ for $k=1,...,2^{N}$ and
$n=1,...,N$.

Using this notation, and defining $\hat{T}_{l,q}[\pm]=\hat{T}^{\pm}_{l,q}$, we
express the action of a single layer as

$\hat{u}_{l}(\vec{x},\vec{\theta})=\sum^{2^{Q}}_{k=1}\left(\bigotimes_{q=1}^{Q}\hat{T}_{l,q}\left[P_{q}^{k}\right]\right)\exp{\left(i\sum_{q=1}^{Q}P_{q}^{k}\vec{\omega}_{l,q}\vec{x}\right)}=\sum_{k=1}^{2^{Q}}\hat{B}^{k}_{l}e^{i\vec{m}^{k}_{l}\cdot\vec{x}},$
(24)

with
$\hat{B}^{k}_{l}=\hat{B}^{k}_{l}(\alpha,\beta)=\bigotimes_{q=1}^{Q}\hat{T}_{l,q}\left[P_{q}^{k}\right]$
and
$\vec{m}^{k}_{l}=\vec{m}^{k}_{l}(\vec{\omega})=(\sum_{q=1}^{Q}P_{q}^{k}\vec{\omega}_{l,q})$.
For illustration purposes, in the case $Q=3$, ignoring the $l$ subindex in
$\hat{T}^{\pm}_{l,q}$ and $\vec{\omega}_{l,q}$ for clarity,

$\displaystyle\hat{u}_{l}(\vec{x},\vec{\theta})=(\hat{T}_{3}^{+}e^{i\vec{\omega}_{3}}+\hat{T}^{-}_{3}e^{-i\vec{\omega}_{3}})\otimes(\hat{T}_{2}^{+}e^{i\vec{\omega}_{2}}+\hat{T}^{-}_{2}e^{-i\vec{\omega}_{2}})\otimes(\hat{T}_{1}^{+}e^{i\vec{\omega}_{1}}+\hat{T}^{-}_{1}e^{-i\vec{\omega}_{1}})=$
$\displaystyle\hat{T}^{+}_{3}\otimes\hat{T}^{+}_{2}\otimes\hat{T}^{+}_{1}e^{i(\vec{\omega}_{3}+\vec{\omega}_{2}+\vec{\omega}_{1})\cdot\vec{x}}+\hat{T}^{+}_{3}\otimes\hat{T}^{+}_{2}\otimes\hat{T}^{-}_{1}e^{i(\vec{\omega}_{3}+\vec{\omega}_{2}-\vec{\omega}_{1})\cdot\vec{x}}+$
$\displaystyle\hat{T}^{+}_{3}\otimes\hat{T}^{-}_{2}\otimes\hat{T}^{-}_{1}e^{i(\vec{\omega}_{3}-\vec{\omega}_{2}-\vec{\omega}_{1})\cdot\vec{x}}+\cdots+\hat{T}^{-}_{3}\hat{\otimes}\hat{T}^{-}_{2}\otimes\hat{T}^{-}_{1}e^{i(-\vec{\omega}_{3}-\vec{\omega}_{2}-\vec{\omega}_{1})\cdot\vec{x}}$
(25)

Introducing entanglement on these architectures leads to a similar
decomposition as Eq. (24), replacing the expansion coefficients
$\hat{B}^{k}_{l}\rightarrow\hat{E}^{l}\hat{B}^{k}_{l}$, while the frequency
spectrum remains the same.

Introducing the expresssion for $\hat{u}_{l}$ obtained in Eq. (24) into Eq.
(22), the action of a the whole transformation produced by the QNN is written
as a sum of harmonic terms,

$\displaystyle\hat{U}(\vec{x},\vec{\theta})=\prod_{l=1}^{L}\hat{u}_{l}=\prod_{l=1}^{L}\sum_{k=1}^{2^{Q}}\hat{B}^{k}_{l}e^{i\vec{m}^{k}_{l}\cdot
x}=$
$\displaystyle\sum_{j=1}^{(2^{Q})^{L}}\left(\prod_{l=1}^{L}\hat{B}_{l}^{W^{j}_{l}}\right)\exp{\left(i\sum_{l=1}^{L}\vec{m}_{l}^{W^{j}_{l}}\cdot\vec{x}\right)}=\sum_{j=1}^{2^{LQ}}\hat{C}_{j}e^{i\vec{M}_{j}\cdot\vec{x}},$
(26)

where we have defined the expansion coefficients
$\hat{C}_{j}=\hat{C}_{j}(\alpha,\beta)=\prod_{l=1}^{L}\hat{B}_{l}^{W^{j}_{l}}$,
the Fourier frquencies
$\vec{M}_{j}=\vec{M}_{j}(\vec{\omega})=\sum_{l}\vec{m}_{l}^{W^{j}_{l}}$ and
the indices $W^{j}_{l}$ are defined in a similar fashion as $P$; we take
$W(L,Q)$ as the set of ordered combinations of $L$ elements which take values
in $\\{1,2,...,2^{Q}\\}$. For $L=3$ and $Q=2$,

$W(3,2)=\\{\\{1,1,1\\},\\{1,1,2\\},\\{1,2,1\\},...,\\{4,4,3\\},\\{4,4,4\\}\\}.$
(27)

We denote with $W^{j}$ the $j-th$ element in $W$, and $W^{j}_{l}$ the $l$-th
element in $W^{j}$. For the previous example $W^{2}=\\{1,1,2\\}$ and
$W^{2}_{1}=1$. Alternatively, defining the element
$X_{l,q}^{j}=P_{q}^{W^{j}_{l}}$, which takes the values $\pm$, we can write
$\hat{C}_{j}=\prod_{l=1}^{L}\bigotimes_{q=1}^{Q}\hat{T}_{l,q}[X_{l,q}^{j}]$
and
$\vec{M}_{j}(\vec{\omega})=\sum_{l=1}^{L}\sum_{q=1}^{Q}X_{l,q}^{j}\vec{\omega}_{l,q}$.
From the expansion in Eq. (26), we deduce that the total number of harmonics
in the expansion of the unitary $\hat{U}_{L,Q}(\vec{x},\vec{\theta})$ is
$\mathcal{N}_{p}=2^{LQ}$.

##  Appendix B Derivation of the output of the network

The mapping of the QNN into the output function is obtained performing a
measurement in the evolved system according to Eq. (4),

$\displaystyle
F(\vec{x},\vec{\theta})=\bra{0}\hat{U}^{\dagger}\hat{\mathcal{M}}\hat{U}\ket{0}=$
$\displaystyle\sum_{i=1}^{2^{LQ}}\sum_{j=1}^{2^{LQ}}\bra{0}\hat{C}^{\dagger}_{i}\hat{\mathcal{M}}\hat{C}_{j}\ket{0}e^{i(\vec{M}_{j}-\vec{M}_{i})\cdot\vec{x}}=\sum_{i=1}^{2^{LQ}}\sum_{j=1}^{2^{LQ}}\bra{0}\hat{C}^{\dagger}_{i}\hat{\mathcal{M}}\hat{C}_{j}\ket{0}e^{i\vec{\Omega}_{j,i}\cdot\vec{x}}.$
(28)

We want to write the double sum in a single sum in order to count the total
number of harmonics in the expansion. Some frequencies may appear more than
once, meaning that $\vec{M}_{j}-\vec{M}_{i}=\vec{M}_{k}-\vec{M}_{n}$ for some
combination of $i,j\neq k,n$. We can group together the $2^{2LQ}$ terms
arising from the double sum of Eq. (28) and analyze all the possible values
that $\vec{M}_{j}(\vec{\omega})-\vec{M}_{i}(\vec{\omega})$ can yield,

$\vec{\Omega}_{j,i}=\vec{M}_{j}-\vec{M}_{i}=\sum_{l=1}^{L}\sum_{q=1}^{Q}(X_{l,q}^{j}-X_{l,q}^{i})\vec{\omega}_{l,q}.$
(29)

Since $X_{l,q}^{j}=\\{1,-1\\}$, the difference $X_{l,q}^{j}-X_{l,q}^{i}$ can
take the values $\\{2,0,-2\\}$. For each frequency in Eq. (29), there exist
three possible coefficients, and since there are $LQ$ frequencies
$\vec{\omega}_{l,q}$, we have $3^{LQ}$ different harmonic frequencies
$\vec{\Omega}_{j,i}$. As previously, let us consider $R(D)$, with $D=LQ$, the
group of the permutations $\\{2,0,-2\\}$ of $R$ elements, such as
$R=\\{\\{2,2,...,2\\},\\{2,2,...,0\\},\\{2,2,...,-2\\},...,\\{-2,-2,...,-2\\}\\}$.
In this sense, $R^{k}$ is the $k$-th element in $R$ and $R_{d}^{k}$ the $d$-th
element in $R^{k}$, with $k=1,...,3^{D}$ and $d=1,...,D$. Therefore, we can
equivalently consider

$\vec{\Omega}_{j,i}=\sum_{q=1}^{Q}\sum_{l=1}^{L}(X_{l,q}^{j}-X_{l,q}^{i})\vec{\omega}_{l,q}\hskip
14.22636pt\longleftrightarrow\hskip
14.22636pt\vec{\Omega}_{k}=\sum_{d=1}^{D}R^{k}_{d}\vec{\omega}^{*}_{d}.$ (30)

Defining $\vec{\omega}^{*}_{d}$ through a bijection with $\vec{\omega}_{l,q}$,
such that $\vec{\omega}^{*}_{1}=\vec{\omega}_{1,1}$,
$\vec{\omega}^{*}_{2}=\vec{\omega}_{1,2}$ …
$\vec{\omega}^{*}_{D}=\vec{\omega}_{L,Q}$. This allows us to rewrite Eq. (28),
as a single sum with $3^{LQ}$ terms,

$F(\vec{x},\vec{\theta})=\sum_{k=1}^{3^{LQ}}\bra{0}\hat{A}(R^{k})\ket{0}\cdot\exp\left(i\sum_{d}R^{k}_{d}\vec{\omega}^{*}_{d}\cdot\vec{x}\right)=\sum_{k=1}^{3^{LQ}}c_{k}e^{i\vec{\Omega}_{k}\cdot\vec{x}},$
(31)

with $\hat{A}(R^{k})=\hat{A}^{k}$ an operator that comes from the combination
of $\hat{C}^{\dagger}_{i}\hat{\mathcal{M}}\hat{C}_{j}$ and is associated to
one of the $3^{LQ}$ frequencies. For example, the coefficient
$\hat{A}({\\{2,0,-2,0\\}})$ is associated with the frequency
$2\vec{\omega}^{*}_{1}-2\vec{\omega}^{*}_{3}$ (or, equivalently,
$2\vec{\omega}_{1,1}-2\vec{\omega}_{2,1}$). Some of the coefficients
$c_{k}=c_{k}(\alpha,\beta)=\bra{0}\hat{A}^{k}\ket{0}$ might be zero, and thus
the sum only contains $N_{h}\leq 3^{LQ}$. Our goal in the next section is to
determine the number of non-vanishing coefficients depending on the operator
that is being used to map the QNN to the output function.

##  Appendix C Counting non-zero expansion coefficients for QNNs using a
single qubit

In this section of the appendix, we start by computing the matrices
$\hat{A}(R^{k})$ from Eq. (31) for a QNN with just one qubit and $L$ layers.
The strategy consists in adding a new layer to a circuit with $L$ layers,
$\hat{U}_{L}\rightarrow\hat{U}_{L+1}=\hat{U}_{L}\hat{u}_{l}$. For a single
qubit network, we simply have
$\hat{u}_{l}=(\hat{T}_{l}^{+}e^{i\vec{\omega}_{l}\cdot\vec{x}}+\hat{T}_{l}^{-}e^{-i\vec{\omega}_{l}\cdot\vec{x}})$.
We analyze the action of adding this extra layer,

$\displaystyle\hat{U}_{L}^{\dagger}\hat{\mathcal{M}}\hat{U}_{L}=\sum_{k=1}^{3^{L}}\hat{A}_{L}^{k}e^{i\vec{\Omega}_{k}\vec{x}}\rightarrow$
$\displaystyle\hat{u}_{l}^{\dagger}\hat{U}_{L}^{\dagger}\hat{\mathcal{M}}\hat{U}_{L}\hat{u}_{l}=\sum_{k=1}^{3^{L}}\left[(\hat{T}_{l}^{+})^{\dagger}e^{-i\vec{\omega}_{l}\cdot\vec{x}}+(\hat{T}_{l}^{-})^{\dagger}e^{i\vec{\omega}_{l}\cdot\vec{x}}\right]\hat{A}_{L}^{k}\left[\hat{T}_{l}^{+}e^{i\vec{\omega}_{l}\cdot\vec{x}}+\hat{T}_{l}^{-}e^{-i\vec{\omega}_{l}\cdot\vec{x}}\right]e^{i\vec{\Omega}_{k}\cdot\vec{x}}=$
$\displaystyle\sum_{k=1}^{3^{L}}\bigg{\\{}\left[(\hat{T}_{l}^{+})^{\dagger}\hat{A}^{k}_{L}\hat{T}^{+}_{l}+(\hat{T}_{l}^{-})^{\dagger}\hat{A}^{k}_{L}\hat{T}^{-}_{l}\right]e^{i\vec{\Omega}_{k}\cdot\vec{x}}+(\hat{T}_{l}^{-})^{\dagger}\hat{A}^{k}_{L}\hat{T}^{+}_{l}e^{i(\vec{\Omega}_{k}+2\vec{\omega}_{l})\cdot\vec{x}}+(\hat{T}_{l}^{+})^{\dagger}\hat{A}^{k}_{L}\hat{T}^{-}_{l}e^{i(\vec{\Omega}_{k}-2\vec{\omega}_{l})\cdot\vec{x}}\bigg{\\}}=$
$\displaystyle\sum_{k=1}^{3^{L}}\sum_{a\in\\{0,2,-2\\}}\mathcal{H}_{l}^{a}[\hat{A}^{k}_{L}]e^{i(\vec{\Omega}_{k}+a\vec{\omega}_{l})\cdot\vec{x}},$
(32)

Where we have introduced the superoperators

$\displaystyle\mathcal{H}_{l}^{0}[\hat{O}]=(\hat{T}^{+}_{l})^{\dagger}\hat{O}\hat{T}^{+}_{l}+(\hat{T}^{-}_{l})^{\dagger}\hat{O}\hat{T}^{-}_{l},$
$\displaystyle\mathcal{H}_{l}^{+2}[\hat{O}]=(\hat{T}^{-}_{l})^{\dagger}\hat{O}\hat{T}^{+}_{l},$
$\displaystyle\mathcal{H}_{l}^{-2}[\hat{O}]=(\hat{T}^{+}_{l})^{\dagger}\hat{O}\hat{T}^{-}_{l}.$
(33)

From the derivation in Eq. (32), we observe that adding a layer to a single
qubit network adds three new frequencies to the Fourier spectrum:
$\\{\Omega_{k}\\}\rightarrow\\{\Omega_{k},\Omega_{k}+2\omega_{l},\Omega_{k}-2\omega_{l}\\}$.
The coefficients of the new network are directly obtained by applying the
operators $\mathcal{H}_{l}^{a}$ to the coefficients $A^{k}_{L}$, resulting in
a new set of coefficients
$\\{A^{k}_{L}\\}\rightarrow\\{A^{k}_{L+1}\\}=\\{\mathcal{H}_{l}^{0}[A^{k}_{L}],\mathcal{H}_{l}^{+2}[A^{k}_{L}],\mathcal{H}_{l}^{-2}[A^{k}_{L}]\\}$.
By defining the operation of adding an element $a\in\\{0,2,-2\\}$ to a list
$R^{k}$ such as $\\{R^{k},a\\}=\\{R^{k}_{1},R^{k}_{2},...,R_{L}^{L},a\\}$, we
can now identify the relationship between the coefficients, when considering
$l=L+1$,

$\hat{A}_{L+1}(\\{R^{k},a\\})=\mathcal{H}_{L+1}^{a}[\hat{A}_{L}(R^{k})]$ (34)

from where we can obtain any term recursively as

$\hat{A}^{k}=\hat{A}(R^{k})=\mathcal{H}_{1}^{R^{k}_{1}}[\mathcal{H}_{2}^{R^{k}_{2}}[...\mathcal{H}_{L-1}^{R^{k}_{L-1}}[[\mathcal{H}^{R^{k}_{L}}_{L}[\hat{\mathcal{M}}]]]...]].$
(35)

The result of Eq. (35) is very powerful as it allows to study the dependence
of the coefficients on the measured operator $\hat{\mathcal{M}}$ for different
layer architectures and find those that are zero. In particular, for the layer
definition of Eq. (3), we show how $\mathcal{H}_{l,q}^{a}$ acts on each
operator of the spin basis,

$\displaystyle\mathcal{H}_{l}^{0}[\hat{\sigma}_{x}]=\sin{(2\alpha_{l})}\hat{\sigma}_{z}+\cos{(2\alpha_{l})}\hat{\sigma}_{x},$
$\displaystyle\mathcal{H}_{l}^{0}[\hat{\sigma}_{y}]=0,$
$\displaystyle\mathcal{H}_{l}^{0}[\hat{\sigma}_{z}]=0,$
$\displaystyle\mathcal{H}_{l}^{0}[\mathbb{I}]=\mathbb{I},$ (36)

$\displaystyle\mathcal{H}_{l}^{2}[\hat{\sigma}_{x}]=0,$
$\displaystyle\mathcal{H}_{l}^{2}[\hat{\sigma}_{y}]=i\frac{e^{2i\beta_{l}}}{2}[\cos(2\alpha_{l})\hat{\sigma}_{z}-\sin(2\alpha_{l})\hat{\sigma}_{x}-i\hat{\sigma}_{y}],$
$\displaystyle\mathcal{H}_{l}^{2}[\hat{\sigma}_{z}]=\frac{e^{2i\beta_{l}}}{2}[\cos(2\alpha_{l})\hat{\sigma}_{z}-\sin(2\alpha_{l})\hat{\sigma}_{x}-i\hat{\sigma}_{y}],$
$\displaystyle\mathcal{H}_{l}^{2}[\mathbb{I}]=0,$ (37)

$\displaystyle\mathcal{H}_{l}^{-2}[\hat{\sigma}_{x}]=0,$
$\displaystyle\mathcal{H}_{l}^{-2}[\hat{\sigma}_{y}]=-i\frac{e^{2i\beta_{l}}}{2}[\cos(2\alpha_{l})\hat{\sigma}_{z}+\sin(2\alpha_{l})\hat{\sigma}_{x}+i\hat{\sigma}_{y}],$
$\displaystyle\mathcal{H}_{l}^{-2}[\hat{\sigma}_{z}]=\frac{e^{2i\beta_{l}}}{2}[\cos(2\alpha_{l})\hat{\sigma}_{z}+\sin(2\alpha_{l})\hat{\sigma}_{x}+i\hat{\sigma}_{y}],$
$\displaystyle\mathcal{H}_{l}^{-2}[\mathbb{I}]=0.$ (38)

If we use the observable $\hat{\mathcal{M}}=\hat{\sigma}_{z}$, then all the
terms associated with sequences such that $R_{L}^{k}=0$, i.e.,
$R^{k}=\\{R^{k}_{1},...,R_{L-1}^{k},0\\}$, return a coefficient equal to zero
due to $\mathcal{H}_{l}^{0}[\hat{\sigma}_{z}]=0$. The recursive relation gives

$\displaystyle\hat{A}(\\{R^{k}_{1},...,R_{L-1}^{k},0\\})=\mathcal{H}_{1}^{R^{k}_{1}}[\mathcal{H}_{2}^{R^{k}_{2}}[...\mathcal{H}_{L-1}^{R^{k}_{L-1}}[\mathcal{H}^{0}_{L}[\hat{\sigma}_{z}]]...]]$
$\displaystyle=\mathcal{H}_{1}^{R^{k}_{1}}[\mathcal{H}_{2}^{R^{k}_{2}}[...\mathcal{H}_{L-1}^{R^{k}_{L-1}}[0]...]]=0$
(39)

In this case, from the $3^{L}$ sequences that we initially had, one third of
them satisfy $R^{k}_{L}=0$. Therefore, the total number of non-zero harmonic
is $N_{h}=2\chi/3$, where we have defined the maximum number of harmonics as
$\chi=\max(N_{h}(L=1))=3^{L}$. Alternatively, if we choose the observable
$\hat{\mathcal{M}}=\hat{\sigma}_{x}$, only sequences with terms $R^{k}_{L}=0$
give non-vanishing terms, reaching $\chi/3$ harmonics, since
$\mathcal{H}_{l}^{2}[\hat{\sigma}_{x}]=\mathcal{H}_{l}^{-2}[\hat{\sigma}_{x}]=0$.
The harmonics that are obtained with $\hat{\mathcal{M}}=\hat{\sigma}_{x}$ are
precisely those that cannot be obtained with
$\hat{\mathcal{M}}=\hat{\sigma}_{z}$. Thus, using a general operator
$\hat{O}=u_{o}\mathbb{I}+u_{x}\hat{\sigma}_{x}+u_{y}\hat{\sigma}_{y}+u_{z}\hat{\sigma}_{z}$,
we can potentially obtain $N_{h}=\chi$ harmonics.

The layer definition also plays an important role on the number of harmonics.
For example, if we define the layer as

$\hat{\mathcal{L}}^{*}_{l}(\vec{x},\vec{\theta}_{l})=R_{x}(2\alpha_{l})R_{y}(2(\vec{\omega}_{l}\cdot\vec{x}+\beta_{l})),$
(40)

the new superoperators $\mathcal{H}^{*}$ act such that

$\mathcal{H}_{l}^{*0}[\hat{\sigma}_{y}]\propto\mathcal{H}_{l}^{*0}[\mathcal{H}_{l}^{*\pm
2}[\hat{O}]]\propto\hat{\sigma}_{y}.$ (41)

Therefore, we can see that any sequence ending with an array of zeros after a
$\pm 2$ term, such as $R^{k}=\\{...,2,0\\}$ or $R^{k}=\\{...,-2,0,0,0\\}$,
will make $\hat{A}(R^{k})\propto\hat{\sigma}_{y}$. This implies that any
sequence ending with a zero necessarily result in a $\hat{\sigma}_{y}$ term,
except the sequence containing all zeros. Since
$\bra{0}\hat{\sigma}_{y}\ket{0}=0$, the only non-zero coefficients are those
for which $R^{k}_{L}$ takes the value $\\{2,-2\\}$. The total number of
harmonics in this case is $N^{h}=2\chi/3+1$, after including the sequence
containing all zeros. Now, with the application of
$\hat{\mathcal{L}}^{*}_{l}(\vec{x},\vec{\theta}_{l})$, the total number of
harmonics is independent of the measured spin operator, which is a totally
different behavior from the case previously explored.

##  Appendix D Counting non-zero expansion coefficients for QNNs using
multiple qubits

We can compute the superoperators $\mathcal{H}_{l}$ for more qubits. When
adding a new layer to an existing network with Q qubits, the number of
harmonics gets multiplied by $3^{Q}$. Thus, the number of required
superoperators scales as $3^{Q}$, one for each new frequency, meaning that,
for more qubits, the problem becomes more complicated. However, the structure
of these superoperators is predictable and allows to introduce the role of
entanglement. We can follow the same procedure used to derive Eq. 32, and we
find

$\hat{u}_{l}^{\dagger}\hat{U}_{L}^{\dagger}\hat{\mathcal{M}}\hat{U}_{L}\hat{u}_{l}=\sum_{k=1}^{3^{L}}\sum_{a_{1}\in\\{0,2,-2\\}}\sum_{a_{2}\in\\{0,2,-2\\}}...\sum_{a_{Q}\in\\{0,2,-2\\}}\mathcal{H}_{l}^{a_{1},a_{2},...,a_{Q}}[\hat{A}^{k}_{L}]e^{i(\vec{\Omega}_{k}+\sum_{q=1}^{Q}a_{q}\vec{\omega}_{l,q})\cdot\vec{x}}.$
(42)

Let us take an example with two qubits

$\displaystyle\mathcal{H}^{0,0}_{l}[\hat{O}]=\sum_{i={\pm}}\sum_{j={\pm}}(\hat{T}_{l,2}^{i}\otimes\hat{T}_{l,1}^{j})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{i}\otimes\hat{T}_{l,1}^{j}),$
$\displaystyle\mathcal{H}^{0,2}_{l}[\hat{O}]=\sum_{i={\pm}}(\hat{T}_{l,2}^{i}\otimes\hat{T}_{l,1}^{-})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{i}\otimes\hat{T}_{l,1}^{+}),$
$\displaystyle\mathcal{H}^{0,-2}_{l}[\hat{O}]=\sum_{i={\pm}}(\hat{T}_{l,2}^{i}\otimes\hat{T}_{l,1}^{+})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{i}\otimes\hat{T}_{l,1}^{-}),$
$\displaystyle\mathcal{H}^{2,0}_{l}[\hat{O}]=\sum_{i={\pm}}(\hat{T}_{l,2}^{-}\otimes\hat{T}_{l,1}^{i})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{+}\otimes\hat{T}_{l,1}^{i}),$
$\displaystyle\mathcal{H}^{2,2}_{l}[\hat{O}]=(\hat{T}_{l,2}^{-}\otimes\hat{T}_{l,1}^{-})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{+}\otimes\hat{T}_{l,1}^{+}),$
$\displaystyle\mathcal{H}^{2,-2}_{l}[\hat{O}]=(\hat{T}_{l,2}^{-}\otimes\hat{T}_{l,1}^{+})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{+}\otimes\hat{T}_{l,1}^{-}),$
$\displaystyle\mathcal{H}^{-2,0}_{l}[\hat{O}]=\sum_{i={\pm}}(\hat{T}_{l,2}^{+}\otimes\hat{T}_{l,1}^{i})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{-}\otimes\hat{T}_{l,1}^{i}),$
$\displaystyle\mathcal{H}^{-2,2}_{l}[\hat{O}]=(\hat{T}_{l,2}^{+}\otimes\hat{T}_{l,1}^{-})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{-}\otimes\hat{T}_{l,1}^{+}),$
$\displaystyle\mathcal{H}^{-2,-2}_{l}[\hat{O}]=(\hat{T}_{l,2}^{+}\otimes\hat{T}_{l,1}^{+})^{\dagger}(\hat{E}_{2}^{l})^{\dagger}\hat{O}\hat{E}_{2}^{l}(\hat{T}_{l,2}^{-}\otimes\hat{T}_{l,1}^{-}).$
(43)

We can build the operator

$\hat{A}(R^{k})=\mathcal{H}_{1}^{R^{k}_{1},R^{k}_{2}}[\mathcal{H}_{2,2}^{R^{k}_{3},R^{k}_{4}}[...\mathcal{H}_{L-1}^{R^{k}_{D-3},R^{k}_{D-2}}[\mathcal{H}_{L}^{R^{k}_{D-1},R^{k}_{D}}[\hat{\mathcal{M}}]]...]].$
(44)

In particular, considering the layer architecture given in Eq. (3), and
introducing $\hat{E}_{2}^{l}$ in all layers, we measure first the operator
$\hat{\mathcal{M}}=\hat{\sigma}_{z}^{tot}=\hat{\sigma}_{z}^{1}+\sigma_{z}^{2}$.
In this case,
$\mathcal{H}^{0,2}_{l}[\hat{\sigma}_{z}^{tot}]=\mathcal{H}^{0,-2}_{l}[\hat{\sigma}_{z}^{tot}]=\mathcal{H}^{2,0}_{l}[\hat{\sigma}_{z}^{tot}]=\mathcal{H}^{-2,0}_{l}[\hat{\sigma}_{z}^{tot}]=0$.
Any sequence ending with one and only one zero in their two last terms, such
as $R^{k}=\\{...,a_{D-3},0,2\\}$ or $R^{k}=\\{...,a_{D-3},-2,0\\}$, yields a
harmonic with zero amplitude, and as this happens for four out of nine cases,
the total number of harmonics is $N_{h}\leq 4\chi^{2}/9$. If we use the
operator $\hat{\mathcal{M}}=\hat{\sigma}_{z}^{1}\otimes\hat{\sigma}_{z}^{2}$
instead, the analysis yields $N_{h}\leq 5\chi^{2}/9$ harmonics. The maximum
scaling can be achieved by measuring a general operator
$\hat{O}_{1}\otimes\hat{O}_{2}$.

We can generalize the expressions from Eq. (43) to multiple qubits by
analyzing the three possible values of $a_{q}$ in
$\mathcal{H}_{l}^{...,a_{q},...}:$:

$\displaystyle\mathcal{H}_{l}^{...,0,...}[\hat{O}]=(...\otimes(\hat{T}_{l,q}^{+})^{\dagger}\otimes...)(\hat{E}^{l})^{\dagger}\hat{O}\hat{E}^{l}(...\otimes\hat{T}_{l,q}^{+}\otimes...)+(...\otimes(\hat{T}_{l,q}^{-})^{\dagger}\otimes...)(\hat{E}^{l})^{\dagger}\hat{O}\hat{E}^{l}(...\otimes\hat{T}_{l,q}^{-}\otimes...)=$
$\displaystyle=\sum_{i=\pm}(...\otimes(\hat{T}_{l,q}^{i})^{\dagger}\otimes...)(\hat{E}^{l})^{\dagger}\hat{O}\hat{E}^{l}(...\otimes\hat{T}_{l,q}^{i}\otimes...)$
$\displaystyle\mathcal{H}_{l}^{...,2,...}[\hat{O}]=(...\otimes(\hat{T}_{l,q}^{-})^{\dagger}\otimes...)(\hat{E}^{l})^{\dagger}\hat{O}\hat{E}^{l}(...\otimes\hat{T}_{l,q}^{+}\otimes...),$
$\displaystyle\mathcal{H}_{l}^{...,-2,...}[\hat{O}]=(...\otimes(\hat{T}_{l,q}^{+})^{\dagger}\otimes...)(\hat{E}^{l})^{\dagger}\hat{O}\hat{E}^{l}(...\otimes\hat{T}_{l,q}^{-}\otimes...).$
(45)

With these definitions, we can derive similar equations to Eq. (44) and
analytically investigate the structure of the coefficients $c_{j}$ for any
architecture.

##  Appendix E The teacher output

###  E.1 A particular example

In this section, we discuss in more detail the results presented in Fig. 5,
where different $\mathcal{A}_{1}^{S}(4,L)$ students are trained for five
different $\mathcal{A}_{1}^{T}(4,8)$ models. In Fig. 6a, we present the
prediction map $F(\vec{x},\vec{\theta})$ of one of these
$\mathcal{A}_{1}^{T}(4,8)$ teachers. The teacher prediction map is used as the
target function $f(\vec{x})$ for the outputs $F(\vec{x},\vec{\theta})$ of the
students. For the case $L=8$, we present in Figs. 6b and 6c the predictions of
two S-models, $\mathcal{A}_{0}^{S}(4,8)$ and $\mathcal{A}_{1}^{S}(4,8)$,
respectively, both initialized with the same set of parameters $\vec{\theta}$.

Alternatively to the values of the cost function presented in Fig. 5, these
prediction maps allow to study the predictions locally, and understand which
patterns and structures the models are able to learn. The student without
entanglement is able to reproduce the overall pattern of the teacher. However,
it fails when replicating smaller details, which we associate to higher
frequencies. On the contrary, $\mathcal{A}_{1}^{S}$ makes a better work at
approximating the teacher. Note that, in this case, both the teacher
$\mathcal{A}_{1}^{T}(4,8)$ and student $\mathcal{A}_{1}^{S}(4,8)$ models are
the same. However, the optimizer gets stuck in a local minimum during the
minimization process, leading to small discrepancies
$f(\vec{x})-F(\vec{x},\vec{\theta})$, see Fig. 6d.

Figure 6: (a) Bi-dimensional output generated by a teacher
$\mathcal{A}_{1}^{T}(4,8)$. The contourplot corresponds to the target function
to be approximated by the different student models. (b) Corresponging
approximation generated by a student $\mathcal{A}_{0}^{S}(4,8)$ without
entanglement and (c) by $\mathcal{A}_{1}^{S}(4,8)$ with entanglement. (d)
Difference $f(\vec{x})-F(\vec{x},\vec{\theta})$ for both student models.

###  E.2 Target function: the role of entanglement

Figure 7: Mean loss error $\bar{\mathcal{C}}$ as a function of the number of
layers $L$ for different student models $\mathcal{A}_{0}^{S}(L,Q)$ (red
dashed-lines) and $\mathcal{A}_{1}^{S}(L,Q)$ (blue solid-lines). The target
function is prepared by T-models $\mathcal{A}_{0}^{T}(L_{Max},Q)$ without and
($\mathcal{A}_{1}^{T}(L_{Max},Q)$ with) entanglement, top (bottom) panels.
From left to right, the teacher QNN contains $Q=2$, $Q=3$, and $Q=4$ qubits.
Same parameters as in Fig. 5.

We have analyzed in the main text the capabilities of both
$\mathcal{A}_{0}^{S}$ and $\mathcal{A}_{1}^{S}$ models to approximate a target
function generated by an entangling QNN. Here, for completeness, we analyze
the case of data generated by non-entangling models,
$\mathcal{A}_{0}^{T}(L_{Max},Q)$, and analyze the differences with the
previous case. According to E.1, the role of entanglement in the T-models is
to provide a more complex target function to be approximated, and thus
increases the demands of the students. For the
$\mathcal{A}_{0}^{T}(L_{Max},Q)$ models, see (a)-(c) in Fig. 7, we observe a
small difference in the global performance of both $\mathcal{A}_{0}^{S}$ and
$\mathcal{A}_{1}^{S}$. In such cases, the overall accuracy of the entangling
S-model beats the non-entangling QNN, except when simple networks (low number
of layers $L$) are considered. This may be attribute to an overestimation of
$f(\vec{x})$ due to the non-tunability of a larger set of coefficients $c_{j}$
in the parial Fourier expansion (13). Note that increasing the complexity of
the QNN (increasing $Q$) displaces the $\bar{\mathcal{C}}$ plots crossing
point to the left, indicating that small entangling networks
$\mathcal{A}_{1}^{S}$ surpass the $\mathcal{A}_{0}^{S}$ accuracy earlier. We
also observe a faster saturation of the accuracy $\bar{\mathcal{C}}$ in
$\mathcal{A}_{0}^{S}$ compared to $\mathcal{A}_{1}^{S}$. Although both QNN
have the same restricted number $N_{p}=LQ(2+n)$ of adjustable parameters and
thus accessibility to the expansion Fourier basis elements
$e^{i\vec{\Omega}_{j}\cdot\vec{x}}$, the larger spectra $\Omega_{j}$ of
$\mathcal{A}_{1}^{S}$ compared to $\mathcal{A}_{0}^{S}$ reduces the
missestimation of $f(\vec{x})$ by the non-tunable expansion coefficients
$c_{j}$. As we see at the bottom panels (d)-(f)

## References

  * [1] Kak S C 1995 Quantum neural computing Quantum Neural Computing (Advances in Imaging and Electron Physics vol 94) ed Hawkes P W (Elsevier) pp 259–313 URL https://www.sciencedirect.com/science/article/pii/S1076567008701472
  * [2] Biamonte J, Wittek P, Pancotti N, Rebentrost P, Wiebe N and Lloyd S 2017 Nature 549 195–202 URL https://doi.org/10.1038/nature23474
  * [3] Schuld M and Killoran N 2019 Phys. Rev. Lett. 122(4) 040504 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.040504
  * [4] Carrasquilla J and Melko R G 2017 Nature Physics 13 431–434
  * [5] Deng D L, Li X and Das Sarma S 2017 Phys. Rev. X 7(2) 021021 URL https://link.aps.org/doi/10.1103/PhysRevX.7.021021
  * [6] Torlai G, Mazzola G, Carrasquilla J and et al 2018 Nature Physics 14 447–450
  * [7] Fösel T, Tighineanu P, Weiss T and Marquardt F 2018 Phys. Rev. X 8(3) 031084 URL https://link.aps.org/doi/10.1103/PhysRevX.8.031084
  * [8] Aharon N, Rotem A, McGuinness L and et al 2019 Sci. Rep. 9 17802
  * [9] Gupta S and Zia R 2001 Journal of Computer and System Sciences 63 355–383 ISSN 0022-0000 URL https://www.sciencedirect.com/science/article/pii/S0022000001917696
  * [10] Paparo G D, Dunjko V, Makmal A, Martin-Delgado M A and Briegel H J 2014 Phys. Rev. X 4(3) 031002 URL https://link.aps.org/doi/10.1103/PhysRevX.4.031002
  * [11] Benedetti M, Realpe-Gómez J, Biswas R and Perdomo-Ortiz A 2017 Phys. Rev. X 7(4) 041052 URL https://link.aps.org/doi/10.1103/PhysRevX.7.041052
  * [12] Sentís G, Guta M and Adesso G 2015 EPJ Quantum Technology 2 17
  * [13] Havlícek V, Corcóles A D, Temme K, Harrow A W and et al 2019 Nature 567(7747) 209
  * [14] Abba A, D Sutter C Z, Lucchi A, Figalli A and Woerner S 2021 Nat Comput Sci 1 403–409
  * [15] Qian Y, Wang X, Du Y, Wu X and Tao D 2021 The dilemma of quantum neural networks URL https://arxiv.org/abs/2106.04975
  * [16] Schuld M and Killoran N 2022 Is quantum advantage the right goal for quantum machine learning? URL https://arxiv.org/abs/2203.01340
  * [17] Funcke L, Hartung T, Jansen K, Kühn S and Stornati P 2021 Quantum 5 422 ISSN 2521-327X URL https://doi.org/10.22331/q-2021-03-29-422
  * [18] Beer K, List D, Müller G, Osborne T J and Struckmann C 2021 Training quantum neural networks on nisq devices URL https://arxiv.org/abs/2104.06081
  * [19] Gratsea A and Humbeli P 2022 Quantum Mach. Intell. 4(2) 15
  * [20] Wilkinson S A and Hartmann M J 2022 Evaluating the performance of sigmoid quantum perceptrons in quantum neural networks URL https://arxiv.org/abs/2208.06198
  * [21] Casas B and Cervera-Lierta A 2023 Phys. Rev. A 107(6) 062612 URL https://link.aps.org/doi/10.1103/PhysRevA.107.062612
  * [22] Cao Y, Guerreschi G G and Aspuru-Guzik A 2017 Quantum neuron: an elementary building block for machine learning on quantum computers URL https://arxiv.org/abs/1711.11240
  * [23] Farhi E and Neven H 2018 Classification with quantum neural networks on near term processors URL https://arxiv.org/abs/1802.06002
  * [24] Torrontegui E and García-Ripoll J J 2019 EPL (Europhysics Letters) 125 30004 URL https://doi.org/10.1209/0295-5075/125/30004
  * [25] Pérez-Salinas A, Cervera-Lierta A, Gil-Fuster E and Latorre J I 2020 Quantum 4 226 ISSN 2521-327X URL https://doi.org/10.22331/q-2020-02-06-226
  * [26] Mangini S, Tacchino F, Gerace D, Macchiavello C and Bajoni D 2020 Machine Learning: Science and Technology 1 045008 URL https://doi.org/10.1088/2632-2153/abaf98
  * [27] Ban Y, Chen X, Torrontegui E, Solano E and Casanova J 2021 Sci. Rep. 11 5783
  * [28] Ban Y, Torrontegui E and Casanova J 2023 Sci. Rep. 13 9096
  * [29] Huber P, Haber J, Barthel P, García-Ripoll J J, Torrontegui E and Wunderlich C 2021 URL https://arxiv.org/abs/2111.08977
  * [30] Dutta T, Pérez-Salinas A, Cheng J P S, Latorre J I and Mukherjee M 2021 Single-qubit universal classifier implemented on an ion-trap quantum device URL https://arxiv.org/abs/2106.14059
  * [31] Havlícek V, Corcóles A D, Temme K, Harrow A W and et al 2019 Nature 567 209–212
  * [32] Pechal M, Roy F, Wilkinson S A, Salis G, Werninghaus M, Hartmann M J and Filipp S 2021 Direct implementation of a perceptron in superconducting circuit quantum hardware URL https://arxiv.org/abs/2111.12669
  * [33] Moreira M S, Guerreschi G G, Vlothuizen W, van Straten J, van Someren H, Marques J, Ali H, Muthusubramanian N, Zachariadis C, Beekman M, Haider N, A Bruno C A, Matsuura A and Dicarlo L 2022 APS March Meeting Abstract T37. 009
  * [34] Cybenko G 1989 Mathematics of Control, Signals, and Systems 2 303–314 ISSN 0932-4194 URL http://link.springer.com/10.1007/BF02551274
  * [35] Nielsen M A and Chuang I L 2010 Quantum Computation and Quantum Information: 10th Anniversary Edition (Cambridge University Press)
  * [36] Benedetti M, Lloyd E, Sack S and Fiorentini M 2019 Quantum Science and Technology 4 043001 URL https://dx.doi.org/10.1088/2058-9565/ab4eb5
  * [37] Pérez-Salinas A, López-Núñez D, García-Sáez A, Forn-Díaz P and Latorre J I 2021 Phys. Rev. A 104(1) 012405 URL https://link.aps.org/doi/10.1103/PhysRevA.104.012405
  * [38] Schuld M, Sweke R and Meyer J J 2021 Physical Review A 103 032430
  * [39] Ben-Aryeh Y 2014 Quantum and classical correlations in bell three and four qubits, related to hilbert-schmidt decomposition URL https://arxiv.org/abs/1411.2720
  * [40] Casas B and Cervera-Lierta A 2023 arXiv preprint arXiv:2302.03389
  * [41] Zimmer M, Viappiani P and Weng P 2014 Teacher-student framework: a reinforcement learning approach AAMAS Workshop Autonomous Robots and Multirobot Systems
  * [42] Hornik K 1990 Neural Networks 4 251–257
  * [43] Barron A 1993 IEEE Transactions on Information Theory 39 930–945
  * [44] Roca-Jerat S, Román-Roche J and Zueco D 2023 Qudit machine learning (Preprint 2308.16230)
  * [45] Gao X and Duan L M 2017 Nat Commun 8
  * [46] Harney C, Pirandola S, Ferraro A and Paternostro M 2020 New Journal of Physics 22 045001 URL https://doi.org/10.1088/1367-2630/ab783d


capbtabboxtable[][]

# Modeling Animal Vocalizations through Synthesizers

Masato Hagiwara
Earth Species Project
<EMAIL_ADDRESS>
Maddie Cusimano
Earth Species Project
<EMAIL_ADDRESS>Jen-Yu Liu
Earth Species Project
<EMAIL_ADDRESS>

Figure 1: Method overview. Given an animal sound, synthesizer parameters are
optimized for reconstruction. The predicted sound can be modified (e.g. pitch-
shift, denoising, envelope modulation).

##  1 Introduction

Modeling real-world sound is a fundamental problem in the creative use of
machine learning and many other fields, including human speech processing and
bioacoustics. Transformer-based generative models [20, 6] are known to produce
realistic sound, although they have limited control and are hard to interpret.
Recently, lighter-weight models that incorporate structured modules and domain
knowledge, notably DDSP [8, 29], have been shown to produce high-quality
musical sound. However, a lack of signal-processing knowledge may hinder users
from effectively manipulating the synthesis parameters, of which there can be
over a hundred per frame.

As an alternative, we aim to use modular synthesizers, i.e., compositional,
parametric electronic musical instruments, for modeling non-music
sounds111Audio samples are at https://earthspecies.github.io/animalsynth/.
Synthesizers are lightweight and designed in part for control, which make them
a plausible candidate model. However, inferring synthesizer parameters given a
target sound, i.e., the parameter inference task [15, 11], is not trivial for
general sounds. Research utilizing modern optimization techniques has
typically focused on musical sound [30, 9, 18, 2, 4]. In this work, we
optimize a differentiable synthesizer from TorchSynth [28] in order to model,
emulate, and creatively generate animal vocalizations. We compare an array of
optimization methods, from gradient-based search to genetic algorithms, for
inferring its parameters, and then demonstrate how one can control and
interpret the parameters for modeling non-music sounds. The current work is
intended as a creative tool for artists to integrate animal sound into their
work, but we hope future work in this direction could also enable new
synthesis methods in bioacoustics [16].

##  2 Parameter Inference

Method |  Loss ($\downarrow$) |  Acc. ($\uparrow$)
---|---|---
Random search | 0.861 | 0.068
Gradient (Adam) | 1.553 | 0.027
Variational optimization | 0.891 | 0.034
Genetic algorithm | 0.633 | 0.114
Differential evolution | 0.738 | 0.068
PGPE [23] | 1.000 | 0.011
CMA-ES [12] | 0.891 | 0.027
Metropolis MCMC | 0.658 | 0.091
Bayesian optimization [31] | 0.737 | 0.034
Encoder | 2.131 | 0.034
Original | — | 0.841

Figure 2: Comparison of optimization methods

Figure 3: Plot of parameters per category. Parameters are normalized to the
[0,1] range.

The default synthesizer implemented in TorchSynth named Voice has a relatively
simple architecture with a total of 78 parameters, consisting of two voltage-
controlled oscillators (VCOs) and a noise generator, whose pitch and amplitude
envelopes are modulated by a mix of low-frequency oscillators (LFOs). For any
non-trivial synthesizer, it is extremely difficult to adjust parameters
manually to reproduce the target sound, necessitating automated optimization
techniques.

We first compare several optimization techniques for parameter inference. As
the target sound dataset, we used the animal sounds (e.g., dog, cat, birds) in
the first fold of ESC-50 [22], an environmental sound dataset. We used the
multi-scale STFT loss implemented in auraloss [25] as the objective. In order
to assess the quality of reconstruction, we fine-tuned the pretrained VGGish
[14] model on the remaining four folds of ESC-50 on the 50-class sound
classification task (which includes both animal and non-animal categories) and
used its accuracy. Intuitively, if an optimization method achieves good
reconstruction, it can also “fool” the classifier trained on real sounds.

Table 3 shows the optimization techniques and their metrics. See the
Supplemental Information for the details of individual methods. Overall,
evolutionary algorithms (the genetic algorithm and differential evolution)
worked relatively well, while gradient-based methods performed poorly, due to
the fact that the objective function is highly nonlinear and complex with
respect to synthesizer parameters. This trend is also observed for other
highly complex optimization problems [26]. Predicting parameters directly with
a neural network encoder was not successful either, possibly due to the
discrepancies between the real sounds and the training data (artificially
generated sound by TorchSynth itself).

One potential way to improve the gradient-based methods is to use a proxy
network that simulates the synthesizer, and replace the original synthesizer
with the proxy network in parameter inference [24]. We trained a MelGAN
[17]-based proxy network that takes the TorchSynth parameters as input and the
corresponding waveform as the target. We found that the reconstructions were
poor. In contrast to audio effects, converting 78 parameters to a 2-second
waveform (88,200 samples) may have its own challenges due to the large
difference in the scales of the input and output.

After inferring parameters for each sample/class, controlling and generating
new sounds are trivial. For example, one can change pitch, de-noise, and/or
modify the amplitude envelope by changing just one or few synthesizer
parameters as shown in the spectrograms in Figure 1.

##  3 Interpretation and Generation

One benefit of synthesizer parameter inference is that the fitted parameters
are interpretable. In Figure 3, we choose three animal categories (cow, pig,
cat) for which the classifier accuracy is relatively high and plot the
fundamental frequency (f0) and the duration of the note inferred by the
genetic algorithm. The sound instances cluster based on their acoustic
properties, i.e., cats make short and high-pitched sounds while cow
vocalizations are longer and lower-pitched.

Finally, as a proof of concept, we generated new cat sounds by first fitting a
Gaussian distribution over the inferred cat parameters (green squares in
Figure 3) in the unnormalized space, then generating 100 new sounds by
sampling from this distribution. The fine-tuned VGGish model correctly
classified 18 samples out of 100, demonstrating that we can generate plausible
new sounds by sampling synthesizer parameters from the fitted distribution.

## Ethical Considerations

With the progress of audio synthesis technology in recent years, synthesized
human voices could be used in harmful activities [21, 7, 19]. Similarly,
playing synthesized animal sounds in the wild could also be harmful to animals
[5, 13], meaning released models could have potential negative effects. One
way to mitigate such risks is to accompany the generative models with
corresponding detection models [3].

With the above considerations, simpler models (like the one used in this work)
could be a good alternative to more powerful generative models to be released
for creative use. Its constrained synthesis makes it difficult to synthesize
authentic sounds, while still capturing evocative aspects of the target sounds
and providing rich variations to be used in creative design.

## Acknowledgement

The authors would like to thank Benjamin Hoffman, Felix Effenberger, and Katie
Zacarian at Earth Species Project for their valuable feedback on this project.

## References

  * [1] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  arXiv, 2019.
  * [2] Oren Barkan, David Tsiris, Ori Katz, and Noam Koenigstein.  Inversynth: Deep estimation of synthesizer parameter configurations from audio signals.  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 27(12):2385–2396, 2019.
  * [3] Zalán Borsos, Raphaël Marinier, Damien Vincent, Eugene Kharitonov, Olivier Pietquin, Matt Sharifi, Olivier Teboul, David Grangier, Marco Tagliasacchi, and Neil Zeghidour.  Audiolm: a language modeling approach to audio generation, 2022.
  * [4] Zui Chen, Yansen Jing, Shengcheng Yuan, Yifei Xu, Jian Wu, and Hang Zhao.  Sound2synth: Interpreting sound via FM synthesizer parameters estimation.  In Proceedings of the Thirty-First International Joint Conference on Artificial Intelligence. International Joint Conferences on Artificial Intelligence Organization, jul 2022.
  * [5] Innes Cuthill.  Field experiments in animal behaviour: methods and ethics.  Animal Behaviour, 42(6):1007–1014, 1991.
  * [6] Prafulla Dhariwal, Heewoo Jun, Christine Payne, Jong Wook Kim, Alec Radford, and Ilya Sutskever.  Jukebox: A generative model for music.  ArXiv, abs/2005.00341, 2020.
  * [7] Nicholas Diakopoulos and Deborah Johnson.  Anticipating and addressing the ethical implications of deepfakes in the context of elections.  New Media & Society, 23(7):2072–2098, 2021.
  * [8] Jesse Engel, Lamtharn (Hanoi) Hantrakul, Chenjie Gu, and Adam Roberts.  Ddsp: Differentiable digital signal processing.  In International Conference on Learning Representations, 2020.
  * [9] Philippe Esling, Naotake Masuda, Adrien Bardet, Romeo Despres, and Axel Chemla-Romeu-Santos.  Flow synthesizer: Universal audio synthesizer control with normalizing flows.  Applied Sciences, 10(1), 2020.
  * [10] Ahmed Fawzy Gad.  Pygad: An intuitive genetic algorithm python library, 2021.
  * [11] Ricardo A. Garcia.  Automatic design of sound synthesis techniques by means of genetic programming.  Journal of the Audio Engineering Society, october 2002.
  * [12] Nikolaus Hansen and Andreas Ostermeier.  Completely derandomized self-adaptation in evolution strategies.  Evolutionary Computation, 9(2):159–195, 2001.
  * [13] J. Berton C. Harris and David G. Haskell.  Simulated birdwatchers’ playback affects the behavior of two tropical birds.  PLOS ONE, 8(10), 10 2013.
  * [14] Shawn Hershey, Sourish Chaudhuri, Daniel P. W. Ellis, Jort F. Gemmeke, Aren Jansen, R. Channing Moore, Manoj Plakal, Devin Platt, Rif A. Saurous, Bryan Seybold, Malcolm Slaney, Ron J. Weiss, and Kevin Wilson.  Cnn architectures for large-scale audio classification.  In 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 131–135, 2017.
  * [15] J. Justice.  Analytic signal processing in music computation.  IEEE Transactions on Acoustics, Speech, and Signal Processing, 27(6):670–684, 1979.
  * [16] Stephanie L. King.  You talkin’ to me? interactive playback is a powerful yet underused tool in animal communication research.  Biology Letters, 11, 2015.
  * [17] Kundan Kumar, Rithesh Kumar, Thibault de Boissiere, Lucas Gestin, Wei Zhen Teoh, Jose Sotelo, Alexandre de Brébisson, Yoshua Bengio, and Aaron C Courville.  Melgan: Generative adversarial networks for conditional waveform synthesis.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [18] Gwendal Le Vaillant, Thierry Dutoit, and Sébastien Dekeyser.  Improving synthesizer programming from variational autoencoders latent space.  In 2021 24th International Conference on Digital Audio Effects (DAFx), pages 276–283. IEEE, 2021.
  * [19] Kyungyun Lee, Gladys Hitt, Emily Terada, and Jin Ha Lee.  Ethics of singing voice synthesis: Perceptions of users and developers.  In Proc. International Society for Music Information Retrieval Conference, 2022.
  * [20] Naihan Li, Shujie Liu, Yanqing Liu, Sheng Zhao, and Ming Liu.  Neural speech synthesis with transformer network.  In Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence and Thirty-First Innovative Applications of Artificial Intelligence Conference and Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, AAAI’19/IAAI’19/EAAI’19. AAAI Press, 2019.
  * [21] Edvinas Meskys, Aidas Liaudanskas, Julija Kalpokiene, and Paulius Jurcys.  Regulating deep fakes: legal and ethical considerations.  Journal of Intellectual Property Law & Practice, 15(1):24–31, 01 2020.
  * [22] Karol J. Piczak.  ESC: Dataset for Environmental Sound Classification.  In Proceedings of the 23rd Annual ACM Conference on Multimedia, pages 1015–1018. ACM Press, 2015.
  * [23] Frank Sehnke, Christian Osendorfer, Thomas Rückstieß, Alex Graves, Jan Peters, and Jürgen Schmidhuber.  Parameter-exploring policy gradients.  Neural Networks, 23(4):551–559, 2010.  The 18th International Conference on Artificial Neural Networks, ICANN 2008.
  * [24] Christian J. Steinmetz.  Learning to mix with neural audio effects in the waveform domain.  PhD thesis, Universitat Pompeu Fabra, September 2020.
  * [25] Christian J. Steinmetz and Joshua D. Reiss.  auraloss: Audio focused loss functions in PyTorch.  In Digital Music Research Network One-day Workshop (DMRN+15), 2020\.
  * [26] Yingtao Tian and David Ha.  Modern evolution strategies for creativity: Fitting concrete images and abstract concepts.  In Artificial Intelligence in Music, Sound, Art and Design: 11th International Conference, EvoMUSART 2022, Held as Part of EvoStar 2022, Madrid, Spain, April 20–22, 2022, Proceedings, page 275–291, Berlin, Heidelberg, 2022. Springer-Verlag.
  * [27] Nihat Engin Toklu, Paweł Liskowski, and Rupesh Kumar Srivastava.  Clipup: A simple and powerful optimizer for distribution-based policy evolution.  In International Conference on Parallel Problem Solving from Nature, pages 515–527. Springer, 2020.
  * [28] Joseph Turian, Jordie Shier, George Tzanetakis, Kirk McNally, and Max Henry.  One billion audio sounds from GPU-enabled modular synthesis.  In Proceedings of the 23rd International Conference on Digital Audio Effects (DAFx2020), September 2021.
  * [29] Yusong Wu, Ethan Manilow, Yi Deng, Rigel Swavely, Kyle Kastner, Tim Cooijmans, Aaron Courville, Cheng-Zhi Anna Huang, and Jesse Engel.  MIDI-DDSP: Detailed control of musical performance via hierarchical modeling.  In International Conference on Learning Representations, 2022.
  * [30] Matthew John Yee-King, Leon Fedden, and Mark d’Inverno.  Automatic programming of vst sound synthesizers using deep networks and other techniques.  IEEE Transactions on Emerging Topics in Computational Intelligence, 2(2):150–159, 2018.
  * [31] Ozaki Yoshihiko, Tanigaki Yuki, Watanabe Shuhei, and Onishi Masaki.  Multiobjective tree-structured parzen estimator for computationally expensive optimization problems.  In The Genetic and Evolutionary Computation Conference (GECCO2020), 2020.

## Supplemental Information

### Parameter Inference

Model. We used the default synthesizer implemented in TorchSynth, named Voice,
with the default ‘voice’ nebula. See
https://torchsynth.readthedocs.io/en/latest/modular-design/modular-
principles.html##synth-architectures for details.

Target dataset. ESC-50 [22] contains ten animal sound classes: dog, rooster,
pig, cow, frog, cat, hen, insects flying (e.g., buzzing), sheep, crow, and
birds. To create the target sound dataset, we took the subset of sounds in the
first fold that were in these ten classes. All sounds were two seconds and
sampled at 44.1 kHz.

Reconstruction loss. We used a mel-scaled multi-scale short-time Fourier
transform (STFT) loss from the auraloss package [25], with equally weighted
linear-scale magnitude and spectral convergence loss terms. Four resolutions
had FFT, hop and window sizes corresponding to (2048, 1024, 512, 256), (256,
128, 64, 32) and (1024, 512, 256, 128) respectively. All resolutions used 45
mel-frequency bins.

### Optimization

Synthesizer parameter were initialized with uniform random samples, unless
otherwise noted. Hyperparameters were manually chosen.

Random search. We took 1000 random uniform samples from all 78 parameters. For
each sound, we selected the sample that minimized the reconstruction loss.

Gradient (Adam). We used 200 iterations of an adaptive gradient descent
algorithm (Adam) on all 78 parameters, with a learning rate of 0.001. For each
sound, we selected the parameters on the final iteration.

Variational optimization. The original 78 parameters of the modular
synthesizer were used in conjunction with an additional 78 variance parameters
to define a Beta distribution over each synthesizer parameter. The
distributions were initialized to uniform and then optimized with Adam for 500
iterations and a learning rate of 0.001, using a batch size of 200 samples.
For each sound, we selected the sample from the last iteration that minimized
the reconstruction loss.

Genetic algorithm. We used the pygad package [10]. For 100 iterations of
steady-state selection, four parents were selected and single-point crossovers
and random mutations (10% of parameters) were used to evolve the parameters.
For each sound, we selected the individual that minimized the reconstruction
loss.

Differential evolution. We used the default differential evolution function
from the scipy package, with a maximum of 20 generations and a population size
of 10. For each sound, we selected the individual that minimized the
reconstruction loss.

Policy gradients with parameter-based exploration (PGPE). We use PGPE
implemented in the pgpelib package [27], a distribution-based evolutionary
algorithm that estimates gradients with samples around a center solution. We
use a population size of 100 evolved over 100 generations using the ClipUp
optimizer. For each sound, we use the center solution of the final generation.

Covariance matrix adaptation Evolutionary Strategy (CMA-ES). We used CMA-ES as
implemented in the pycma package [12], with a maximum of 200 iterations. For
each sound, we use the mean solution.

Metropolis algorithm. We implemented Markov Chain Monte Carlo (MCMC) using the
Metropolis algorithm. The synthesizer was supplemented with a uniform prior
for each parameter and a Gaussian noise model ($\sigma=0.1$) for the
likelihood. On each iteration, we randomly selected parameters to be resampled
($p=0.1$) with a uniform random proposal function. We collected 10000
posterior samples. For each sound, we selected the posterior sample that
minimized the reconstruction loss.

Bayesian optimization. We used the default hyperparameter optimization in the
optuna package [1] (Tree-Structured Parzen Estimator [31]), with 1000 trials.

Encoder. We fine-tuned the pretrained VGGish [14] with a fully-connected
regression layer on top to predict the synthesizer parameters from the
spectrogram. The model is trained with an MSE loss on a dataset of 60,000
sounds generated from TorchSynth.

### Generation

The Gaussian distribution over the inferred cat parameters was fit in the
unnormalized space (keyboard.duration in sec, keyboard.midi_f0 in midi-
values).