[1]Quantinuum, Partnership House, Carlisle Place, London SW1P 1BX, United
Kingdom

[2]Quantinuum, Terrington House, 13-15 Hills Road, Cambridge CB2 1NL, United
Kingdom

[3]Quantinuum K.K., Otemachi Financial City Grand Cube 3F, 1-9-2 Otemachi,
Chiyoda-ku, Tokyo, Japan

[4]Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
RIKEN, Wako, Saitama 351-0198, Japan

# Quantum state preparation for multivariate functions

Matthias Rosenkranz<EMAIL_ADDRESS>Eric Brunner  Gabriel
Marin-Sanchez  Nathan Fitzpatrick  Silas Dilkes  Yao Tang  Yuta Kikuchi
Marcello Benedetti<EMAIL_ADDRESS>
(31 May, 2024)

###### Abstract

A fundamental step of any quantum algorithm is the preparation of qubit
registers in a suitable initial state. Often qubit registers represent a
discretization of continuous variables and the initial state is defined by a
multivariate function. We develop protocols for preparing quantum states whose
amplitudes encode multivariate functions by linearly combining block-encodings
of Fourier and Chebyshev basis functions. Without relying on arithmetic
circuits, quantum Fourier transforms, or multivariate quantum signal
processing, our algorithms are simpler and more effective than previous
proposals. We analyze requirements both asymptotically and pragmatically in
terms of near/medium-term resources. Numerically, we prepare bivariate
Student’s t-distributions, 2D Ricker wavelets and electron wavefunctions in a
3D Coulomb potential, which are initial states with potential applications in
finance, physics and chemistry simulations. Finally, we prepare bivariate
Gaussian distributions on the Quantinuum H2-1 trapped-ion quantum processor
using 24 qubits and up to 237 two-qubit gates.

##  1 Introduction

State preparation is a computational problem that appears as a sub-routine in
many quantum algorithms. In Hamiltonian simulation, we approximate the real-
time evolution $\mathrm{e}^{-\mathrm{i}tH}\ket{\psi}$ of an initially prepared
state $\ket{\psi}$. Quantum algorithms for linear systems of equations
$A\ket{x}=\ket{b}$ require the preparation of $\ket{b}$ in the first place.
Other examples are found in the literature of quantum algorithms for partial
differential equations [1, 2], financial instrument pricing [3, 4, 5], quantum
chemistry simulations [6, 7, 8, 9, 10], high-energy or nuclear physics [11,
12, 13], and other applications [14].

The computational complexity of state preparation is well understood [15, 16,
17, 18, 19, 20] and it is connected to that of synthesizing circuits for
arbitrary unitary matrices. Algorithms for preparing generic $n$-qubit states
incur an unavoidable exponential scaling, either in terms of circuit depth
[21] or number of ancilla qubits [22]. Under certain assumptions, however, the
scaling can be improved. For the special case of $k$-sparse states, i.e.
states with $k$ non-zero amplitudes, there exists an efficient preparation
using $\log(nk)$ circuit depth and $\mathcal{O}(nk\log k)$ ancilla qubits
[22]. For more structured states, such as when the amplitudes are defined by a
function, there exist a number of efficient algorithms. For example, [23, 24,
25, 26] consider certain classes of probability density functions; [27, 28,
29] consider real-valued functions approximated by polynomials; [30] considers
complex-valued continuous functions. Quantum arithmetic circuits are among the
most expensive components in state preparation algorithms. Recent work has
aimed at limiting their use, potentially leading to a reduction of the number
of qubits and gates. When the amplitudes are unknown and provided by a black-
box, the use of arithmetic circuits can be limited to simple comparisons
against some initial reference state [31, 32, 33], or replaced by linear
combination of unitaries circuits [34]. When amplitudes are given by a known
function instead, quantum signal processing (QSP) [35, 36, 37, 38, 39] can be
used to avoid arithmetic circuits altogether [40].

In this work, we contribute to the state preparation problem by considering
the case where the amplitudes are defined by multivariate functions. The
multivariate case is seldom discussed in the literature and presents its
unique challenges. For example, multivariate versions of QSP appear to be much
more complicated and highly constrained [41, 42, 43]. We attack the
multivariate state preparation problem using a simple set of primitives and
without using arithmetic circuits. We make use of (i) function approximation
by truncated Fourier or Chebyshev series, (ii) block-encoding of the
corresponding basis functions, and (iii) linear combination of unitaries. The
resources required by our algorithms depend on three hyper-parameters: the
number of dimensions $D$, the number of bits $n$ used for the discretization
of each dimension, and the degree $d$ used for function approximation in each
dimension. We obtain a number of two-qubit gates that scales as
$\mathcal{O}(d^{D}+Dn\log d)$ for Fourier series and as
$\mathcal{O}\left(d^{D}+Ddn\log n\right)$ for Chebyshev series. The total
number of qubits scales as $\mathcal{O}\left(Dn+D\log d\right)$ for both
methods.

For the special case of Fourier series approximation, the closest work to ours
is the Fourier series loader (FSL) [44]. This method loads the Fourier
coefficients in the main qubit register, hence avoiding the use of ancilla
qubits. It then employs the quantum Fourier transform (QFT) to prepare the
target state using additional $\mathcal{O}(n^{2})$ gates. Related ideas of
using the QFT for Fourier series interpolation of a target function from a
coarse grid to a finer grid were earlier explored in [45]. In comparison, our
method uses ancilla qubits, but it is more general, and cuts down the number
of two-qubit gates by avoiding QFT. A related method called the Walsh series
loader (WSL) [46] approximates the target function by a Walsh series, which is
loaded on the quantum computer using techniques from [47]. While WSL is
shallower than FSL, it introduces a further approximation error which also
affects the probability of success. Our method does not have this additional
error.

##  2 Methods

We wish to prepare a quantum state with amplitudes proportional to a
multivariate target function. We consider complex target functions,
$f\colon[-1,1]^{D}\rightarrow\mathbb{C}$. For simplicity of exposition, let us
first discuss the one-dimensional case, $D=1$. We assume that $f$ has a
Fourier or Chebyshev series representation of the form
$f(x)=\sum_{k=-\infty}^{\infty}\hat{c}_{k}\mathrm{e}^{\mathrm{i}\pi kx}$ or
$f(x)=\sum_{k=0}^{\infty}\hat{c}_{k}T_{k}(x)$, respectively, with $T_{k}(x)$
the $k$-th Chebyshev polynomial of the first kind. For the Fourier case, this
requires the function to be periodic with period $2$ on the whole interval
$[-1,1]$. However, we also work with non-periodic target functions as Fourier
series by focusing on the interval $[0,1]$ and using the interval $[-1,0)$ for
its periodic extension (see Appendix B).

Next, we approximate the target function with a finite Fourier or Chebyshev
series of the form

$\displaystyle f_{d}^{F}(x)$
$\displaystyle=\sum_{k=-d}^{d}c_{k}\mathrm{e}^{\mathrm{i}\pi kx}$ (Fourier)
(1) and $\displaystyle f_{d}^{C}(x)$
$\displaystyle=\sum_{k=0}^{d}c_{k}T_{k}(x)$ $\displaystyle\quad\text{
(Chebyshev)},$ (2)

respectively. The largest index in the sum, $d$, is the _degree_ of the
polynomial (for which $c_{d}$ or $c_{-d}$ does not vanish). We use $f_{d}$ to
indicate either $f_{d}^{F}$ or $f_{d}^{C}$, which should be clear from
context. If the coefficients of the infinite series are known, one
approximation is truncation by setting $c_{k}=\hat{c}_{k}$ for all $|k|\leq
d$. Another approximation is interpolation such that $f_{d}$ matches $f$ in a
number of interpolation points. Computing accurate coefficients for this
interpolant can be done with the fast Fourier transform (FFT) with classical
preprocessing cost $\mathcal{O}(d\log d)$ (see Appendix A–B). A rule-of-thumb
is that the smoother the target function the faster its finite series
approximations converges with increasing $d$ [48, 49].

For the multivariate case, we consider approximations with products of the
corresponding basis functions. For example, in two dimensions

$\displaystyle f_{d_{x},d_{y}}^{F}(x,y)$
$\displaystyle=\sum_{k=-d_{x}}^{d_{x}}\sum_{l=-d_{y}}^{d_{y}}c_{k,l}\mathrm{e}^{\mathrm{i}\pi
kx}\mathrm{e}^{\mathrm{i}\pi ly}\quad\text{(Fourier)}$ (3) and $\displaystyle
f_{d_{x},d_{y}}^{C}(x,y)$
$\displaystyle=\sum_{k=0}^{d_{x}}\sum_{l=0}^{d_{y}}c_{k,l}T_{k}(x)T_{l}(y)\quad\text{(Chebyshev)}.$
(4)

This generalizes to $D$ dimensions with degrees
$\bm{d}=(d_{1},d_{2},\dots,d_{D})$. The following protocols prepare a quantum
state $\sum_{\bm{x}}g_{\bm{d}}(\bm{x})\ket{\bm{x}}$ with
$g_{\bm{d}}(\bm{x})=f_{\bm{d}}(\bm{x})/\sqrt{\sum_{\bm{x}}|f_{\bm{d}}(\bm{x})|^{2}}$
for all $\bm{x}$ on a grid of size $2^{n_{1}}\times
2^{n_{2}}\times\cdots\times 2^{n_{D}}$.

###  2.1 Discretization

In this work, each variable is assigned to its own ‘main’ register. We use the
notation $x$ for the variable and $n_{x}$ for the number of qubits in the
register. The variable is therefore discretized on a grid of $2^{n_{x}}$
points. In the univariate setting, we simplify the notation to $n=n_{x}$.

Recall that when working with Fourier series approximations, we focus on the
non-negative part of the domain, $x\in[0,1]$. In this case the uniform grid of
points is represented by the $n$-qubit diagonal matrix

$H^{\mathrm{F}}=\frac{1}{2^{n}-1}\matrixquantity(0&&&\\\ &1&&\\\ &&\ddots&\\\
&&&2^{n}-1).$ (5)

For example, with $n=2$ we have grid points $0,\tfrac{1}{3},\tfrac{2}{3},1$.

When working with Chebyshev series approximations we don’t need periodic
extensions. Thus we use the negative part of the domain as well, $x\in[-1,1]$.
The uniform grid of point is then represented by the $n$-qubit diagonal matrix

$H^{\mathrm{C}}=2H^{\mathrm{F}}-I^{\otimes{n}},$ (6)

where $I=\matrixquantity(1&0\\\ 0&1)$ is the identity matrix. For example,
with $n=2$ we have grid points $-1,-\tfrac{1}{3},\tfrac{1}{3},1$. It is
important to note that this grid is uniform and should not be confused with
the grid of Chebyshev nodes commonly used in polynomial interpolation.

In the following we often use $H$ to indicate either $H^{\mathrm{F}}$ or
$H^{\mathrm{C}}$, which should be clear from the context. Note that both
matrices (5) and (6) have the operator norm $\|H\|=1$. Grids for multivariate
functions are built by taking tensor products $H\otimes\cdots\otimes H$, and
this does not change the operator norm.

###  2.2 Block-encoding

To construct the Fourier and Chebyshev approximations we shall block-encode
the corresponding basis functions. We begin with the standard definition of
block-encoding. An $(n+m)$-qubit unitary $W_{L}$ is an
$(\lambda,m,\epsilon)$-block-encoding of $L\in\mathbb{C}^{2^{n}\times 2^{n}}$
if

$\left\|L-\lambda\left(\bra{0}^{\otimes m}\otimes I^{\otimes
n}\right)W_{L}\left(\ket{0}^{\otimes m}\otimes I^{\otimes
n}\right)\right\|\leq\epsilon,$ (7)

where $m$ is the number of ancilla qubits, $\lambda$ is a normalization
constant such that $\lambda\geq\|L\|-\epsilon$ and $\epsilon$ is the error.

We make use of the linear combination of unitaries (LCU) [50] to construct
block-encoding circuits. Assume the input matrix is given as a weighted sum of
$K$ unitary operators $\\{U_{k}\\}_{k}$, i.e.
$L=\sum_{k=0}^{K-1}\lambda_{k}U_{k}$. To deal with negative and complex
coefficients, we write them in polar form,
$\lambda_{k}=\mathrm{e}^{\mathrm{i}\gamma_{k}}|\lambda_{k}|$, and compute the
normalization $\lambda=\sum_{k=0}^{K-1}|\lambda_{k}|$. LCU uses
$m=\lceil\log_{2}K\rceil$ ancilla qubits and the following unitary operators

$\displaystyle A$
$\displaystyle=\sum_{k=0}^{K-1}\sqrt{\frac{|\lambda_{k}|}{\lambda}}\outerproduct{k}{0}^{\otimes
m}\otimes I^{\otimes n}+\text{ u.c. },$ (8) $\displaystyle B$
$\displaystyle=\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes U_{k},$ (9)
$\displaystyle C$
$\displaystyle=\sum_{k=0}^{K-1}\mathrm{e}^{\mathrm{i}\gamma_{k}}\outerproduct{k}{k}\otimes
I^{\otimes n},$ (10)

where u.c. stands for unitary completion. With these definitions, the operator
$A^{\dagger}CBA$ yields

$\left(\bra{0}^{\otimes m}\otimes I^{\otimes
n}\right)\;A^{\dagger}CBA\;\left(\ket{0}^{\otimes m}\otimes I^{\otimes
n}\right)=\\\ \frac{1}{\lambda}\sum_{k=0}^{K-1}\lambda_{k}U_{k},$ (11)

which is a $(\lambda,m,0)$-block-encoding of $L$.

Let us discuss the important special case where $U_{k}=U^{k}$ are powers of
some unitary $U$. When $k\geq 0$ is encoded to the state as a binary string
$\ket{k}=\ket{k_{m-1}}\otimes\dots\otimes\ket{k_{0}}$, corresponding to the
binary representation $k=2^{0}k_{0}+2^{1}k_{1}+\cdots+2^{m-1}k_{m-1}$, the $B$
operator in Eq. (9) takes the form

$\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes
U^{k}=\prod_{i=0}^{m-1}\left[\outerproduct{0}{0}_{i}\otimes I^{\otimes
n}+\outerproduct{1}{1}_{i}\otimes U^{2^{i}}\right]$ (12)

where we highlight the relevant qubits with a subscript and assume identity
gates on all other qubits. This simple identity is illustrated in Fig. 1. It
significantly reduces the resources required for LCU in our use case.

Note that the approximate functions we wish to encode, Eqs. (1) and (2), are
weighted sums of basis functions. In the following sections we implement the
$k$-th basis functions as the operator $U_{k}$ and then use LCU to block-
encode the whole weighted sum. The circuit construction for the operators $A$
and $C$ does not change depending on whether the Fourier or Chebyshev basis is
used, while the circuit construction for the operator $B$ does.

[row sep=0.6cm, between origins, column sep=0.2cm] & 1 1 1 1
1 1 1 1
[3]U^0 [3]U^1 [3]U^2 [3]U^3

= [row sep=0.6cm, between origins, column sep=0.2cm] & 2
1
[3]U^2^0 [3]U^2^1

Figure 1: Circuit identity Eq. (12) for $\sum_{k}\outerproduct{k}{k}\otimes
U^{k}$ with $K=4$ and $m=2$. [row sep=0.7cm, between origins, column
sep=0.2cm] & 1 1 1 1
1 1 1 1
[3]e^iπ0 H^F [3]e^iπ1 H^F [3]e^iπ2 H^F [3]e^iπ3 H^F

