††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work

# Linear-depth quantum circuits for loading Fourier approximations of
arbitrary functions

Mudassir Moosa<EMAIL_ADDRESS>Department of Physics, Cornell University,
Ithaca, NY 14853, USA  Department of Physics and Astronomy, Purdue University,
West Lafayette, IN 47907, USA  Thomas W. Watts<EMAIL_ADDRESS>School of
Applied and Engineering Physics, Cornell University, Ithaca, NY 14853, USA
Yiyou Chen  School of Applied and Engineering Physics, Cornell University,
Ithaca, NY 14853, USA  Department of Electrical and Computer Engineering,
Princeton University, Princeton, NJ 08544, USA  Abhijat Sarma  Department of
Physics, Cornell University, Ithaca, NY 14853, USA  Peter L. McMahon
<EMAIL_ADDRESS>School of Applied and Engineering Physics, Cornell
University, Ithaca, NY 14853, USA

###### Abstract

The ability to efficiently load functions on quantum computers with high
fidelity is essential for many quantum algorithms, including those for solving
partial differential equations and Monte Carlo estimation. In this work, we
introduce the Fourier Series Loader (FSL) method for preparing quantum states
that exactly encode multi-dimensional Fourier series using linear-depth
quantum circuits. Specifically, the FSL method prepares a ($Dn$)-qubit state
encoding the $2^{Dn}$-point uniform discretization of a $D$-dimensional
function specified by a $D$-dimensional Fourier series. A free parameter, $m$,
which must be less than $n$, determines the number of Fourier coefficients,
$2^{D(m+1)}$, used to represent the function. The FSL method uses a quantum
circuit of depth at most
$2(n-2)+\lceil\log_{2}(n-m)\rceil+2^{D(m+1)+2}-2D(m+1)$, which is linear in
the number of Fourier coefficients, and linear in the number of qubits ($Dn$)
despite the fact that the loaded function’s discretization is over
exponentially many ($2^{Dn}$) points. The FSL circuit consists of at most
$Dn+2^{D(m+1)+1}-1$ single-qubit and $Dn(n+1)/2+2^{D(m+1)+1}-3D(m+1)-2$ two-
qubit gates; we present a classical compilation algorithm with runtime
$O(2^{3D(m+1)})$ to determine the FSL circuit for a given Fourier series. The
FSL method allows for the highly accurate loading of complex-valued functions
that are well-approximated by a Fourier series with finitely many terms. We
report results from noiseless quantum circuit simulations, illustrating the
capability of the FSL method to load various continuous 1D functions, and a
discontinuous 1D function, on 20 qubits with infidelities of less than
$10^{-6}$ and $10^{-3}$, respectively. We also demonstrate the practicality of
the FSL method for near-term quantum computers by presenting experiments
performed on the Quantinuum H$1$-$1$ and H$1$-$2$ trapped-ion quantum
computers: we loaded a complex-valued function on 3 qubits with a fidelity of
over $95\%$, as well as various 1D real-valued functions on up to 6 qubits
with classical fidelities $\approx 99\%$, and a 2D function on 10 qubits with
a classical fidelity $\approx 94\%$.

##  I Introduction

Efficiently loading classical data on quantum computers is an important
ingredient in many quantum algorithms. For example, quantum-linear-solver
algorithms require an efficient preparation of a state $\ket{b}$ in order to
compute the solution $\ket{x}\propto A^{-1}\ket{b}$ to a large linear system
$A$ [1, 2, 3, 4]. In general, exactly loading $2^{n}$ dimensional data into a
state of $n$ qubits requires a quantum circuit of $O(2^{n})$ quantum
operations [5, 6, 7, 8]. Exponential space and time complexities for the read-
in component of a quantum algorithm will compromise any potential for an
exponential speed-up over comparable classical algorithms [9]. Hence, finding
an efficient method for accurately loading classical input data into quantum
computers is a problem with vast applications.

In many practical applications, classical input data takes the form of values
of a function or a distribution on a uniform discretized grid (or mesh) [10,
11, 12, 13, 14, 15, 16, 17, 18, 19]. It is well-known that a large class of
functions and distributions can be accurately approximated by a Fourier series
of only a few terms. In this work, we exploit this fact in order to load a
truncated Fourier-series approximation to a general multivariate target
function $f:[0,1]^{D}\to\mathbb{C}$. More precisely, given a truncated Fourier
series consisting of $2^{D(m+1)}$ Fourier modes,
$f_{(m)}(x_{1},\cdots,x_{D})$, that approximates the target function of $D$
variables, $f(x_{1},\cdots,x_{D})$, we present a method to exactly prepare a
quantum state of $Dn$ qubits of the form
$\ket{f_{(m)}}\,=\,\sum_{k_{1}=0}^{2^{n}-1}\cdots\sum_{k_{D}=0}^{2^{n}-1}f_{(m)}(k_{1}/2^{n},\cdots,k_{D}/2^{n})\ket{k_{1}}\otimes\cdots\otimes\ket{k_{D}}\,$.
For single variable functions (i.e., $D=1$), the target state becomes
$\ket{f_{(m)}}\,=\,\sum_{k=0}^{2^{n}-1}f_{(m)}(k/2^{n})\ket{k}\,$. This allows
us to load the target function with practically arbitrarily high fidelity so
long as a sufficient number of Fourier modes are used. We refer to our method
as the Fourier Series Loader (FSL).

The FSL method consists of three main steps: First, we find the dominant
$2^{m+1}$ Fourier coefficients for a given function. This can be done
classically with a time complexity of $O(m2^{m})$ using sparse Fourier
transform algorithms [20, 21, 22, 23]. The second step is to load the
$2^{m+1}\ll 2^{n}$ Fourier coefficients into a sparse quantum state of $n$
qubits. We achieve this using a quantum circuit of depth $O(\log(n)+2^{m})$,
which can be found in $O(8^{m})$ classical runtime, and consists of
$O(n+2^{m})$ and $O(2^{m})$ two-qubit and single-qubit gates, respectively.
This is an improvement over a generic sparse state preparation algorithm that
sports a classical runtime of $O(mn4^{m})$ and produces a circuit of depth
$O(n2^{m})$ which consists of $O(n2^{m})$ two-qubit gates and $O(m2^{m}+n)$
single-qubit gates [24]. Finally, the third step of the FSL method is to apply
the inverse quantum Fourier transform to generate a ‘real-space’
representation of the desired function from the loaded Fourier coefficients.
The full quantum circuit implementing our method is shown in Fig. (1a) and has
a depth of $O(n+2^{m})$ and consists of $O(n^{2}+2^{m})$ two-qubit and
$O(n+2^{m})$ single-qubit gates. The only error introduced by the FSL method
is due to the approximation of the given function by a truncated Fourier
series. This truncation error is determined by the number of Fourier modes
which dictates the requisite value of $m$. In particular, the infidelity
between the target state and the state prepared by the FSL method decays
exponentially with $m$.

It is worth mentioning that the FSL method is analogous to a state preparation
method based on Fourier interpolation [25, 26]. The interpolation-based method
loads the values of the target function at $2^{m+1}$ discretized points and
then interpolates to $2^{n}$ points using the quantum Fourier transform. The
FSL method, as discussed above, instead loads $2^{m+1}$ Fourier coefficients
and then produces a Fourier series representation of the target function using
the quantum Fourier transform. Despite their similarities, the FSL method has
various advantages over the interpolation-based method including, but not
limited to, superior accuracy and reduced gate count.

Figure 1: (a) The Fourier Series Loader (FSL) method loads a Fourier series
approximation of the target function into a quantum state of $n$ qubits.
Implementing the FSL method requires approximating the target function with a
Fourier series of $2^{m+1}$ terms and finding a unitary circuit $U_{c}$ to
load the $2^{m+1}$ Fourier coefficients on $m+1$ qubits, both of which can be
done on a classical computer. The integer $m$ should be chosen such that the
truncation error of the Fourier series, and hence, the infidelity of the
prepared state is within a desired threshold. Once a suitable unitary circuit
$U_{c}$ is found, the Fourier series approximation of the target function is
loaded on $n$ qubits by running the FSL circuit on a quantum computer. The FSL
circuit consists of three parts. First, the unitary $U_{c}$ loads $2^{m+1}$
Fourier coefficients in $m+1$ qubits and prepares the state given in Eq. (1).
Then, a cascade of CNOT gates entangles the remaining $n-m-1$ qubits with the
$m+1$ qubits storing the Fourier coefficients which efficiently implements the
technique of zero padding (see Eq. (2)). Finally, the inverse QFT at the end
of the circuit transforms the loaded Fourier coefficients into a state
encoding the Fourier series approximation of the target function. (b) One
possible implementation of the unitary $U_{c}$ is using a cascade of uniformly
controlled rotations, which can be used to prepare an arbitrary quantum state.
A uniformly $k$-controlled rotation operator, shown in the shaded region,
decomposes as a series of controlled rotations conditioned on each
computational basis state of $k$ control qubits. (c) Another possible
implementation of the unitary $U_{c}$ is using a circuit that can prepare an
arbitrary quantum state once the Schmidt decomposition of the target state is
computed. The unitary $A$ in the circuit loads the Schmidt coefficients
whereas unitaries $U$ and $V$ rotate the computational basis into the Schmidt
basis. (d) The FSL method loads a non-periodic function $f$ by first loading a
periodic extension $F$ on $n$ qubits and $1$ ancilla qubit. Then the desired
state is prepared either through disentangling the ancilla qubit (left) or
through measurement-based operations (right).

The FSL has several attractive features that make it appealing for loading
classical functions on quantum computers. First, we have complete control over
the accuracy of our results as it is determined by the number of Fourier
coefficients used to represent the target function. Secondly, no time-
consuming classical optimization of a parameterized quantum circuit or matrix
product state is needed unlike other comparable methods for function loading
in the literature [27, 28, 29, 25]. In fact, all the angles of rotation in our
quantum circuits can be determined from the Fourier coefficients of the
function of interest as we explain in Sec. (II). Thirdly, the FSL works
equally well for loading real-valued and complex-valued functions. More
remarkably, the FSL can be easily generalized to load a $D$-dimensional
Fourier series into a quantum state of $Dn$ qubits using a quantum circuit of
$O(n)$ depth and $O(Dn^{2})$ quantum gates. Finally, due to the low depth of
the circuit in Fig. (1a), the FSL method is capable of loading functions on
near-term noisy quantum computers as we demonstrate in Sec. (III).

The FSL method is particularly useful for quantum algorithms for solving
differential equations that require the preparation of states encoding initial
conditions [13, 11, 14, 10, 12]. Other contexts in which the FSL method is
useful include quantum image processing [30], Monte Carlo methods based on
quantum amplitude estimation algorithms [18, 17], and computing generalized
inner products [19]. These last two contexts find applications in the domain
of option pricing and risk analysis [18, 17, 19]. Though this is not the goal
of the present work, we elaborate in the supplementary Sec. (F) on how the FSL
method can be utilized to load an image into a quantum state.

The rest of this paper is organized as follows: In Sec. (II), we present the
details of the FSL method and the accompanying quantum circuit
implementations. Starting with the simplest use cases, we demonstrate how our
FSL method can be used to load periodic functions of a single variable. Then
we discuss how the FSL method can be easily modified to load non-periodic,
piece-wise discontinuous, and multivariate functions. For each case, we
demonstrate high fidelity loading of various functions using the FSL method by
quantum simulations performed using Qiskit [31]. Next, in Sec. (III), we
discuss several experiments of loading functions of one and two variables on
the Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers [32]. The results of
these experiments demonstrate that the FSL method can load functions on up to
$10$ qubits on a present-day noisy quantum computer. Finally, we conclude with
a comparison of the FSL method with other state preparation methods for
loading functions to quantum states and some possible future directions in
Sec. (IV).

##  II Methods

Our goal is to load the values of any arbitrary complex-valued function to the
amplitudes of a quantum state of $n$ qubits. For example, given a one-
dimensional function $f:[0,1]\to\mathbb{C}$, we want to prepare an $n$-qubit
state $\ket{f}=\sum_{k=0}^{2^{n}-1}\,f_{k}\ket{k}$, where $\ket{k}$ is the
$k^{\text{th}}$ computational basis state and $f_{k}=f(k/2^{n})$ is the value
of the function at $k^{\text{th}}$ grid point. Without loss of generality, we
assume that the function $f$ is such that the corresponding quantum state
representation $\ket{f}$ is normalized.

We begin by considering the simplest case of a function $f$ given by a finite
Fourier series of the form $f(x)\,=\,\sum_{k=-M}^{M}c_{k}e^{-i2\pi kx}\,$,
where $M:=2^{m}-1<2^{n-1}$. We claim that the quantum state $\ket{f}$
corresponding to this Fourier series can be prepared exactly using the quantum
circuit shown in Fig. (1a) in which $U_{c}$ denotes an arbitrary
implementation of the unitary that maps $\ket{0}^{\otimes(m+1)}$ state to
$\ket{\tilde{c}}$, where

$\displaystyle\ket{\tilde{c}}\,:=\,2^{n/2}\sum_{k=0}^{M}c_{k}\ket{k}+2^{n/2}\sum_{k=1}^{M}c_{-k}\ket{2^{m+1}-k}\,.$
(1)

The cascade of $n-m-1$ CNOT gates that follow $U_{c}$ entangle all $n$ qubits,
mapping $\ket{0}^{\otimes(n-m-1)}\otimes\ket{\tilde{c}}$ to an $n$-qubit state
$\ket{c}$, where

$\displaystyle\ket{c}\,:=\,2^{n/2}\sum_{k=0}^{M}c_{k}\ket{k}+2^{n/2}\sum_{k=1}^{M}c_{-k}\ket{2^{n}-k}\,.$
(2)

Finally, the inverse Quantum Fourier Transform (QFT) maps $\ket{c}$ to the
desired state $\ket{f}$. A more detailed derivation is provided in the
Supplementary Sec. (A).

The QFT and its inverse can be implemented with a quantum circuit of depth
$O(n)$ and with $O(n^{2})$ quantum gates. The cascade of CNOT gates in Fig.
(1a) can be implemented in $O(\log(n))$ depth. In the worst case, the
implementation of $(m+1)$-qubit unitary $U_{c}$ will require $O(2^{m})$
quantum gates and a circuit of depth $O(2^{m})$. However, note that $m$ is
entirely fixed by the number of Fourier modes in our function and does not
scale with the total number of qubits, $n$. Hence, for a Fourier series with a
fixed number of terms, the circuit in Fig. (1a) has $O(n)$ depth and has
$O(n^{2})$ quantum gates.

The problem of preparing the $n$-qubit state $\ket{f}$ has now reduced to the
problem of preparing a $(m+1)$-qubit state $\ket{\tilde{c}}$. When the target
Fourier series consists of only a few terms, the construction of an efficient
circuit that prepares $\ket{\tilde{c}}$ is usually straightforward. However,
when loading more complex functions that require many Fourier coefficients, we
propose two general methods for constructing the state $\ket{\tilde{c}}$.

### Implementations of $\boldsymbol{U_{c}}$

Our first proposal is to prepare the state $\ket{\tilde{c}}$ using a cascade
of uniformly controlled rotations as shown in Fig. (1b), which can be used to
prepare any arbitrary quantum state as was shown by Möttönen et al. [7].
Möttönen et al. proposed this circuit as an extension of a circuit of Zalka,
Grover, and Rudolph [5, 6] which only involves a cascade of uniformly
controlled $R_{y}$ gates, and hence, can only prepare states with positive,
real-valued amplitudes. Moreover, Möttönen et al. made two interesting
observations that render this circuit a practical choice for implementing
$U_{c}$. First, this circuit can be constructed using $2^{m+2}$ single-qubit
rotations and $2^{m+2}$ CNOT gates. Secondly, there exists a set of formulae
for calculating the angles of rotation for all single-qubit gates in the
circuit when provided access to the amplitudes of the target state
$\ket{\tilde{c}}$ [7]. These formulae can be evaluated on a classical computer
in $O(8^{m})$ time. (See Supplementary Sec. (B) for more details.)

Our second proposal is based on the Schmidt decomposition. We can treat
$\ket{\tilde{c}}$ as the overall quantum state of a bipartite system by
partitioning $m+1$ qubits into $(m+1)/2$ and $(m+1)/2$ qubits if $m$ is odd or
into $m/2$ and $m/2+1$ qubits if $m$ is even. Then, we can classically compute
the Schmidt decomposition of $\ket{\tilde{c}}$ and express it as
$\ket{\tilde{c}}=U\otimes V\,\sum_{k=0}^{r}\alpha_{k}\ket{k}\otimes\ket{k}\,$,
where $\alpha_{k}$ are the Schmidt coefficients (singular values), $U$ and $V$
are unitary operators, and $r=2^{(m+1)/2}$ ($r=2^{m/2}$) if $m$ is odd (even).
Given this decomposition, we can prepare the state $\ket{\tilde{c}}$ using the
circuit shown in Fig. (1c) [33]. The gate $A$ in Fig. (1c) is taken to be the
uniformly controlled rotations that encode the Schmidt coefficients in the
state $\ket{\alpha}\,=\,\sum_{k=0}^{r}\alpha_{k}\ket{k}\,,$ and the gates $U$
and $V$ are generated using methods for unitary synthesis methods such as
those described in [34]. The classical computational cost to implement this
method grows exponentially in $m$, which is due to the cost of finding the
Schmidt decomposition of the target state $\ket{\tilde{c}}$ and the cost of
computing the gate decompositions of unitaries $A$, $U$, and $V$.

These two methods complement each other nicely. On the one hand, the first
method has a lower classical pre-processing cost than the second method. On
the other hand, the first method leads to a quantum circuit with a higher gate
count than the second method for generic states. Since a lower gate count is
more desirable when working with noisy quantum computers, we used the second
method for most of the experiments we performed on real-world quantum
computers as we discuss in Sec. (III). However, we used the first method when
performing quantum simulations which we discuss in the next subsection.