= [row sep=0.7cm, between origins, column sep=0.2cm] & 2 3 4
1 2 3
P(π20+023-1) P(π20+123-1)
P(π21+023-1) P(π21+123-1)
P(π22+023-1) P(π22+123-1)  Figure 2: Example of circuit synthesis for the
Fourier basis with Hamiltonian $H^{F}$, $a=2$ ancilla and $n=3$ main qubits.
The special structure of this Hamiltonian leads to an efficient circuit made
of $a\times n$ controlled-phase gates.

###  2.3 Fourier basis functions

We seek a circuit construction for the Fourier basis $B$ operator
$B^{F}=\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes U_{k}$ with $U_{k}$ a block-
encoding of the $k$-th Fourier basis function. As the univariate Fourier basis
in Eq. (1) is defined as $\\{\mathrm{e}^{\mathrm{i}\pi
k^{\prime}x}\\}_{k^{\prime}=-d}^{d}$, the $k^{\prime}$-th Fourier basis
function applied to $H^{F}$ is the real-time evolution operator
$U_{k^{\prime}}=U^{k^{\prime}}=\mathrm{e}^{\mathrm{i}\pi k^{\prime}H^{F}}$. By
the definition in Eq. (7), a unitary operator is a trivial $(1,0,0)$-block-
encoding of itself. However, we need to show that there exists an efficient
circuit to implement the operator for any value of $k^{\prime}$. It can be
verified by inspection that

$\displaystyle\mathrm{e}^{\mathrm{i}\pi
k^{\prime}H^{F}}=\bigotimes_{j=0}^{n-1}P_{j}\left(\frac{\pi
k^{\prime}2^{j}}{2^{n}-1}\right),$ (13)

where

$\displaystyle P(\theta)=\begin{pmatrix}1&0\\\
0&\mathrm{e}^{\mathrm{i}\theta}\end{pmatrix},$ (14)

is the phase shift gate and subscript $j$ indicates application of the gate to
the $j$-th qubit. This is a block-encoding circuit for any Fourier basis
function and uses only $n$ single-qubit phase shift gates.

To map this block-encoding to the special case discussed at the end of Sec.
2.2, we first write $U_{k}=U^{-d}U^{k}$ with $k=0,1,\dots,2d$. Next we
implement the operator
$\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes\mathrm{e}^{\mathrm{i}\pi kH^{F}}$
using Eq. (12) with $U=\mathrm{e}^{\mathrm{i}\pi H^{F}}$ and an ancilla
register of size $m=a=\log_{2}(2d+1)$, and synthesize the circuit efficiently
as

$\displaystyle\begin{split}&\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes\mathrm{e}^{\mathrm{i}\pi
kH^{F}}=\\\
&\prod_{i=0}^{a-1}\prod_{j=0}^{n-1}\left[\outerproduct{0}{0}_{i}\otimes
I_{j}+\outerproduct{1}{1}_{i}\otimes P_{j}\left(\frac{\pi
2^{j+i}}{2^{n}-1}\right)\right].\end{split}$ (15)

This can be verified using the factorization $\mathrm{e}^{\mathrm{i}\pi
kH^{F}}=\prod_{i=0}^{a-1}\mathrm{e}^{\mathrm{i}\pi k_{i}2^{i}H^{F}}$, followed
by Eq. (13), and using $P_{j}(0)=I_{j}$. An example circuit diagram is
provided in Fig. 2. Finally, we construct

$\displaystyle
B^{F}=U^{-d}\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes\mathrm{e}^{\mathrm{i}\pi
kH^{F}}$ (16)

by applying the circuit Eq. (15) followed by the circuit for $U^{-d}$ using
Eq. (13) one more time. We have thus obtained a controlled $(1,0,0)$-block-
encoding circuit for the Fourier basis functions. It requires only $a\times n$
controlled-phase gates and additional $n$ single-qubit gates from the circuit
for $U^{-d}$. A controlled-phase gate can be implemented with two controlled-
not (C$X$) gates plus single-qubit gates.

###  2.4 Chebyshev basis functions

We seek a circuit construction for the Chebyshev basis $B$ operator
$B^{C}=\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes U_{k}$ with $U_{k}$ a block-
encoding of the $k$-th Chebyshev polynomial. As the univariate Chebyshev basis
in Eq. (2) is defined as $\\{T_{k}(x)=\cos(k\cos^{-1}(x))\\}_{k=0}^{d}$, the
$k$-th Chebyshev basis function applied to $H^{\mathrm{C}}$ is not a unitary
operator. We proceed by first block-encoding $H^{\mathrm{C}}$, and then using
the algebra of block-encoding circuits to obtain $T_{k}(H^{C})$.

Using Eqs. (5) and (6), we obtain

$\displaystyle
H^{C}=\sum_{j=0}^{n-1}\left(-\frac{2^{j}}{2^{n}-1}\right)Z_{j}=\sum_{j=0}^{n-1}\left(\frac{2^{j}}{2^{n}-1}\right)X_{j}Z_{j}X_{j}.$
(17)

Here $Z_{j}$ and $X_{j}$ indicate the Pauli operators $Z=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}$ and $X=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}$ applied to
the $j$-th qubit, respectively. The result is a linear combination of $n$
unitaries with which we can define a $(1,b,0)$-block-encoding of
$H^{\mathrm{C}}$, $V=A^{\dagger}_{V}C_{V}B_{V}A_{V}$, where $C_{V}=I$ and

$\displaystyle A_{V}$
$\displaystyle=\sum_{j=0}^{n-1}\sqrt{\frac{2^{j}}{2^{n}-1}}\outerproduct{j}{0}^{\otimes
m}\otimes I^{\otimes n}+\text{ u.c. },$ (18) $\displaystyle B_{V}$
$\displaystyle=\sum_{j=0}^{n-1}\outerproduct{j}{j}\otimes X_{j}Z_{j}X_{j}.$
(19)

Such a block-encoding requires $m=b=\lceil\log_{2}n\rceil$ ancilla qubits.

(a) [row sep=0.7cm, between origins, column sep=0.2cm] & 1
[2]U
= [row sep=0.7cm, between origins, column sep=0.19cm] & [style=NormalRed,wire
style=dotted]1 1 1 [style=NormalRed,wire style=dotted]1 1
A [2]B C A^† R

(b) [row sep=0.7cm, between origins, column sep=0.2cm] & 1
[2]U^2
= [row sep=0.75cm, between origins, column sep=0.19cm] &
[style=NormalBlue,wire style=dotted]1 [style=NormalBlue,wire style=dotted]1
[style=NormalYellow,wire style=dotted]1 [style=NormalBlue,wire style=dotted]1
[style=NormalBlue,wire style=dotted]1 [style=NormalYellow,wire style=dotted]1
A [2]B C A^† R A [2]B C A^† R

Figure 3: Simplifications used in the block-encoding circuit of Chebyshev
basis functions. (a) For controlled-$U$, where $U$ is any LCU circuit, one can
remove the controls depicted in red. (b) For the specific case where $B^{2}=I$
and $C^{2}=I$ it is possible to simplify controlled-$U^{2}$ as well. One can
remove either the blue or the yellow colored controls, depending on which
option is more convenient. Our method uses controlled-$U^{2^{i}}$ and thus
benefits significantly from these techniques.

Next, Lemma 8 in Ref. [51] shows that a block-encoding $V$ can be used to
derive another block-encoding $U_{V}$ that satisfies the qubitization
condition

$\displaystyle U_{V}=\bigoplus_{\xi}\begin{pmatrix}\xi&-\sqrt{1-\xi^{2}}\\\
\sqrt{1-\xi^{2}}&\xi\end{pmatrix}_{\xi},$ (20)

where $\xi$’s are the eigenvalues of the block-encoded marix and each term is
represented in the basis
$\left\\{|G_{\xi}\rangle,|G_{\xi}^{\bot}\rangle\right\\}$ defined in Ref.
[51]. Moreover, Corollary 9 in Ref. [51] shows that in the special case where
$V^{2}=I$ (which holds true here given that $B_{V}^{2}=I$ and $C_{V}=I$) the
qubitized block-encoding takes a simple form

$\displaystyle U_{V}=\left(\big{(}2\outerproduct{0}{0}^{\otimes b}-I^{\otimes
b}\big{)}\otimes I^{\otimes n}\right)V.$ (21)

We now combine two results from the literature. The first result is Lemma 16
in Ref. [52]. It states that for positive integer $k$

$\begin{pmatrix}\xi&-\sqrt{1-\xi^{2}}\\\
\sqrt{1-\xi^{2}}&\xi\end{pmatrix}^{k}=\\\
\begin{pmatrix}T_{k}(\xi)&-\sqrt{1-\xi^{2}}\;u_{k-1}(\xi)\\\
\sqrt{1-\xi^{2}}\;u_{k-1}(\xi)&T_{k}(\xi)\end{pmatrix},$ (22)

where $u_{k}$ are Chebyshev polynomials of the second kind. For $k=0$ we
define the above expression as the identity matrix. The second result we use
is Lemma 54 in Ref. [36] which deals with products of block-encodings. In the
case where the encoded matrices are unitaries and their normalization constant
is $\lambda=1$, the product of $(1,b,0)$-block-encodings yields the
$(1,b,0)$-block-encoding of the product. Moreover, since $U_{V}$ is the direct
sum in Eq. (20), powers of $U_{V}$ preserve the qubitization condition. It
follows that

$\displaystyle U^{k}_{V}$
$\displaystyle=\bigoplus_{\xi}\begin{pmatrix}T_{k}(\xi)&-\sqrt{1-\xi^{2}}\;u_{k-1}(\xi)\\\
\sqrt{1-\xi^{2}}\;u_{k-1}(\xi)&T_{k}(\xi)\end{pmatrix}_{\xi},$ (23)

and $U^{0}_{V}=I$. This is a $(1,b,0)$-block-encoding of the $k$-th Chebyshev
polynomial applied to $H^{C}$, i.e. $T_{k}(H^{C})=(\bra{0}^{\otimes b}\otimes
I^{\otimes n})U_{V}^{k}(\ket{0}^{\otimes b}\otimes I^{\otimes n})$. Note that
the algebra of qubitized block-encodings is also used in Ref. [53] to
implement Krylov basis vectors and Lanczos’ method on a quantum computer. We
emphasize that our method works because $\lambda=\|H^{C}\|=1$. For
$\lambda\neq 1$ we would not obtain the correct result after rescaling the
eigenvalues of $H$ since $T_{k}(H/\lambda)\neq T_{k}(H)/\lambda$ (but see
Section 5.1 of Ref. [54] for approximate solutions when this issue arises).

To construct the finite Chebyshev series Eq. (2) from the block-encoding
circuit $U^{k}_{V}$ we insert $U^{k}_{V}$ into Eq. (9). To implement
$U^{k}_{V}$ controlled on the value of $k$ we again use the identity in Eq.
(12) with $m=a$ and obtain

$\begin{split}B^{C}&=\sum_{k=0}^{K-1}\outerproduct{k}{k}\otimes U^{k}_{V}\\\
&=\prod_{i=0}^{a-1}\left[\outerproduct{0}{0}_{i}\otimes I^{\otimes
n+b}+\outerproduct{1}{1}_{i}\otimes U^{2^{i}}_{V}\right].\end{split}$ (24)

We can save some resources when controlling $U_{V}$. First, we can remove the
control on $A_{V}$ and $A_{V}^{\dagger}$ since they cancel when no operation
is performed in between. Second, when controlling $U_{V}^{2}$, we can either
remove the control on $B_{V}$ or on $R=2\outerproduct{0}{0}^{\otimes
b}-I^{\otimes b}$. These improvements are shown in the circuit diagrams in
Fig. 3 and for our purposes we remove the control on $B_{V}$.

To provide a gate count for $B^{C}$ in Eq. (24), we start with a gate count
for $\outerproduct{1}{1}_{i}\otimes U_{V}^{2}$ recalling
$U_{V}=RA_{V}^{\dagger}C_{V}B_{V}A_{V}$. The operators $A_{V}$ and
$A_{V}^{\dagger}$ in $U_{V}^{2}$ can be implemented with a $b$-qubit state
preparation circuit with real coefficients, requiring $2^{b}$ C$X$ gates each
[17]. $B_{V}$ costs $n$ $b$-qubit controlled $Z$ gates, which we denote as
C${}^{b}Z$. Controlled $R$ costs one C${}^{b}Z$ gate. In total
$\outerproduct{1}{1}_{i}\otimes U_{V}^{2}$ costs $2(n+1)$ C${}^{b}Z$ gates and
$2\cdot 2^{b}$ C$X$ gates. For Eq. (24) we need to implement $2^{i-1}$ copies
of $\outerproduct{1}{1}_{i}\otimes U_{V}^{2}$. By summing up the gate counts
for each copy we get the following totals for $B^{C}$:
$\sum_{i=0}^{a-1}2^{i}(n+1)=(2^{a}-1)(n+1)$ C${}^{b}Z$ gates and
$\sum_{i=0}^{a-1}2^{i}\cdot 2^{b}=(2^{a}-1)2^{b}$ C$X$ gates.

We wish to reduce gate counts to two-qubit gates only. Each C${}^{b}Z$ gate
can be implemented with $\mathcal{O}(b)$ gates [55], giving a total for
$B^{C}$ of
$\mathcal{O}\left[\left(2^{a}-1\right)\left((n+1)b+2^{b}\right)\right]=\mathcal{O}(2^{a}nb)$
C$X$ gates. Here we kept only dominant terms in $d$ and $n$.

[row sep=0.7cm,between origins] [3]ancillas for
coefficients $\ket{0}^{\otimes a}$ & [3]A 3
[wires=9,steps=4,style=dotted,cap=round,inner sep=2pt]$B$ … [3]C [3]A^†
[3]$\bra{0}^{\otimes a}$
2 …
… 1
[2]ancillas for
block-encoding $\ket{0}^{\otimes b}$
(Chebyshev only) [6]U^2^0 [6]U^2^1 … [6]U^2^a-1 [2]$\bra{0}^{\otimes b}$
…
[4]main register $\ket{+}^{\otimes n}$ … [4]$\propto\sum_{x}f_{d}(x)\ket{x}$
…
…
…

Figure 4: Quantum circuit for state preparation in one dimension. The
amplitudes are defined by a finite series, such as Fourier or Chebyshev.
Circuit $U$ is the block-encoding of the corresponding basis functions and may
require $b$ additional ancilla qubits. Each bit in the coefficient register
drives the application of powers of the block-encoding circuit. For example,
in the subspace where the coefficient ancillas are in state $\ket{k}$, we
apply $U^{k}$ to the other two registers. The algorithm succeeds if we measure
the zero state on all $a+b$ ancilla qubits.

###  2.5 State preparation in 1 and 2 dimensions

We now utilize our building-blocks for the construction of quantum states
defined by truncated series expansions. In this section we drop the
superscripts $F$ for Fourier and $C$ for Chebyshev, since the protocols are
almost identical. Then, the $n$-qubit discretization of the domain ($2^{n}$
points) is denoted as $H$, and the block-encoding of the $k$-th basis function
is denoted as $U_{k}$.

We begin with one-dimensional functions $f$ approximated by either Eq. (1) or
Eq. (2), here both denoted by $f_{d}$. Without further assumptions, the number
of terms is $K=2d+1$ for Fourier series and $K=d+1$ for Chebyshev series. Each
basis function is associated to a complex coefficient $c_{k}$ which we rewrite
as $c_{k}=\mathrm{e}^{\mathrm{i}\phi_{k}}|c_{k}|$. Then

$\displaystyle
f_{d}(H)=\sum_{k=0}^{K-1}c_{k}U^{k}=\mathcal{N}\sum_{k}\mathrm{e}^{\mathrm{i}\phi_{k}}\frac{|c_{k}|}{\mathcal{N}}U^{k},$
(25)

where we normalize the coefficients with $\mathcal{N}=\sum_{k}|c_{k}|$. We
encode the coefficients using operators $A$ and $C$ (Eqs. (8) and (10)), and
an ancilla register of size $a=\lceil\log_{2}K\rceil$. We then use the
operator $B$ defined in Eq. (16) for Fourier and Eq. (24) for Chebyshev
noting, however, that the Chebyshev method requires additional $b$ ancilla
qubits. We thus obtain a $(\mathcal{N},a+b,0)$-block-encoding $A^{\dagger}CBA$
with

$\left(\bra{0}^{\otimes(a+b)}\otimes I^{\otimes
n}\right)\;A^{\dagger}CBA\;\left(\ket{0}^{\otimes(a+b)}\otimes I^{\otimes
n}\right)=\\\ \frac{1}{\mathcal{N}}f_{d}(H).$ (26)

We apply $A^{\dagger}CBA$ to $\ket{0}^{\otimes(a+b)}\otimes\ket{+}^{\otimes
n}$ with the uniform superposition $\ket{+}^{\otimes
n}=\frac{1}{\sqrt{2^{n}}}\sum_{x}\ket{x}$ on the main register. Upon measuring
the zero state for all the ancilla qubits, the protocol generates the
unnormalized state $\frac{1}{\mathcal{N}\sqrt{2^{n}}}\sum_{x}f_{d}(x)\ket{x}$.
The probability of success is given by
$p_{\mathrm{success}}=\frac{\sum_{x}|f_{d}(x)|^{2}}{\mathcal{N}^{2}2^{n}}$.
Renormalizing the state (dividing by $\sqrt{p_{\mathrm{success}}}$) we obtain
the desired result
$\sum_{x}g_{d}(x)\ket{x}:=\frac{1}{\sqrt{\sum_{x}|f_{d}(x)|^{2}}}\sum_{x}f_{d}(x)\ket{x}$.
The full circuit is sketched in Fig. 4.

[row sep=0.75cm, between origins] $\ket{0}$ & a_x [2]A 2
[wires=4,steps=2,style=dotted,cap=round,inner sep=0pt,label
style=yshift=-1pt]$B_{X}$ … [2]C [2]A^†
$\ket{0}$ a_y 3 [wires=5,steps=2,style=dotted,cap=round,inner sep=0pt,label
style=yshift=-2pt]$B_{Y}$ …
$\ket{0}$ b_x [2]U …
$\ket{+}$ n_x …
$\ket{0}$ b_y [2]U …
$\ket{+}$ n_y …

Figure 5: Quantum circuit for state preparation in two dimensions.

With the one-dimensional warm up, we have all the elements to construct states
defined by bivariate functions $f(x,y)$. Again we start by discretizing the
variables using $n_{x}$ and $n_{y}$ qubits in the main registers, which gives
us a two-dimensional grid of $2^{n_{x}}\times 2^{n_{y}}$ points. We then
approximate $f$ as a finite series of degrees $d_{x}$ and $d_{y}$. The
coefficients are rewritten as
$c_{k,l}=\mathrm{e}^{\mathrm{i}\phi_{k,l}}|c_{k,l}|$ where $k=0,\dots,K_{x}-1$
and $l=0,\dots,K_{y}-1$, and $K_{x}$, $K_{y}$ depend on the series type and
degree as before. The normalization constant becomes
$\mathcal{N}=\sum_{k}\sum_{l}|c_{k,l}|$. The LCU operators are now defined as:

$\displaystyle A$
$\displaystyle=\sum_{k=0}^{K_{x}-1}\sum_{l=0}^{K_{y}-1}\sqrt{\frac{|c_{k,l}|}{\mathcal{N}}}\outerproduct{k}{0}^{\otimes{a_{x}}}\otimes\outerproduct{l}{0}^{\otimes{a_{y}}}\otimes
I^{\otimes(b_{x}+b_{y}+n_{x}+n_{y})}+\text{ u.c. },$ (27) $\displaystyle
B_{X}$ $\displaystyle=\sum_{k=0}^{K_{x}-1}\outerproduct{k}{k}\otimes
I^{\otimes a_{y}}\otimes U^{k}\otimes I^{\otimes(b_{y}+n_{y})},$ (28)
$\displaystyle B_{Y}$ $\displaystyle=\sum_{l=0}^{K_{y}-1}I^{\otimes
a_{x}}\otimes\outerproduct{l}{l}\otimes I^{\otimes(b_{x}+n_{x})}\otimes
U^{l},$ (29) $\displaystyle C$
$\displaystyle=\sum_{k=0}^{K_{x}-1}\sum_{l=0}^{K_{y}-1}\mathrm{e}^{\mathrm{i}\phi_{k,l}}\outerproduct{k}{k}\otimes\outerproduct{l}{l}\otimes
I^{\otimes(b_{x}+b_{y}+n_{x}+n_{y})}.$ (30)

Applying the circuit to the uniform superposition of both main registers and
post-selecting on the zero state of all the ancillas yields the desired state

$\begin{split}\left(\bra{0}^{\otimes(a_{x}+a_{y}+b_{x}+b_{y})}\otimes
I^{\otimes(n_{x}+n_{y})}\right)A^{\dagger}CB_{Y}B_{X}A\left(\ket{0}^{\otimes(a_{x}+a_{y}+b_{x}+b_{y})}\otimes\ket{+}^{\otimes(n_{x}+n_{y})}\right)=\\\
\frac{1}{\mathcal{N}\sqrt{2^{n_{x}+n_{y}}}}\sum_{x}\sum_{y}f_{d}(x,y)\ket{x}\otimes\ket{y}.\end{split}$
(31)

The protocol succeeds with probability

$\displaystyle
p_{\mathrm{success}}=\frac{\sum_{x}\sum_{y}|f_{d}(x,y)|^{2}}{\mathcal{N}^{2}2^{n_{x}+n_{y}}}.$
(32)

To obtain the desired result we renormalize Eq. (31) by dividing by
$\sqrt{p_{\text{success}}}$. The full circuit for bivariate state preparation
is sketched in Fig. 5. The generalization of this protocol to arbitrary number
of variables is straightforward.

###  2.6 Resource scaling for state preparation in arbitrary dimension

We can extrapolate the resource costs given in the previous sections for the
one-dimensional case to estimate the scaling of state preparation in arbitrary
dimensions. We focus on the number of qubits and C$X$ gate count as a proxy
for the two-qubit gate count, which is the main resource cost for our
algorithm on near-term hardware. We neglect single-qubit gates.

Let us recall the resource requirements for state preparation with one
variable ($D=1$). The operator $B$ costs $\mathcal{O}(an)$ C$X$ gates for the
Fourier approach (Sec. 2.3) and $\mathcal{O}\left(2^{a}nb\right)$ C$X$ gates
for the Chebyshev approach (Sec. 2.4). The size of the ancilla register $a$ is
determined by the number of coefficients in the approximating series $f_{d}$,
namely, $a=\lceil\log_{2}(2d+1)\rceil$ for Fourier series and
$a=\lceil\log_{2}(d+1)\rceil$ for Chebyshev series. The block-encodings of the
Chebyshev basis functions require an additional ancilla register of size
$b=\lceil\log_{2}n\rceil$. $A$, $A^{\dagger}$, Eq. (8), and $C$, Eq. (10), can
be implemented with $a$-qubit circuits requiring
$\mathcal{O}\left(2^{a}\right)$ C$X$ gates each.

| Fourier | Chebyshev
---|---|---
Target function | $[0,1]^{D}\rightarrow\mathbb{C}$ | $[-1,1]^{D}\rightarrow\mathbb{C}$
Classical preprocessing | Cost of computing coefficients, e.g. $\mathcal{O}(Dd^{D}\log d)$ for FFT
Main qubits | $Dn$ | $Dn$
Ancilla qubits | $D\lceil\log_{2}(2d+1)\rceil$ | $D\lceil\log_{2}(d+1)\rceil+D\lceil\log_{2}n\rceil$
Two-qubit gates | $\mathcal{O}(d^{D}+Dn\log d)$ | $\mathcal{O}(d^{D}+Ddn\log n)$
Success probability | $\mathcal{N}^{-2}2^{-Dn}\sum_{\bm{x}}|f_{d}(\bm{x})|^{2}$ | $\mathcal{N}^{-2}2^{-Dn}\sum_{\bm{x}}|f_{d}(\bm{x})|^{2}$
Table 1: Resources of our algorithm for preparing a quantum state with
amplitudes proportional to a $D$-dimensional function represented on an
equidistant grid of total size $2^{Dn}$. The method uses a Fourier or
Chebyshev polynomial approximation $f_{d}$ with maximal degree $d$ in each
dimension. $\mathcal{N}=\sum_{\bm{k}}|c_{\bm{k}}|$ is a normalization
dependent on the expansion coefficients $c_{\bm{k}}$.

We now provide the resource requirements for state preparation in arbitrary
dimension. First, we assume that there are $D$ variables, each discretized
with $n$ qubits, i.e. a total of $Dn$ qubits in the main registers. Second, we
assume the same degree $d$ for each variable. The coefficients are encoded in
ancilla registers $a_{1},a_{2},\dots,a_{D}$ each with size $a$. Then the total
number of ancilla qubits is $Da=D\lceil\log_{2}(2d+1)\rceil$ and
$D(a+b)=D(\lceil\log_{2}(d+1)\rceil+\lceil\log_{2}n\rceil)$ for the Fourier
and Chebyshev approaches, respectively, where we inserted the respective
expressions from the 1D case for $a$ and $b$. The qubit counts are summarized
in Table 1.

In general, the coefficients in a multivariate approximation $f_{d}$ may take
arbitrary values. Then the operators $A$, $A^{\dagger}$ and $C$ can be
implemented with $Da$-qubit circuits requiring
$\mathcal{O}\left(2^{Da}\right)$ C$X$ gates each. We have $D$ operators
$B_{1},B_{2},\dots,B_{D}$. Each $B_{i}$ is controlled only on the
corresponding coefficient ancilla register $a_{i}$ and acts only on the
corresponding main register $n_{i}$. For the Chebyshev case each $B_{i}$ also
acts on the corresponding block-encoding ancilla register $b_{i}$, cf. Fig. 5.
Hence, implementing all $B_{i}$ is $D$ times the 1D implementation cost of
$B$: $\mathcal{O}(Dan)$ C$X$ gates for the Fourier approach and
$\mathcal{O}\left(D2^{a}nb\right)$ C$X$ gates for the Chebyshev approach.
Summing the costs for $A$, $A^{\dagger}$, $C$ and all $B_{i}$ and inserting
the values for $a$ and $b$ yields
$\mathcal{O}\left(2^{D\lceil\log_{2}(2d+1)\rceil}+Dn\lceil\log_{2}(2d+1)\rceil\right)$
C$X$ gates for the Fourier approach and
$\mathcal{O}\left(2^{D\lceil\log_{2}(d+1)\rceil}+Dn\lceil\log_{2}n\rceil
2^{\lceil\log_{2}(d+1)\rceil}\right)$ C$X$ gates for the Chebyshev approach.
The simplified expressions $\mathcal{O}\left(d^{D}+Dn\log d\right)$ and
$\mathcal{O}\left(d^{D}+Ddn\log n\right)$, respectively, are provided as the
two-qubit gate count in Table 1.

Some remarks are in order. First, if we wanted to boost the probability of
success to unity we could use amplitude estimation. That would increase the
number of gates by a multiplicative factor of
$\mathcal{O}\left(p_{\text{success}}^{-1/2}\right)$. Second, the degree $d$
for each dimension must be chosen depending on the desired accuracy
$\varepsilon$. The scaling of $d(\varepsilon)$ is problem-dependent and not
analyzed in this work. However, the protocol is efficient only if
$d\in\mathcal{O}(\mathrm{poly}(n))$, which limits the accuracy that can be
achieved in general.

##  3 Results

With our first numerical simulations we visually inspect the block-encodings
of Fourier and Chebyshev basis functions. In Fig. 6(a) we plot the diagonal of
Eq. (13) for $k=0,\dots,4$. In Fig. 6(b) we show the diagonal of
$(\bra{0}^{\otimes b}\otimes I^{\otimes n})U^{k}(\ket{0}^{\otimes b}\otimes
I^{\otimes n})$ where $U^{k}$ is given in Eq. (23). In both cases we use
$n=5$, leading to a $x$-axis discretization with $32$ points. As expected the
data points closely track $\mathrm{e}^{\mathrm{i}k\pi H}$ in Fig. 6(a), and
$T_{k}(H)$ in Fig. 6(b).

Figure 6: Block-encodings of the Fourier and Chebyshev basis functions.
Markers are the diagonal of the Fourier block encoding, Eq. (13) in (a), and
$(\bra{0}^{\otimes b}\otimes I^{\otimes n})U^{k}(\ket{0}^{\otimes b}\otimes
I^{\otimes n})$ with Chebyshev block encoding $U^{k}$, Eq. (23) in (b), for
$n=5$ qubits (32 grid points). Lines are the corresponding continuous basis
functions $\mathrm{e}^{\mathrm{i}\pi kx}$ and $T_{k}(x)$.

###  3.1 2D Ricker wavelet via Chebyshev series

To illustrate the Chebyshev approach, we consider a 2D Ricker wavelet centered
at $(0,0)$,

$f(x,y)=\frac{1}{\pi\sigma^{4}}\left(1-\frac{x^{2}+y^{2}}{2\sigma^{2}}\right)\mathrm{e}^{-\frac{x^{2}+y^{2}}{2\sigma^{2}}}.$
(33)

Wavelets of this form are used, e.g. in image or seismic data applications
[56, 57, 58].

Figure 7: Resources for preparing an approximation to the 2D Ricker wavelet
with products of Chebyshev polynomials of degree $d$ in each dimension. (a)
Discretized target function and the absolute difference to the rescaled
amplitudes from a noiseless simulation of the algorithm for $d=7$. (b) Maximum
absolute error of the prepared state evaluated at the grid points for
different grid sizes. The difference between the prepared state and the
Chebyshev approximation is at machine precision (not shown). Hence, the
maximum error of the noiseless simulation against the target function is upper
bounded by the maximum error of the degree-$d$ Chebyshev interpolant over the
(non-discretized) domain (grey line – floored at machine precision). The
dashed line illustrates a super-exponential asymptotic error bound
$\mathcal{O}\left(\mathrm{e}^{-d\log d}\right)$. (c) Two-qubit gate count for
different grid sizes compiled with pytket optimization level 2 for the H2-1
default gate set. The line markers are a least-square fit of the data (round
markers) to the worst-case scaling with $\alpha_{0}=-29.3,-47.3,-49.3$,
$\alpha_{1}=3.14,4.47,4.39$, $\alpha_{2}=3.00,3.00,3.00$ for $n=4,5,6$ and
$D=2$. (d) Empirical (markers) and asymptotic success probability for
$n\rightarrow\infty$ and $d\rightarrow\infty$ (line).