Finally, even though the asymptotic scaling of the computational cost of
finding $U_{c}$ using the aforementioned two methods is exponential in $m$, we
find that the computational time to find implement $U_{c}$ is reasonable in
practice. For example, for $m=10$ (i.e., around $2000$ Fourier coefficients),
it takes less than $6$ and $3$ seconds to find $U_{c}$ using the first and the
second method respectively. More details about the benchmarking of the CPU
wall clock time can be found in Supplementary Sec. (C).

### Function Loading

Generally, only a few terms in the truncated Fourier series are needed to
well-approximate any periodic function $f$. Now, given a Fourier series
approximation $f_{(m)}=\sum_{k=-M}^{M}c_{k}e^{-i2\pi kx}$ of a periodic
function $f$, we can use the FSL method to prepare the state
$\ket{f_{(m)}}\approx\ket{f}$. The only source of error between the prepared
state $\ket{f_{(m)}}$ and the target state $\ket{f}$ is the truncation error
from the Fourier series approximation. The error as measured by infidelity,
$\epsilon_{(m)}:=1-\left|\braket{f}{f_{(m)}}\right|^{2}\,$, decays
exponentially with $m$ in the limit of large $n$ as we show in the
supplementary Sec. (D). Moreover, the rate of exponential decay of infidelity
depends on the smoothness of the function $f$. The bounds on infidelity
derived in the supplementary Sec. (D) can be used to determine the number of
Fourier modes needed to prepare the state $\ket{f}$ with a specified
infidelity.

In order to approximate a function by a Fourier series, we first need to find
its Fourier coefficients, which leads to additional classical pre-processing
costs. The fast Fourier transform (FFT) has a computational cost of
$O(n2^{n})$, which takes away any hope of exponential advantage. Fortunately,
there exists various sparse Fourier transform algorithms that can be used to
find $M\ll 2^{n}$ dominant Fourier coefficients efficiently [20]. In
particular, algorithms developed in [21, 22, 23] compute the $M$ non-zero
coefficients in time $O(M\log M)$. However, for our implementations of the FSL
method, we used the conventional FFT.

To demonstrate the usefulness of the FSL method, we performed the quantum
simulation of loading various periodic functions into a $20$ qubit state using
the ‘statevector_simulator’ provided by Qiskit [31]. For concreteness, we
chose $m=6$ (i.e., $128$ Fourier modes) and loaded the Fourier coefficients
using the uniformly controlled rotations shown in Fig. (1b). We were able to
load non-trivial functions such as $x^{x}$, a sinc function, the reflected put
option function of Ref. [35], and the wavefunction of the first excited state
of a quantum harmonic oscillator with an infidelity of less than $10^{-6}$,
$10^{-7}$, $10^{-6}$, and $10^{-12}$, respectively. The results of the
simulation are shown in Fig. (2a)-(2d).

Figure 2: Quantum circuit simulations of loading each of the following 1D
functions into a quantum state of $n=20$ qubits using the FSL method: (a)
$x^{x}$, (b) sinc, (c) reflected put option, (d) first excited state of a
quantum oscillator, (e) hyperbolic tangent, and (f) a piece-wise discontinuous
function. Continuous functions in (a)-(e) were approximated with $128$ Fourier
coefficients (i.e., $m=6$) whereas the piece-wise discontinuous function in
(f) was approximated with $512$ Fourier coefficients (i.e., $m=8$). These
functions were loaded with infidelity of less than (a) $10^{-6}$, (b)
$10^{-7}$, (c) $10^{-6}$, (d) $10^{-12}$, (e) $10^{-8}$, and (f) $10^{-3}$.
All of the simulations are performed using the ‘statevector_simulator’
provided by Qiskit [31] and the Fourier coefficients were loaded using the
uniformly controlled rotations shown in Fig. (1b). The results are presented
up to normalization.

The FSL method can be modified to load non-periodic functions as well. Given a
non-periodic function $f$ on the domain $[0,1]$, we can define a function $F$
on the extended domain $[0,2]$ such that $F(x)=f(x)/\sqrt{2}$ for $0\leq x\leq
1$ and $F(x)=f(2-x)/\sqrt{2}$ for $1\leq x\leq 2$. By construction, $F$ is
periodic on the domain $[0,2]$, and hence, can be loaded with high fidelity
into a quantum state of $(n+1)$ qubits using the FSL method. We introduce an
ancilla qubit so that there will be twice as many grid points on the extended
domain which ensures that the grid spacing remains $1/2^{n}$. The state
$\ket{F}$ of $n+1$ qubits is of the form

$\displaystyle\ket{F}\,=\,\frac{1}{\sqrt{2}}\ket{0}\otimes\sum_{k=0}^{2^{n}-1}f_{k}\ket{k}+\frac{1}{\sqrt{2}}\ket{1}\otimes\sum_{k=0}^{2^{n}-1}f_{2^{n}-1-k}\ket{k}\,.$
(3)

With the observation that $X^{\otimes n}\ket{k}\,=\,\ket{2^{n}-1-k}$, where
$X$ is the Pauli-$X$ operator, we find that the state $\ket{F}$ can be written
as

$\displaystyle\ket{F}\,=\,\frac{1}{\sqrt{2}}\ket{0}\otimes\ket{f}+\frac{1}{\sqrt{2}}\ket{1}\otimes
X^{\otimes n}\ket{f}\,.$ (4)

Once we prepare $\ket{F}$ using the FSL, we measure the ancilla qubit in the
computational basis. We do nothing to the $n$ qubits if the measurement output
is $0$, but we apply $X$ to each of the $n$ qubits if the measurement output
is $1$. In either case, the state of the $n$ qubits is $\ket{f}$ as desired.
Alternatively, we can prepare the state $\ket{F}$ and apply CNOT gates on each
of the $n$ qubits and a Hadamard gate on the ancilla qubit. This disentangles
the ancilla qubit from the other $n$ qubits and maps $\ket{F}$ to
$\ket{0}\otimes\ket{f}\,$. The circuits implementing both methods are provided
in Fig. (1d).

In order to demonstrate the FSL’s ability to load non-periodic functions, we
performed the simulation of loading a hyperbolic tangent function into a $20$
qubit state. Once again, we chose $m=6$ and used the uniformly controlled
rotations to load the Fourier coefficients. We were able to load the
hyperbolic tangent with infidelity of less than $10^{-8}$ and the result of
this simulation is shown in Fig. (2e).

Generalizing further, Fourier series approximations of piece-wise
discontinuous functions typically exhibit large oscillations at the
discontinuities resulting in overshooting; this is known as the Gibbs
phenomenon. Depending on the application, these large oscillations may not be
desirable. One way to suppress these oscillations is to apply Fourier filters,
such as a Lanczos $\sigma$-factor, to the Fourier coefficients:
$c_{k}\to\,\left(\text{sinc}(\pi k/M)\right)^{a}c_{k}\,$, where $M$ is the
order of the partial Fourier series and $a\in\mathbb{R}^{+}\,$ [36]. The
Lanczos filter is just one of many Fourier filters available and is by no
means the best choice for every discontinuous function of interest [37, 38,
39]. The FSL method can easily incorporate these Fourier filters with a small
additional classical pre-processing cost. Instead of directly loading in the
Fourier coefficients using the unitary $U_{c}$ in Fig. (1a), one can load the
filtered Fourier coefficients. This ability of the FSL method to suppress the
Gibbs phenomenon gives it an edge over other QFT-based state preparation
methods [25, 29]. (See Sec. (IV) for more details.) As an example, we
performed the quantum simulation of loading a piece-wise discontinuous shown
in Fig. (2f) into a $20$ qubit state. We found that with $m=8$, we were able
to prepare the desired state with an infidelity of around $10^{-3}$.

Remarkably, the FSL method can be generalized to functions of more than one
variable. The approach is the same: we first find a multi-variable Fourier
series approximation of the target function, prepare a sparse state of Fourier
coefficients using either the uniformly controlled rotation or the Schmidt-
decomposition circuit, and then apply inverse QFT operators. Further details,
including the circuit diagrams for loading functions of two variables, are
presented in the Supplementary Sec. (A). Here, we simply demonstrate the
practicality of the FSL method by presenting the simulation results of loading
a two-dimensional sinc function into a $20$ qubits state, $10$ qubits for each
dimension. We approximated this two-dimensional sinc function with $256$
Fourier coefficients, which we loaded to $8$ qubits using the cascade of
uniformly controlled rotations method for implementing $U_{c}$. With only
$256$ Fourier modes, we were able to load the sinc function with infidelity of
$5.0\times 10^{-4}$. The results of the simulation performed using Qiskit’s
‘statevector_simulator’ are presented in Fig. (3).

(a) Exact two-dimensional sinc function(b) Simulation of the two-dimensional
FSL
method$x$$y$$0.0$$0.5$$1.0$$0.0$$0.5$$1.0$$0.000$$0.002$$0.004$$0.006$$f(x,y)$$x$$y$$0.0$$0.5$$1.0$$0.0$$0.5$$1.0$$0.000$$0.002$$0.004$$0.006$$f(x,y)$$-0.001$$0.000$$0.001$$0.002$$0.003$$0.004$$0.005$
Figure 3: Quantum circuit simulation demonstrating the FSL method being used
to load a 2D function. Specifically, a two-dimensional sinc function was
loaded into a quantum state of $20$ qubits (i.e., $n=10$) with infidelity of
$5.0\times 10^{-4}$. In this simulation, $256$ Fourier coefficients (i.e.,
$m=3$) were loaded using the uniformly controlled rotations shown in Fig.
(1b). This simulation was performed using the ‘statevector_simulator’ provided
by Qiskit [31].

We now conclude this introduction to the mechanisms of the FSL method. In the
next section, we present the experimental results from running the FSL method
on the Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers. Our results
indicate that the FSL method is able to load functions with high fidelity even
in the presence of noise in near-term quantum computers.

##  III Experimental Results

In the previous section, we demonstrated that the FSL method is capable of
preparing function-encoding states with high fidelity in an ideal noiseless
setting. However, present-day quantum computers are very noisy and exhibit
complex errors when running highly structured quantum programs [40]. To test
how the FSL method performs in the presence of noise, we performed experiments
on the Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers which we accessed
remotely via an Open QASM-based API [32]. System Model H$1$ quantum computers
have, on average, single-qubit gate infidelities of $4\times 10^{-5}$, two-
qubit gate infidelities of $3\times 10^{-3}$, state preparation and
measurement (SPAM) error of $3\times 10^{-3}$, and a measurement cross-talk
error of $2\times 10^{-5}$. More impressively, System Model H$1$ quantum
computers support full connectivity between qubits and allows for mid-circuit
measurements, the latter of which we used extensively as we discuss below.
More details about the specifications for the System Model H$1$ quantum
computers, H1-1 and H1-2, are provided in the supplementary Sec. (E).

$-0.25$$0$$0.25$(a)(b)$\text{Re}\big{\\{}\langle k|f\rangle\big{\\}}$ vs
$\text{Re}\big{\\{}\langle k|f_{\text{tomo}}\rangle\big{\\}}$
$\text{Im}\big{\\{}\langle k|f\rangle\big{\\}}$ vs $\text{Im}\big{\\{}\langle
k|f_{\text{tomo}}\rangle\big{\\}}$ $\langle{k}|{f}\rangle$$\langle
k|{f_{\text{tomo}}}\rangle$$0$$1$$2$$3$$k$$4$$5$$6$$7$$0$$1$$2$$3$$k$$4$$5$$6$$7$$-0.3$$-0.2$$-0.1$$0.0$$0.1$$0.2$$0.3$$-0.4$$-0.2$$0.0$$0.2$$0.4$

Figure 4: Experimental results of loading a complex-valued function,
$f(x)\,=\,\left(\cos(2\pi x)-1.5i\cos(6\pi x)\right)/\sqrt{13}\,$, into a
three qubit state on the Quantinuum H$1$-$1$ quantum computer using the
quantum circuit shown in Fig. (12). (a) The loaded quantum state was
reconstructed through quantum state tomography based on maximum likelihood
estimation, and has a fidelity of $95.5\%$ with the target state and a purity
of $0.918$. (b) Comparison of the elements of the target function state with
the eigenvector of $\rho_{\text{tomo}}$ corresponding to the dominant
eigenvalue. As evident from these plots, the reconstructed state,
$\rho_{\text{tomo}}$, is in good agreement with the target state,
$\rho_{\text{target}}\,=\,|f\rangle\langle f|$. Hence, the FSL method can
accurately load both the phases and the amplitudes of a complex-valued target
state.

For our first experiment, we considered a complex-valued function
$f(x)\,=\,\left(\cos(2\pi x)-1.5i\cos(6\pi x)\right)/\sqrt{13}$. We loaded
this function into a quantum state of three qubits using the quantum circuit
shown in Fig. (12) in the supplementary Sec. (E). To check how close the
prepared state is to the target state, we performed the state tomography where
we measured each qubit in the eigenbasis of Pauli-$Z$, Pauli-$X$, and
Pauli-$Y$ operators by applying $\\{I,H,R_{x}(\pi/2)\\}$ on each qubit before
measuring it. From the results of $800$ measurements for each of the $27$
experiments, we estimated the probabilities of projection,
$P_{\mu_{1},\mu_{2},\mu_{3}}$, onto states
$\ket{\psi_{\mu_{1},\mu_{2},\mu_{3}}}\,\equiv\,\ket{\psi_{\mu_{1}}}\otimes\ket{\psi_{\mu_{2}}}\otimes\ket{\psi_{\mu_{3}}}\,$,
where $\mu_{i}\in\\{0,1,2,3\\}$ and $\ket{\psi_{0}}=\ket{0}\,$,
$\ket{\psi_{1}}=\ket{1}\,$,
$\ket{\psi_{2}}=\left(\ket{0}+\ket{1}\right)/\sqrt{2}\,$, and
$\ket{\psi_{3}}=\left(\ket{0}+i\ket{1}\right)/\sqrt{2}\,$. We assumed the
density matrix to be of the form
$\rho(T)\,=\,\left(T^{\dagger}T\right)/\text{tr}\left(T^{\dagger}T\right)\,$,
where $T$ is a $8\times 8$ lower triangular matrix with real-valued diagonal
elements. The matrix $\rho(T)$ which best fits the measured probabilities of
projection onto states $\ket{\psi_{\mu_{1},\mu_{2},\mu_{3}}}$ was determined
by minimizing the cost function, [41, 42]

$\displaystyle
C(T)\,=\,\sum_{\mu_{1}=0}^{3}\sum_{\mu_{2}=0}^{3}\sum_{\mu_{3}=0}^{3}\,\frac{\Big{(}\bra{\psi_{\mu_{1},\mu_{2},\mu_{3}}}\rho(T)\ket{\psi_{\mu_{1},\mu_{2},\mu_{3}}}\,-\,P_{\mu_{1},\mu_{2},\mu_{3}}\Big{)}^{2}}{\bra{\psi_{\mu_{1},\mu_{2},\mu_{3}}}\rho(T)\ket{\psi_{\mu_{1},\mu_{2},\mu_{3}}}},$
(5)

with respect to $T$. The reconstructed state, $\rho_{\text{tomo}}$ through
this tomography method has a fidelity of $95.5\%$ with the target state,
$\rho_{\text{target}}=|f\rangle\langle f|$, and a purity of
$\text{tr}(\rho_{\text{tomo}}^{2})=0.918$. The matrix elements of the target
density matrix, $\rho_{\text{target}}$, and the measured density matrix,
$\rho_{\text{tomo}}$, are shown in Fig. (4a). As evident from this figure, the
matrix elements of $\rho_{\text{target}}$ are in good agreement with those of
$\rho_{\text{tomo}}$. Moreover, the comparison of the target function with the
‘measured’ function is shown in Fig. (4b). Since the measured density matrix
$\rho_{\text{tomo}}$ is not a pure state, the notion of the measured function
is not uniquely defined. For concreteness, we defined the measured function
such that its values at $8$ discretized points are the elements of the state
$\ket{f_{\text{tomo}}}$ in the computational basis, where
$\ket{f_{\text{tomo}}}$ is the eigenvector of $\rho_{\text{tomo}}$
corresponding to the dominant eigenvalue 111Of course, $\ket{f_{\text{tomo}}}$
is only defined up to a global phase. We fix this phase by demanding that
$\langle f|f_{\text{tomo}}\rangle$ is real and positive..

The above approach, based on the maximum likelihood estimation, is especially
useful as it restricts the resulting density matrix, $\rho(T)$, to being
Hermitian, positive, and normalized. In the supplementary Sec. (E), we discuss
an alternate ‘direct reconstruction’ approach where we measured the density
matrix without imposing any constraints. When using the ‘direct
reconstruction’ of the density matrix, we calculated a fidelity of around
$94.0\%$ between the measured and the target state. However, as is usually the
case with such unconstrained approaches [41], the reconstructed matrix has
negative eigenvalues, and hence, is not a valid density matrix.

We now discuss experiments we performed to test the FSL method’s potential for
loading functions into a quantum state of five or more qubits. Note that due
to the high cost of running a quantum computer, performing similar quantum
state tomography experiments to verify the accuracy of the FSL method is not a
viable option for more than three qubits. Even though there are economical
algorithms to perform state tomography using randomized measurements [44, 45]
or machine learning techniques [46, 47, 48, 49], they are beyond the scope of
this work.

To circumvent the issue of the verification of the prepared quantum state, we
simplified our task in the rest of the experiments and only verified the
amplitudes of the prepared quantum state by performing measurements in the
computational basis. In particular, we first loaded $\sqrt{f(x)}$ for a given
function $f(x)$ into a quantum state
$\ket{\sqrt{f}}=\sum_{k=0}^{2^{n}-1}\sqrt{f_{k}}\ket{k}$ using the FSL method
and then estimated the measurement probabilities, $|\langle
k|\sqrt{f}\rangle|^{2}$, by measuring all the qubits in the computational
basis sufficiently many times. This allowed us to compare the measurement
probabilities with the value of the target function $f(x)$ at $x=k/2^{n}$. It
is worth mentioning that even though we are unable to verify the accuracy of
the local phases in the quantum state prepared using the FSL method, we still
expect the experiments described above to test the performance of the FSL
method on a present-day noisy quantum computer. Note that we already have
enough evidence that the FSL method is, in principle, capable of loading
correct amplitudes and local phases based on the simulation results in Fig.
(2) and Fig. (3) and the proof of the FSL method provided in the supplementary
Sec. (A). Therefore, the only thing that needs to be checked is if the quantum
state prepared using the FSL circuit can survive the hardware noise present in
noisy quantum computers. For example, if the depth of the FSL quantum circuit,
despite being $O(n)$, is not small enough, the resultant noisy state will be
far from the target state, in which case even the measured amplitudes will be
far from the actual amplitudes. Thus, measuring amplitudes of the prepared
state and comparing them with those of the target state provides a nice
alternate to more expensive state tomography to check the practicality of the
FSL method for present-day noisy quantum computers.

Before we proceed, let us also point out that since we are only comparing the
amplitudes of the prepared state to the target function, it seems more
reasonable to quantify the accuracy of the FSL method in terms of the
classical fidelity between the measurement probabilities (in the computational
basis) and the target function. Recall that the classical fidelity between
probability distributions $p$ and $q$ is defined as
$F(p,q)\,=\,\left(\sum_{i}\sqrt{p_{i}}\sqrt{q_{i}}\,\right)^{2}\,$.

For our first experiment involving five or more qubits, we considered a bi-
modal Gaussian function given by a superposition of two Gaussian functions;
see Eq. (51) in the supplementary Sec. (E). We loaded this function into a
$5$-qubit state and ran $5000$ shots of measurements. The measurement
probabilities are in good agreement with the target function as shown in Fig.
(5a). In fact, we obtained a classical fidelity of $99.8\%$ between the
measured probabilities and the target function. In this experiment, we used
the Schmidt-decomposition circuit shown in Fig. (1c) to load the Fourier
coefficients. Further details about this circuit can be found in the
supplementary Sec. (E).

(a) Bi-modal Gaussian function on $5$ qubitsClassically-controlled QFTTextbook
QFT$x$$0$$5$$10$$15$$20$$25$$30$$f(x)$$0.00$$0.02$$0.04$$0.06$$0.08$$x$$0$$5$$10$$15$$20$$25$$30$$f(x)$$0.00$$0.02$$0.04$$0.06$$0.08$$x$$0$$10$$20$$30$$40$$50$$60$$f(x)$$0.00$$0.01$$0.02$$0.03$$0.04$$0.05$$0.06$$0.07$$x$$0$$5$$10$$15$$20$$25$$30$$f(x)$$0.00$$0.02$$0.04$$0.06$$0.08$$0.10$$x$$0$$10$$20$$30$$40$$50$$60$$f(x)$$0.00$$0.01$$0.02$$0.03$$0.04$$0.05$(c)
Lorentzian function on $6$ qubits(b) Log-normal distribution on $6$ qubits(d)
‘Spiky’ function on $5$ qubits Figure 5: Experimental results of loading
functions on the Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers using the
FSL method. The histograms (blue bars) represent the measurement probabilities
(in the computational basis) determined from $5,000$ measurements for (a) and
(d) and $10,000$ measurements for (b) and (c), whereas the orange dots
represent the values of target functions. The classical fidelities between the
empirical measurement probabilities and the target function are (a) more than
$99.7\%$, (b) $98.8\%$, (c) $99.7\%$, and (d) $98.9\%$. The Schmidt-
decomposition circuit shown in Fig. (1c) was used to load the Fourier
coefficients in experiments (a), (b), and (c), whereas the quantum circuit in
Fig. (13) was used to load the ‘spiky’ function in experiment (d).

When all the qubits are measured right after the application of the inverse
QFT, as is the case in the above experiment, it is possible to replace all the
controlled phase gates in the inverse QFT with classically-controlled single-
qubit rotation gates [50]. Utilizing this classically-controlled
implementation of the inverse QFT, we repeated the above bi-modal Gaussian
function experiment and obtained a fidelity of $99.7\%$ between the measured
probabilities and the target function. The measurement probabilities for this
experiment are also shown in Fig. (5a).

As evident from the measured fidelities and the experimental results in Fig.
(5a), the FSL method works equally well with the ‘textbook’ implementation and
the classically-controlled implementation of the inverse QFT. In other words,
the effect of hardware noise is similar for both implementations of the
inverse QFT. However, the implementation of the classically-controlled QFT was
more economical as the experiment using the classically-controlled QFT
required almost three times fewer compute credits than the experiment with the
textbook QFT. For this reason, we used the classically-controlled QFT to
perform the rest of the function-loading experiments that we discuss.

Next, we loaded a log-normal distribution and a Lorentzian function given in
Eq. (52) and Eq. (53), respectively, to six qubits. From the results of
$10,000$ measurements, we obtained a classical fidelity of $98.8\%$ for the
log-normal distribution and a classical fidelity of $99.7\%$ for the
Lorentzian function. The measured probabilities for these experiments are
shown in Fig. (5b) and Fig. (5c) respectively. For these experiments, we used
the Schmidt-decomposition circuit shown in Fig. (1c) to load the Fourier
coefficients. Further details about these circuits are provided in the
supplementary Sec. (E).

After testing the FSL method’s ability to load smooth functions in the
presence of noise, we next considered a ‘spiky’ function given by a
superposition of low-frequency ($\omega=4\pi$) and high-frequency
($\omega=20\pi$) cosine waves; see Eq. (54). We loaded this function into a
$5$-qubit state using the quantum circuit shown in Fig. (13) and performed
$5000$ measurements. From these measurements, we determined the measurement
probabilities which are shown in Fig. (5d) and obtained a classical fidelity
of $98.9\%$ between the measurement probabilities and the target function.

Encouraged by the success of the FSL method to load functions of a single
variable on five and six qubits states, we wanted to verify how the FSL method
performs when loading a function of two variables on $10$ qubits of a noisy
quantum computer. Specifically, we considered a function
$f(x,y)\,=\,\exp\left(-(x-\mu_{1})^{2}/\sigma_{11}^{2}\,-\,(y-\mu_{1})^{2}/\sigma_{12}^{2}\right)\,+\,\lambda\,\exp\left(-(x-\mu_{2})^{2}/\sigma_{21}^{2}\,-\,(y-\mu_{2})^{2}/\sigma_{22}^{2}\right)\,$
with $\mu_{1}=0.65$, $\mu_{2}=0.35$, $\sigma_{11}=\sigma_{22}=\sqrt{1/50}$,
$\sigma_{12}=\sqrt{1/40}$, $\sigma_{21}=\sqrt{1/30}$, and $\lambda=0.5\,$. We
then loaded $\sqrt{f(x,y)}$ into a $10$-qubit state (five qubits for each
variable). A graph of the function $f(x,y)$ is shown in Fig. (6a). For
concreteness, we used the Schmidt-decomposition circuit shown in Fig. (1c) to
load the Fourier coefficients of this function. Performing this experiment is
expensive due to the large number of measurements required to convincingly
compare the amplitudes of the prepared state with the target function,
$f(x,y)$.

Figure 6: (a) The function
$f(x,y)=\exp\left(-(x-\mu_{1})^{2}/\sigma_{11}^{2}-(y-\mu_{1})^{2}/\sigma_{12}^{2}\right)+\lambda\exp\left(-(x-\mu_{2})^{2}/\sigma_{21}^{2}-(y-\mu_{2})^{2}/\sigma_{22}^{2}\right)$
with $\mu_{1}=0.65$, $\mu_{2}=0.35$, $\sigma_{11}=\sigma_{22}=\sqrt{1/50}$,
$\sigma_{12}=\sqrt{1/40}$, $\sigma_{21}=\sqrt{1/30}$, and $\lambda=0.5\,$ we
use to test the performance of the FSL method to load functions of two
variables on a present-day noisy quantum computer. (b) Measurement
probabilities are determined by first loading $\sqrt{f(x,y)}$ into a quantum
state of $10$ qubits of the Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers
and then measuring all the qubits in the computational basis. The Schmidt-
decomposition circuit shown in Fig. (1c) was used to load the Fourier
coefficients of this function. Even with $14,031$ measurements, we observe
that the general shape of the target function has emerged and we obtain a
fidelity of $94.0\%$ between the measured probabilities and the target
function. (c) To see how the experimental results in (b) would improve with
the number of measurements, we perform the simulation of $100,000$
measurements on the Quantinuum H$1$-$1$ and H$1$-$2$ emulators whose noise
models are calibrated to simulate the noise present in the Quantinuum H$1$-$1$
and H$1$-$2$ quantum computers. These results demonstrate that the FSL method
can load multi-variable functions on up to $10$ qubits on a near-term noisy
quantum computer.

Due to limited compute resources, we only managed to perform $14,031$
measurements on the Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers. Even
with this small number of measurements, we observe the correct overall shape
of the target function as shown in Fig. (6b). Furthermore, we obtained a
fidelity of $94.0\%$ between the measured probabilities and the target
function. We expect these results to improve further with the number of
measurements. To justify this expectation, we performed a simulation of the
experiment on the Quantinuum H$1$-$1$ and H$1$-$2$ emulators whose noise
models and parameters match those of the Quantinuum H$1$-$1$ and H$1$-$2$
quantum computers [51]. The resulting probabilities determined from $100,000$
simulated measurements are shown in Fig. (6c) and have a fidelity of $94.1\%$
with the target function.

The results presented in this section clearly demonstrate that the FSL method
can load functions into a quantum state of up to $10$ qubits on present-day
noisy quantum computers. We now conclude with a comparison of the FSL method
with previously known state preparation methods and some possible future
directions.

##  IV Discussion

###  IV.1 Comparison with Previous Work

Many function-loading algorithms have been proposed in the literature. In this
section, we discuss a variety of algorithms that are comparable with the FSL
method. A comparative summary of the algorithms we consider is provided in
Table (1).

We start by discussing other function-loading methods that also use the
quantum Fourier transform (QFT) and have a close resemblance with the FSL
method. One such method is based on Fourier interpolation or trigonometric
interpolation, which uses a Fourier transform to approximate the values of a
function everywhere given the values of the function on a subset of points in
its domain. A quantum circuit implementing this interpolation technique was
proposed in Refs. [25, 26], and was generalized for loading images on a
quantum computer in Ref. [52]. This circuit is similar to the circuit that we
have proposed in Fig. (1a) but with $U_{c}$ replaced with
QFT$\mbox{}_{m^{\prime}}U_{f}$, where $U_{f}$ is any unitary that loads the
values of the function $f$ sampled at $2^{m^{\prime}}$ equidistant grid points
to an $m^{\prime}$-qubit state, and QFT$\mbox{}_{m^{\prime}}$ is a QFT
operator acting on $m^{\prime}$ qubits. Hence, the circuit implementing
Fourier interpolation also has $O(n^{2})$ gate count and $O(n)$ depth.
However, despite their similarity, the FSL method has several advantages over
the interpolation method. Firstly, depending on how the unitary $U_{f}$ is
implemented, the circuit implementing the Fourier interpolation may have a
higher gate count than the FSL circuit in Fig. (1a) due to an additional QFT
operator. Secondly, it is known that the error associated with Fourier
interpolation is always greater than the truncation error due to aliasing
[53]. Therefore, the FSL method is capable of higher accuracy than an
interpolation method. Thirdly, as we demonstrated in Sec. (II), it is easy to
incorporate Fourier filters when using the FSL method, which can suppress the
Gibbs phenomenon present in the Fourier-series approximation of piece-wise
discontinuous functions. It is not clear how to suppress the Gibbs phenomenon
using the interpolation method. Finally, we do not expect the interpolation
method to perform well with functions having sharp peaks as these peaks may
not be captured by the values of the functions at some subset of points.
Consider the $5$-qubit spiky function state in Fig. (4d) as an example. Since
this spiky function has a Fourier mode with frequency $\omega=10\pi$, the
Nyquist theorem dictates that we need the value of the function at
$2^{m^{\prime}}\geq 20$ points to get an accurate result using the
interpolation method. This implies that the interpolation method with
$m^{\prime}<n=5$ cannot accurately load this spiky function into a $5$-qubit
state. The FSL method, on the other hand, can accurately load this function
into a $5$-qubit state as we have demonstrated.

|  Type |  Requires |  Classical |  Classical |  Gate |  Circuit |  Ancilla |
---|---|---|---|---|---|---|---|---
Methods |  of |  Classical |  Time |  Space |  Count |  Depth |  Qubits |  NISQ-y?
|  State |  Optimization |  Complexity |  Complexity |  |  |  |
Fourier Interpolation [25] |  Smooth functions |  No |  $O(\log\omega_{\text{max}})$ |  $O(\log\omega_{\text{max}})$ |  $O(n^{2})$ |  $O(n)$ |  $0$ to $1$ |  Yes
QFT Sampler [29] |  Periodic distributions |  Yes |  $\Omega(p)$ |  $O(p)$ |  $O(n^{2})$ |  $O(n)$ |  $0$ |  Yes
Quantum GAN  [27] |  Arbitrary distributions |  Yes |  $\Omega(p)$ |  $O(p)$ |  $O(nL)$ |  $O(nL)$ |  $0$ |  Yes
Matrix Product State [28] |  Smooth functions |  Yes |  $O(n)$ |  $O(n)$ |  $O(n)$ |  $O(n)$ |  0 |  Yes
Variational ZGR [54] |  Smooth functions |  Yes |  $O(2^{n})$ |  $O(2^{n})$ |  $\Omega(n^{2})$ |  $\Omega(n)$ |  $0$ |  Yes
Hamiltonian Simulation  [55] |  Arbitrary functions |  No |  $O(1)$ |  $O(1)$ |  $O(T_{\text{oracle}})$ |  $O(T_{\text{oracle}})$ |  $O(n+n^{\prime})$ |  No
QSVT  [56] |  Arbitrary functions |  No |  $O(d_{\text{poly}}^{3})$ |  $O(d_{\text{poly}})$ |  $O(nd_{\text{poly}}R)$ |  $O(nd_{\text{poly}}R)$ |  $3$ to $4$ |  No
Kitaev-Webb [57] |  Gaussian distribution |  No |  $O(1)$ |  $O(1)$ |  $O(\text{poly}(n))$ |  $O(\text{poly}(n))$ |  $\Omega(n^{\prime})$ |  No
Discrete Random Walks [58] |  Gaussian distribution |  No |  $O(1)$ |  $O(1)$ |  $O(n\log n)$ |  $O(n\log n)$ |  $O(n)$ |  Yes
Powers of Cosine Approximation [59] |  Gaussian distribution |  No |  $O(1)$ |  $O(1)$ |  $O(n^{2})$ |  $O(n)$ |  $0$ |  Yes
FSL (This work) |  Arbitrary functions |  No |  $O(8^{m})$ |  $O(2^{m})$ |  $O(n^{2})$ |  $O(n)$ |  $0$ to $1$ |  Yes
Table 1: A comparison between various function-loading algorithms. The
classical cost of the Fourier interpolation method is determined by the
maximum frequency, denoted by $\omega_{\max}$, in the spectrum of the target
function as per the Shannon-Nyquist theorem. In variational algorithms, such
as Variational QFT and Quantum GAN, the number of parameters and the number of
layers in a variational circuit are denoted by $p$ and $L$ respectively. The
degree of the polynomial used to approximate the target function and the
number of rounds of amplitude amplification needed in the QSVT method are
denoted by $d_{\text{poly}}$ and $R$ respectively. In fully quantum
algorithms, such as Hamiltonian simulation and the Kitaev-Webb method, the
function/Hamiltonian/angles of rotation are stored using quantum registers.
The number of qubits used to store this information is denoted by
$n^{\prime}$. Furthermore, $T_{\text{oracle}}$ denotes the circuit depth of
the oracle used in the implementation of the corresponding method. Finally,
the ‘NISQ-y’ category classifies those methods whose circuits are relatively
shallow [e.g., circuit depth of less than $1000$ for states of less than $100$
qubits] thereby allowing for the preparation of states with high fidelity
[e.g., $>90\%$] on noisy, intermediate-scale quantum (NISQ) devices.

The Fourier transform has also been used in Ref. [29] to load probability
distributions as the amplitudes of a quantum state. In fact, the circuit
proposed in Ref. [29] to load a distribution consists of a unitary operator
acting on a small subset of qubits followed by a QFT operator acting on all of
the qubits, and hence, is similar to the circuit that we have proposed in Fig.
(1a). However, there are several key differences. First, the unitary operator
that precedes the QFT operator is taken to be a function of variational
parameters which needs to be determined through classical optimization.
Secondly, and more importantly, the subset of qubits on which the variational
unitary acted on were not entangled with the rest of the qubits before the
application of the QFT operator. In other words, there are no CNOT gates
similar to those in Fig. (1a). This has interesting consequences as it can be
shown that this circuit can only load functions of the form
$f(x)=\sum_{k=0}^{2^{m}-1}c_{k}e^{-i2\pi kx}$, where $c_{k}$ depend on the
variational parameters. This is only half of the Fourier series as it is
missing the complex conjugate terms $c_{-k}=c^{*}_{k}$. Therefore, even though
this circuit can be trained to sample a target probability distribution, it
can not be used to load distributions as a subroutine of another quantum
algorithm.