We interpolate the function on the domain $[-1,1]^{2}$ with the Chebyshev
series Eq. (4) and compute its coefficients via the fast Fourier transform
(see Appendix A). We choose the same degree and same grid size in each
variable, i.e. $d_{x}=d_{y}=d$ and $n_{x}=n_{y}=n$. Figure 7 shows results for
approximating the 2D Ricker wavelet with $\sigma=0.5$ by a Chebyshev series of
varying degree $d$ and for different grid sizes $2^{n}\times 2^{n}$. Figure
7(a) illustrates the discretized target function and its absolute difference
to the amplitudes of the prepared quantum state for degree $d=7$ and grid size
$2^{4}\times 2^{4}$. For this comparison we rescale the amplitudes from a
noiseless simulation of the algorithm appropriately,

$\tilde{g}_{d}(x,y)=g_{d}(x,y)\sqrt{\sum_{i=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}|f(x_{i},y_{i})|^{2}},$
(34)

where $x_{i},y_{i}$ are the evaluated grid points. Figure 7(b) shows the
maximum error $\max_{(x,y)}|f(x,y)-\tilde{g}_{d}(x,y)|$ where the maximum is
over grid points. Increasing the degree leads to rapid convergence of the
amplitudes to the target function. We observe that the error between the
rescaled amplitudes and the Chebyshev approximation $f_{d}$ is at machine
precision (not shown). Hence, the convergence of our method to the target
function is determined by the convergence of the Chebyshev approximation. The
grey line shows the maximum error of the Chebyshev approximation,
$\max_{(x,y)\in[-1,1]^{2}}|f(x,y)-f_{d}(x,y)|$, evaluated numerically over the
whole domain (not just the grid) via global optimization. The line is
consistent with supergeometric convergence to the target function. The
numerics are floored at machine precision from around degree 30. The data from
the algorithm’s simulation is indeed upper bounded by this maximum error
demonstrating that the accuracy tracks the Chebyshev approximation error. The
dashed line illustrates a supergeometric bound
$\mathcal{O}\left(\mathrm{e}^{-d\log d}\right)$. This indicates that the
accuracy $\varepsilon$ is bounded by $\mathcal{O}\left(\mathrm{e}^{-d\log
d}\right)$. Hence, we can choose degree
$d=\mathcal{O}\left(\mathrm{e}^{W_{0}[\log(1/\varepsilon)]}\right)=\mathcal{O}\left[\log(1/\varepsilon)\right]$
for a fixed target error $0<\varepsilon<1$ with $W_{0}$ the principal branch
of the Lambert W function. As a rule-of-thumb we expect at least geometric
convergence to hold for analytic functions.

Figure 7(c) shows the number of two-qubit gates for varying degrees and
different grid sizes. Circuits are compiled with pytket for the default native
gate set of the H2-1 quantum computer at optimization level $2$ [59]. The
default native two-qubit gate is
$ZZ\text{Phase}(\theta)=\mathrm{e}^{-\mathrm{i}\frac{\theta}{2}Z\otimes Z}$
with arbitrary angle $\theta$. We perform a least-square fit of the data to
the scaling of two-qubit gates derived in Sec. 2.6,
$\alpha_{0}+\alpha_{1}D2^{\lceil\log_{2}(d+1)\rceil}n\lceil\log_{2}n\rceil+\alpha_{2}2^{D\lceil\log_{2}(d+1)\rceil}$,
where $\alpha_{0},\alpha_{1},\alpha_{2}$ are the fitted parameters and we
added a constant offset. The fit in Fig. 7(c) shows good agreement. For large
$n$ and $d$ this scaling law is upper bounded by the simplified expression in
Table 1. The best parameters from the fit (see caption of Fig. 7) are in line
with expectations from a worst-case analysis. For example, the term
$\alpha_{2}2^{D\lceil\log_{2}(d+1)\rceil}$ comes from the implementation of
the operators $A$, $C$, $A^{\dagger}$ in Fig. 5 and the value $\alpha_{2}=3$
indicates that each of the operators is implemented with the expected worst-
case cost $2^{D\lceil\log_{2}(d+1)\rceil}$. We note some residual variance in
the coefficients $\alpha_{0}$ and $\alpha_{1}$ depending on $n$, which we
attribute to suppressed, non-dominant terms and compiler optimizations. The
total number of qubits can be computed from Table 1. For example, the largest
circuit in Fig. 7(c) at $d=63$, $n=6$ requires 30 qubits (12 for the main
register and 18 for ancilla registers). The success probability shown in Fig.
7(d) saturates at around $8\%$, which is close to the asymptotic success
probability for $n\rightarrow\infty$, $d\rightarrow\infty$. We compute the
asymptotic success probability as
$p_{\text{success}}^{*}=\int_{[-1,1]^{2}}\mathrm{d}x\mathrm{d}y|f(x,y)|^{2}/4\big{(}\sum_{k,l=0}^{30}|c_{k,l}|\big{)}^{2}$,
which is the continuous version of Eq. (32). The sum in the numerator runs to
$k,l=30$ because higher-order coefficients do not contribute noticeably and
fall off super-exponentially. The factor $1/4$ is from the size of the
integration domain.

###  3.2 Bivariate Student’s t-distribution via Fourier series

To illustrate the Fourier approach, we consider the bivariate non-standardized
Student’s t-distribution with one degree of freedom ($\nu=1$, this is also
known as the bivariate Cauchy distribution)

$\displaystyle
f(\bm{x})=\frac{1}{2\pi\sqrt{\det(\Sigma)}}\left[1+(\bm{x}-\bm{\mu})^{T}\Sigma^{-1}(\bm{x}-\bm{\mu})\right]^{-\frac{3}{2}},$
(35)

where $\bm{x}=(x,y)^{T}$ and $\bm{\mu}$ and $\Sigma$ are location and scale
parameters, respectively. Student’s t-distribution is characterized by heavy
tails with applications in finance to model heavy-tailed return distributions
[60] among others.

Figure 8: Resources for preparing an approximation to the bivariate Student’s
t-distribution with a 2D Fourier series approximation of degree $d$ in each
dimension. (a) Discretized target function and the difference to the rescaled
amplitudes from a noiseless simulation of the algorithm for $d=7$. (b) Maximum
absolute error of the prepared state evaluated at the grid points for
different grid sizes. The difference between the prepared state and the
Fourier approximation is at machine precision (not shown). Hence, the maximum
error of the noiseless simulation against the target function is upper bounded
by the maximum error of the degree-$d$ Fourier interpolant over the (non-
discretized) domain (grey line – floored at machine precision). The dashed
line illustrates the asymptotic error bound $\mathcal{O}(1/d)$. (c) Two-qubit
gate count of the algorithm for different grid sizes compiled with pytket
optimization level 2 for the H2-1 gate set. The line markers are a least-
square fit of the data (round markers) to the worst-case scaling with
$\alpha_{0}=-4.91,-4.91,-4.91$, $\alpha_{1}=0.78,0.82,0.85$,
$\alpha_{2}=2.99,2.99,2.99$ for $n=4,5,6$ and $D=2$. (d) Empirical (markers)
and asymptotic success probability for $n\rightarrow\infty$ and
$d\rightarrow\infty$ (line).

For the Fourier approach we focus on the domain $\bm{x}\in[0,1]^{2}$ and
choose $\bm{\mu}=(0.5,0.5)^{T}$ and, $\Sigma=\begin{pmatrix}0.05&0\\\
0&0.05\end{pmatrix}$. Note that the bivariate Student’s t-distribution does
not factorize into two univariate Student’s t-distributions. First, we
periodically extend the function from $[0,1]^{2}$ to $[-1,1]^{2}$, interpolate
it with the Fourier series in Eq. (3) and compute coefficients with the fast
Fourier transform following Appendix B. We choose the same degree and same
grid size in each variable, i.e. $d_{x}=d_{y}=d$ and $n_{x}=n_{y}=n$. Figure 8
shows the results of approximating the periodically extended Student’s
t-distribution by a Fourier series of varying degree $d$ and for different
grid sizes $2^{n}\times 2^{n}$. The maximum error in Fig. 8(b) scales
asymptotically as $\mathcal{O}(1/d)$. This is caused by the periodic extension
introducing a discontinuity in the first derivative of the extended function
at the boundaries $x=0,1$ and $y=0,1$. In the univariate case a rule-of-thumb
is that for $k-1$ times continuously differentiable functions and a $k$-th
derivative with bounded variation, the asymptotic interpolation error is
$\mathcal{O}(1/d^{k})$ [49, see Theorem 7.2 for an equivalent statement for
Chebyshev interpolants]. Our numerics indicate this also holds for the
bivariate case of the considered function. To reach higher convergence rate we
could e.g. adjust the periodic extension method to smoothen the sharp corners
at the boundaries.

Figure 8(c) shows the number of two-qubit gates after compilation to the
Quantinuum native gate set with optimization level 2. A least-square fit of
the data to the scaling of two-qubits gates derived in Sec. 2.6,
$\alpha_{0}+\alpha_{1}Dn\lceil\log_{2}(2d+1)\rceil+\alpha_{2}2^{D\lceil\log_{2}(2d+1)\rceil}$,
shows excellent agreement. The small coefficients $\alpha_{1}<1$ (see caption
of Fig. 8) indicate savings from the compilation step. As in the Chebyshev
case, $\alpha_{2}\approx 3$ indicates worst-case synthesis of the operators
$A$, $C$, $A^{\dagger}$ in the LCU on the coefficient registers. The largest
circuit in Fig. 8(c) at $d=63$, $n=6$ requires 26 qubits – 12 for the main
register and 14 for the ancilla register, cf. Table 1. Similar to the
Chebyshev case the success probability in Fig. 8(d) saturates close to its
asymptotic value of $p_{\text{success}}^{*}=7.52\%$. This is computed as
$p_{\text{success}}^{*}=\int_{[0,1]^{2}}\mathrm{d}x\mathrm{d}y|f(x,y)|^{2}/\big{(}\sum_{k,l=-63}^{63}|c_{k,l}|\big{)}^{2}$.
The sum runs to degree $63$ because the coefficients beyond this degree did
not contribute noticeably (tested up to $k,l=128$) and fall off with
increasing degree.

###  3.3 Single electron in periodic 3D Coulomb potential

Our protocols are applicable even when the exact target function is unknown,
but we have access to the coefficients of an approximation by other means,
e.g. through the solution of a differential equation. To demonstrate this, in
this section we use our Fourier approach for the preparation of wavefunctions
of a single electron in a 3D translationally invariant Coulomb potential.

The prepared quantum state can serve as an initial state for the quantum
simulation of materials in a unit cell of atoms with nuclear electronic
Coulomb interactions. The unit cell wavefunction with cell volume equal to 1,
is expressed as a linear combination of $N^{3}$ plane wave basis functions via

$\psi(\bm{r})=\sum_{\bm{k}\in
G}c_{\bm{k}}\mathrm{e}^{\mathrm{i}\pi\bm{k}^{T}\bm{r}}.$ (36)

Here, $\bm{k}=(k_{x},k_{y},k_{z})^{T}$ runs over
$G=\left[-\frac{N}{2},\frac{N-1}{2}\right]^{3}\subset\mathbb{Z}^{3}$,
$\bm{r}=(x,y,z)^{T}\in[0,1]^{3}$ is the real-space position of the electron in
the unit cell. Equation (36) is a trigonometric polynomial in 3 variables with
maximal degree $N/2$ in each dimension, to which we apply our Fourier
approach.

In order to obtain the coefficients $c_{\bm{k}}$, we insert the expansion (36)
into the time-independent Schrödinger equation and solve the resulting set of
linear equations. In the plane wave basis the kinetic energy matrix is
diagonal with entries given by

$T_{\bm{k},\bm{k}}=\frac{(\hbar\pi\|\bm{k}\|)^{2}}{2m}\outerproduct{\bm{k}}{\bm{k}}.$
(37)

The electron-nuclear attraction matrix elements have the following convenient
analytical form in the basis of plane waves

$U_{\bm{k},\bm{k}^{\prime}}=\frac{4}{\pi}\sum_{\ell=1}^{L}w_{\ell}\frac{\mathrm{e}^{\mathrm{i}\pi(\bm{k}^{\prime}-\bm{k})^{T}\mathbf{R}_{\ell}}}{\|\bm{k}-\bm{k}^{\prime}\|^{2}}\outerproduct{\bm{k}}{\bm{k}^{\prime}}.$
(38)

$\mathbf{R}_{\ell}$ is the position of nucleus $\ell$ in real space and
$w_{\ell}$ is the corresponding nuclear weight and $L$ is the number of
nuclei. The $\pi$ factors arise from the definition of the state in Eq. (36).

The electronic Hamiltonian matrix elements are then given by
$H_{\bm{k},\bm{k}^{\prime}}^{\mathrm{(el)}}=-T_{\bm{k},\bm{k}}\delta_{\bm{k},\bm{k}^{\prime}}-U_{\bm{k},\bm{k}^{\prime}}$
and the time-independent Schrödinger equation becomes
$H^{\mathrm{(el)}}\bm{c}=E\bm{c}$ with $\bm{c}$ the vector of coefficients in
Eq. (36). Solutions $\bm{c}$ for the first and second lowest eigenenergies
result in the upper bounds to the ground and excited state wavefunction of the
non-interacting single electronic system. We use those coefficients in our
multivariate Fourier approach to prepare the quantum states for the ground and
first excited state wavefunctions.

We choose the simplest 3D periodic unit cell for demonstration, which is a
single Coulomb potential of nuclear weight and charge $1$ at the center of a
$(1\times 1\times 1)$ unit cell at position $(0.5,0.5,0.5)$. Other more
advanced unit cells can be used, such as face-centered cubic (FCC). In the
left panels of Fig. 9 we present the real part of the ground (top) and first
excited state wavefunctions (bottom) computed directly with Eq. (36) and
coefficients from the solution of the time-independent Schrödinger equation
above. We then use the same set of coefficients to construct the circuits of
our Fourier-based protocol for preparing those wavefunctions. In the right
panels of Fig. 9 we present the real part of the resulting wavefunctions
obtained from a statevector simulation of those circuits. The simulated
circuits replicate to numerical precision the results from the direct
numerical computation.

Figure 9: Real parts of the ground state (top panels) and first excited state (bottom) electronic wavefunctions for a unit cell with a single nucleus at the center at $(0.5,0.5,0.5)$. Left panels: direct evaluation of Eq. (36) with a total of $N^{3}=512$ coefficients from the numerical solution of the time-independent Schrödinger equation. Right panels: equivalent solutions from the Fourier approach of our method evaluated with a statevector simulator using $n=7$ qubits per spatial dimension and $a=3$ ancilla qubits per dimension for a total of 30 qubits. $n$ | $a$ | $D$ | Grid points ($2^{Dn}$) | Coefficients ($2^{Da}$) | Total qubits ($Dn+Da)$ | Two-qubit gates
---|---|---|---|---|---|---
4 | 3 | 3 | $4096$ | $512$ | 21 | $1602$
5 | 3 | 3 | $32\,768$ | $512$ | 24 | $1620$
6 | 3 | 3 | $262\,144$ | $512$ | 27 | $1638$
7 | 3 | 3 | $2\,097\,152$ | $512$ | 30 | $1656$
7 | 4 | 3 | $2\,097\,152$ | $4096$ | 33 | $12\,450$
Table 2: Circuit compilation results for the solutions of the Schrödinger
equation in a 3D periodic Coulomb potential with varying numbers of grid
points and Fourier coefficients. With $n$ the number of qubits per spatial
dimension and $a=\lceil\log_{2}N\rceil$ the number of ancillas per dimension
the total number of qubits is $3n+3a$. The number of two-qubit gates is
determined by compiling the circuits to the gate set given in the main text.
The row marked with * corresponds to the setting of Fig. 9.