We now discuss other state preparation methods that are not based on Fourier
transforms and compare them with the FSL method. Consider the variational
algorithm based on the Generative Adversarial Network (GAN) which involves a
variational quantum circuit (generator) and a classical neural network
(discriminator) [27]. The generator and the discriminator are simultaneously
optimized until the discriminator cannot distinguish the samples from the
generator with the given samples from the target distribution. The variational
circuit that was used in Ref. [27] consists of $L$ consecutive layers, each of
depth $O(n)$ and containing $O(n)$ quantum gates. However, the computational
cost of optimizing a variational quantum circuit scales at least linearly in
the number of parameters $p$, which in this case is $O(nL)$. Despite the
usefulness of a quantum GAN for state preparations, there are limitations. For
example, by construction, a quantum GAN can only learn to load distributions
and not arbitrary complex-valued functions into a quantum state. Moreover, the
performance of a quantum GAN-based method highly depends on the choice of the
initial state on which the variational circuit acts [27]. In general, it can
take a significant amount of time to train a quantum GAN, further hindering
the overall time complexity of this state preparation method. Therefore, we
expect that the FSL method, which does not suffer from these limitations,
offers a good alternative for loading functions and distributions into a
quantum state.

Another interesting state preparation method is based on matrix product states
(MPS) [28, 12]. The main insight that was used in Ref. [28] was that smooth,
differentiable, real-valued (SDR) functions can be well-approximated by piece-
wise polynomials. This led to an approximation of SDR function states using
MPS of bond dimension $K(p+1)$, where $K$ is the number of subdomains and $p$
is the order of the polynomials [60]. This MPS was further approximated by a
MPS of bond dimension $2$ using the MPS compression algorithm which has a time
complexity of $O(n)$ [61]. The resulting compressed MPS was converted into a
quantum circuit of $O(n)$ depth and gate count [62]. The linear depth circuit
and gate count make this approach exceedingly efficient for loading SDR
functions into a quantum state. However, the approximation of the SDR
functions by piece-wise polynomials and the compression of the MPS lead to a
build-up of error. Indeed, the fidelities reported in [28] are lower than the
fidelities that we obtained in Sec. (II) from the FSL method for which the
only source of error is the truncation error in the Fourier series.
Furthermore, the FSL method can load a larger class of functions than the MPS
method. Specifically, the FSL method can load complex-valued and non-smooth
functions.

In in Sec. (II), we proposed a method for encoding the target Fourier
coefficients via a unitary $U_{c}$ where $U_{c}$ is taken to be the circuit of
Zalka, Grover, and Rudolph (ZGR) which consists of a cascade of uniformly
controlled rotations. Recently, another state preparation algorithm also based
on the ZGR circuit was proposed in Ref. [54]. In fact, it was observed in Ref.
[54] that for a class of positive integrable functions for which
$|\partial_{x}^{2}\log f^{2}(x)|\leq 8\pi$ on the domain $[0,1]$, most of the
angles appearing in the ZGR circuits are close to each other. Hence, it is
possible to replace uniformly $k$-controlled rotation gates with single
rotation gates for all $k>k_{0}$ for some $k_{0}>0$. This approximation
results in a reduction of the gate count and the depth of the circuit from
$O(2^{n})$ to $O(2^{k_{0}})$. However, for functions that are non-positive or
piece-wise defined, not all of the $k$-controlled rotation gates for $k>k_{0}$
can be replaced with single qubit gates. Furthermore, the angles of these
gates have to be determined through classical optimization. Hence, the total
gate count for this variational algorithm is $\Omega(n^{2})$ whereas the
classical time and space complexity is exponential in $n$ [54].

The methods that we have discussed until now have all required some sort of
classical pre-processing. A state preparation method that leverages adiabatic
time evolution was recently proposed in Ref. [55] which does not require
classical pre-processing at the cost of significantly deeper circuits. The
idea behind this method is to encode the target state as a rank-$1$ projector
Hamiltonian $\hat{H}=|f\rangle\langle f|$, then perform the adiabatic
evolution using low-rank, $1$-sparse Hamiltonian simulation techniques [63,
64, 65, 66]. The time evolution is simulated using Trotterization-like time
steps that suppress adiabatic and simulation errors. Each time step requires
$O(1)$ call to an oracle which implements the Hamiltonian using the arithmetic
operations. Interestingly, it was proven in Ref. [55] that the query
complexity of this algorithm approaches a constant asymptotically in the
number of qubits. However, note that the simulation results reported in Ref.
[55] show that on the order of $10^{5}$ Trotterization steps are required to
achieve an error of $10^{-5}$ for $n>5$ qubits. Given that the implementation
of the oracle requires quadratic gate count and circuit depth [67], the total
resources needed to run this fully-quantum state preparation algorithm are
orders of magnitude higher than the present-day technology. This is not
surprising as the authors intended the method would be executed on a fault-
tolerant quantum computer. It is worth appreciating that this method is also
capable of preparing arbitrary quantum states (i.e., loading arbitrary
discrete functions), however, the associated complexity can differ greatly
from the function loading use case.

Recently, a function-loading algorithm based on the quantum singular value
transformation (QSVT), was proposed [56]. The basic idea of the QSVT is to
find a block encoding of $h(A)$ given a block-encoding $U_{A}$ of a Hermitian
operator $A=\sum_{k}a_{k}|k\rangle\langle k|$, where $h(x)$ is a polynomial of
degree $d_{poly}$. By cleverly choosing $A=\sum_{k}\sin(k/N)|k\rangle\langle
k|$ and choosing $h(x)$ to be a polynomial approximation of $f(\arcsin(x))$,
the authors were able to prepare a block-encoding of
$\sum_{k}f(k/N)|k\rangle\langle k|$. The desired state could then be prepared
with high probability by applying this block-encoding to $\ket{+}^{\otimes
n}\otimes\ket{0}_{\text{ancilla}}$. In many cases, several rounds of amplitude
amplification are required to boost the probability of success closer to 1.
The number of rounds of amplitude estimation, $R$, needed is inversely related
to the “$L_{2}$-norm filling ratio” of the target function which is
approximately the ratio $\|f\|_{2}/\|f\|_{\infty}$. Given that quantum circuit
that implements the QSVT requires $d_{poly}/2$ applications of $U_{A}$ and
$U_{A}^{\dagger}$ [56], and $U_{A}$ for $A=\sum_{k}\sin(k/N)|k\rangle\langle
k|$ can be implemented with a circuit of depth $n$, the full quantum circuit
needed to prepare the desired state has a depth of $O(nd_{poly}R)$. Moreover,
as discussed in Ref. [56], the full quantum circuit consists of
$O(nd_{poly}R)$ quantum gates and requires at most $4$ ancilla qubits. Just
like the Hamiltonian simulation method discussed above, this method is only
suitable for fault-tolerant quantum computers. However, this method replaces
the requirement for an oracle needed in the Hamiltonian simulation method with
a classical pre-processing step in which all the angles of rotations in the
quantum circuit can be determined using algorithms such as those described in
Refs. [68, 69, 70].

There are many other state preparation algorithms that are tailored to a
specific distribution. The most notable example is the algorithm of Kitaev and
Webb that loads a Gaussian distribution using a recursive algorithm [57]. At
each recursive step, the algorithm calculates the angles of rotations directly
on an ancilla quantum register and rotates the non-ancillary qubits by that
angle. Even though this algorithm asymptotically only requires polynomial
resources, a recent study argued that this algorithm is more costly than a
generic exponentially scaling algorithm for $n\leq 15$ [71]. This is due to
the large complexity of the quantum arithmetic circuits [67]. The method of
Kitaev and Webb can also be generalized to load $D$-dimensional correlated
Gaussian distributions. The algorithm first loads $D$ uncorrelated one-
dimensional Gaussian distributions and then applies a ‘shearing’
transformation which implements the change of basis required to introduce the
correlations between different dimensions. The implementation of the shearing
transformation requires a circuit involving $(2n+1)+2\log_{2}(D+1)$ ancillae
qubits and $O(n^{2}D^{2})$ CNOT gates [71]. In contrast, the FSL method does
not require the shearing transformation as it can load the correlated
distribution directly and has a gate count of $O(n^{2}D)$.

There are a few other proposals to load a Gaussian distribution that are,
unlike the Kitaev-Webb algorithm, suitable for near-term quantum devices For
example, Rattew et al. have proposed a Gaussian distribution loading method
that is provably robust to bit-flip and phase-flip errors [58]. This method is
based on a discrete random walk inspired by Galton machines whose walk
operator maps a computational basis state $\ket{j}$ to a superposition state
$\ket{j}+\ket{j+1}$ up to a normalization factor. Rattew et al. propose first
loading a low-resolution Gaussian distribution on a small number of qubits
using a method such as the cascade of controlled rotations depicted in Fig.
(1b), then iteratively qubits are added in the
$\ket{+}=\left(\ket{0}+\ket{1}\right)/\sqrt{2}$ state followed by a small
number of applications of a quantum walk operator until the target number of
qubit is reached. The computational complexity of this approach depends upon
the specific implementation of the walk operator which requires a choice of
quantum adder circuit. Hence, if the logarithmic depth adder presented by
Draper et al. [72] is used, then the overall quantum gate complexity of this
Gaussian loading method is $O(n\log n)$ and $O(n)$ ancilla qubits are
required. Alternatively, a QFT-based adder can be used which leads to an
overall quantum gate complexity of $O(n^{2})$ and only a single ancilla is
required.

Another method to approximately load Gaussian distributions into a quantum
state was introduced by Markov et al. in Ref. [59]. This method is based on
approximations of Gaussian distributions by powers of trigonometric functions.
In particular, these approximations can be derived from the limit
$\lim_{r\to\infty}\cos^{r}(x/\sqrt{r})=e^{-x^{2}/2}$. To illustrate their
method, consider the task of loading the 2nd power approximation
$p_{2}(x)=\sigma^{-1}\cos^{2}(\pi(x-\mu)/2\sigma)$ for $\mu-\sigma\leq
x\leq\mu+\sigma$. The main insight of Ref. [59] is that a quantum state
$\ket{\sqrt{p_{2}}}$ can be prepared by applying the inverse QFT operator
applied to the state $(\ket{0}-\ket{2^{n}-1})/\sqrt{2}$. This method is nearly
a special case of the FSL method introduced in this work with the exception
that their method neglects to control for complex phases introduced by the QFT
operator.

###  IV.2 Future Directions

The FSL method inherits the majority of its circuit complexity from the
quantum Fourier transform. Therefore, a natural question to ask is how well
the FSL method performs when using an approximate QFT sporting which would
lead to reduced circuit complexity. Specifically, it has been shown that an
approximate QFT (AQFT) with error bounded by $\varepsilon$ can be implemented
by a circuit with depth bounded above by $O(\log n+\log\log(1/\varepsilon))$
and bounded below by $\Omega(\log n)$ [73]. Therefore, using an AQFT to
implement an approximate FSL would yield a logarithmic-depth function loading
circuit. In order to justify the usage of an AQFT for practical applications,
a rigorous study of the performance of an approximate FSL is needed.

In this work, we have only focused on amplitude embedding where the values of
the functions are loaded as the amplitudes of the quantum state. There exists
another type of embedding known as the key-value embedding, in which the
function $f:\\{0,\ldots,2^{n}-1\\}\to[0,2^{m}-1]$ is encoded as an
$(n+m)$-qubit state,
$\frac{1}{\sqrt{2^{n}}}\sum_{k}\ket{k}_{n}\otimes\ket{f_{k}}_{m}$, where
$\ket{k}_{n}$ is the $n$-qubit ‘key’ register and $\ket{f_{k}}$ is the
$m$-qubit ‘value’ register. The difficulty in preparing such a key-value pair
state is that the circuit preparing the state on the value register has to be
controlled by the key register, which leads to a large gate count. Various
efficient algorithms have been proposed for the key-value embedding of the
functions/distributions [74, 75, 76]. An interesting future direction is to
see if the FSL method introduced in this work can be generalized to
efficiently generate key-value embeddings.

It is also worth emphasizing that the Fourier series is just one possible
orthonormal series representation of functions of interest. There are many
other basis functions – such as wavelets, Chebyshev polynomials, etc. – that
could be more suitable than a Fourier series expansion for some classes of
functions. In this work, we only considered the Fourier-series representation
given that the implementation of the Fourier transform on a quantum computer
requires fewer gates than the implementation of the wavelet transform [77, 78]
and discrete cosine transform [79]. Nevertheless, it could be worthwhile to
investigate the performance of the truncated wavelet series, and other
possible representations, for function loading on quantum computers.

## Data and code availability

A user-friendly implementation of the FSL method, along with various examples,
is available at https://github.com/mcmahon-lab/Fourier-Series-Loader. All the
experimental data gathered from the experiments performed on Quantinuum System
Model H$1$ and the code to analyze the experimental results and to produce the
figures presented in this work are available at
https://doi.org/10.5281/zenodo.8331675.

## Acknowledgements

This research used resources of the Oak Ridge Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC$05$-$00$OR$22725$. We thank Travis Humble (Oak Ridge National Laboratory)
for his support of our research. We also thank Ryan Landfield (ORNL/OLCF User
Support) and Brian Neyenhuis (Quantinuum) for technical assistance and for
providing the hardware specifications included in this work. It is also a
pleasure to thank Javier Gonzalez-Conde, Keisuke Fujii, Constantin Gonciulea,
Vanio Markov, Juan José Garcia-Rípoll, Arthur Rattew, Mikel Sanz, Charlee
Stefanski, Naoki Yamamoto, and Christa Zoufal for useful feedback on a draft
of this manuscript. The work of M.M. was supported by the US Department of
Energy under grant numbers DE-SC$0014123$ and DE-SC$0007884$, and the
QuantiSED Fermilab consortium. P.L.M. gratefully acknowledges support from a
Sloan Foundation Fellowship and a David and Lucile Packard Foundation
Fellowship, and acknowledges membership in the CIFAR Quantum Information
Science Program as an Azrieli Global Scholar.

Software used: The circuit diagrams in this paper and the supplementary
sections were prepared using quantikz package [80]. All the plots were
prepared using Matplotlib [81]. To perform the maximum likelihood estimation
for the tomography experiment, the cost function was minimized using the
generic minimize function in SciPy [82].

## References

  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Physical Review Letters 103, 10.1103/physrevlett.103.150502 (2009).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017).
  * Tong  _et al._ [2021] Y. Tong, D. An, N. Wiebe, and L. Lin, Fast inversion, preconditioned quantum linear system solvers, fast green's-function computation, and fast evaluation of matrix functions, Physical Review A 104, 10.1103/physreva.104.032422 (2021).
  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm, Physical Review Letters 110, 10.1103/physrevlett.110.250504 (2013).
  * Zalka [1998] C. Zalka, Simulating quantum systems on a quantum computer, Proceedings of the Royal Society of London Series A 454, 313 (1998), arXiv:quant-ph/9603026 [quant-ph] .
  * Grover and Rudolph [2002] L. Grover and T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions, arXiv e-prints , quant-ph/0208112 (2002), arXiv:quant-ph/0208112 [quant-ph] .
  * Möttönen  _et al._ [2005] M. Möttönen, J. Vartiainen, V. Bergholm, and M. Salomaa, Transformation of quantum states using uniformly controlled rotations, Quantum Information & Computation 5, 467 (2005).
  * Plesch and Brukner [2011] M. Plesch and Č . Brukner, Quantum-state preparation with universal gate decompositions, Physical Review A 83, 10.1103/physreva.83.032302 (2011).
  * Aaronson [2015] S. Aaronson, Read the fine print, Nature Physics 11, 291 (2015).
  * Benenti and Strini [2008] G. Benenti and G. Strini, Quantum simulation of the single-particle schrödinger equation, American Journal of Physics 76, 657 (2008).
  * Berry  _et al._ [2017] D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, Quantum algorithm for linear differential equations with exponentially improved dependence on precision, Communications in Mathematical Physics 356, 1057 (2017).
  * Lubasch  _et al._ [2020] M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, Variational quantum algorithms for nonlinear problems, Phys. Rev. A 101, 010301 (2020), arXiv:1907.09032 [quant-ph] .
  * Childs  _et al._ [2021] A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021).
  * Liu  _et al._ [2021] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2026805118 (2021).
  * Lemieux  _et al._ [2020] J. Lemieux, B. Heim, D. Poulin, K. Svore, and M. Troyer, Efficient quantum walk circuits for metropolis-hastings algorithm, Quantum 4, 287 (2020).
  * Wiebe  _et al._ [2012] N. Wiebe, D. Braun, and S. Lloyd, Quantum algorithm for data fitting, Physical Review Letters 109, 10.1103/physrevlett.109.050505 (2012).
  * Stamatopoulos  _et al._ [2020] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner, Option pricing using quantum computers, Quantum 4, 291 (2020).
  * Woerner and Egger [2019] S. Woerner and D. J. Egger, Quantum risk analysis, npj Quantum Information 5, 15 (2019).
  * Markov  _et al._ [2022] V. Markov, C. Stefanski, A. Rao, and C. Gonciulea, A generalized quantum inner product and applications to financial engineering (2022).
  * Gilbert  _et al._ [2014] A. C. Gilbert, P. Indyk, M. Iwen, and L. Schmidt, Recent developments in the sparse fourier transform: A compressed fourier transform for big data, IEEE Signal Processing Magazine 31, 91 (2014).
  * Lawlor  _et al._ [2013] D. Lawlor, Y. Wang, and A. Christlieb, Adaptive sub-linear time fourier algorithms, Advances in Adaptive Data Analysis 05, 1350003 (2013), https://doi.org/10.1142/S1793536913500039 .
  * Ghazi  _et al._ [2013] B. Ghazi, H. Hassanieh, P. Indyk, D. Katabi, E. Price, and L. Shi, Sample-optimal average-case sparse fourier transform in two dimensions, 2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton) , 1258 (2013).
  * Pawar and Ramchandran [2013] S. A. Pawar and K. Ramchandran, Computing a k-sparse n-length discrete fourier transform using at most 4k samples and o(k log k) complexity, 2013 IEEE International Symposium on Information Theory , 464 (2013).
  * Gleinig and Hoefler [2021] N. Gleinig and T. Hoefler, An efficient algorithm for sparse quantum state preparation, in  _2021 58th ACM/IEEE Design Automation Conference (DAC)_ (2021) pp. 433–438.
  * García-Ripoll [2021] J. J. García-Ripoll, Quantum-inspired algorithms for multivariate analysis: from interpolation to partial differential equations, Quantum 5, 431 (2021).
  * García-Molina  _et al._ [2022] P. García-Molina, J. Rodríguez-Mediavilla, and J. J. García-Ripoll, Quantum fourier analysis for multivariate functions and applications to a class of schrödinger-type partial differential equations, Physical Review A (2022).
  * Zoufal  _et al._ [2019] C. Zoufal, A. Lucchi, and S. Woerner, Quantum Generative Adversarial Networks for learning and loading random distributions, npj Quantum Information 5, 103 (2019), arXiv:1904.00043 [quant-ph] .
  * Holmes and Matsuura [2020] A. Holmes and A. Y. Matsuura, Efficient Quantum Circuits for Accurate State Preparation of Smooth, Differentiable Functions, arXiv e-prints , arXiv:2005.04351 (2020), arXiv:2005.04351 [quant-ph] .
  * Endo  _et al._ [2020] K. Endo, T. Nakamura, K. Fujii, and N. Yamamoto, Quantum self-learning Monte Carlo and quantum-inspired Fourier transform sampler, Physical Review Research 2, 043442 (2020), arXiv:2005.14075 [quant-ph] .
  * Le  _et al._ [2011a] P. Q. Le, F. Dong, and K. Hirota, A flexible representation of quantum images for polynomial preparation, image compression, and processing operations, Quantum Information Processing 10, 63 (2011a).
  * Anis  _et al._ [2021] M. S. Anis, Abby-Mitchell, H. Abraham, AduOffei, R. Agarwal, G. Agliardi, M. Aharoni, I. Y. Akhalwaya, G. Aleksandrowicz, T. Alexander, M. Amy, S. Anagolum, Anthony-Gandon, E. Arbel, A. Asfaw, A. Athalye, A. Avkhadiev, C. Azaustre, P. BHOLE, A. Banerjee, S. Banerjee, W. Bang, A. Bansal, P. Barkoutsos, A. Barnawal, G. Barron, G. S. Barron, L. Bello, Y. Ben-Haim, M. C. Bennett, D. Bevenius, D. Bhatnagar, A. Bhobe, P. Bianchini, L. S. Bishop, C. Blank, S. Bolos, S. Bopardikar, S. Bosch, S. Brandhofer, Brandon, S. Bravyi, N. Bronn, Bryce-Fuller, D. Bucher, A. Burov, F. Cabrera, P. Calpin, L. Capelluto, J. Carballo, G. Carrascal, A. Carriker, I. Carvalho, A. Chen, C.-F. Chen, E. Chen, J. C. Chen, R. Chen, F. Chevallier, K. Chinda, R. Cholarajan, J. M. Chow, S. Churchill, CisterMoke, C. Claus, C. Clauss, C. Clothier, R. Cocking, R. Cocuzzo, J. Connor, F. Correa, Z. Crockett, A. J. Cross, A. W. Cross, S. Cross, J. Cruz-Benito, C. Culver, A. D. Córcoles-Gonzales, N. D, S. Dague, T. E. Dandachi, A. N. Dangwal, J. Daniel, M. Daniels, M. Dartiailh, A. R. Davila, F. Debouni, A. Dekusar, A. Deshmukh, M. Deshpande, D. Ding, J. Doi, E. M. Dow, E. Drechsler, E. Dumitrescu, K. Dumon, I. Duran, K. EL-Safty, E. Eastman, G. Eberle, A. Ebrahimi, P. Eendebak, D. Egger, ElePT, Emilio, A. Espiricueta, M. Everitt, D. Facoetti, Farida, P. M. Fernández, S. Ferracin, D. Ferrari, A. H. Ferrera, R. Fouilland, A. Frisch, A. Fuhrer, B. Fuller, M. GEORGE, J. Gacon, B. G. Gago, C. Gambella, J. M. Gambetta, A. Gammanpila, L. Garcia, T. Garg, S. Garion, J. R. Garrison, J. Garrison, T. Gates, L. Gil, A. Gilliam, A. Giridharan, J. Gomez-Mosquera, Gonzalo, S. de la Puente González, J. Gorzinski, I. Gould, D. Greenberg, D. Grinko, W. Guan, D. Guijo, J. A. Gunnels, H. Gupta, N. Gupta, J. M. Günther, M. Haglund, I. Haide, I. Hamamura, O. C. Hamido, F. Harkins, K. Hartman, A. Hasan, V. Havlicek, J. Hellmers, Ł. Herok, S. Hillmich, H. Horii, C. Howington, S. Hu, W. Hu, J. Huang, R. Huisman, H. Imai, T. Imamichi, K. Ishizaki, Ishwor, R. Iten, T. Itoko, A. Ivrii, A. Javadi, A. Javadi-Abhari, W. Javed, Q. Jianhua, M. Jivrajani, K. Johns, S. Johnstun, Jonathan-Shoemaker, JosDenmark, JoshDumo, J. Judge, T. Kachmann, A. Kale, N. Kanazawa, J. Kane, Kang-Bae, A. Kapila, A. Karazeev, P. Kassebaum, T. Kehrer, J. Kelso, S. Kelso, V. Khanderao, S. King, Y. Kobayashi, Kovi11Day, A. Kovyrshin, R. Krishnakumar, V. Krishnan, K. Krsulich, P. Kumkar, G. Kus, R. LaRose, E. Lacal, R. Lambert, H. Landa, J. Lapeyre, J. Latone, S. Lawrence, C. Lee, G. Li, J. Lishman, D. Liu, P. Liu, Lolcroc, A. K. M, L. Madden, Y. Maeng, S. Maheshkar, K. Majmudar, A. Malyshev, M. E. Mandouh, J. Manela, Manjula, J. Marecek, M. Marques, K. Marwaha, D. Maslov, P. Maszota, D. Mathews, A. Matsuo, F. Mazhandu, D. McClure, M. McElaney, C. McGarry, D. McKay, D. McPherson, S. Meesala, D. Meirom, C. Mendell, T. Metcalfe, M. Mevissen, A. Meyer, A. Mezzacapo, R. Midha, D. Miller, Z. Minev, A. Mitchell, N. Moll, A. Montanez, G. Monteiro, M. D. Mooring, R. Morales, N. Moran, D. Morcuende, S. Mostafa, M. Motta, R. Moyard, P. Murali, J. Müggenburg, T. NEMOZ, D. Nadlinger, K. Nakanishi, G. Nannicini, P. Nation, E. Navarro, Y. Naveh, S. W. Neagle, P. Neuweiler, A. Ngoueya, T. Nguyen, J. Nicander, Nick-Singstock, P. Niroula, H. Norlen, NuoWenLei, L. J. O’Riordan, O. Ogunbayo, P. Ollitrault, T. Onodera, R. Otaolea, S. Oud, D. Padilha, H. Paik, S. Pal, Y. Pang, A. Panigrahi, V. R. Pascuzzi, S. Perriello, E. Peterson, A. Phan, K. Pilch, F. Piro, M. Pistoia, C. Piveteau, J. Plewa, P. Pocreau, A. Pozas-Kerstjens, R. Pracht, M. Prokop, V. Prutyanov, S. Puri, D. Puzzuoli, J. Pérez, Quant02, Quintiii, R. I. Rahman, A. Raja, R. Rajeev, I. Rajput, N. Ramagiri, A. Rao, R. Raymond, O. Reardon-Smith, R. M.-C. Redondo, M. Reuter, J. Rice, M. Riedemann, Rietesh, D. Risinger, M. L. Rocca, D. M. Rodríguez, RohithKarur, B. Rosand, M. Rossmannek, M. Ryu, T. SAPV, N. R. C. Sa, A. Saha, A. Ash-Saki, S. Sanand, M. Sandberg, H. Sandesara, R. Sapra, H. Sargsyan, A. Sarkar, N. Sathaye, B. Schmitt, C. Schnabel, Z. Schoenfeld, T. L. Scholten, E. Schoute, M. Schulterbrandt, J. Schwarm, J. Seaward, Sergi, I. F. Sertage, K. Setia, F. Shah, N. Shammah, R. Sharma, Y. Shi, J. Shoemaker, A. Silva, A. Simonetto, D. Singh, D. Singh, P. Singh, P. Singkanipa, Y. Siraichi, Siri, J. Sistos, I. Sitdikov, S. Sivarajah, M. B. Sletfjerding, J. A. Smolin, M. Soeken, I. O. Sokolov, I. Sokolov, V. P. Soloviev, SooluThomas, Starfish, D. Steenken, M. Stypulkoski, A. Suau, S. Sun, K. J. Sung, M. Suwama, O. Słowik, H. Takahashi, T. Takawale, I. Tavernelli, C. Taylor, P. Taylour, S. Thomas, K. Tian, M. Tillet, M. Tod, M. Tomasik, C. Tornow, E. de la Torre, J. L. S. Toural, K. Trabing, M. Treinish, D. Trenev, TrishaPe, F. Truger, G. Tsilimigkounakis, D. Tulsi, W. Turner, Y. Vaknin, C. R. Valcarce, F. Varchon, A. Vartak, A. C. Vazquez, P. Vijaywargiya, V. Villar, B. Vishnu, D. Vogt-Lee, C. Vuillot, J. Weaver, J. Weidenfeller, R. Wieczorek, J. A. Wildstrom, J. Wilson, E. Winston, WinterSoldier, J. J. Woehr, S. Woerner, R. Woo, C. J. Wood, R. Wood, S. Wood, J. Wootton, M. Wright, L. Xing, J. YU, B. Yang, U. Yang, J. Yao, D. Yeralin, R. Yonekura, D. Yonge-Mallo, R. Yoshida, R. Young, J. Yu, L. Yu, C. Zachow, L. Zdanski, H. Zhang, I. Zidaru, C. Zoufal, aeddins ibm, alexzhang13, b63, bartek bartlomiej, bcamorrison, brandhsn, charmerDark, deeplokhande, dekel.meirom, dime10, dlasecki, ehchen, fanizzamarco, fs1132429, gadial, galeinston, georgezhou20, georgios ts, gruu, hhorii, hykavitha, itoko, jeppevinkel, jessica angel7, jezerjojo14, jliu45, jscott2, klinvill, krutik2966, ma5x, michelle4654, msuwama, nico lgrs, ntgiwsvp, ordmoj, sagar pahwa, pritamsinha2304, ryancocuzzo, saktar unr, saswati qiskit, septembrr, sethmerkel, sg495, shaashwat, smturro2, sternparky, strickroman, tigerjack, tsura crisaldo, upsideon, vadebayo49, welien, willhbang, wmurphy collabstar, yang.luh, and M. Čepulkovskis, Qiskit: An open-source framework for quantum computing (2021).
  * Quantinuum H1- [1] Quantinuum H1-1, H1-2., https://www.quantinuum.com/, December 04, 2021 - September 23, 2022.
  * Abhijith  _et al._ [2018] J. Abhijith, A. Adedoyin, J. Ambrosiano, P. Anisimov, A. Bärtschi, W. Casper, G. Chennupati, C. Coffrin, H. Djidjev, D. Gunter, S. Karra, N. Lemons, S. Lin, A. Malyzhenkov, D. Mascarenas, S. Mniszewski, B. Nadiga, D. O’Malley, D. Oyen, S. Pakin, L. Prasad, R. Roberts, P. Romero, N. Santhi, N. Sinitsyn, P. J. Swart, J. G. Wendelberger, B. Yoon, R. Zamora, W. Zhu, S. Eidenbenz, P. J. Coles, M. Vuffray, and A. Y. Lokhov, Quantum Algorithm Implementations for Beginners, arXiv e-prints , arXiv:1804.03719 (2018), arXiv:1804.03719 [cs.ET] .
  * Krol  _et al._ [2021] A. M. Krol, A. Sarkar, I. Ashraf, Z. Al-Ars, and K. Bertels, Efficient decomposition of unitary matrices in quantum circuit compilers (2021).
  * Gonzalez-Conde  _et al._ [2021] J. Gonzalez-Conde, Á. Rodríguez-Rozas, E. Solano, and M. Sanz, Simulating option price dynamics with exponential quantum speedup, arXiv e-prints , arXiv:2101.04023 (2021), arXiv:2101.04023 [quant-ph] .
  * Hamming [1986] R. W. Hamming,  _Numerical Methods for Scientists and Engineers, 2nd edition_ (Dover, New York, 1986) pp. 534–536.
  * Gottlieb  _et al._ [1992] D. Gottlieb, C.-W. Shu, A. Solomonoff, and H. Vandeven, On the gibbs phenomenon i: recovering exponential accuracy from the fourier partial sum of a nonperiodic analytic function, Journal of Computational and Applied Mathematics 43, 81 (1992).
  * Cai  _et al._ [1992] W. Cai, D. Gottlieb, and C.-W. Shu, On one-sided filters for spectral fourier approximations of discontinuous functions, SIAM Journal on Numerical Analysis 29, 905 (1992).
  * Kopriva [1987] D. A. Kopriva, A practical assessment of spectral accuracy for hyperbolic problems with discontinuities, Journal of Scientific Computing 2, 249 (1987).
  * Proctor  _et al._ [2021] T. Proctor, S. Seritan, K. Rudinger, E. Nielsen, R. Blume-Kohout, and K. Young, Scalable randomized benchmarking of quantum computers using mirror circuits (2021).
  * James  _et al._ [2001] D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White, Measurement of qubits, Physical Review A 64, 052312 (2001).
  * Takeda  _et al._ [2021] K. Takeda, A. Noiri, T. Nakajima, J. Yoneda, T. Kobayashi, and S. Tarucha, Quantum tomography of an entangled three-qubit state in silicon, Nature Nanotechnology 16, 965 (2021), arXiv:2010.10316 [cond-mat.mes-hall] .
  * Note [1] Of course, $\mathinner{|{f_{\text{tomo}}}\rangle}$ is only defined up to a global phase. We fix this phase by demanding that $\langle f|f_{\text{tomo}}\rangle$ is real and positive.
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics , 1 (2020).
  * Elben  _et al._ [2022] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox (2022).
  * Torlai  _et al._ [2018] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-network quantum state tomography, Nature Physics 14, 447 (2018).
  * Neugebauer  _et al._ [2020] M. Neugebauer, L. Fischer, A. Jäger, S. Czischek, S. Jochim, M. Weidemüller, and M. Gärttner, Neural-network quantum state tomography in a two-qubit experiment, Physical Review A 102 (2020).
  * Cha  _et al._ [2021] P. Cha, P. H. Ginsparg, F. Wu, J. F. Carrasquilla, P. L. McMahon, and E.-A. Kim, Attention-based quantum tomography, Machine Learning: Science and Technology 3 (2021).
  * Zuo  _et al._ [2021] Y. Zuo, C. Cao, N. Cao, X. Lai, B. Zeng, and S. Du, All-optical neural network quantum state tomography (2021).
  * Griffiths and Niu [1996] R. B. Griffiths and C.-S. Niu, Semiclassical fourier transform for quantum computation, Phys. Rev. Lett. 76, 3228 (1996).
  * h1- [2022a] Quantinuum system model h1 emulator (2022a), December 01, 2022.
  * Ramos-Calderer [2022] S. Ramos-Calderer, Efficient quantum interpolation of natural data, arXiv e-prints , arXiv:2203.06196 (2022), arXiv:2203.06196 [quant-ph] .
  * Kompenhans [2016] M. Kompenhans, Adaptation strategies for discontinuous galerkin spectral element methods by means of truncation error estimations (2016).
  * Marin-Sanchez  _et al._ [2021] G. Marin-Sanchez, J. Gonzalez-Conde, and M. Sanz, Quantum algorithms for approximate function loading, arXiv e-prints , arXiv:2111.07933 (2021), arXiv:2111.07933 [quant-ph] .
  * Rattew and Koczor [2022] A. G. Rattew and B. Koczor, Preparing Arbitrary Continuous Functions in Quantum Registers With Logarithmic Complexity, arXiv e-prints , arXiv:2205.00519 (2022), arXiv:2205.00519 [quant-ph] .
  * McArdle  _et al._ [2022] S. McArdle, A. Gilyén, and M. Berta, Quantum state preparation without coherent arithmetic (2022).
  * Kitaev and Webb [2008] A. Kitaev and W. A. Webb, Wavefunction preparation and resampling using a quantum computer, arXiv e-prints , arXiv:0801.0342 (2008), arXiv:0801.0342 [quant-ph] .
  * Rattew  _et al._ [2021] A. G. Rattew, Y. Sun, P. Minssen, and M. Pistoia, The efficient preparation of normal distributions in quantum registers, Quantum 5, 609 (2021).
  * Markov  _et al._ [2022] V. Markov, C. Stefanski, A. Rao, and C. Gonciulea, A Generalized Quantum Inner Product and Applications to Financial Engineering, arXiv e-prints , arXiv:2201.09845 (2022), arXiv:2201.09845 [quant-ph] .
  * Oseledets [2013] I. Oseledets, Constructive representation of functions in low-rank tensor formats, Constructive Approximation 37, 1 (2013).
  * Schollwoeck [2011] U. Schollwoeck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Ran [2020] S.-J. Ran, Encoding of matrix product states into quantum circuits of one- and two-qubit gates, Physical Review A (2020).
  * Rebentrost  _et al._ [2018] P. Rebentrost, A. Steffens, I. Marvian, and S. Lloyd, Quantum singular-value decomposition of nonsparse low-rank matrices, Phys. Rev. A 97, 012327 (2018).
  * Aharonov and Ta-Shma [2003] D. Aharonov and A. Ta-Shma, Adiabatic Quantum State Generation and Statistical Zero Knowledge, arXiv e-prints , quant-ph/0301023 (2003), arXiv:quant-ph/0301023 [quant-ph] .
  * Childs  _et al._ [2002] A. M. Childs, R. Cleve, E. Deotto, E. Farhi, S. Gutmann, and D. A. Spielman, Exponential algorithmic speedup by quantum walk, arXiv e-prints , quant-ph/0209131 (2002), arXiv:quant-ph/0209131 [quant-ph] .
  * Berry  _et al._ [2007] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Efficient Quantum Algorithms for Simulating Sparse Hamiltonians, Communications in Mathematical Physics 270, 359 (2007), arXiv:quant-ph/0508139 [quant-ph] .
  * Häner  _et al._ [2018] T. Häner, M. Roetteler, and K. M. Svore, Optimizing Quantum Circuits for Arithmetic, arXiv e-prints , arXiv:1805.12445 (2018), arXiv:1805.12445 [quant-ph] .
  * Haah [2019] J. Haah, Product Decomposition of Periodic Functions in Quantum Signal Processing, Quantum 3, 190 (2019), arXiv:1806.10236 [quant-ph] .
  * Dong  _et al._ [2021] Y. Dong, X. Meng, K. B. Whaley, and L. Lin, Efficient phase-factor evaluation in quantum signal processing, Phys. Rev. A 103, 042419 (2021), arXiv:2002.11649 [quant-ph] .
  * Chao  _et al._ [2020] R. Chao, D. Ding, A. Gilyen, C. Huang, and M. Szegedy, Finding Angles for Quantum Signal Processing with Machine Precision, arXiv e-prints , arXiv:2003.02831 (2020), arXiv:2003.02831 [quant-ph] .
  * Bauer  _et al._ [2021] C. W. Bauer, P. Deliyannis, M. Freytsis, and B. Nachman, Practical considerations for the preparation of multivariate gaussian states on quantum computers (2021).
  * Draper  _et al._ [2006] T. G. Draper, S. Kutin, E. M. Rains, and K. M. Svore, A logarithmic-depth quantum carry-lookahead adder, Quantum Inf. Comput. 6, 351 (2006).
  * Cleve and Watrous [2000] R. Cleve and J. Watrous, Fast parallel circuits for the quantum fourier transform (2000).
  * Araujo  _et al._ [2021] I. F. Araujo, D. K. Park, F. Petruccione, and A. J. da Silva, A divide-and-conquer algorithm for quantum state preparation, Scientific Reports 11 (2021).
  * Zhang  _et al._ [2022] X.-M. Zhang, T. Li, and X. Yuan, Quantum state preparation with optimal circuit depth: Implementations and applications (2022).
  * Stefanski  _et al._ [2022] C. Stefanski, V. Markov, and C. Gonciulea, Quantum Amplitude Interpolation, arXiv e-prints , arXiv:2203.08758 (2022), arXiv:2203.08758 [quant-ph] .
  * Fijany and Williams [1998] A. Fijany and C. P. Williams, Quantum Wavelet Transforms: Fast Algorithms and Complete Circuits, arXiv e-prints , quant-ph/9809004 (1998), arXiv:quant-ph/9809004 [quant-ph] .
  * Li  _et al._ [2018] H. Li, P. Fan, H. Xia, S. Song, and X. He, The multi-level and multi-dimensional quantum wavelet packet transforms, Scientific Reports 8 (2018).
  * Klappenecker and Roetteler [2001] A. Klappenecker and M. Roetteler, Discrete Cosine Transforms on Quantum Computers, arXiv e-prints , quant-ph/0111038 (2001), arXiv:quant-ph/0111038 [quant-ph] .
  * Kay [2018] A. Kay, Tutorial on the quantikz package, arXiv: Quantum Physics (2018).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9, 90 (2007).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Giardina and Chirlian [1972] C. Giardina and P. Chirlian, Bounds on the truncation error of periodic signals, IEEE Transactions on Circuit Theory 19, 206 (1972).
  * h1- [2022b] Quantinuum system model h1 data sheet (2022b), October 03, 2022.
  * Note [2] These specifications were provided by Quantinuum upon request.
  * Zhang  _et al._ [2014] Y. Zhang, K. Lu, and Y. Gao, Qsobel: A novel quantum image edge extraction algorithm, Science China Information Sciences 58, 1 (2014).
  * Yao  _et al._ [2017] X. Yao, H. Wang, Z. Liao, M.-C. Chen, J.-W. Pan, J. Y. Li, K. Zhang, X. Lin, Z. Wang, Z. Luo, W. Zheng, J. Li, M. sheng Zhao, X. Peng, and D. Suter, Quantum image processing and its application to edge detection: Theory and experiment, arXiv: Quantum Physics (2017).
  * Le  _et al._ [2011b] P. Q. Le, F. Dong, and K. Hirota, A flexible representation of quantum images for polynomial preparation, image compression, and processing operations, Quantum Information Processing 10, 63 (2011b).
  * Owen [2020] M. Owen, Photograph of the moon by michael on unsplash (2020).