Table 2 shows the resource requirements from compilation of the circuits of
our method for preparing the wavefunctions with different numbers of grid
points and coefficients. pytket [59] is used to construct the circuits, using
the gateset {C$X$, $R_{Z}$, $R_{X}$}. This confirms that the dominant cost of
the algorithm comes from the addition of the ancilla qubits that contain the
Fourier coefficients, which is expected from Table 1. Increasing the spatial
resolution (given by $2^{3n}$, where $n$ is the number of qubits per spatial
dimension), leads to little increase in gate count due to the efficient
control structure presented in Fig 2. Using a total of 9 and 12 ancilla qubits
(i.e. registers of size $a=3,4$ per dimension) allows one to encode $512$ and
$4096$ Fourier coefficients, respectively. Thus, the large increase in the
number of gates is due to the 9 vs. 12 qubit circuits for $A$, $C$ and
$A^{\dagger}$ (cf. the term $\mathcal{O}(d^{D})$ for two-qubit gates in Table
1). For the 30 qubit case marked with *, the success probability of generating
the ground ($\Psi_{0}$) and first excited ($\Psi_{1}$) states are $0.7047$ and
$0.2392$, respectively.

Although these are only single electron wavefunctions, they can be used as a
starting point for many-electron wavefunctions by anti-symmetrizing the
spatial registers of multiple single electron states following the procedure
proposed in Ref. [61], or as the initial state for real or imaginary time
evolution under an interacting electron Hamiltonian in the first quantized
picture.

Figure 10:  Experimental results for preparing a 2D Fourier series approximation of the bivariate Gaussian distribution. We show two setting, an uncorrelated Gaussian (upper panel) and a correlated Gaussian (lower panel); parameters are given in the main text. (a, d) Histograms of the raw two-dimensional measurement counts with 50 bins along each axis (upper plots), and the resulting estimates $\hat{\mathfrak{f}}_{d}(x,y)$ (lower plots, obtained from kernel density estimates) of the target functions $\mathfrak{f}_{d}(x,y)$. (b, e) Marginal probabilities $w_{x}(x),w_{y}(y)$ of the measurement data (thin solid lines) and the corresponding kernel density estimates $k_{x}(x),k_{y}(y)$ (bold solid lines). (c, f) Direct comparison between marginal estimates $\hat{\mathfrak{f}}_{d,x}(x),\hat{\mathfrak{f}}_{d,y}(y)$ (solid lines) and the target functions’ marginals $\mathfrak{f}_{d,x}(x),\mathfrak{f}_{d,y}(y)$ (dashed lines).  |  | $\mu_{d,x},\hat{\mu}_{d,x}$ | $\mu_{d,y},\hat{\mu}_{d,y}$ | $\sigma_{d,x}^{2},\hat{\sigma}_{d,x}^{2}$ | $\sigma_{d,y}^{2},\hat{\sigma}_{d,y}^{2}$ | $\rho_{d},\hat{\rho}_{d}$ | $F_{x}$ | $F_{y}$ | $F$
---|---|---|---|---|---|---|---|---|---
Uncorrelated | $\mathfrak{f}_{d}(x,y)$ | $0.5$ | $0.5$ | $0.0413$ | $0.0299$ | $0.0$ | $0.998$ | $0.984$ | $0.986$
Gaussian | $\hat{\mathfrak{f}}_{d}(x,y)$ | $0.494$ | $0.495$ | $0.0468$ | $0.0396$ | $-0.015$
Correlated | $\mathfrak{f}_{d}(x,y)$ | $0.5$ | $0.5$ | $0.0415$ | $0.0321$ | $0.353$ | $0.991$ | $0.963$ | $0.918$
Gaussian | $\hat{\mathfrak{f}}_{d}(x,y)$ | $0.493$ | $0.488$ | $0.0519$ | $0.0510$ | $0.088$
Table 3:  Experimental results for preparing a 2D Fourier series approximation
of the bivariate Gaussian distribution, see Fig. 10. We consider two settings,
the preparation of an uncorrelated and a correlated Gaussian. The table
collects, for both settings, mean values, variances and correlation
coefficients [cf. Eq. (40)] of the target densities $\mathfrak{f}_{d}(x,y)$
and their estimates $\hat{\mathfrak{f}}_{d}(x,y)$, based on kernel density
estimation. Moreover, we show the fidelities $F_{x},F_{y}$ [cf. Eq. (41)]
between marginal estimates,
$\hat{\mathfrak{f}}_{d,x},\hat{\mathfrak{f}}_{d,y}$, and target functions’
marginals, $\mathfrak{f}_{d,x},\mathfrak{f}_{d,y}$, and the full fidelity $F$
between $\hat{\mathfrak{f}}_{d}$ and $\mathfrak{f}_{d}$ [cf. Eq. (41)].

###  3.4 Bivariate Gaussian on quantum hardware

This section demonstrates our state preparation protocol on the H2-1 trapped-
ion quantum computer and illustrates the effect of noise. We focus on
analyzing the output signal without any error mitigation. Error mitigation
protocols, if needed, would likely be performed in combination with a
downstream quantum algorithm, not after state preparation, which is only the
first step of the computation.

We choose a two-dimensional Gaussian distribution as a target function. The
low overheads of our method allow us to prepare the Gaussian on a fine grid
with 24 qubits using up to $237$ two-qubit gates with meaningful output
signal. The Gaussian distribution is given by

$\displaystyle
f(\bm{x})=\frac{1}{2\pi\sqrt{\det(\Sigma)}}\exp[-\tfrac{1}{2}(\bm{x}-\bm{\mu})^{T}\Sigma^{-1}(\bm{x}-\bm{\mu})],$
(39)

where $\bm{x}=(x,y)^{T}$, $\bm{\mu}$ is the mean and $\Sigma$ is the positive
semi-definite covariance matrix. Let $Z=X\times Y$ be a 2D random variable
distributed according to $f(x,y)$. The covariance matrix can be parametrized
by the marginal variances $\sigma_{x}^{2}=\mathrm{Var}(X)$ and
$\sigma_{y}^{2}=\mathrm{Var}(Y)$ (diagonal entries of $\Sigma$), and the
(Pearson) correlation coefficient
$\rho=\mathrm{Cov}(X,Y)/\sigma_{x}\sigma_{y}$ (setting off-diagonal entries
$\rho\sigma_{x}\sigma_{y}$). For this demonstration, we use
$\bm{\mu}=(\mu_{x},\mu_{y})^{T}=(0.5,0.5)^{T}$,
$\sigma_{x}=0.22,\sigma_{y}=0.18$, and the two choices $\rho=0.0$ and
$\rho=0.4$. The first setting defines an uncorrelated Gaussian, i.e. the
bivariate density $f(x,y)$ is the product of two univariate densities. The
second setting defines a correlated Gaussian, which does not factorize.

For both settings we use $n_{x}=n_{y}=9$ qubits for spatial discretization of
the domain $[0,1]^{2}$ (cf. Fig. 5) into $$512$\times$512$=$262\,144$$ grid
points (for comparison, Ref. [44] prepared a bivariate function on a grid of
size $32\times 32$). We use a Fourier series approximation $f_{d}(x,y)$ of
degrees $d=d_{x}=d_{y}=3$, which leads to a total of $(2d+1)^{2}=49$
coefficients. Note that we use a different method to compute the Fourier
coefficients than the one used in Sec. 3.2 because the approximation at the
chosen degree is better for the specific target, see Appendix C.2. To encode
the Fourier coefficients we require ancilla registers with $a_{x}=a_{y}=3$
qubits each, cf. Fig. 5. Overall, this results in a circuit with $24$ qubits.
We test our approach on the H2-1 trapped-ion quantum computer [62, 63]. At the
time of this demonstration in April and May 2024, it provides $32$ fully-
connected qubits with average single- and two-qubit gate infidelities of
$3\text{\times}{10}^{-5}$ and $2\text{\times}{10}^{-3}$, respectively, and
average state preparation and measurement error of $2\text{\times}{10}^{-3}$.
Detailed specifications are given in Appendix C.1.

The Fourier approximation $f_{d=3}(x,y)$ of the (positive) Gaussian density in
Eq. (39) can have noticeable negative values in $[0,1]^{2}$. This is visible
in the small kinks of the red dashed line (showing a marginal over
$|f_{d}(x,y)|$) in Fig. 10(c) close to the boundaries of the domain $[0,1]$.
To avoid expensive state-tomographic methods, we only read out the magnitudes
$|f_{d}(x,y)|$ of the Fourier series via measurements in the computational
basis and compare our experimental data to the positive target function
$\mathfrak{f}_{d}(x,y)=|f_{d}(x,y)|$. We run each experiment with $20\,000$
shots. The success probability in our setting is around $10\%$ (for details
see Appendix C.2). Let $w(x_{i},y_{j})$ be the probability to measure grid
point $(x_{i},y_{j})\in[0,1]^{2},i,j=1,\dots,512$, after post-selection on the
zero state in the ancilla register (cf. Fig. 4), and let
$w_{x}(x_{i})=\sum_{j=1}^{512}w(x_{i},y_{j})$,
$w_{y}(y_{j})=\sum_{i=1}^{512}w(x_{i},y_{j})$ be its marginals. The upper
panels of Fig. 10(a) and (d) show histograms of the joint measurement
distribution $w(x,y)$ for the uncorrelated and correlated settings,
respectively. Figures 10(b) and (e) show the corresponding $x$\- and
$y$-marginals as thin lines. Notably, the measured signals suffer from
substantial fluctuations due to shot noise.

To compare the empirical probabilities to the target $\mathfrak{f}_{d}(x,y)$
in a meaningful way, we smoothen the experimental data with kernel density
estimation (KDE) [64, 65, 66]. The KDE yields a continuous approximation of
the density underlying a finite data sample. This is obtained by replacing the
raw counts with, in our case, Gaussian kernels of a given bandwidth $h$ (see
Appendix C.3). The KDE becomes smoother with increasing bandwidth $h$. On the
other hand, for too large $h$ the KDE loses information about the underlying
data. In the limit $h\rightarrow\infty$, the KDE approaches a perfectly flat
distribution. The “optimal” bandwidth can be computed via minimization of an
appropriate cost function. A naive choice would be, for example, to consider a
distance measure or statistical divergence (such as the relative entropy)
between KDE and the target function $\mathfrak{f}_{d}(x,y)$. However, this
choice introduces an obvious bias. Instead we choose to optimize $h$ via cross
validation from the measurement data alone (see Appendix C).

Let $k(x,y)$, $k_{x}(x)$ and $k_{y}(y)$ be the KDE [cf. Eq. (57)]
corresponding to the empirical probabilities $w(x_{i},y_{j})$, $w_{x}(x_{i})$
and $w_{y}(y_{j})$ ($i,j=1,\dots,512$), rescaled such that
$\sum_{i,j=1}^{512}k(x_{i},y_{j})=\sum_{i=1}^{512}k_{x}(x_{i})=\sum_{j=1}^{512}k_{y}(y_{j})=1$.
The thick solid lines in Figs. 10(b) and (e) are the marginal KDE
$k_{x}(x_{i})$ and $k_{y}(y_{j})$. To compare to $\mathfrak{f}_{d}(x,y)$, we
need to take the square root, which defines our experimental estimates
$\hat{\mathfrak{f}}_{d}(x,y)=\sqrt{k(x,y)}$,
$\hat{\mathfrak{f}}_{d,x}(x)=\sqrt{k_{x}(x)}$ and
$\hat{\mathfrak{f}}_{d,y}(y)=\sqrt{k_{y}(y)}$. The lower panels of Fig. 10(a)
and (d) display $\hat{\mathfrak{f}}_{d}(x,y)$. In Fig. 10(c) and (f), we
compare the estimated marginals (solid lines) to the target marginals
$\mathfrak{f}_{d,x},\mathfrak{f}_{d,y}$ of $\mathfrak{f}_{d}$ (dashed lines).
For direct comparison, we rescale all functions
$\hat{\mathfrak{f}}_{d},\hat{\mathfrak{f}}_{d,x},\hat{\mathfrak{f}}_{d,y},\mathfrak{f}_{d,x},\mathfrak{f}_{d,y}$
such that they integrate to one on their corresponding domains, i.e. that they
represent proper probability densities. Visually, we observe a good agreement.
We clearly recognize an elongated shape in the lower plot of Fig. 10(d), as
expected from a correlated 2D Gaussian.

For a quantitative comparison we determine the mean values and the covariance
matrix of the target function $\mathfrak{f}_{d}(x,y)$ and its estimate
$\hat{\mathfrak{f}}_{d}(x,y)$ (obtained from the experimental data):

$\begin{split}\mu_{d,x}&=\int_{[0,1]^{2}}\mathrm{d}x\mathrm{d}y\,x\mathfrak{f}_{d}(x,y),\\\
\mu_{d,y}&=\int_{[0,1]^{2}}\mathrm{d}x\mathrm{d}y\,y\mathfrak{f}_{d}(x,y),\\\
\sigma_{d,x}^{2}&=\int_{[0,1]^{2}}\mathrm{d}x\mathrm{d}y\,(x-\mu_{d,x})^{2}\mathfrak{f}_{d}(x,y),\\\
\sigma_{d,y}^{2}&=\int_{[0,1]^{2}}\mathrm{d}x\mathrm{d}y\,(y-\mu_{d,y})^{2}\mathfrak{f}_{d}(x,y),\\\
\rho_{d}&=\frac{\int_{[0,1]^{2}}\mathrm{d}x\mathrm{d}y\,(x-\mu_{d,x})(y-\mu_{d,y})\mathfrak{f}_{d}(x,y)}{\sigma_{d,x}\sigma_{d,y}},\end{split}$
(40)

and
$\hat{\mu}_{d,x},\hat{\mu}_{d,y},\hat{\sigma}_{d,x}^{2},\hat{\sigma}_{d,x}^{2},\hat{\rho}_{d}$
computed equivalently by replacing all quantities with their “hat” versions
(e.g. replace $\mathfrak{f}_{d}(x,y)$ with $\hat{\mathfrak{f}}_{d}(x,y)$).
Furthermore, we compute the classical fidelity between the target and
measurements

$F=\frac{\left(\sum_{i,j=1}^{512}\sqrt{\mathfrak{f}_{d}(x_{i},y_{j})\hat{\mathfrak{f}}_{d}(x_{i},y_{j})}\right)^{2}}{\sum_{i,j=1}^{512}\mathfrak{f}_{d}(x_{i},y_{j})\sum_{i,j=1}^{512}\hat{\mathfrak{f}}_{d}(x_{i},y_{j})},$
(41)