##  Appendix A Derivation of the Fourier Series Loader method

In Sec. (II), we claimed that the FSL method can load a Fourier series of one
or more variables into a quantum state. We demonstrated this capability by
loading a function of one and two variables on the Quantinuum H$1$-$1$ and
H$1$-$2$ quantum computers. In this appendix, we provide a detailed derivation
of the FSL method.

###  A.1 FSL method for a Fourier series of single variable

Consider a truncated Fourier series of the form

$\displaystyle f(x)\,=\,\sum_{p=-M}^{M}\,c_{p}\,e^{-i2\pi px}\,,$ (6)

where $M=2^{m}-1$ denotes the number of Fourier coefficients. Our goal is then
to prepare the following state of $n$ qubits:

$\displaystyle\ket{f}_{n}\,=\,\sum_{k=0}^{2^{n}-1}\,f_{k}\,\ket{k}_{n}\,,$ (7)

where $f_{k}\,=\,f\left(x=k/2^{n}\right)\,=\,\sum_{p=-M}^{M}\,c_{p}\,e^{-i2\pi
pk/2^{n}}\,$. Here, the subscript $n$ in $\ket{\hskip 2.84526pt\cdot\hskip
2.84526pt}_{n}$ specifies the number of qubits.

Let us assume the unitary $U_{c}$ in the quantum circuit in Fig. (1a) acts on
$(m+1)$ qubits and maps $\ket{0}^{\otimes(m+1)}$ state to $\ket{\tilde{c}}$,
where

$\displaystyle\ket{\tilde{c}}_{m+1}\,:=\,2^{n/2}\sum_{p=0}^{M}c_{p}\ket{p}_{m+1}+2^{n/2}\sum_{p=1}^{M}c_{-p}\ket{2^{m+1}-p}_{m+1}\,.$
(8)

Then we claim that the circuit in Fig. (1a) maps initial $\ket{0}^{\otimes n}$
state to the target $\ket{f}_{n}$ state in Eq. (7). To verify this claim, we
traverse through the circuit in Fig. (1a) and analyze the state after each
step. After the action on $U_{c}$ on the lower $(m+1)$ qubits, the initial
$\ket{0}^{\otimes n}$ state becomes

$\displaystyle
2^{n/2}\sum_{p=0}^{M}c_{p}\ket{0}^{\otimes(n-m-1)}\otimes\ket{p}_{m+1}+2^{n/2}\sum_{p=1}^{M}c_{-p}\ket{0}^{\otimes(n-m-1)}\otimes\ket{2^{m+1}-p}_{m+1}\,,$
(9)

which can equivalently be written as

$\displaystyle
2^{n/2}\sum_{p=0}^{M}c_{p}\ket{0}^{\otimes(n-m-1)}\otimes\ket{0}\otimes\ket{p}_{m}+2^{n/2}\sum_{p=1}^{M}c_{-p}\ket{0}^{\otimes(n-m-1)}\otimes\ket{1}\otimes\ket{2^{m}-p}_{m}\,.$
(10)

The next step following the action of $U_{c}$ is the action of the cascade of
$n-m-1$ CNOT gates as shown in Fig. (1a). After the action of these CNOT
gates, the state of $n$ qubits becomes

$\displaystyle
2^{n/2}\sum_{p=0}^{M}c_{p}\ket{0}^{\otimes(n-m)}\otimes\ket{p}_{m}+2^{n/2}\sum_{p=1}^{M}c_{-p}\ket{1}^{\otimes(n-m)}\otimes\ket{2^{m}-p}_{m}\,,$
(11)

which can equivalently be written as

$\displaystyle
2^{n/2}\sum_{p=0}^{M}c_{p}\ket{p}_{n}+2^{n/2}\sum_{p=1}^{M}c_{-p}\ket{2^{n}-p}_{n}\,.$
(12)

The last step in the FSL circuit shown in Fig. (1a) is the action of the
inverse Quantum Fourier Transform (QFT). Recall that the inverse QFT acts on a
computational basis state as

$\displaystyle\left(\text{QFT}\right)^{\dagger}\ket{p}_{n}\,=\,\frac{1}{2^{n/2}}\,\sum_{k=0}^{2^{n}-1}\,e^{-i2\pi
pk/2^{n}}\,\ket{k}_{n}\,.$ (13)

This implies that the state of $n$ qubits after the action of the inverse QFT
becomes

$\displaystyle 2^{n/2}$
$\displaystyle\sum_{p=0}^{M}c_{p}\left(\text{QFT}\right)^{\dagger}\ket{p}_{n}+2^{n/2}\sum_{p=1}^{M}c_{-p}\left(\text{QFT}\right)^{\dagger}\ket{2^{n}-p}_{n}\,,$
$\displaystyle=$ $\displaystyle\sum_{p=0}^{M}c_{p}\sum_{k=0}^{2^{n}}e^{-i2\pi
kp/2^{n}}\ket{k}_{n}+\sum_{p=1}^{M}c_{-p}\sum_{k=0}^{2^{n}}e^{-i2\pi
k(2^{n}-p)/2^{n}}\ket{k}_{n}\,,$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{2^{n}}\,\sum_{p=-M}^{M}c_{p}e^{-i2\pi
kp/2^{n}}\,\ket{k}_{n}\,,$ (14)

which is precisely the target state in Eq. (7). This finishes our derivation
of the FSL method for single-variable functions.

###  A.2 FSL method for a multi-variate Fourier series

Here, we provide the details for extending the FSL method to loading higher-
dimensional functions and provide the corresponding quantum circuit. For
simplicity, we consider the case of a Fourier series of two variables.
Although, our construction can easily be generalized to Fourier series of more
than two variables.

Any arbitrary truncated Fourier Series of two variables can be written as

$\displaystyle f(x,y)\,=\,\sum_{p=-M}^{M}\sum_{q=-M}^{M}\,c_{p,q}\,e^{-i2\pi
px}\,e^{-i2\pi qy}\,,$ (15)

where $M=2^{m}-1$. Our goal is then to prepare the following state of $2n$
qubits ($n$ qubits for each variable):

$\displaystyle\ket{f}_{2n}\,=\,\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}\,f_{k,j}\,\ket{k}_{n}\otimes\ket{j}_{n}\,,$
(16)

where
$f_{k,j}\,=\,f\left(k/2^{n},j/2^{n}\right)\,=\,\sum_{p=-M}^{M}\sum_{q=-M}^{M}\,c_{p,q}\,e^{-i2\pi
pk/2^{n}}\,e^{-i2\pi qj/2^{n}}\,$. Here, the subscript $n$ in $\ket{\hskip
2.84526pt\cdot\hskip 2.84526pt}_{n}$ specifies the number of qubits.

Figure 7: The quantum circuit implementing the two-dimensional FSL method.
Similarly to the one-dimensional FSL circuit in Fig. (1), this circuit
consists of an unitary $U_{c}$ which loads the Fourier coefficients, a cascade
of CNOT gates, and inverse QFTs.

We claim that the state in Eq. (16) can be prepared using the circuit shown in
Fig. (7). In this circuit, the unitary operator $U_{c}$ is any unitary
operator that maps $\ket{0}^{\otimes 2(m+1)}$ to state
$\ket{\tilde{c}}_{2(m+1)}$, where

$\displaystyle\ket{\tilde{c}}_{2(m+1)}\,=$
$\displaystyle\,\,2^{n}\,\sum_{p=0}^{M}\sum_{q=0}^{M}\,c_{p,q}\,\ket{p}_{m+1}\otimes\ket{q}_{m+1}\,+\,2^{n}\,\sum_{p=0}^{M}\sum_{q=1}^{M}\,c_{p,-q}\,\ket{p}_{m+1}\otimes\ket{2^{m+1}-q}_{m+1}$
(17) $\displaystyle+$
$\displaystyle\,\,\,2^{n}\,\sum_{p=1}^{M}\sum_{q=0}^{M}\,c_{-p,q}\,\ket{2^{m+1}-p}_{m+1}\otimes\ket{q}_{m+1}\,+\,2^{n}\,\sum_{p=1}^{M}\sum_{q=1}^{M}\,c_{-p,-q}\,\ket{2^{m+1}-p}_{m+1}\otimes\ket{2^{m+1}-q}_{m+1}\,.$

Note that this state is analogous to the state $\ket{\tilde{c}}$ in Eq. (8).

To verify our claim, we traverse through the circuit in Fig. (7) and analyze
the state after each step: the initial state of the $2n$ qubits is
$\ket{0}^{\otimes 2n}$. Then, after applying the $U_{c}$ operator, the state
of $2n$ qubits is

$\displaystyle\,2^{n}\,\sum_{p=0}^{M}\sum_{q=0}^{M}\,c_{p,q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{p}_{m+1}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{q}_{m+1}\,$
(18) $\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=0}^{M}\sum_{q=1}^{M}\,c_{p,-q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{p}_{m+1}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{2^{m+1}-q}_{m+1}$
$\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=1}^{M}\sum_{q=0}^{M}\,c_{-p,q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{2^{m+1}-p}_{m+1}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{q}_{m+1}\,$
$\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=1}^{M}\sum_{q=1}^{M}\,c_{-p,-q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{2^{m+1}-p}_{m+1}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{2^{m+1}-q}_{m+1}\,.$

The above state can equivalently be written as

$\displaystyle\,2^{n}\,\sum_{p=0}^{M}\sum_{q=0}^{M}\,c_{p,q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{0}\otimes\ket{p}_{m}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{0}\otimes\ket{q}_{m}\,$
(19) $\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=0}^{M}\sum_{q=1}^{M}\,c_{p,-q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{0}\otimes\ket{p}_{m}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{1}\otimes\ket{2^{m}-q}_{m}$
$\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=1}^{M}\sum_{q=0}^{M}\,c_{-p,q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{1}\otimes\ket{2^{m}-p}_{m}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{0}\otimes\ket{q}_{m}\,$
$\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=1}^{M}\sum_{q=1}^{M}\,c_{-p,-q}\,\ket{0}^{\otimes(n-m-1)}\otimes\ket{1}\otimes\ket{2^{m}-p}_{m}\otimes\ket{0}^{\otimes(n-m-1)}\otimes\ket{1}\otimes\ket{2^{m}-q}_{m}\,.$

Following the application of the controlled $X^{\otimes(n-m-1)}$ gates in Fig.
(7), the above state is mapped to

$\displaystyle\,2^{n}\,\sum_{p=0}^{M}\sum_{q=0}^{M}\,c_{p,q}\,\ket{0}^{\otimes(n-m)}\otimes\ket{p}_{m}\otimes\ket{0}^{\otimes(n-m)}\otimes\ket{q}_{m}\,$
(20) $\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=0}^{M}\sum_{q=1}^{M}\,c_{p,-q}\,\ket{0}^{\otimes(n-m)}\otimes\ket{p}_{m}\otimes\ket{1}^{\otimes(n-m)}\otimes\ket{2^{m}-q}_{m}$
$\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=1}^{M}\sum_{q=0}^{M}\,c_{-p,q}\,\ket{1}^{\otimes(n-m)}\otimes\ket{2^{m}-p}_{m}\otimes\ket{0}^{\otimes(n-m)}\otimes\ket{q}_{m}\,$
$\displaystyle+\,$
$\displaystyle\,2^{n}\,\sum_{p=1}^{M}\sum_{q=1}^{M}\,c_{-p,-q}\,\ket{1}^{\otimes(n-m)}\otimes\ket{2^{m}-p}_{m}\otimes\ket{1}^{\otimes(n-m)}\otimes\ket{2^{m}-q}_{m}\,.$

We now denote this state by $\ket{c}_{2n}$ and write it equivalently as

$\displaystyle\ket{c}_{2n}\,=$
$\displaystyle\,\,2^{n}\,\sum_{p=0}^{M}\sum_{q=0}^{M}\,c_{p,q}\,\ket{p}_{n}\otimes\ket{q}_{n}\,+\,2^{n}\,\sum_{p=0}^{M}\sum_{q=1}^{M}\,c_{p,-q}\,\ket{p}_{n}\otimes\ket{2^{n}-q}_{n}$
(21) $\displaystyle+$
$\displaystyle\,\,\,2^{n}\,\sum_{p=1}^{M}\sum_{q=0}^{M}\,c_{-p,q}\,\ket{2^{n}-p}_{n}\otimes\ket{q}_{n}\,+\,2^{n}\,\sum_{p=1}^{M}\sum_{q=1}^{M}\,c_{-p,-q}\,\ket{2^{n}-p}_{n}\otimes\ket{2^{n}-q}_{n}\,.$

Finally, we consider the action of the inverse QFTs. The state of $2n$ qubits
following the application of the inverse QFTs is given by

$\displaystyle\left(\text{QFT}\right)^{\dagger}\otimes\left(\text{QFT}\right)^{\dagger}\,\ket{c}_{2n}\,=$
$\displaystyle\,\,\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}\sum_{p=0}^{M}\sum_{q=0}^{M}\,c_{p,q}\,e^{-i2\pi
pk/2^{n}}e^{-i2\pi qj/2^{n}}\,\ket{k}_{n}\otimes\ket{j}_{n}\,$ (22)
$\displaystyle+$
$\displaystyle\,\,\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}\sum_{p=0}^{M}\sum_{q=1}^{M}\,c_{p,-q}\,e^{-i2\pi
pk/2^{n}}e^{i2\pi qj/2^{n}}\,\ket{k}_{n}\otimes\ket{j}_{n}$ $\displaystyle+$
$\displaystyle\,\,\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}\sum_{p=1}^{M}\sum_{q=0}^{M}\,c_{-p,q}\,e^{i2\pi
pk/2^{n}}e^{-i2\pi qj/2^{n}}\,\ket{k}_{n}\otimes\ket{j}_{n}\,$
$\displaystyle+$
$\displaystyle\,\,\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}\sum_{p=1}^{M}\sum_{q=1}^{M}\,c_{-p,-q}\,e^{i2\pi
pk/2^{n}}e^{i2\pi qj/2^{n}}\,\ket{k}_{n}\otimes\ket{j}_{n}\,.$

By combining terms, we can write this state as

$\displaystyle\left(\text{QFT}\right)^{\dagger}\otimes\left(\text{QFT}\right)^{\dagger}\,\ket{c}_{2n}\,=$
$\displaystyle\,\,\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}\left(\,\sum_{p=-M}^{M}\sum_{q=-M}^{M}\,c_{p,q}\,e^{-i2\pi
pk/2^{n}}e^{-i2\pi qj/2^{n}}\right)\,\ket{k}_{n}\otimes\ket{j}_{n}\,,$ (23)

which is exactly the desired state $\ket{f}_{2n}$ in Eq. (16). This verifies
our claim that the circuit in Fig. (7) can be used for the exact loading of a
Fourier series of two or more variables.

##  Appendix B Uniformly controlled rotations

In the main text, we proposed a cascade of uniformly controlled rotations
shown in Fig. (1b) as one possible method to load the Fourier coefficients to
the amplitudes of a quantum state $\ket{\tilde{c}}=U_{c}\ket{0}^{\otimes m}$.
This method has many desirable features. For example, an efficient
implementation of this circuit $U_{c}$ in terms of CNOT gates and rotation
gates is known [7]. Moreover, all the angles of rotation can be classically
determined given the target state [7]. Since this method plays a significant
role in the FSL method, we briefly review the results of [7] in this appendix.

As discussed in Sec. (II), the $m$-qubit circuit shown in Fig. (1b) can be
used to prepare an arbitrary $m$-qubit state $\ket{\psi}$. The circuit in Fig.
(1b) has three distinct parts, (i) a $R_{z}$ gate, (ii) uniformly controlled
$R_{y}$ gates, and (iii) uniformly controlled $R_{z}$ gates. Hence, the action
of this circuit can be succinctly written as

$\displaystyle\ket{\psi}\,=$
$\displaystyle\,\,\prod_{j=0}^{m-1}\,F^{(z)}_{j}\left[\boldsymbol{\alpha}_{m-1-j}^{(z)}\right]\,\cdot\,\prod_{j=0}^{m-1}\,F^{(y)}_{j}\left[\boldsymbol{\alpha}_{m-1-j}^{(y)}\right]\,\cdot\,R_{z}(-\phi_{\psi})\,\ket{0}^{\otimes
n}\,.$ (24)

In this equation,
$\boldsymbol{\alpha}^{(y/z)}_{j}=\left(\alpha^{(y/z)}_{j,0},\alpha^{(y/z)}_{j,1},\dots,\alpha^{(y/z)}_{j,2^{m-1-j}-1}\right)$
is a vector of size $2^{m-1-j}$, and
$F^{(y/z)}_{j}\left[\boldsymbol{\alpha}^{(y/z)}_{m-1-j}\right]$ is a uniformly
$j$-controlled rotation operator which acts on the $(j+1)^{\text{th}}$ qubit
as $R_{y/z}\big{(}\alpha^{(y/z)}_{m-1-j,k}\big{)}$ if the state of the first
$j$ qubits is $\ket{k}$ for $k=\\{0,1,\dots,2^{j}-1\\}$ [5, 6, 7]. A pictorial
representation of $F^{(y)}_{j}$ for $j=3$ is given in Fig. (8).

Given the $m$-qubit target state
$\ket{\psi}\,=\,\sum_{i=0}^{2^{m}-1}|\psi_{i}|\,e^{i\omega_{i}}\,\ket{i}\,$,
the angles $\boldsymbol{\alpha}^{(y)}_{j}$ for $j=\\{0,1,\dots,m-1\\}$ in Eq.
(24) can be determined using [5, 6, 7]

$\displaystyle\alpha_{j,k}^{(y)}\,=\,2\arcsin\left[\frac{\sum_{\ell=0}^{2^{j}-1}\,|\psi_{(2k+1)2^{j}+\ell}|^{2}}{\sum_{\ell=0}^{2^{j+1}-1}\,|\psi_{k2^{j+1}+\ell}|^{2}}\right]\,.$
(25)

Moreover, the angles $\boldsymbol{\alpha}^{(z)}_{j}$ for
$j=\\{0,1,\dots,m-1\\}$ and $\phi_{\psi}$ can be determined using [7]

$\displaystyle\alpha_{j,k}^{(z)}\,=\,2^{-j}\,\sum_{\ell=0}^{2^{j}}\,\left(\omega_{(2k+1)2^{j}+\ell}-\omega_{k2^{j+1}+\ell}\right)\,,$
(26)

and

$\displaystyle\phi_{\psi}\,=\,2^{1-m}\,\sum_{j=0}^{2^{m}-1}\,\omega_{j}\,,$
(27)

respectively. Hence, given the target state $\ket{\psi}$, all the angles in
Eq. (24) can be evaluated on a classical computer.

The next step after finding the angles of uniformly controlled rotations, is
to efficiently implement the uniformly controlled gates in terms of simple
gates such as rotations and CNOT gates. It was shown in [7] that
$F_{j}^{(y/z)}$ for $j>0$ can be implemented using $2^{j}$ $R_{y/z}$ gates
acting on the $(j+1)^{\text{th}}$ qubit. The angle of the $k^{\text{th}}$
rotation for $k=\\{0,1,\dots,2^{j}-1\\}$ is given by [7]

$\displaystyle\theta_{m-1-j,k}^{(y/z)}\,=\,\sum_{\ell=0}^{2^{j}-1}\,M_{k\ell}\,\alpha_{m-1-j,\ell}^{(y/z)}\,;\quad\quad\quad
M_{k\ell}\,=\,2^{-j}\,\left(-1\right)^{b_{\ell}\cdot g_{k}}\,,$ (28)

where $b_{k}$ and $g_{k}$ denote the binary code and binary Gray code of the
integer $k$ respectively. In addition to the $2^{j}$ rotation gates, the
efficient implementation of $F_{j}^{(y/z)}$ also has $2^{j}$ CNOT gates acting
on the $(j+1)^{\text{th}}$ qubit as shown in Fig. (8). The last CNOT gate is
controlled by the first qubit. Whereas, the control qubit for the
$k^{\text{th}}$ CNOT, excluding the last CNOT, has to be determined based on
which bit of the binary Gray code of integer $k$ differs from that of $k-1$.
For example, the control qubit for the first CNOT in Fig. (8) is the
$3^{\text{rd}}$ qubit because the third bit in the Gray code of $1$ ($001$)
and $0$ ($000$) are different. Whereas, the control qubit for the second CNOT
is the $2^{\text{nd}}$ qubit because the second bit in the Gray code of $2$
($011$) and $1$ ($001$) are different.

Figure 8: Equivalent implementations of uniformly $3$-controlled rotations.
The implementation in terms of rotation gates and CNOT gates was derived in
Ref. [7].

The Python code used to implement the cascade of uniformly controlled
rotations and to perform the pre-processing step discussed above is provided
in the code repository accompanying this paper. This finishes our review of
Ref. [7].

##  Appendix C Benchmarking classical pre-processing time

Figure 9: Plots for the CPU wall clock time, $T_{\text{cl}}$, passed during
the classical pre-processing step of finding a circuit implementation of the
unitary $U_{c}$ using the uniformly controlled rotations (left) and using the
Schmidt decomposition method (right). We find that even though $T_{\text{cl}}$
scales polynomially with the number of Fourier coefficients $M=2^{m+1}$, the
values of $T_{\text{cl}}$ are reasonable for large enough values of $M$. For
example, $T_{cl}<6$ seconds for uniformly controlled rotations and $T_{cl}<3$
seconds for the Schmidt decomposition method for $M\sim 2000$ ($m=10)$.

An attractive feature of the FSL method is that the classical pre-processing
cost needed to implement the FSL circuit is constant with respect to the total
number of qubits, $n$. However, an important step during the classical pre-
processing involves compiling a circuit implementation of a unitary $U_{c}$
that loads the $M=2^{(m+1)}$ Fourier coefficients on $m+1$ qubits. This cost
scales polynomially with the number of Fourier coefficients. In fact, the
asymptotic cost of implementing $U_{c}$ either using the cascade of uniformly
controlled rotations as in Fig. (1b) or using the Schmidt decomposition method
as in Fig. (1c) scales as $O(M^{3})$. If this classical pre-processing cost
were too high, then our method would not be suitable for functions for which a
large number of Fourier coefficients are needed to approximate the target
function with the desired error. For this reason, we benchmark the classical
pre-processing time needed to implement the unitary $U_{c}$. More precisely,
we recorded the CPU wall clock time, $T_{\text{cl}}$, required to find a
circuit implementation of the unitary $U_{c}$ for randomly generated $2^{m+1}$
Fourier coefficients for $6\leq m\leq 14$, and the results are presented in
Fig. (9). We found that the classical pre-processing time, $T_{\text{cl}}$,
for $m=12$ is less than two minutes using the uniformly controlled rotations
and is less than a minute using the Schmidt decomposition. The benchmark was
performed on a $2020$ $13$-inch Macbook Pro with $3.2$ GHz M1 CPU and $8$ GB
main memory. This benchmark demonstrates that the classical pre-processing
time needed to implement the FSL method is computationally tractable even for
sufficiently large values of $m$.

Figure 10: Quantum circuit simulation of loading a piece-wise discontinuous
function into a quantum state of $n=20$ qubits using the FSL method with a
quarter million Fourier coefficients. The total CPU wall clock time for this
simulation was around $15$ minutes which included the classical pre-processing
time needed to find a circuit implementation of loading Fourier coefficients
and the time needed to simulate the execution of the circuit.

Furthermore, we recorded the total compute time needed to perform a classical
simulation of the FSL method’s quantum circuit loading the discontinuous
function we presented in Fig. (2f) with $m=17$ (i.e., approximately a quarter
million Fourier coefficients) and $n=20$. We found that it took a total of
around $15$ minutes ($13$ minutes for finding $U_{c}$ using the Schmidt
decomposition method and $2$ minutes for executing the circuit in Qiskit). The
results of the simulation are presented in Fig. (10).

##  Appendix D Error analysis

In the main text, we stated that the infidelity of the FSL method decays
exponentially with $m$ at a rate that depends on the smoothness of the target
function. In Fig. (11), we explicitly verify this exponential decay of
infidelity for all of the functions considered in Sec. (II). In particular, we
show that the infidelity decays as $O(2^{-m})$ for the discontinuous function
in Fig. (2f), whereas it decays as $O(2^{-3m})$ for continuous functions such
as $x^{x}$ in Fig. (2a), sinc function in Fig. (2b), reflected put option
function in Fig. (2c), hyperbolic tangent function in Fig. (2e), and the two-
dimensional sinc function considered in Fig. (3). Moreover, for an infinitely
differentiable function such as the wavefunction of the first excited state of
the quantum harmonic oscillator considered in Fig. (2d), the infidelity decays
doubly exponential in $m$.

In the rest of this appendix, we argue that the exponential decay of
infidelity is true in general. We do this by deriving a bound on infidelity as
a function of $m$. For simplicity, we only consider functions of a single
variable in our analysis, but the generalization to multivariate functions is
straightforward. Our derivation is similar to the analysis in Ref. [83] in
which the bounds on the truncation error, as a function of the number of terms
in the Fourier series, were obtained.

Note that any function encoding state
$\ket{f}\,=\,\sum_{k=0}^{2^{n}-1}\,f_{k}\ket{k}$ can exactly be written as
$\ket{f}\,=\,\left(\text{QFT}\right)^{\dagger}\ket{c}\,$, where

$\displaystyle\ket{c}\,=\,\sum_{k=0}^{2^{n-1}}\,c_{k}\ket{k}\,+\,\sum_{k=1}^{2^{n-1}-1}\,c_{-k}\ket{2^{n}-k}\,,$
(29)

and $c_{k}$ are the discrete Fourier coefficients of $f(x)$:

$\displaystyle
c_{k}\,=\,\frac{1}{2^{n/2}}\,\sum_{\ell=0}^{2^{n}-1}\,f_{\ell}\,e^{i2\pi
k\ell/2^{n}}\,.$ (30)

The normalization of the state $\ket{f}$ implies that the state $\ket{c}$ is
normalized. That is,

$\displaystyle\sum_{k=-2^{n-1}+1}^{2^{n-1}}\,|c_{k}|^{2}\,=\,1\,.$ (31)

The FSL method approximates the state $\ket{f}$ by the state
$\ket{f_{(m)}}\,=\,\left(\text{QFT}\right)^{\dagger}\ket{c_{(m)}}\,$, where

$\displaystyle\ket{c_{(m)}}\,=\,\mathcal{N}^{-1/2}\,\sum_{k=0}^{2^{m}-1}\,c_{k}\ket{k}\,+\,\mathcal{N}^{-1/2}\,\sum_{k=1}^{2^{m}-1}\,c_{-k}\ket{2^{n}-k}\,,$
(32)

and $\mathcal{N}$ is normalization constant:

$\displaystyle\mathcal{N}\,=\,\sum_{k=-2^{m}+1}^{2^{m}-1}\,|c_{k}|^{2}\,.$
(33)

The infidelity between the state $\ket{f}$ and $\ket{f_{(m)}}$ can be
calculated in terms of the Fourier coefficients as follows:

$\displaystyle\epsilon_{(m)}\,=$
$\displaystyle\,\,1-\left|\braket{f}{f_{(m)}}\right|^{2}\,,$ $\displaystyle=$
$\displaystyle\,\,1-\left|\braket{c}{c_{(m)}}\right|^{2}\,,$ $\displaystyle=$
$\displaystyle\,\,1-\mathcal{N}^{-1}\,\left(\sum_{k=-2^{m}+1}^{2^{m}-1}\,|c_{k}|^{2}\right)^{2}\,,$
$\displaystyle=$
$\displaystyle\,\,1-\sum_{k=-2^{m}+1}^{2^{m}-1}\,|c_{k}|^{2}\,,$
$\displaystyle=$
$\displaystyle\,\,\sum_{k=2^{m}}^{2^{n-1}}\,|c_{k}|^{2}\,+\,\sum_{k=2^{m}}^{2^{n-1}-1}\,|c_{-k}|^{2}\,,$
(34)

where we have used Eq. (31) in the last step and Eq. (33) in the step before
that.

Figure 11: Plots for infidelity, $\epsilon_{(m)}$, as a function of $m$ for
all of the functions considered in Sec. (II). The colored lines are used to
plot the values of the infidelities, whereas the colored dots denote the
linear function determined by linear fitting the values of the infidelities.
The slope of the linear fit function describing the infidelity for the
function in Fig. (2f) is approximately $1$, whereas the slope is approximately
$3$ for the functions in Fig. (2a), Fig. (2b), Fig. (2c), Fig. (2e), and Fig.
(3). Plot (b) illustrates the doubly exponential decay of infidelity for the
function in Fig. (2d).

In order to derive a bound on the infidelity, we first need to derive a bound
on the Fourier coefficients. To do this, we define a difference operator
$\Delta$ such that

$\displaystyle\Delta f_{\ell}\,=:\,f_{\ell+1}-f_{\ell}\,.$ (35)

From this definition, we deduce that the action of this difference operator on
a product of functions is given by

$\displaystyle\Delta(f_{\ell}g_{\ell})\,=\,(\Delta f_{\ell})\cdot
g_{\ell+1}\,+\,f_{\ell}\cdot(\Delta g_{\ell})\,.$ (36)

This is a discrete analog of the action of the derivative operator on a
product of functions. Now, if we take $g_{\ell}\,=\,e^{i2\pi k\ell/2^{n}}$ in
the above expression and rearrange the terms, we get

$\displaystyle f_{\ell}e^{i2\pi k\ell/2^{n}}\,=\,\frac{1}{e^{i2\pi
k/2^{n}}-1}\,\Delta(f_{\ell}e^{i2\pi k\ell/2^{n}})\,-\,\frac{1}{1-e^{-i2\pi
k/2^{n}}}\,(\Delta f_{\ell})\cdot e^{i2\pi k\ell/2^{n}}\,.$ (37)

Note we have applied the identity $\Delta e^{i2\pi
k\ell/2^{n}}\,=\,\left(e^{i2\pi k/2^{n}}-1\right)\,\cdot\,e^{i2\pi
k\ell/2^{n}}\,$ to obtain the above expression. Now using Eq. (30), we find
that the Fourier coefficient $c_{k}$ can be written as