where $x_{i},y_{i}$ runs over the grid points. Table 3 compares the results
for both settings. Overall, we observe a good agreement between the targets
and the experimentally prepared functions, both in terms of statistical
quantifiers and fidelities. For the uncorrelated 2D Gaussian we obtain an
overall fidelity of $F=0.986$, for the correlated one we obtain $F=0.918$.

We observe a larger deviation between target and prepared function for the
correlated Gaussian. For example, the observed variances
$\hat{\sigma}_{d,x}^{2},\hat{\sigma}_{d,x}^{2}$ overestimate the target
variances $\sigma_{d,x}^{2},\sigma_{d,x}^{2}$ in this setting. Although we
clearly identify a non-zero correlation coefficient in the setting of the
correlated Gaussian, the estimated $\hat{\rho}_{d}=0.088$ is significantly
smaller, roughly by a factor of four, than the target value $\rho_{d}=0.353$.
The pronounced discrepancy between $\mathfrak{f}_{d}$ and
$\hat{\mathfrak{f}}_{d}$ in the correlated setting, in comparison to the
uncorrelated Gaussian, is expected. This is because, in the latter scenario,
we are able to exploit the factorization of $\mathfrak{f}_{d}(x,y)$ into two
univariate functions also at the level of circuit synthesis. Indeed we obtain
much shallower circuits for the uncorrelated setting, which are much less
prone to noise, as discussed in more detail in App. C.2.

##  4 Conclusions

This work introduces methods for the preparation of quantum states that encode
multivariate functions, and analyzes their resource requirements. As a proof-
of-concept, we demonstrate the successful preparation of representative
initial states (with potential applications e.g. in quantum chemistry, physics
and finance simulations) on comparatively high-resolution grids. We execute a
24-qubit experiment on Quantinuum’s H2-1 trapped-ion quantum computer and
analyse the impact of noise on the prepared states.

Many functions of interest cannot be naturally encoded in a quantum state
(e.g. $1/x$ due to a singularity), and function approximation becomes a
necessary step to obtain practical protocols. Our methods are based on
multivariate versions of Fourier and Chebyshev series approximations. We note
that power series can be mapped to Chebyshev series [67, 68] where our method
applies as-is. Walsh series used e.g. in [46, 47] can be implemented as a
special case of our Fourier method. It is also possible to use different basis
functions for each dimension in a hybrid Fourier-Chebyshev method. Moreover,
the building blocks of our method can be applied recursively to create more
expressive basis functions, e.g. wavelets and kernels, which are then linearly
combined to form more complex functions. Ultimately, the approximation method
and error budget will depend on the use case (e.g. whether the target function
is periodic) and specifics of the available quantum hardware (e.g. two-qubit
gate fidelity).

QSP [35, 36, 37, 38, 39] has attracted significant interest as a framework for
developing quantum algorithms via polynomial transformations of block-encoded
matrices. Multivariate versions of QSP have been proposed and analyzed,
although questions about expressive power remain open [41, 42, 43]. Our work
shows that linear combination of unitaries [50], with its simplicity and
effectiveness, is a valid alternative for the task of preparing quantum states
that encode multivariate functions in their amplitudes. An interesting
research direction consists of using (univariate) QSP to further process the
matrices produced by our method. For example, suppose we wanted to prepare a
state that encodes a hypersurface – a set of dimension $(D-1)$ defined by a
polynomial equation of the form $f(x_{1},\dots,x_{D})=0$. After discretizing
each dimension via $H^{\otimes D}$ and using our method to construct the
diagonal operator $f(H^{\otimes D})$, we use QSP to filter its eigenvalues. In
particular, we want a polynomial $p$ such that $p(0)=1$ and $|p(y)|$ is small
for $y\in[-1,-\epsilon]\cup[\epsilon,1]$. An optimal filter is given by Lin
and Tong in [69]. Applying the resulting operator to the uniform superposition
we obtain a state proportional to
$\sum_{x_{1}}\cdots\sum_{x_{D}}p(f(x_{1},\dots,x_{D}))\ket{x_{1},\dots,x_{D}}$,
approximately encoding the hypersurface. Remarkably, this proposal uses
univariate QSP for a fundamentally multivariate problem.

## Acknowledgements

We thank Alexandre Krajenbrink and Enrico Rinaldi for helpful feedback on the
manuscript.

## References

  * [1] Sarah K. Leyton and Tobias J. Osborne.  “A quantum algorithm to solve nonlinear differential equations” (2008).  arXiv:0812.4423.
  * [2] Andrew M. Childs, Jin-Peng Liu, and Aaron Ostrander.  “High-precision quantum algorithms for partial differential equations”.  Quantum 5, 574 (2021).
  * [3] Dylan Herman, Cody Googin, Xiaoyuan Liu, Alexey Galda, Ilya Safro, Yue Sun, Marco Pistoia, and Yuri Alexeev.  “A survey of quantum computing for finance” (2022).  arXiv:2201.02773.
  * [4] Nikitas Stamatopoulos and William J. Zeng.  “Derivative Pricing using Quantum Signal Processing”.  Quantum 8, 1322 (2024).
  * [5] Ismail Yunus Akhalwaya, Adam Connolly, Roland Guichard, Steven Herbert, Cahit Kargi, Alexandre Krajenbrink, Michael Lubasch, Conor Mc Keever, Julien Sorci, Michael Spranger, and Ifan Williams.  “A Modular Engine for Quantum Monte Carlo Integration” (2023).  arXiv:2308.06081.
  * [6] Alán Aspuru-Guzik, Anthony D. Dutoi, Peter J. Love, and Martin Head-Gordon.  “Simulated quantum computation of molecular energies”.  Science 309, 1704–1707 (2005).
  * [7] Ivan Kassal, Stephen P. Jordan, Peter J. Love, Masoud Mohseni, and Alán Aspuru-Guzik.  “Polynomial-time quantum algorithm for the simulation of chemical dynamics”.  Proceedings of the National Academy of Sciences 105, 18681–18686 (2008).
  * [8] Yudong Cao, Jonathan Romero, Jonathan P Olson, Matthias Degroote, Peter D Johnson, Mária Kieferová, Ian D Kivlichan, Tim Menke, Borja Peropadre, Nicolas P D Sawaya, Sukin Sim, Libor Veis, and Alán Aspuru-Guzik.  “Quantum chemistry in the age of quantum computing”.  Chem. Rev. 119, 10856–10915 (2019).
  * [9] Sam McArdle, Suguru Endo, Alán Aspuru-Guzik, Simon C Benjamin, and Xiao Yuan.  “Quantum computational chemistry”.  Rev. Mod. Phys. 92, 015003 (2020).
  * [10] Jie Liu, Yi Fan, Zhenyu Li, and Jinlong Yang.  “Quantum algorithms for electronic structures: basis sets and boundary conditions”.  Chem. Soc. Rev. 51, 3263–3279 (2022).
  * [11] Stephen P. Jordan, Keith S. M. Lee, and John Preskill.  “Quantum algorithms for quantum field theories”.  Science 336, 1130–1133 (2012).
  * [12] Mari Carmen Bañuls, Rainer Blatt, Jacopo Catani, Alessio Celi, Juan Ignacio Cirac, Marcello Dalmonte, Leonardo Fallani, Karl Jansen, Maciej Lewenstein, Simone Montangero, Christine A. Muschik, Benni Reznik, Enrique Rico, Luca Tagliacozzo, Karel Van Acoleyen, Frank Verstraete, Uwe-Jens Wiese, Matthew Wingate, Jakub Zakrzewski, and Peter Zoller.  “Simulating lattice gauge theories within quantum technologies”.  The European Physical Journal D 74, 165 (2020).
  * [13] Christian W. Bauer, Zohreh Davoudi, A. Baha Balantekin, Tanmoy Bhattacharya, Marcela Carena, Wibe A. de Jong, Patrick Draper, Aida El-Khadra, Nate Gemelke, Masanori Hanada, Dmitri Kharzeev, Henry Lamm, Ying-Ying Li, Junyu Liu, Mikhail Lukin, Yannick Meurice, Christopher Monroe, Benjamin Nachman, Guido Pagano, John Preskill, Enrico Rinaldi, Alessandro Roggero, David I. Santiago, Martin J. Savage, Irfan Siddiqi, George Siopsis, David Van Zanten, Nathan Wiebe, Yukari Yamauchi, Kübra Yeter-Aydeniz, and Silvia Zorzetti.  “Quantum simulation for high-energy physics”.  PRX Quantum 4, 027001 (2023).
  * [14] Alexander M. Dalzell, Sam McArdle, Mario Berta, Przemyslaw Bienias, Chi-Fang Chen, András Gilyén, Connor T. Hann, Michael J. Kastoryano, Emil T. Khabiboulline, Aleksander Kubica, Grant Salton, Samson Wang, and Fernando G. S. L. Brandão.  “Quantum algorithms: A survey of applications and end-to-end complexities” (2023).  arXiv:2310.03011.
  * [15] Juha J. Vartiainen, Mikko Möttönen, and Martti M. Salomaa.  “Efficient decomposition of quantum gates”.  Phys. Rev. Lett. 92, 177902 (2004).
  * [16] Mikko Möttönen, Juha J. Vartiainen, Ville Bergholm, and Martti M. Salomaa.  “Transformation of quantum states using uniformly controlled rotations”.  Quant. Inf. Comp. 5, 467 (2005).  arXiv:quant-ph/0407010.
  * [17] Ville Bergholm, Juha J. Vartiainen, Mikko Möttönen, and Martti M. Salomaa.  “Quantum circuits with uniformly controlled one-qubit gates”.  Phys. Rev. A 71, 052330 (2005).
  * [18] Vivek V. Shende, Stephen S. Bullock, and Igor L. Markov.  “Synthesis of quantum-logic circuits”.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25, 1000–1010 (2006).  arXiv:quant-ph/0406176.
  * [19] Martin Plesch and Časlav Brukner.  “Quantum-state preparation with universal gate decompositions”.  Phys. Rev. A 83, 032302 (2011).
  * [20] Guang Hao Low, Vadym Kliuchnikov, and Luke Schaeffer.  “Trading T-gates for dirty qubits in state preparation and unitary synthesis” (2018).  arXiv:1812.00954.
  * [21] Xiaoming Sun, Guojing Tian, Shuai Yang, Pei Yuan, and Shengyu Zhang.  “Asymptotically optimal circuit depth for quantum state preparation and general unitary synthesis”.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 42, 3301–3314 (2023).
  * [22] Xiao-Ming Zhang, Tongyang Li, and Xiao Yuan.  “Quantum state preparation with optimal circuit depth: Implementations and applications”.  Phys. Rev. Lett. 129, 230504 (2022).
  * [23] Lov Grover and Terry Rudolph.  “Creating superpositions that correspond to efficiently integrable probability distributions” (2002).  arXiv:quant-ph/0208112.
  * [24] Alexei Kitaev and William A. Webb.  “Wavefunction preparation and resampling using a quantum computer” (2009).  arXiv:0801.0342.
  * [25] Arthur G. Rattew, Yue Sun, Pierre Minssen, and Marco Pistoia.  “The Efficient Preparation of Normal Distributions in Quantum Registers”.  Quantum 5, 609 (2021).
  * [26] Gabriel Marin-Sanchez, Javier Gonzalez-Conde, and Mikel Sanz.  “Quantum algorithms for approximate function loading”.  Phys. Rev. Res. 5, 033114 (2023).
  * [27] Adam Holmes and A. Y. Matsuura.  “Efficient quantum circuits for accurate state preparation of smooth, differentiable functions” (2020).  arXiv:2005.04351.
  * [28] Javier Gonzalez-Conde, Thomas W. Watts, Pablo Rodriguez-Grasa, and Mikel Sanz.  “Efficient quantum amplitude encoding of polynomial functions”.  Quantum 8, 1297 (2024).
  * [29] Jason Iaconis, Sonika Johri, and Elton Yechao Zhu.  “Quantum state preparation of normal distributions using matrix product states”.  npj Quantum Information 10, 15 (2024).
  * [30] Arthur G. Rattew and Bálint Koczor.  “Preparing arbitrary continuous functions in quantum registers with logarithmic complexity” (2022).  arXiv:2205.00519.
  * [31] Yuval R. Sanders, Guang Hao Low, Artur Scherer, and Dominic W. Berry.  “Black-box quantum state preparation without arithmetic”.  Phys. Rev. Lett. 122, 020502 (2019).
  * [32] Johannes Bausch.  “Fast Black-Box Quantum State Preparation”.  Quantum 6, 773 (2022).
  * [33] Jessica Lemieux, Matteo Lostaglio, Sam Pallister, William Pol, Karthik Seetharam, Sukin Sim, and Burak Şahinoğlu.  “Quantum sampling algorithms for quantum state preparation and matrix block-encoding” (2024).  arXiv:2405.11436.
  * [34] Shengbin Wang, Zhimin Wang, Guolong Cui, Shangshang Shi, Ruimin Shang, Lixin Fan, Wendong Li, Zhiqiang Wei, and Yongjian Gu.  “Fast black-box quantum state preparation based on linear combination of unitaries”.  Quantum Information Processing 20, 270 (2021).
  * [35] Guang Hao Low and Isaac L. Chuang.  “Optimal Hamiltonian Simulation by Quantum Signal Processing”.  Phys. Rev. Lett. 118, 010501 (2017).
  * [36] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  “Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics”.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing.  ACM (2019).  arXiv:1806.01838.
  * [37] Thais L. Silva, Márcio M. Taddei, Stefano Carrazza, and Leandro Aolita.  “Fragmented imaginary-time evolution for early-stage quantum signal processors”.  Scientific Reports 13, 18258 (2023).
  * [38] John M. Martyn, Zane M. Rossi, Andrew K. Tan, and Isaac L. Chuang.  “Grand unification of quantum algorithms”.  PRX Quantum 2, 040203 (2021).
  * [39] Yuta Kikuchi, Conor Mc Keever, Luuk Coopmans, Michael Lubasch, and Marcello Benedetti.  “Realization of quantum signal processing on a noisy quantum computer”.  npj Quantum Information9 (2023).
  * [40] Sam McArdle, András Gilyén, and Mario Berta.  “Quantum state preparation without coherent arithmetic” (2022).  arXiv:2210.14892.
  * [41] Zane M. Rossi and Isaac L. Chuang.  “Multivariable quantum signal processing (M-QSP): prophecies of the two-headed oracle”.  Quantum 6, 811 (2022).  arXiv:2205.06261.
  * [42] Hitomi Mori, Keisuke Fujii, and Kaoru Mizuta.  “Comment on "Multivariable quantum signal processing (M-QSP): prophecies of the two-headed oracle"” (2023).  arXiv:2310.00918.
  * [43] Balázs Németh, Blanka Kövér, Boglárka Kulcsár, Roland Botond Miklósi, and András Gilyén.  “On variants of multivariate quantum signal processing and their characterizations” (2023).  arXiv:2312.09072.
  * [44] Mudassir Moosa, Thomas W Watts, Yiyou Chen, Abhijat Sarma, and Peter L McMahon.  “Linear-depth quantum circuits for loading fourier approximations of arbitrary functions”.  Quantum Science and Technology 9, 015002 (2023).
  * [45] Juan José García-Ripoll.  “Quantum-inspired algorithms for multivariate analysis: From interpolation to partial differential equations”.  Quantum 5, 431 (2021).
  * [46] Julien Zylberman and Fabrice Debbasch.  “Efficient Quantum State Preparation with Walsh Series” (2023).  arXiv:2307.08384.
  * [47] Jonathan Welch, Daniel Greenbaum, Sarah Mostame, and Alan Aspuru-Guzik.  “Efficient quantum circuits for diagonal unitaries without ancillas”.  New Journal of Physics 16, 033040 (2014).
  * [48] John P. Boyd.  “Chebyshev and Fourier Spectral Methods”.  Dover Publications. Mineola, N.Y. (2001).  2nd revised edition.
  * [49] Lloyd N. Trefethen.  “Approximation Theory and Approximation Practice, Extended Edition”.  Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics. Philadelphia, PA (2019).
  * [50] Andrew M. Childs and Nathan Wiebe.  “Hamiltonian simulation using linear combinations of unitary operations”.  Quantum Info. Comput. 12, 901–924 (2012).  arXiv:1202.5822.
  * [51] Guang Hao Low and Isaac L. Chuang.  “Hamiltonian simulation by qubitization”.  Quantum 3, 163 (2019).
  * [52] Andrew M. Childs, Robin Kothari, and Rolando D. Somma.  “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision”.  SIAM Journal on Computing 46, 1920–1950 (2017).
  * [53] William Kirby, Mario Motta, and Antonio Mezzacapo.  “Exact and efficient Lanczos method on a quantum computer”.  Quantum 7, 1018 (2023).
  * [54] Daan Camps, Lin Lin, Roel Van Beeumen, and Chao Yang.  “Explicit quantum circuits for block encodings of certain sparse matrices”.  SIAM Journal on Matrix Analysis and Applications 45, 801–827 (2024).  arXiv:2203.10236.
  * [55] Xiaoming Sun Junhong Nie, Wei Zi.  “Quantum circuit for multi-qubit toffoli gate with optimal resource” (2024).  arXiv:2402.05053.
  * [56] Yanghua Wang.  “Generalized seismic wavelets”.  Geophysical Journal International 203, 1172–1178 (2015).
  * [57] J. P. Antoine.  “Wavelet analysis of signals and images, a grand tour”.  Ciencias Matemáticas 18, 113–144 (2000).
  * [58] Lewis Wright, Conor Mc Keever, Jeremy T. First, Rory Johnston, Jeremy Tillay, Skylar Chaney, Matthias Rosenkranz, and Michael Lubasch.  “Noisy intermediate-scale quantum simulation of the one-dimensional wave equation” (2024).  arxiv:2402.19247.
  * [59] Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington, and Ross Duncan.  “T$|$ket$\rangle$: A retargetable compiler for NISQ devices”.  Quantum Science and Technology 6, 014003 (2020).  arxiv:2003.10611.
  * [60] Ching-Wai (Jeremy) Chiu, Haroon Mumtaz, and Gábor Pintér.  “Forecasting with VAR models: Fat tails and stochastic volatility”.  International Journal of Forecasting 33, 1124–1143 (2017).
  * [61] Hans Hon Sang Chan, Richard Meister, Tyson Jones, David P. Tew, and Simon C. Benjamin.  “Grid-based methods for chemistry simulations on a quantum computer”.  Science Advances 9, eabo7484 (2023).
  * [62] S. A. Moses, C. H. Baldwin, M. S. Allman, R. Ancona, L. Ascarrunz, C. Barnes, J. Bartolotta, B. Bjork, P. Blanchard, M. Bohn, J. G. Bohnet, N. C. Brown, N. Q. Burdick, W. C. Burton, S. L. Campbell, J. P. Campora, C. Carron, J. Chambers, J. W. Chan, Y. H. Chen, A. Chernoguzov, E. Chertkov, J. Colina, J. P. Curtis, R. Daniel, M. DeCross, D. Deen, C. Delaney, J. M. Dreiling, C. T. Ertsgaard, J. Esposito, B. Estey, M. Fabrikant, C. Figgatt, C. Foltz, M. Foss-Feig, D. Francois, J. P. Gaebler, T. M. Gatterman, C. N. Gilbreth, J. Giles, E. Glynn, A. Hall, A. M. Hankin, A. Hansen, D. Hayes, B. Higashi, I. M. Hoffman, B. Horning, J. J. Hout, R. Jacobs, J. Johansen, L. Jones, J. Karcz, T. Klein, P. Lauria, P. Lee, D. Liefer, S. T. Lu, D. Lucchetti, C. Lytle, A. Malm, M. Matheny, B. Mathewson, K. Mayer, D. B. Miller, M. Mills, B. Neyenhuis, L. Nugent, S. Olson, J. Parks, G. N. Price, Z. Price, M. Pugh, A. Ransford, A. P. Reed, C. Roman, M. Rowe, C. Ryan-Anderson, S. Sanders, J. Sedlacek, P. Shevchuk, P. Siegfried, T. Skripka, B. Spaun, R. T. Sprenkle, R. P. Stutz, M. Swallows, R. I. Tobey, A. Tran, T. Tran, E. Vogt, C. Volin, J. Walker, A. M. Zolot, and J. M. Pino.  “A race-track trapped-ion quantum processor”.  Phys. Rev. X 13, 041052 (2023).
  * [63] “Quantinuum System Model H2 Powered by Honeywell”.  url: https://www.quantinuum.com/hardware/h2.  (accessed: 2024-03-25).
  * [64] Duin.  “On the Choice of Smoothing Parameters for Parzen Estimators of Probability Density Functions”.  IEEE Transactions on Computers C-25, 1175–1179 (1976).
  * [65] Mats Rudemo.  “Empirical choice of histograms and kernel density estimators”.  Scandinavian Journal of Statistics 9, 65–78 (1982).  url: http://www.jstor.org/stable/4615859.
  * [66] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  “The Elements of Statistical Learning”.  Springer Series in Statistics. Springer New York. New York, NY (2009).
  * [67] Henry C. Thacher.  “Conversion of a power to a series of Chebyshev polynomials”.  Commun. ACM 7, 181–182 (1964).
  * [68] Hippolyte Nyengeri, Rénovat Nizigiyimana, Jean-Pierre Mutankana, Henry Bayaga, and Ferdinand Bayubahe.  “Power and Chebyshev Series Transformation Formulas with Applications to Solving Ordinary Differential Equations via the Fröbenius and Taylor’s Methods”.  Open Access Library Journal 8, 1–19 (2021).
  * [69] Lin Lin and Yu Tong.  “Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems”.  Quantum 4, 361 (2020).
  * [70] D. Elliott, D.F. Paget, G.M. Phillips, and P.J. Taylor.  “Error of truncated chebyshev series and other near minimax polynomial approximations”.  Journal of Approximation Theory 50, 49–57 (1987).
  * [71] Alex Townsend and Lloyd N. Trefethen.  “An Extension of Chebfun to Two Dimensions”.  SIAM Journal on Scientific Computing 35, C495–C518 (2013).
  * [72] Lloyd N. Trefethen.  “Multivariate polynomial approximation in the hypercube”.  Proceedings of the American Mathematical Society 145, 4837–4844 (2017).
  * [73] Charles Fefferman.  “On the convergence of multiple Fourier series”.  Bulletin of the American Mathematical Society 77, 744–745 (1971).
  * [74] Charles Fefferman.  “On the divergence of multiple Fourier series”.  Bulletin of the American Mathematical Society 77, 191–195 (1971).
  * [75] Adrian W. Bowman.  “An Alternative Method of Cross-Validation for the Smoothing of Density Estimates”.  Biometrika 71, 353–360 (1984).
  * [76] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  “Scikit-learn: Machine learning in Python”.  Journal of Machine Learning Research 12, 2825–2830 (2011).  url: http://jmlr.org/papers/v12/pedregosa11a.html.