$\displaystyle c_{k}\,=$
$\displaystyle\,\,\frac{1}{2^{n/2}}\,\frac{1}{e^{i2\pi
k/2^{n}}-1}\,\sum_{\ell=0}^{2^{n}-1}\Delta(f_{\ell}e^{i2\pi
k\ell/2^{n}})\,-\frac{1}{2^{n/2}}\,\frac{1}{1-e^{-i2\pi
k/2^{n}}}\,\sum_{\ell=0}^{2^{n}-1}\,(\Delta f_{\ell})\cdot e^{i2\pi
k\ell/2^{n}}\,,$ $\displaystyle=$
$\displaystyle\,\,\frac{1}{2^{n/2}}\,\frac{1}{e^{i2\pi
k/2^{n}}-1}\,\Big{(}f_{2^{n}}-f_{0}\Big{)}\,-\frac{1}{2^{n/2}}\,\frac{1}{1-e^{-i2\pi
k/2^{n}}}\,\sum_{\ell=0}^{2^{n}-1}\,(\Delta f_{\ell})\cdot e^{i2\pi
k\ell/2^{n}}\,,$ $\displaystyle=$
$\displaystyle\,\,-\frac{1}{2^{n/2}}\,\frac{1}{1-e^{-i2\pi
k/2^{n}}}\,\sum_{\ell=0}^{2^{n}-1}\,(\Delta f_{\ell})\cdot e^{i2\pi
k\ell/2^{n}}\,,$ (38)

where we have used the periodicity of $f$ to cancel the first term. Now, using
the triangle inequality, we get

$\displaystyle|c_{k}|\,=$
$\displaystyle\,\,\frac{1}{2^{n/2}}\,\frac{1}{2\sin(\pi
k/2^{n})}\,\left|\sum_{\ell=0}^{2^{n}-1}\,(\Delta f_{\ell})\cdot e^{i2\pi
k\ell/2^{n}}\right|\,,$ $\displaystyle\leq$
$\displaystyle\,\,\frac{1}{2^{n/2}}\,\frac{1}{2\sin\big{(}\pi
k/2^{n}\big{)}}\,\sum_{\ell=0}^{2^{n}-1}\,\big{|}(\Delta f_{\ell})\big{|}\,,$
$\displaystyle=$ $\displaystyle\,\,\frac{1}{2^{n/2}}\,\frac{1}{2\sin\big{(}\pi
k/2^{n}\big{)}}\,\Big{\|}\ket{\Delta f}\Big{\|}_{1}\,,$ (39)

where $\Big{\|}\ket{\Delta f}\Big{\|}_{1}$ is the $1$-norm of the state
$\ket{\Delta f}\,=\,\sum_{\ell=0}^{2^{n}-1}\,\Delta f_{\ell}\ket{\ell}\,$.

Due to this bound on the Fourier coefficients, the infidelity in Eq. (34) is
bounded by

$\displaystyle\epsilon_{(m)}\,\leq\,\frac{\,\Big{\|}\ket{\Delta
f}\Big{\|}_{1}^{2}}{2^{n+1}}\,\sum_{k=2^{m}}^{2^{n-1}}\frac{1}{\Big{(}\sin\big{(}\pi
k/2^{n}\big{)}\Big{)}^{2}}\,.$ (40)

Bounding the sum by an integral provides the bound

$\displaystyle\epsilon_{(m)}\,\leq$
$\displaystyle\,\,\frac{\Big{\|}\ket{\Delta
f}\Big{\|}_{1}^{2}}{2^{n+1}}\,\int_{k=2^{m}}^{2^{n-1}}dk\,\frac{1}{\Big{(}\sin\big{(}\pi
k/2^{n}\big{)}\Big{)}^{2}}\,,$ $\displaystyle=$
$\displaystyle\,\,\frac{1}{2\pi}\,\Big{\|}\ket{\Delta
f}\Big{\|}_{1}^{2}\,\cot\left(\pi 2^{m}/2^{n}\right)\,.$ (41)

The bound derived above holds for an arbitrary number of qubits, $n$. To find
the asymptotic limit of this bound, we first need to determine how the
$1$-norm of the state $\ket{\Delta f}$ scales with large values of $n$. Note
that the $1$-norm of $\ket{\Delta f}$ is determined by the values of
$f_{\ell}$ at local maximas, local minimas, and the end-points of the domain.
That is,

$\displaystyle\Big{\|}\ket{\Delta f}\Big{\|}_{1}\,=$
$\displaystyle\,\,\sum_{\ell=0}^{2^{n}-1}\,\big{|}(\Delta f_{\ell})\big{|}\,,$
$\displaystyle=$ $\displaystyle\,\,-\big{(}\text{sgn}(\Delta
f_{0})-\text{sgn}(\Delta
f_{2^{n}-1})\big{)}\,f_{0}\,+\,2\sum_{\ell\in\ell_{\text{max}}}f_{\ell}\,-\,2\sum_{\ell\in\ell_{\text{min}}}f_{\ell}\,,$
(42)

where $\ell_{\text{max}}$ is the set of points when $f_{\ell}$ is locally
maximum (i.e., $f_{\ell}>f_{\ell\pm 1}$) and $\ell_{\text{min}}$ is the set of
points when $f_{\ell}$ is locally minimum (i.e., $f_{\ell}<f_{\ell\pm 1}\,$).
The normalization condition for the state $\ket{f}$ implies that
$f_{\ell}\,\sim 1/\sqrt{2^{n}}$ for large values of $n$. This means that the
$1$-norm also scales as

$\displaystyle\Big{\|}\ket{\Delta f}\Big{\|}_{1}\sim\frac{C}{\sqrt{2^{n}}}\,,$
(43)

where $C>0$ is a constant. With this scaling result, we find that the bound on
$\epsilon_{(m)}$ approaches

$\displaystyle\epsilon_{(m)}\,\leq\,\,\frac{C^{2}}{\pi^{2}2^{m+1}}$ (44)

in the large limit as $n\to\infty$.

The bound in Eq. (44) is valid for any periodic target function, $f$. We can
obtain stronger bounds if we make certain assumptions about the continuity of
$\Delta^{p}f_{k}$, where

$\displaystyle\Delta^{p}f_{\ell}\,=:\,\Delta^{p-1}f_{\ell+1}-\Delta^{p-1}f_{\ell}\,,$
(45)

where $p$ is a positive integer, and $f_{\ell}$ is periodically extended
beyond $\ell\in\\{0,1,\cdots,2^{n}-1\\}$ according to
$f_{2^{n}+\ell}=f_{\ell}$. Let us assume that the $\Delta^{q}f_{\ell}$ for all
$q\leq p$ is continuous by which we mean that
$\Delta^{q}f_{\ell}\sim(1/2^{n})^{q+1/2}$ for $q\leq p$ and that
$\Delta^{q}f_{\ell}\sim(1/2^{n})^{p+1/2}$ for $q\geq p$. This implies that

$\displaystyle\Big{\|}\ket{\Delta^{p+1}f}\Big{\|}_{1}\sim\frac{C_{p}}{(2^{n})^{p+1/2}}\,.$
(46)

To see why this leads to a stronger bound, we repeat the analysis in Eq. (38)
$p$ times to get

$\displaystyle c_{k}\,=\,\frac{1}{2^{n/2}}\,\left(\frac{-1}{1-e^{-i2\pi
k/2^{n}}}\right)^{p+1}\,\,\sum_{\ell=0}^{2^{n}-1}\,(\Delta^{p+1}f_{\ell})\cdot
e^{i2\pi k\ell/2^{n}}\,.$ (47)

This leads to a stronger bound of the Fourier coefficients,

$\displaystyle|c_{k}|\,\leq\,\frac{1}{2^{n/2}}\,\frac{1}{\left(2\sin\big{(}\pi
k/2^{n}\big{)}\right)^{p+1}}\,\Big{\|}\ket{\Delta^{p+1}f}\Big{\|}_{1}\,,$ (48)

which leads to a stronger bound on infidelity,

$\displaystyle\epsilon_{(m)}\,\leq$
$\displaystyle\,\,\frac{\Big{\|}\ket{\Delta^{p+1}f}\Big{\|}_{1}^{2}\,}{2^{n+1+2p}}\,\sum_{k=2^{m}}^{2^{n-1}}\frac{1}{\Big{(}\sin\big{(}\pi
k/2^{n}\big{)}\Big{)}^{2p+2}}\,,$ $\displaystyle\leq$
$\displaystyle\,\,\frac{\Big{\|}\ket{\Delta^{p+1}f}\Big{\|}_{1}^{2}\,}{2^{n+1+2p}}\,\int_{k=2^{m}}^{2^{n-1}}dk\,\frac{1}{\Big{(}\sin\big{(}\pi
k/2^{n}\big{)}\Big{)}^{2p+2}}\,,$ $\displaystyle=$
$\displaystyle\,\,\frac{\Big{\|}\ket{\Delta^{p+1}f}\Big{\|}_{1}^{2}\,}{2^{2p+2}\pi}\,\left[B\Big{(}1;-\frac{p+1}{2},\frac{1}{2}\Big{)}\,-\,B\Big{(}\sin^{2}\left(\frac{\pi
2^{m}}{2^{n}}\right);-\frac{p+1}{2},\frac{1}{2}\Big{)}\right]\,,$ (49)

where $B(x;a,b)$ is the incomplete Beta function. Asymptotically, this bound
approaches

$\displaystyle\epsilon_{(m)}\,\leq\,\,\frac{C_{p}^{2}}{(2p+1)\pi^{2p+2}2^{(m+1)(2p+1)}}\,.$
(50)

This concludes our discussion of the upper bound on the infidelity.

##  Appendix E Details of experiments on the Quantiuum System Model H1 Quantum
Computers

In Sec. (III), we presented the results of the experiments we performed on the
Quantinuum H$1$-$1$ and H$1$-$2$ quantum computers. Here, we provide further
details about those experiments and their implementation.

Figure 12: The FSL circuit used to load the complex-valued function
$f(x)=\left(\cos(2\pi x)-1.5i\cos(6\pi x)\right)/\sqrt{13}$ on $3$-qubit
state. In order to verify the prepared state, we apply single qubit gates from
$\\{I,H,R_{x}(\pi/2)\\}$ to each qubit before measuring it in the
computational basis. We repeat this process for all operations in the set
$\\{I,H,R_{x}(\pi/2)\\}^{\otimes 3}$ as depicted in the above figure. The
results of these $27$ experiments are processed to reconstruct the density
matrix through maximum likelihood estimation.

The first experiment we performed involved loading a complex-valued function
$f(x)=\left(\cos(2\pi x)-1.5i\cos(6\pi x)\right)/\sqrt{13}\,$. Since this
function is already in the form of a Fourier series, the pre-processing step
of approximating a function by a truncated Fourier series is not needed.
Moreover, since there is no truncation error in this case, the FSL method can
load this function exactly. The FSL circuit used to load this function into a
$3$-qubit state is shown in Fig. (12). It is worthwhile to point out that this
circuit is a special case of the circuit shown in Fig. (1a) as the unitary
$U_{c}$ is acting on all of the $n=3$ qubits. This, of course, is only
possible for functions that are sufficiently sparse in Fourier space so that
all the non-zero Fourier coefficients can be loaded with a shallow circuit.

To verify the state prepared by the circuit in Fig. (12), we performed quantum
state tomography where we applied single-qubit gates from the set
$\\{I,H,R_{x}(\pi/2)\\}$ on each of the qubits before measuring them. The
reason for using $H$ rather than $R_{y}(\pi/2)$ for rotating the measurement
basis was the observation that the inverse QFT, used in the FSL circuit, ends
with a $H$ gate. Therefore, in some of the $27$ experiments, we were able to
ignore two successive $H$ gates due to their cancellation.

As we discussed in Sec. (III), we performed maximum likelihood estimation to
reconstruct the density matrix from the measurement results. Here, we discuss
the results of an alternate reconstruction of the state. Recall that a
$3$-qubit density matrix can be written as:
$\rho\,=\,\frac{1}{8}\,\sum_{\mu_{1}=0}^{3}\sum_{\mu_{2}=0}^{3}\sum_{\mu_{3}=0}^{3}\,s_{\mu_{1}\mu_{2}\mu_{3}}\,\sigma_{\mu_{1}}\otimes\sigma_{\mu_{2}}\otimes\sigma_{\mu_{3}}\,$,
where
$\\{\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3}\\}=\\{I,\sigma_{x},\sigma_{y},\sigma_{z}\\}$
are the standard Pauli matrices. Thus, the reconstruction of the state is
reduced to the reconstruction of the coefficients $s_{\mu_{1}\mu_{2}\mu_{3}}$.
Since
$s_{\mu_{1}\mu_{2}\mu_{3}}\,=\,\text{tr}\left(\rho\cdot\sigma_{\mu_{1}}\otimes\sigma_{\mu_{2}}\otimes\sigma_{\mu_{3}}\right)\,$
are just the expectation values of tensor products of Pauli operators, they
can be estimated from the results of the measurements in $27$ different bases.
We found that the matrix reconstructed in this way has a fidelity of $94.0\%$
which is almost the same as the fidelity we obtained through maximum
likelihood estimation. However, this matrix is not a valid density matrix as
it possesses some negative eigenvalues. Nevertheless, we find it encouraging
that the fidelities obtained using two different methods are almost in
agreement with each other.

In the rest of the experiments we performed, we used the FSL method to load
$\sqrt{f(x)}$ into a quantum state for a given function $f(x)$ and then
compared the amplitudes of the prepared state to the value of the target
function. The precise expressions for the functions we considered are:

$\displaystyle\text{Bi-modal Gaussian:}\quad\quad f(x)\,\sim$
$\displaystyle\,\,\left(1-\lambda\right)e^{-(x-1/4)^{2}/2\sigma^{2}}\,+\,\lambda
e^{-(x-3/4)^{2}/2\sigma^{2}}\,\quad\quad\quad\text{with}\,(\lambda,\sigma)=(0.3,0.1)\,,$
(51) $\displaystyle\text{Log-normal:}\,\,\,\,\,\,\,\quad\quad\quad\quad
f(x)\,\sim$
$\displaystyle\,\,\frac{1}{x}e^{-\left(\log\left(x/q\right)\right)^{2}/2\sigma^{2}}\,\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\text{with}\,(q,\sigma)=(0.2,0.5)\,\,,$
(52) $\displaystyle\quad\text{Lorentzian:}\quad\,\,\,\quad\quad\quad\quad
f(x)\,\sim$
$\displaystyle\,\,\frac{1}{1+(x-0.5)^{2}/\sigma^{2}}\,\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\text{with}\,\,\sigma=0.1\,\,,$
(53)
$\displaystyle\text{`Spiky':}\,\,\,\,\,\,\,\,\quad\quad\quad\quad\quad\quad
f(x)\,\sim$ $\displaystyle\,\,\left(\cos(4\pi x)\,+\lambda\,\cos(20\pi
x)\right)^{2}\quad\quad\quad\quad\quad\quad\quad\quad\quad\text{with}\,\,\lambda=2.5\,\,,$
(54)

where $\sim$ denotes equality up to normalization. In addition to these
functions of a single variable, we also considered the following function of
two variables:

$\displaystyle
f(x,y)\,\sim\,\exp\left(-(x-\mu_{1})^{2}/\sigma_{11}^{2}\,-\,(y-\mu_{1})^{2}/\sigma_{12}^{2}\right)\,+\,\lambda\,\exp\left(-(x-\mu_{2})^{2}/\sigma_{21}^{2}\,-\,(y-\mu_{2})^{2}/\sigma_{22}^{2}\right)\,,$
(55)

with $\mu_{1}=0.65$, $\mu_{2}=0.35$, $\sigma_{11}=\sigma_{22}=\sqrt{1/50}$,
$\sigma_{12}=\sqrt{1/40}$, $\sigma_{21}=\sqrt{1/30}$, and $\lambda=0.5\,$.

To load these functions on a quantum computer using the FSL method, we need to
find the first few dominant Fourier coefficients of the square root of these
functions. For bi-modal Gaussian function, log-normal distribution, and
Lorentzian function, we used the first $8$, $16$, and $16$ Fourier
coefficients, respectively. For the two-variable function, we used the first
$64$ Fourier coefficients. Once the Fourier coefficients are determined, we
need to find a suitable unitary $U_{c}$ which can load these Fourier
coefficients. For all the functions except the ‘spiky’ function in Eq. (54),
we implemented $U_{c}$ using the Schmidt-decomposition circuit shown in Fig.
(1c). We first determined the unitary operators $U$ and $V$ in Fig. (1c) by
performing the Schmidt decomposition of the state $\ket{\tilde{c}}$ from Eq.
(1). Then, we decomposed the unitary operators $U$ and $V$ into single and
two-qubit gates using the built-in transpiler available in Qiskit [31].

On the other hand, when loading the ‘spiky’ function, we do not need to
perform any such classical pre-processing given that only two Fourier modes
are active. This is another instance in which the target state is suitably
sparse in Fourier space such that we are able to construct a shallow quantum
circuit to load all the non-zero Fourier coefficients. The full FSL circuit
that we used to load the ‘spiky’ function is shown in Fig. (13).

Before we executed the FSL circuit on the Quantinuum H$1$-$1$ and H$1$-$2$
quantum computers, we reduced the gate count of the circuits as much as
possible in order to reduce the effect of hardware noise. One obvious but
crucial simplification that we exploited was to ignore all the SWAP gates that
appear in the implementation of the inverse QFT. Instead of applying these
SWAP gates, we manually changed the ordering of the qubits prior to applying
the inverse QFT without SWAPs.

Figure 13: The quantum circuit to load the ‘spiky’ function in Eq. (54) into a
quantum state of $5$-qubit. Since there were only two Fourier modes present in
this case, we were able to construct an efficient to load the Fourier
coefficients.