##  Appendix A Classical preprocessing for the Chebyshev method

###  A.1 One dimension

A Lipshitz continuous function $f$ defined on $[-1,1]$ has a unique expansion
as a Chebyshev series,

$f(x)=\sum_{k=0}^{\infty}\hat{c}_{k}T_{k}(x)$ (42)

with coefficients

$\displaystyle\hat{c}_{k}$
$\displaystyle=\frac{2-\delta_{k,0}}{\pi}\int_{-1}^{1}\frac{\mathrm{d}x}{\sqrt{1-x^{2}}}f(x)T_{k}(x),$
(43)

where $\delta_{k,j}$ is the Kronecker delta (see e.g. [49, Theorem 3.1]). The
Chebyshev polynomials $T_{k}(x)$ satisfy the orthogonality condition

$\displaystyle\int_{-1}^{1}\frac{\mathrm{d}x}{\sqrt{1-x^{2}}}T_{k}(x)T_{j}(x)=\begin{cases}0&\text{if
}k\neq j,\\\ \frac{\pi}{2}&\text{if }k=j\neq 0,\\\ \pi&\text{if
}k=j=0\end{cases}$ (44)

as well as the discrete orthogonality condition

$\sum_{m=0}^{N-1}T_{k}\left(x_{m}^{(N)}\right)T_{j}\left(x_{m}^{(N)}\right)=\begin{cases}0&\text{if
}k\neq j,\\\ \frac{N}{2}&\text{if }k=j\neq 0,\\\ N&\text{if
}k=j=0,\end{cases}$ (45)

where $x_{m}^{(N)}=\cos\big(\frac{\pi(2m+1)}{2N}\big{missing})$ for
$m=0,1,\dots,N-1$ are the roots of the Chebyshev polynomial $T_{N}(x)$.

The analysis can be extended to complex functions. The Chebyshev series
remains convergent within an ellipse in the complex plane with foci $\pm 1$
that does not contain any poles, branch points or other singularities of $f$
[48, Theorem 7].

For the numerical results in this paper we construct the degree-$d$ polynomial
approximation Eq. (2) via interpolation. The interpolant $f_{d}$ is uniquely
determined by finding coefficients $c_{k}$ such that $f_{d}(x_{m})=f(x_{m})$
on a set of $d+1$ points $\left\\{x_{m}\right\\}_{m=0}^{d}$. We choose the
roots of the degree-$(d+1)$ Chebyshev polynomial as interpolation points,
$\left\\{x_{m}^{(d+1)}\right\\}_{m=0}^{d}$, as this entails favorable
convergence properties to the target function.

First, let us consider the interpolation error [70]. If the Chebyshev series
of a function $f$ exists, the error is bounded by
$\max_{x\in[-1,1]}|f(x)-f_{d}(x)|\leq 2\sum_{k=d+1}^{\infty}|\hat{c}_{k}|$
[48, Theorem 21]. Typically, the magnitude of coefficients decreases rapidly
in $k$ with a rate depending on the smoothness of $f$. Given a target function
$f\in C^{(d+1)}[-1,1]$ (i.e. $f$ is a complex function whose $(d+1)$-st
derivative exists and is continuous on the interval $[-1,1]$), the degree-$d$
Chebyshev interpolation incurs an error

$\displaystyle\varepsilon_{d}:=\max_{x\in[-1,1]}|f(x)-f_{d}(x)|=\frac{|f^{(d+1)}(\xi)|}{2^{d}(d+1)!}$
(46)

for some $\xi\in(-1,1)$. The same formula holds for truncation instead of
interpolation, possibly for some other $\xi$ [70]. If $M$ is an upper bound
for $f^{(d+1)}$ on $[-1,1]$ (this exists because $f^{(d+1)}$ is continuous on
the bounded interval $[-1,1]$), then the truncation and interpolant error is
upper-bounded by

$\varepsilon_{d}\leq\frac{M}{2^{d}(d+1)!}.$ (47)

The take-away message is that the Chebyshev interpolant rapidly converges to
the target function if it is sufficiently smooth. The function does not need
to be periodic as is the case for the Fourier expansion.

The coefficients $c_{k}$ of the interpolant can be computed with the fast
Fourier transform. The interpolant $f_{d}$ is a Chebyshev series for a
polynomial of degree at most $d$. Its coefficients are given exactly by
applying the $(d+1)$-point trapezoidal rule to the integral in Eq. (43)
evaluated at the $d+1$ Chebyshev roots [48, Theorem 21]

$c_{k}^{\prime}=\frac{2}{d+1}\sum_{m=0}^{d}f_{d}\left(x_{m}^{(d+1)}\right)T_{k}\left(x_{m}^{(d+1)}\right)=\frac{2}{d+1}\sum_{m=0}^{d}f_{d}\left(x_{m}^{(d+1)}\right)\cos\left(k\pi\frac{2m+1}{2(d+1)}\right).$
(48)

This is the discrete cosine transform, which we evaluate via the fast Fourier
transform with complexity $\mathcal{O}(d\log d)$. The notation
$c_{k}^{\prime}$ means we need to divide the $k=0$ term by $2$,
$c_{0}=c_{0}^{\prime}/2$ while $c_{k}=c_{k}^{\prime}$ for $k>0$. The
coefficients of the interpolant and infinite Chebyshev series are related via
a simple formula [48, Theorem 21].

###  A.2 Two and higher dimensions

To compute the 2D interpolant $f_{d_{x},d_{y}}(x,y)$, Eq. (4), we match the
target function at the grid points
$\left\\{\left(x_{m}^{(d_{x}+1)},y_{n}^{(d_{y}+1)}\right)\right\\}_{m=0,n=0}^{d_{x},d_{y}}$
given by the Chebyshev roots of degrees $(d_{x}+1)$ and $(d_{y}+1)$. In
analogy to the 1D case, the coefficients $c_{k,l}$ of this interpolant are
given by the 2D cosine transform on this grid

$c_{k,l}^{\prime}=\frac{4}{(d_{x}+1)(d_{y}+1)}\sum_{m=0}^{d_{x}}\sum_{n=0}^{d_{y}}f(x_{m}^{(d_{x}+1)},y_{n}^{(d_{y}+1)})\cos\left(k\pi\frac{2m+1}{2(d_{x}+1)}\right)\cos\left(l\pi\frac{2n+1}{2(d_{y}+1)}\right).$
(49)

This can be evaluated with the fast Fourier transform with cost
$\mathcal{O}(d_{x}d_{y}\log(d_{x}d_{y}))$. The notation $c_{k,l}^{\prime}$
means we need to divide the $k=0,l=0$ terms appropriately,
$c_{k,l}=c_{k,l}^{\prime}/(1+\delta_{k,0})(1+\delta_{l,0})$.

The above procedure readily generalizes to higher dimensions. However, for
high degrees or high dimensions it may become computationally more efficient
to compute a low-rank approximation of $f$ and approximate the univariate,
constituent functions as Chebyshev polynomials [71].

Error bounds are more difficult to establish rigorously for the multivariate
case. Ref. [72] shows that multivariate functions with a certain property
display similar behavior of the maximum error as the univariate case, namely
exponential convergence to the target function with increasing (total or
maximal) degree.

##  Appendix B Classical preprocessing for the Fourier method

Here we discuss two preprocessing steps that may be required for implementing
the Fourier approach. We used those methods in Sec. 3.2. First, we describe
one way to periodically extend functions which are not periodic with period 2.
Second, we show how to compute coefficients in a degree-$d$ Fourier series
approximation via the fast Fourier transform.

Without loss of generality consider a _bivariate_ target function
$f:[0,1]^{2}\rightarrow\mathbb{C}$, which is not periodic with period 2. To
obtain a converging Fourier series, we turn the function into a periodic one.
We begin by extending the domain to $[-1,1]^{2}$ by considering the function
$\tilde{f}(x,y)=f(|x|,|y|)$. This step guarantees that the function is
mirrored with respect to both coordinate axes and for both the real and
imaginary part. There are no jump discontinuities in the resulting function,
although the first derivative may be discontinuous. Next, we extend the domain
to $\mathbb{R}^{2}$ by considering the function

$\displaystyle\tilde{f}(x+2n,y+2m)=\tilde{f}(x,y)\qquad\forall
n,m\in\mathbb{N},$ (50)

which is a periodic function of period $2$ in both variables. Then we use this
extended target function for finding Fourier series coefficients.

Next we discuss a method for computing the coefficients of a degree-$d$
trigonometric polynomial approximation to a function via interpolation and the
fast Fourier transform. For simplicity let us consider the one-dimensional
case, i.e. a periodic function $f(x)$ with period 2 on the interval $[-1,1]$
(after applying a periodic extension if necessary). We assume $f$ admits a
Fourier series
$f(x)=\sum_{k=-\infty}^{\infty}\hat{c}_{k}\mathrm{e}^{\mathrm{i}\pi kx}$. Its
coefficients are given by

$\hat{c}_{k}=\frac{1}{2}\int_{-1}^{1}\mathrm{d}xf(x)\mathrm{e}^{-\mathrm{i}\pi
kx}.$ (51)

If those coefficients are known, a valid approximation is truncation of the
infinite Fourier series at degree $d$. If they are unknown, it is often
numerically favourable to seek an interpolant $f_{d}$ of the form Eq. (1) such
that $f_{d}(x_{m})=f(x_{m})$ on the set of points $x_{m}=2m/(2d+1)$ for
$m=0,\dots,2d$. The coefficients of $f_{d}$ are given by

$c_{k}=\frac{1}{2d+1}\sum_{m=0}^{2d}f(x_{m})\mathrm{e}^{-\mathrm{i}\tfrac{2\pi
km}{2d+1}}.$ (52)

This follows from multiplying $f_{d}(x_{m})$ by $\mathrm{e}^{-\mathrm{i}\pi
kx_{m}}$, summing over $m$ and using discrete orthogonality of the Fourier
basis functions:

$\displaystyle\sum_{m=0}^{2d}f_{d}(x_{m})\mathrm{e}^{-\mathrm{i}\pi
kx_{m}}=\sum_{m=0}^{2d}\sum_{l=-d}^{d}c_{l}\mathrm{e}^{\mathrm{i}\pi
x_{m}(l-k)}=\sum_{l=-d}^{d}c_{l}(2d+1)\delta_{l,k}=(2d+1)c_{k}.$ (53)

Inserting $x_{m}$ yields Eq. (52).

Note that Eq. (52) is the discrete Fourier transform of
$\\{f(x_{m})\\}_{m=0}^{2d}$, which can be evaluated with cost
$\mathcal{O}(d\log d)$ via the fast Fourier transform. This method generalizes
to $D$ dimensions and use of the $D$-dimensional fast Fourier transform with
cost $\mathcal{O}[d_{1}d_{2}\cdots d_{D}\log(d_{1}d_{2}\cdots d_{D})]$.

We do not discuss the convergence of multivariate Fourier series in general.
We shall mention however that the examples presented in this paper do enjoy
this property. For bivariate functions in $L^{p}$, with $p>1$, and with a
careful definition of the partial sums, the Fourier series converges as shown
by Fefferman [73, 74].

##  Appendix C Hardware experiment

###  C.1 Hardware specifications

For the experimental demonstration of our state preparation algorithm, we use
the H2-1 trapped-ion quantum computer. Here we summarize the relevant
specifications of H2-1 at the time of the experiments performed in April and
May 2024. For details, see Ref. [63]. H2-1 operates with 32 qubits implemented
with $S_{1/2}$ hyperfine clock states of ${}^{171}\mathrm{Yb}^{+}$ ions. The
native gate set consists of single-qubit rotation gates
$\text{Phased}X(\alpha,\beta)=R_{Z}(\beta)R_{X}(\alpha)R_{Z}(-\beta)$ and two-
qubit gates
$ZZ\text{Phase}(\theta)=\mathrm{e}^{-\mathrm{i}\frac{\theta}{2}Z\otimes Z}$,
parametrized by real angles $\alpha,\beta$ and $\theta$. All the gate
operations as well as initializations and measurements are performed in one of
four gate zones, each operating in parallel. The native two-qubit gates can be
applied to an arbitrary pair of qubits by shuttling ions to one of the gate
zones, enabling all-to-all connectivity. Average single- and two-qubit gate
infidelities are typically $0.003\%$ and $0.2\%$, respectively. State
preparation and measurement error is typically $0.2\%$ on average.

###  C.2 Implementation details

As discussed in Sec. 3.4, we consider the preparation of an uncorrelated and a
correlated 2D Gaussian distribution $f(x,y)$, Eq. (39), with the same mean
value $\bm{\mu}=(\mu_{x},\mu_{y})^{T}$, covariance matrix $\Sigma$
parameterized by the same marginal variances $\sigma_{x}^{2}$,
$\sigma_{y}^{2}$ and two values of the correlation coefficient $\rho$. We
prepare both functions on a 2D grid with $512\times 512$ grid points,
requirering $n_{x}=n_{y}=9$ qubits for the spatial discretization of both
axes.

For both settings, we compute the Fourier series coefficients in $f_{d}(x,y)$
with degree $d=d_{x}=d_{y}=3$ using the following method. First, we note that
both Gaussians are close to $0$ at the boundary of the domain $[0,1]^{2}$. We,
thus, extend $f(x,y)$ to $[-1,1]^{2}$ by setting the function to $0$ in the
extended domain. Next, assuming this function admits a Fourier series, its
coefficients can be approximated via

$\displaystyle\hat{c}_{k,\ell}=\frac{1}{4}\int_{[-1,1]^{2}}\mathrm{d}x\mathrm{d}y\,f(x,y)\mathrm{e}^{-\mathrm{i}\pi(xk+y\ell)}\approx\frac{1}{4}\int\mathrm{d}x\mathrm{d}y\,f(x,y)\mathrm{e}^{\mathrm{i}(x,y)^{T}\bm{\omega}}=\frac{1}{4}\mathrm{e}^{\mathrm{i}\bm{\mu}^{T}\bm{\omega}-\frac{1}{2}\bm{\omega}^{T}\bm{\Sigma}\bm{\omega}}\,.$
(54)

The first equality is a generalisation of Eq. (51) to the bivariate case. The
second approximate equality follows from extending the integration domain to
the full space $\mathbb{R}^{2}$, and setting $\bm{\omega}=-\pi(k,\ell)^{T}$.
This integral is the characterisic function of the bivariate Gaussian, for
which there exists an analytic formula, given by the last expression in Eq.
(54). The coefficients go to $0$ exponentially in $|k|$ and $|l|$, allowing us
to truncate the series at relatively small degrees $d_{x}$ and $d_{y}$. We
numerically observe that this method leads to a faster convergence for our
specific target function compared to the method used in Sec. 3.2. This allows
us to use a lower degree at a similar approximation quality. A lower degree
reduces circuit depth which is crucial for the experimental implementation.

Note that in the setting of the uncorrelated Gaussian, the density $f(x,y)$
[cf. Eq. (39)] and its Fourier series approximation $f_{d}(x,y)$ factorize
into a product of univariate Gaussian densities. This factorization can be
exploited in the circuit construction as well. More precisely, we prepare the
two univariate functions with two separate circuits as given in Fig. 4 (each
with $a=3,n=9$ qubits for the ancilla and main registers, respectively) and
concatenate those. This considerably simplifies the circuit and reduces the
gate complexity, in comparison to the full 2D state preparation circuit, cf.
Fig. 5. For the uncorrelated Gaussian we obtain $123$ single-qubit gates
($\text{Phased}X$), $80$ two-qubit gates ($ZZ\text{Phase}$) and a depth of
$42$ for the full state preparation circuit. For the correlated Gaussian we
obtain $266$ single-qubit gates, $237$ two-qubit gates and circuit depth
$352$.

The success probability of the algorithm is predominantly controlled by the
weight $\sum_{x}\sum_{y}|f_{d}(x,y)|^{2}$ of the prepared function in the
domain $[0,1]^{2}$, see Eq. (32). Based on a noiseless simulation of the
circuit, the ideal success probabilities are given by
$p_{\text{success}}=0.134$ and $p_{\text{success}}=0.139$ for the case of the
correlated and uncorrelated Gaussian, respectively. The interplay between
noise and circuit depth also affects the success probability. Deeper circuits
increase overall noise when executed on noisy hardware, which, in general,
reduces the overlap of the prepared (mixed) state with the target state
$\ket{0}^{\otimes(a_{x}+a_{y})}$ in the ancilla registers, cf. Fig. 5.
Generally, this reduces the success probability. Experimentally, we obtain a
success probability $p_{\text{success}}=$0.1299$$ for the uncorrelated setting
and a success probability $p_{\text{success}}=$0.0877$$ for the correlated
setting.

###  C.3 Kernel density estimation

Let
$\bm{X}=\\{\bm{x}_{1},\dots,\bm{x}_{N_{\mathrm{data}}}\\},\,\bm{x}_{i}\in\mathbb{R}^{D}$
be the observed data points in a given $D$-dimensional domain. Assume the data
is independently drawn from a density $f(\bm{x})$. The kernel density estimate
(KDE) of $f$ based on dataset $\bm{X}$ is given by

$k^{\bm{X}}_{h}(\bm{x})=\frac{1}{N_{\mathrm{data}}}\sum_{i=1}^{N_{\mathrm{data}}}K\left(\frac{\bm{x}-\bm{x}_{i}}{h}\right)\,,$
(55)

where $K$ is the kernel, i.e. a smooth normalized function (which integrates
to one over the given domain), which is often chosen to be Gaussian. The
parameter $h$ is the bandwidth. The KDE yields a smooth approximation of the
density $f$ underlying the observations $\bm{x}_{i}$. The bandwidth $h$ needs
to be chosen based on the competition between the smoothness of the estimate
and the information it contains about the dataset. For a small bandwidth $h$
the KDE tends to overfit the data, and in the limit of $h\rightarrow 0$ the
KDE just represents the original dataset without smoothing. For
$h\rightarrow\infty$, on the other hand, the KDE approaches a flat
distribution, which does not contain any information about the dataset
$\bm{X}$.

We optimize $h$ from the data, i.e. we do not assume access to the underlying
density $f$, using cross validation. For this let $\bm{X}_{i}$ be the dataset
obtained from $\bm{X}$ by excluding the point $\bm{x}_{i}$, and
$k^{\bm{X}_{i}}_{h}$ the corresponding KDE. The probability to observe the
excluded point is given by $k^{\bm{X}_{i}}_{h}(\bm{x}_{i})$. We want to
maximize this probability, averaged over all choices of $i$. In practice one
often works with log-probabilities. Note that, since the logarithm is a
monotonically increasing function, maximizing $k^{\bm{X}_{i}}_{h}(\bm{x}_{i})$
is equivalent to maximizing $\log k^{\bm{X}_{i}}_{h}(\bm{x}_{i})$. This
defines our final optimizer, the mean of log-probabilities taken over $i$:

$q(h)=\frac{1}{N_{\mathrm{data}}}\sum_{i=1}^{N_{\mathrm{data}}}\log
k^{\bm{X}_{i}}_{h}(\bm{x}_{i})\,.$ (56)

The optimal $\hat{h}$ is obtained by maximizing $q$ over $h$. The KDE estimate
of the target density $f$ is then given by

$k^{\bm{X}}=k^{\bm{X}}_{\hat{h}}\,.$ (57)

For a given $h$, the log-probability $\log k^{\bm{X}_{i}}_{h}(\bm{x}_{i})$ is
an unbiased estimator (under the distribution of $\bm{X}$) of the Kullback-
Leibler divergence (up to a sign and a constant additive factor depending on
$f$) between $f$ and $k^{\bm{X}}_{h}$. What this means is that
$\mathbb{E}[\log k^{\bm{X}_{i}}_{h}(\bm{x}_{i})]$ asymptotically approaches
$\mathbb{E}[\int\mathrm{d}\bm{x}f(\bm{x})\log k^{\bm{X}}_{h}(\bm{x})]$ for
$N_{\mathrm{data}}\rightarrow\infty$, where the expectation
$\mathbb{E}[\cdot]$ is taken with respect to the distribution of $\bm{X}$, see
e.g. [64, 65]. There exist alternatives to the log-probabilities which are,
for example, more closely related to the mean integrated squared error (MISE)
between the target density $f$ and the estimator [65], see also [75] for a
comparison of both.

Figure 11:  Cross validation of the kernel density estimates for
$\mathfrak{f}_{d}(x,y)=|f_{d}(x,y)|$ and its marginals. The results are shown
for the uncorrelated and the correlated Gaussian in the left and right plots,
respectively. The maximal value of the mean of log-probabilities $q(h)$, Eq.
(56), and the corresponding bandwidth $h$ are highlighted via vertical and
horizontal dashed lines. For the uncorrelated Gaussian, the maximum of $q(h)$
is achieved for $\hat{h}\approx 0.052$ for the full 2D distribution, and for
$\hat{h}_{x}\approx 0.034$ and $\hat{h}_{y}\approx 0.029$ for its marginals.
For the correlated Gaussian, we obtain optimal bandwidths $\hat{h}\approx
0.045$, $\hat{h}_{x}\approx 0.029$ and $\hat{h}_{y}\approx 0.022$.

In Fig. 11 we apply KDE to the data obtained from our state preparation
hardware experiments, cf. Sec. C.2. To optimize $q(h)$ [cf. Eq. (56)] we use
the scikit-learn [76] implementations of KDE and cross validation. We first
consider the datasets of the full 2D state preparation circuit. The dataset
$\bm{X}$ is given by the two-dimensional points $\bm{x}_{i}=(x_{i},y_{i})$ (as
represented via bitstrings, see Sec. 2.1) measured at the output of our state
preparation algorithm, see Appendix C.1 and C.2. In case of the uncorrelated
Gaussian, we have $N_{\mathrm{data}}=$2598$$ data points, for the correlated
one we have $N_{\mathrm{data}}=$1754$$. The yellow solid curves in both plots
of Fig. 11 show $q$ as a function of the bandwidth $h$, for both settings. We
obtain optimal bandwidths (vertical yellow dashed lines) of $\hat{h}\approx
0.052$ (uncorrelated Gaussian) and $\hat{h}\approx 0.045$ (correlated
Gaussian). We repeat the same procedure for the marginals. The blue and red
solid curves in Fig. 11 show $q(h)$ when we use
$\bm{X}=\\{x_{1},\dots,x_{N_{\mathrm{data}}}\\}$ (the x-marginal) and
$\bm{X}=\\{y_{1},\dots,y_{N_{\mathrm{data}}}\\}$ (the y-marginal),
respectively. Optimal marginal KDE are obtained for bandwidths
$\hat{h}_{x}\approx 0.034,\hat{h}_{y}\approx 0.029$ for the setting of the
uncorrelated Gaussian, and $\hat{h}_{x}\approx 0.029,\hat{h}_{y}\approx 0.022$
for the correlated Gaussian, as indicated by the vertical dashed blue and red
lines in the plots, respectively.