Quantum algorithms for approximate function loading

Gabriel Marin-Sanchez
Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain

Javier Gonzalez-Conde
[Corresponding author: ]<EMAIL_ADDRESS>Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain

Mikel Sanz
[Corresponding author: ]<EMAIL_ADDRESS>Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain
IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5, 48009, Bilbao, Spain

Loading classical data into quantum computers represents an essential stage in many relevant quantum algorithms, especially in the field of quantum machine learning. Therefore, the inefficiency of this loading process means a major bottleneck for the application of these algorithms. Here, we introduce two approximate quantum-state preparation methods inspired by the Grover-Rudolph algorithm, which partially solve the problem of loading real functions. Indeed, by allowing for an infidelity $\epsilon$ and under certain smoothness conditions, we prove that the complexity of Grover-Rudolph algorithm can be reduced from $\mathcal{O}(2^{n})$ to $\mathcal{O}(2^{k_0(\epsilon)})$, with $n$ the number of qubits and $k_0(\epsilon)$ asymptotically independent of $n$. This leads to a dramatic reduction in the number of required two-qubit gates. Aroused by this result, we also propose a variational algorithm capable of loading functions beyond the aforementioned smoothness conditions. Our variational ansatz is explicitly tailored to the landscape of the function, leading to a quasi-optimized number of hyperparameters. This allows us to  achieve high fidelity in the loaded state with high speed convergence for the studied examples.

Introduction.— Quantum computing has triggered a great interest in the last decades due to its theoretical capability to outperform classical information processing. Even though noise and decoherence are major drawbacks for the computational capacity of current quantum computers, quantum advantage has been experimentally achieved  [1, 2, 3]. Unfortunately, these accomplishments lack any industrial or scientific relevance, so the search of a useful application still remains. In this sense, the realistic experimental implementation of many promising quantum algorithms in several fields like solving systems of linear equations [4, 5], performing data fitting [6], computing scattering cross sections [7, 8], pricing financial derivatives [9, 10, 11] or initial conditions in differential equations [12, 13, 14], is constrained by the assumption that data can be efficiently loaded into a quantum device. In this context, the efficient loading of classical data into quantum computers is a particularly important problem, and represents a major bottleneck of the practical application of quantum computation in the NISQ-Era, especially with the emergence of the quantum machine learning field [20, 21, 17, 18, 19, 22].

There exits different quantum embedding techniques transforming classical data into quantum information  [20, 23]. In particular, we can distinguish two main embedding protocols depending on how the information is encoded. On the one hand, the basis embedding, in which each bit value $``0"$ or $``1"$ is map to a computational qubit state $|0\rangle$ or $|1\rangle$, respectively  [24]. In this way, the embedded quantum state corresponds to a uniform superposition of the bit-wise translations of binary strings. On the other hand, the amplitude-embedding technique encodes the normalized vector of classical data, which is now not necessarily binary, into the amplitudes of a quantum state [25, 26, 27, 32, 31, 29, 28, 30, 33, 34, 35, 36, 37]. In particular, these feature maps have been proposed to load discretized real valued functions [28, 29, 26, 27] with relevant applications in loading initial conditions for solving partial derivatives equations [11, 12, 13, 14], computing Monte-Carlo integrations [10, 9, 38] and quantum field theory [39, 40]. However, a practical implementation of these approaches generally incurs into an overhead of resources, which can be reflect into either an exponential number of entangling gates [28, 31, 29, 41, 42, 30], or the employment of a huge number of ancillary qubits [33, 34, 35]. A rather different approach sustained by the Solovay–Kitaev theorem [43, 44] is based on the application of quantum generative models to efficiently accomplish an approximate amplitude encoding of discretized real valued functions [45, 46]. Nonetheless, in these generic ansätze, increasing the number of hyperparameters does not necessarily reflect in improving the expressability of the ansatz to capture the function details [47]. Additionally, these variational methods usually suffer from training problems such as local minima and barren plateaus [48].

In this article, we present two approximate quantum algorithms to load real functions into quantum computers. Our first protocol, inspired by the Grover-Rudolph algorithm [28], reduces its complexity from $\mathcal{O}(2^{n})$ to $\mathcal{O}(2^{k_0(\epsilon)})$,
where $n$ is the number of qubits, $\epsilon$ the infeldelity respect to the exact state, and $k_0(\epsilon)$ is asymptotically independent of $n$. This leads to a dramatic reduction in the number of required two-qubit gates. Inspired in this result, we also introduce and benchmark a variational quantum circuit with applications in a broader family of functions. Our proposed ansatz is adapted to the structure of the function, which intuitively correlates hyperparameters and expressability. Moreover, by taking the angles provided by the Grover-Rudolph protocol, we can define a suitable initial training angle set, which considerably improves the training process, avoiding barren plateaus and local minima.

The article is structured as follows. First, we briefly review the Grover-Rudolph algorithm. Then, we propose and analyze step by step our first protocol. Next, we present our variational quantum circuit and, finally, we show the results of both methods.

The Grover-Rudolph Algorithm.— The method of Grover and Rudolph, originally proposed in Ref. [28], describes a constructive protocol to load into an $n$-qubit quantum state the discretized version $\{f_i\}$ of certain integrable density function $f\colon [x_\text{min}, x_\text{max}] \subset $IR$  \to $IR$^{+}$ [39] as

\begin{equation}
\label{eq:grover}
\ket{\Psi(f)}_n = \sum_{i=0}^{2^n-1}\sqrt{f_i}\ket{i}.
\end{equation}

The protocol provides a constructive algorithm which applies a sequence of operation blocks, $F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})$, to the initial state  $\ket{0}^{\otimes n}$. For $k\in \{1,\dots,n\}$, each $k$-qubit block, $F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})= \sum_{l = 0}^{2^{k-1}-1}\ket{l}\bra{l}\otimes R_{y}(\theta_{l}^{(k-1)})$, corresponds to a uniformly-controlled $y$-axis rotation (see Appendix [ap:ucr] A), where the $l$-th component of the angle vector $\bm{\theta}^{(k-1)}$ is calculated as

\begin{equation}
\label{eqn:angles}
\theta^{(k-1)}_{l} = 2\arccos\left(\sqrt{\frac{\int_{x_\text{min}+ l\delta_k}^{x_\text{min} + (l+1/2)\delta_k}f(x)dx}{\int_{x_\text{min}+ l\delta_k}^{x_\text{min} + (l+1)\delta_k}f(x)dx}}\right).
\end{equation}

Here, $\delta_k = \frac{x_\text{max}-x_\text{min}}{2^{k-1}}$ and $l\in \{0, \dots, 2^{k-1}-1\}$ is the index corresponding to the $(l+1)$-th subinterval of the $2^{k-1}$ partition of the interval $[x_\text{min}, x_\text{max}]$. Each multi-controlled gate comprising $F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})$ is denoted by $\wedge_{k-1}\left(R_y(\theta^{(k-1)}_{l})\right)$ and bisects the $(l+1)$-th partition interval of the function by using the conditional probability of being in right or left side of the interval, as depicted in Fig. <ref>.

operator/.append style=fill=mycolor, opacity=0.3
[column sep=0.2cm, row sep=0.2cm]
[wires=5]$\ket{0}^{\otimes n}$      …          1[3,steps=4,style=rounded corners,fill=gray!25, inner xsep=2pt, opacity=0.4,background]$F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})$      1[3,steps=1,style=dashed,rounded corners,fill=gray!50, inner xsep=0.3pt, opacity=0.4,background]     1     1          …

…         1      1     1     1          …

…         R_y(θ^(2)_0)     R_y(θ^(2)_1)     R_y(θ^(2)_2)     R_y(θ^(2)_3)          …

⋮           ⋮      ⋮      ⋮      ⋮           ⋮

…                              …
[inner sep=0pt] at (0.5, -2.7);

Effect of the uniformly-controlled $y$-axis rotation $F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})$ for $k=3$ and $n$ qubits. Each multi-controlled gate comprising $F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})$ bisects the $(l+1)$-th partition interval of the function by using the conditional probability of being in right or left side of the interval.

Moreover, Möttönen et al. [29] presented a generalization of this algorithm showing how to prepare an arbitrary $n$-qubit quantum state. This protocol requires $2^{n+2}-4n-4$ CNOT gates. Additionally, in [41, 42], the same authors proved that each block $F_{k}^{(k-1)}(\bm{y}, \bm{\theta}^{(k-1)})$ can be implemented with $2^{k-1}$ CNOTs and $2^{k-1}$ single-qubit gates. Therefore, we can conclude that this protocol and its generalizations are theoretically capable to load a discretized density function at the cost of an exponential overhead of resources to prepare the state in Eq. (<ref>).

First Algorithm.— Inspired by the Grover-Rudolph algorithm [28], we present an efficient method to encode discretized density functions into quantum states. By permitting an error in the final state and assuming certain smoothness conditions, an angle clustering significantly reduces the required entangling gates.

Let us firstly define the $n$-qubit representative state of a density function $f(x)$.

Let $f\colon [x_\text{min}, x_\text{max}]\to $IR$^{+}$ be a positive function in $L^2([x_\text{min}, x_\text{max}])$. We define the $n$-qubit normalized representative state of $f(x)$ as the $n$-qubit state $|f(x)\rangle_n=\sum_{j=0}^{2^n-1} f(x_\text{min}+j\delta_n)|j\rangle$, with $\delta_n=\frac{x_\text{max}-x_\text{min}}{2^n-1}$ and $\sum_{j=0}^{2^n-1}f^2(x_\text{min}+j\delta_n) = 1$.

Note that, according to this definition, we encode the discretized function into the amplitude of a quantum state. This is different to Grover-Rudolph algorithm, which does it in the probability, as shown in Eq. (<ref>).

In order to standardize the method, we first transform the interval of the target function to $[0, 1]$, so the units become dimensionless. This change of variables is given by
$ x' = x_{\text{min}} + xL$, with $x'\in[x_{\text{min}}, x_{\text{max}}]$, $x\in[0, 1]$ and $L = x_{\text{max}}-x_{\text{min}}$. From now on, we consider our function standardized, which enables us to state the main theorem that stands this first algorithm.

Theorem 1.— Let $f:[0,1]\rightarrow$IR$^+$ be a positive integrable function in $L^2([0,1])$ and $0\leq\eta\leq 8\pi$ a constant such that $\eta=\sup_{x\in[0,1]}\left|\partial^2_x \log f^2(x) \right|$. Then, it is posible to approximate the $n$-qubit representative state of $f(x)$, $|f(x)\rangle_n$, by a quantum state $|\Psi(f)\rangle_n$ such that the fidelity $|\langle \Psi(f)|f(x)\rangle_n|^2 \geq 1-\epsilon$ with at most $2^{k_0(\epsilon)}-1$ two-qubit gates, with

\begin{equation}\label{eq:k0}
k_0(\epsilon) = \max \left \{ \lceil -\frac{1}{2}\log_2 (4^{-n}-\frac{96}{\eta^2}\log (1-\epsilon))\rceil,2 \right \},
\end{equation}

and the circuit to perform it is provided in Fig. <ref>[fig:circuit-heatmap] (i).

(i) Quantum circuit performing the protocol presented in Theorem 1, based on the Grover-Rudolph method for a system of $n$ qubits. We cluster the angles of the blocks for $k\geq k_0+1>2$, leading to a drastic reduction in the number of gates needed. (ii) Values of $k_0$ for $\eta\in[0, 8\pi]$, $\epsilon\in[0.0001, 0.01]$ and $n\to\infty$, following Eq. (<ref>). We can appreciate that $k_0\leq 10$ for most of the cases.

In the following, let us analyze each part of the algorithm. Firstly, we will prove a bound for the difference between consecutive angles and the algebraic properties of the density function.  The following lemma provides a sufficient condition on the target density function to guarantees an upper bound over the difference between two contiguous angles.

Lemma 1.— Let $f$ be a continuous function such that $f:~[0, 1]\to$IR$^+$ and consider a block comprising a uniformly-controlled rotation of $k$ qubits. Then, the difference between two contiguous angles is bounded by the second derivative of the logarithm of $f$ in the following way:

\begin{equation}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| \leq \frac{\delta_k^2}{4}\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\left.\left(\partial^2_{y}\log f(y)\right)\right\rvert_{y=y'}\right|,
\label{eq4}
\end{equation}

where $\delta_k=\frac{1}{2^{k-1}}$, $l\in\{0, \dots, 2^{k-1}-2\}$.

The proof of this lemma is provided in the Appendix [ap:proofs] B.
As a consequence of this first lemma, if $\left|\partial^2_{y}\log f(y)\right| \leq \eta$, then for each block comprising a uniformly-controlled rotation of $k$ qubits, the difference between any two angles is bounded by

\begin{equation}
\label{eq:cluster}
\left|\theta^{(k-1)}_{l}-\theta^{(k-1)}_{l'}\right| \leq \frac{\delta_k}{4}\eta \ \ \ \   \forall l,l' \in \{0, \dots, 2^k-1\}.
\end{equation}

This result allows us to cluster angles of each block in which Eq. (<ref>) is fulfilled. Thus, we can define a cluster representative angle, $\tilde{\theta}^{(k-1)}$,  as

\begin{equation}
\label{eq:representative}
\left|\tilde{\theta}^{(k-1)}-\theta^{(k-1)}_{l}\right| \leq \frac{\delta_k}{8}\eta\vcentcolon=\eta_k \ \ \ \   \forall l \in \{0, \dots, 2^k-1\},
\end{equation}

$\forall k\geq k_0$, where $k_0+1$ is the index of the first block in which Eq. (<ref>) is fulfilled (successive blocks will also verify it). Note the fact that $\left|\partial^2_{y}\log f(y)\right| \leq \eta$ is a sufficient condition for the bound in Eq. (<ref>), however if in a singularity point it grew slower than $\frac{1}{\delta^2_k}$, the difference between consecutive angles would still vanish. This extension of the theorem is studied in Appendix [ap:general] C.

(i) Simulations of $\ket{f(x)}_8$ and $\ket{\Psi(f)}_8$ resulting of loading the Normal Distribution according to Theorem 1, $n=8$, $\epsilon=0.05$ and different values of $\sigma$: (a) $\sigma=1$, (b) $\sigma=0.6$, (c) $\sigma=0.4$ and (d) $\sigma=0.3$. The numerical results for each experiment is given in Table <ref>. (ii) Simulations of $\ket{\Psi\left(\sin(x)\right)}_5$ and $\ket{\Psi\left(\bm{\theta}\right)}_5$ resulting of loading the normalized sine function by training the variational ansatz, learning rate $\gamma=1.5$, $n=5$, $k_0=2$, and different values of $p(k)$: (a) $p(k)=1$, (b) $p(k)=2$, (c) $p(k)=3$ and (d) $p(k)=k$.

Fidelity of loading $f(x)=\sin{x}$, $x\in[0,3\pi/2]$, for different combinations of $n\in\{5, \dots, 10\}$, learning rate $\gamma=1.5$ and $p(k)$, with $k_0=3$. The axis corresponding to the total number of angles is in logarithmic scale.

Let us now analyze how clustering the angles according to Eq. (<ref>) affects on the final error of the process, measured by means of the fidelity with respect to the exact discretized state $|f(x)\rangle_n$. Considering an $n$-qubit system, the unitary gate to prepare the quantum state representing the target density function can be written in terms of the blocks as $\mathfrak{U}_n =~\mathcal{U}_{n-1}\left(\bm{\theta}^{(n-1)}\right)\cdots\mathcal{U}_{0}\left(\bm{\theta}^{(0)}\right)$, where we define $\mathcal{U}_{k-1}\left(\bm{\theta}^{(k-1)}\right) \coloneqq~ F_{k}^{k-1}(\bm{y}, \bm{\theta}^{(k-1)})\otimes\iden{(n-k)}.$

Let $\tilde{\mathfrak{U}}_n$ denote the operation $\mathfrak{U}_n$ when the rotation angles corresponding to each block $k$ are replaced by a representative, $\tilde{\theta}^{(k-1)}$, such that its difference with any angle of the block is at most $\eta_k$, i.e. $|\theta_l^{(k-1)} -\tilde{\theta}^{(k-1)}|\leq \eta_k$ for $l=0, \dots, 2^{k-1}-1$ and $k=1, \dots, n$. Then, the following lemma can be proven:

Lemma 2.— Consider a system of $n$ qubits and an error $\eta_k$ between any angle of the $k$-th block and its representative such that $\eta_k\leq\pi$, with $k=1, \dots, n$. Then, the fidelity between the final states with and without clustering, $F=|\brazero{n}\mathfrak{U}_n^\dagger \tilde{\mathfrak{U}}_n\ketzero{n}|^2$, satisfies

\begin{equation}
F \geq \prod_{k=1}^n\cos^{2}\left(\eta_k/2\right).
\end{equation}

The proof of this lemma is provided in Appendix [ap:proofs] B. If we now assume that we cluster angles of the blocks for $k\geq k_0+1>2$, then

\begin{equation}
F \geq \prod_{k=k_0+1}^n\cos^{2}\left(\frac{\eta_k \delta_{k>k_0}}{2}\right)\geq e^{-\frac{\eta^2}{24}(4^{-k_0}-4^{-n})}\vcentcolon= F_{k_0}.
\end{equation}

Therefore, if an infidelity $\epsilon=1-F_{k_0}$ is allowed and the angles of all blocks comprised of more than $k_0(\epsilon)$ qubits are clustered, with $k_0$ given by Eq. (<ref>), then the fidelity satisfies $F\geq 1-\epsilon$.

In the asymptotic limit of $n\to\infty$, we have that $k_0(\epsilon)$ tends to

\begin{equation}\label{eq:k0-limit}
k_0(\epsilon) \to \max \left \{ \lceil -\frac{1}{2}\log_2 (-\frac{96}{\eta^2}\log (1-\epsilon))\rceil,2 \right \},
\end{equation}

which is independent of the system size, $n$. In Fig. (<ref>)[fig:results] (ii), we have depicted the values of $k_0$ at this limit, for $\eta\in[0, 8\pi]$ and $\epsilon\in[0.0001, 0.01]$.

Let us finally study the implementation cost of the proposed protocol. Collodo et al. [50] proved that both the error and execution time corresponding to the implementation of single-qubit gates is several orders of magnitude smaller than for two-qubit gates. Thus, we only take into account the latter type of gates and ignore single-qubit operations. Using the result of [41, 42], which illustrates how a uniformly-controlled rotation of $k$ qubits can be implemented with $2^{k-1}$ CNOTs, the complexity of the circuit described in Theorem 1 is $\mathcal{O}(2^{k_0(\epsilon)})$, instead of the $\mathcal{O}(2^n)$ required by the Grover-Rudolph algorithm.

Variational Quantum Circuit.— Based on the previous protocol, we propose a variational ansatz for loading functions beyond the conditions required by Theorem 1, i.e. functions containing zeros or singular points, which lead to singularities in the second logarithmic derivative.

We propose a variational ansatz inspired by Grover-Rudolph, in which each $k$-qubit block satisfying $k>k_0$ comprises $(z+s+1)p(k)$ hyperparameters. Here, $z$ and $s$ are respectively the number of zeros and singular points in the function, $p(k)$ is polynomial number in $k$ denoting the number of hyperparameters allowed per singular point/zero for the $k$-th block, and $k_0=\text{max}\{\ k \  | \ z+s+1\geq 2^k\}$. Assuming $p(k)\geq1$, the minimum number of hyperparameters needed to capture the singular behavior of the function is $z+s+1$. Therefore, for each $k$-qubit block $F_{k}^{(k-1)}(\bm{\hat{y}}, \bm{\theta})$ comprising more than $z+s+1$ parameters, we cluster the angles which do not correspond to the position of zeros or singularities. More details of this proccedure are detailed in Appendix [ap:variational] D.
In this form, this proposal establishes an intuitive correlation between the number of hyperparameters $p(k)$ and the expressability of the circuit to capture the details of the target functions in the most relevant points. By using the decomposition of the multi-controlled rotations [41, 42], the number of two-qubit gates which must be included in the variational circuit sums up to

\begin{equation*}
\# \text{TQGs} = \sum_{k=1}^{k_0}2^{k-1}+(s+z)\sum_{k=k_0+1}^{n}p(k)(80k-398)
\end{equation*}

\begin{equation}
=2^{k_0}-1+r_{k_0+1}(n, s+z),
\end{equation}

where $m$ is the number of zeros/singularities and $r_{k_0+1}(n, s+z)$ grows polynomially with the system size $n$. A remarckable advantage of this ansatz is the training procedure, since Grover-Rudolph algorithm provides a suitable set of initial training angles which considerably enhances the convergence of the protocol with respect to a random initialization. This fact, together with a scaling in the number of hyperparameters $p(k)$ substantially slower than the system size, allow us to avoid that the training procedure to get stuck into both local minima and barren plateaus. A detailed numerical discussion about this is provided in Appendix [ap:variational] D.

Results.— Let us now illustrate the behavior of our techniques with different density functions.

(i) Normal Distribution.—  First, we apply the algorithm given by Theorem 1 to a normal distribution with a mean value $\mu=0.5$ and for different values of the variance $\sigma$. We numerically benchmark the fidelity attained by our first protocol using the value of $k_0$ resulting of Eq.(<ref>) when we assume a maximum infidelity of $\epsilon=0.05$ and a system of $8$ qubits. In Fig. <ref>[fig:results] (i) and Table <ref>, we have depicted the results from the simulations of $\ket{f(x)}_8$ and $\ket{\Psi(f)}_8$, for different values of $\sigma$. We can appreciate that the condition of $F\geq0.95=1-\epsilon$ is not only satisfied in all cases, but the fidelities obtained are considerably better. Furthermore, it is noteworthy the significant reduction in the number of two-qubit gates required to achieve these results. In the worst case, for $\sigma=0.3$, the quantity of gates needed represents the $12.16\%$ of the original set, while the fidelity of the experiment reaches a remarkable $0.99963$.

Numerical data for the simulations depicted in Fig. <ref>[fig:results](i). The values of $k_0$ have been computed using Eq. (<ref>). In addition, the number of two-qubit gates is given by $2^{k_0}-1$, and the last column is the percentage between the required gates and the total given by the Grover-Rudolph algorithm, which for $n=8$ are 255.

$\sigma$   $\eta$   $k_0$   Fidelity   #TQG   %TQG

1.0     2.00      2     0.99961     3      1.18

0.6     5.56      3     0.99931     7      2.75

0.4     12.50     4     0.99931     15     5.88

0.3     22.22     5     0.99961     31     12.16

(ii) Sine Function.— As a second example, we have tested the variational circuit with the normalized sine function $\sin{x}$ in the domain $[0, \frac{3}{2}\pi]$. This example contains zeros in $x=0$ and $x=\pi$ and clearly does not satisfy the conditions in Theorem 1 or a possible generalization, since it is not positive. In order to train the parametrized quantum circuit, we use the mean squared error loss function and the gradient descent as optimization method, with a learning rate $\gamma=1.5$. In Fig. <ref>[fig:results] (ii), we illustrate the arising loaded states of our trained circuits for different values of $p(k)=1$, 2, 3 and $k$ (the last case means the introduction of $k$ hyperparameters per zero of the function) and a system comprising $n=5$ qubits, with $k_0=2$. For the studied cases, our trained ansatz is able to load the sine function with a fidelity larger than $0.97$. Additionally, in Fig. <ref>, we show the resulting fidelity of the different combinations of the number of qubits and $p(k)$, for $n\in\{5, \dots, 10\}$ and $k_0=3$. See Appendix [ap:variational] D for further details on the training and additional examples.

Conclusions.— In this article, we have considered the problem of loading real valued functions into a quantum computer, which is a major bottleneck for solving partial derivatives equations [11, 12, 13, 14], computing Monte-Carlo integrations [10, 9, 38] and quantum field theory [39, 40] and quantum machine learning [20, 21, 17, 19, 22]. Firstly, inspired by the Grover-Rudolph algorithm [28], we have analytically proven a reduction of its complexity in an $n$-qubit system from $\mathcal{O}(2^{n})$ to $\mathcal{O}(2^{k_0(\epsilon)})$, with $\epsilon$ the infidelity respect to the exact state and $k_0(\epsilon)$ asymptotically independent of $n$. This reduction of two-qubit gates leads to a significant speedup, which allows us to implement quantum protocols involving data embeddings in large qubit systems. Additionally, we have generalized this method for functions containing a certain type of singularities, obtaining promising results for density functions with singularities that satisfy the expanded theorem conditions (see Appendix [ap:general] C). Furthermore, we have proposed a variational ansatz inspired in our previous protocol. We have observed that it can efficiently and accurately load functions with zeros and singularities. Our proposed ansatz is tailored to the landscape of the function, providing an intuitive correlation between hyperparameters and expressability. Moreover,  our previous protocol allows us to define a suitable initial training angle set, which considerably improves the training process, avoiding barren plateaus and local minima. As a future work, tensor networks could be used to prove the quasi-optimality in the minimal number of hyperparameters in the variational ansatz.

We thank  J.J. García-Ripoll for the useful discussions regarding the quasi-optimality of the variational ansatz. The authors acknowledge financial support from Spanish Government PGC2018-095113-B-I00 (MCIU/AEI/FEDER, UE), Basque Government IT986-16 and the QUANTEK project from ELKARTEK program (KK-2021/00070), as well as from QMiCS (820505) and OpenSuperQ (820363) of the EU Flagship on Quantum Technologies, and the EU FET-Open projects Quromorphic (828826) and EPIQUS (899368).

[1] F. Arute, K. Arya, R. Babbush et al., Quantum supremacy using a programmable superconducting processor. Nature 574, 505 (2019).

[2] F. Zhong, Y. Deng, J. Qin et al., Strong quantum computational advantage using a superconducting quantum processor. ArXiv:2106.14734 (2021)

[3] F. Zhong, H. Wang, Y. Deng et al., Quantum computational advantage using photons. Science 370, 1460 (2020).

[4] A. W. Harrow, A. Hassidim and S. Lloyd, Quantum algorithm for solving linear systems of equations. Phys. Rev. Lett. vol. 15, no. 103, pp. 150502 (2009).
[5] A. M. Childs, R. Kothari and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.
SIAM Journal on Computing 46, 1920-1950 (2017).

[6] N. Wiebe, D. Braun, and S. Lloyd, Quantum Algorithm for Data Fitting. Phys. Rev. Lett. 109, 050505 (2012).

[7] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned Quantum Linear System Algorithm.  Rev. Lett. 110, 250504 (2013).

[8] A. Scherer, B. Valiron, S.-C. Mau, S. Alexander, E. van den Berg, and T. E. Chapuran, Concrete resource analysis of the quantum linear-system algorithm used to compute the electromagnetic scattering cross section of a 2D target..  Quantum Information Processing 16 (2017).
[9] P. Rebentrost, B. Gupt, T. R. Bromley Quantum computational finance: Monte Carlo pricing of financial derivatives. Phys. Rev. A 98, 022321 (2018).

[10] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner, Option Pricing using Quantum Computers.  Quantum 4, 291 (2020).

[11] J. Gonzalez-Conde, A. Rodríguez-Rozas, E. Solano, and M. Sanz, Pricing Financial Derivatives with Exponential Quantum Speedup
.  arXiv:2101.04023 (2021).

[12] J. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations.  Proceedings of the National Academy of Sciences 118, e2026805118 (2021)

[13] B. Zanger, C. B. Mendl, M. Schulz, and M. Schreiber, Quantum Algorithms for Solving Ordinary Differential Equations via Classical Integration Methods. Quantum 5, 502 (2021).

[14] J. J. García-Ripoll, Quantum-inspired algorithms for multivariate analysis: from interpolation to partial differential equations. Int. J. Mod. Phys. C, 22,  Quantum 5, 431 (2021).

[17]  S. Lloyd, M. Mohseni and P. Rebentrost, Quantum algorithms for supervised and unsupervised machine learning. ArXiv:1307.0411v2 (2013).

[18]  J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, S. Lloyd,  Quantum Machine Learning. Nature 549, 195-202 (2017).

[19]  H. Huang, M. Broughton, M. Mohseni, R. Babbush, S. Boixo, H. Neven, and J. R. McClean,  Power of data in quantum machine learning. Nature Communications volume 12,  2631 (2021)

[20] M. Schuld and N. Killoran, Quantum machine learning in feature Hilbert spaces. Phys. Rev. Lett. 122, 040504 (2019).

[21] S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, and N. Killoran, Quantum embeddings for machine learning. ArXiv:2001.03622 (2020).

[22] N. Wiebe, Key questions for the quantum machine learner to ask themselves. New J. Phys. 22, 091001 (2020).

[23] V. Havlicek et al.,  Supervised learning with quantum enhanced feature spaces. Nature. vol. 567, 209 (2019).

[24] J. A. Cortese and T. M. Braje, Loading Classical Data into a Quantum Computer. arXiv:1803.01958 (2018).

[25] M. Schuld and F. Petruccione, Supervised Learning with Quantum Computers Springer Publishing Company, Incorporated  (2018).

[26] L. K. Grover, Synthesis of Quantum Superpositions by Quantum Computation. Phys. Rev. Lett. 85, 1334 (2000).

[27] Y. R. Sanders, G. H. Low, A. Scherer, and D.W. Berry, Black-box quantum state preparation without arithmetic. Phys. Rev. Lett. 122, 020502 (2019).

[28] L. Grover and T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions. ArXiv:quant-ph/0208112.

[29] M. Mottonen, J. J. Vartiainen, V. Bergholm, and M.M. Salomaa, Transformation of quantum states using uniformly controlled rotations. Quant. Inf. Comp. 5, 467 (2005).

[30] G. Kumar, Encoding classical data into quantum data. 10.13140/RG.2.2.28237.46563. (2020).

[31] M. Plesch and Č. Brukner,  Quantum-state preparation with universal gate decompositions Phys. Rev. A 83, 032302 (2011).

[32] J. Bausch, Fast Black-Box Quantum State Preparation. arXiv:2009.10709 (2020).
[33] I.F. Araujo, D. K. Park, F. Petruccione, and A. J. da Silva ,  A divide-and-conquer algorithm for quantum state preparation. Scientific Reports 11, 6329 (2021).

[34] X. Zhang, M. Yung, X. Yuan ,  Low-depth Quantum State Preparation. arXiv:2102.07533 (2021).

[35] J. Zhao, Y. Wu, G. Guo, and G. Guo, State preparation based on quantum phase estimation. arXiv:1912.05335 (2019).

[36] C. W. Bauer, P. Deliyannis, M. Freytsis, B. Nachman Practical considerations for the preparation of multivariate Gaussian states on quantum computers. 	arXiv:2109.10918 (2021).

[37] S. Anwer, A. Younes, I. Elkabani, A. Elsayed Preparation of quantum superposition using partial negation. 		arXiv:2109.14369 (2021).

[38] A. Montanaro, Quantum speedup of Monte Carlo methods. 	Proc. Roy. Soc. Ser. A, 471, 2181, 20150301, (2015).

[39] N. Klco and M. J. Savage, Digitization of Scalar Fields for Quantum Computing. Phys. Rev. A 99, 052335 (2019).

[40] S. P. Jordan, K. S. M. Lee, J. Preskill, Quantum Algorithms for Quantum Field Theories. Science 336, 1130-1133 (2012).

[41] V. Bergholm, J. J. Vartiainen, M. Mottonen, and M. M. Salomaa, Quantum circuits with uniformly controlled one-qubit gates . Phys. Rev. A 71, 052330 (2005).

[42] M. Mottonen, J. J. Vartiainen, V. Bergholm, and M. M. Salomaa Quantum circuits for general multi-qubit gates. Phys. Rev. Lett. 93, 130502 (2004).

[43] A. Y. Kitaev Quantum computations: algorithms and error correction. Russian Mathematical Surveys. 52 (6): 1191–1249. doi:10.1070/rm1997v052n06abeh002155. ISSN 0036-0279.

[44] Dawson, Christopher M.; Nielsen, Michael (2006-01-01). The Solovay-Kitaev algorithm. Quantum Information  Computation. 6: 81–95. arXiv:quant-ph/0505030. doi:10.26421/QIC6.1-6.

[45] C. Zoufal,   A. Lucchi, and  S. Woerner,  Quantum generative adversarial networks for learning and loading random distributions. Npj Quantum Information 5, 1–9 (2019).

[46] K. Nakaji et al,  Approximate amplitude encoding in shallow parameterized quantum circuits and its application to financial market indicator. arXiv:2103.13211 (2021).

[47] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, P. J. Coles  Variational Quantum Algorithms. Nature Reviews Physics 3, 625-644 (2021).

[48] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush and H. Neven.  Barren plateaus in quantum neural network training landscapes. Nature communications, 9, 1–6, (2018).

[50] M. C. Collodo, J. Herrmann, N. Lacroix, C. K. Andersen, A. Remm, S. Lazar, J.-C. Besse, T. Walter, A. Wallraff, and C. Eichler, Implementation of Conditional Phase Gates Based on Tunable ZZ Interactions. Physical Review Letters 125 (2020).

[51] A. Barenco et al Elementary gates for quantum computation. Phys. Rev. A52 3457  (1995).
[52] A. Cauchy, Methode generale pour la resolution des systemes d'equations simultanees. Comptes rendus de l'Académie des Sciences 25, 536 (1847).

[53]  F. Black and M. Scholes, The Pricing of Options and Corporate Liabilities. Journal of Political Economy 81, 637 (1973).

§ A  GATES DECOMPOSITION AND COMPLEXITY

§.§ Uniformly-Controlled Rotation

Let us consider $n-1$ control qubits, a target $n$-th qubit, an axis of rotation $\bm{u}$, and a vector of $2^{n-1}$ angles $\bm{\theta}$. Then, a uniformly controlled rotation $F_n^{n-1}(\bm{u}, \bm{\theta})$, depicted in Fig. <ref>, is a sequence of multi-controlled rotations $\wedge_{n-1}(R_{\bm{u}}(\theta_i))$ comprising the $2^{n-1}$ combinations of control bits.

Circuit corresponding to the uniformly controlled rotation gate $F_n^{n-1}(\bm{u}, \bm{\theta})$, with $n-1$ control qubits, a target at the $n$-th qubit, an axis of rotation $\bm{u}$, and a set of angles $\bm{\theta}$.

Analytically this gate can be expressed as

\begin{align}
F_{n}^{n-1}(\bm{u}, \bm{\theta}) &=\prod_{i = 0}^{2^{n-1}-1}\left(\ket{i}\bra{i}\otimes R_{\bm{u}}(\theta_{i})+\left(\mathbbm{1}_{n-1\times n-1} - \ket{i}\bra{i}\right)\otimes\mathbbm{1}\right),
\end{align}

and since $\sum_{i = 0}^{2^n-1}\ket{i}\bra{i} = \mathbbm{1}_{n\times n}$ and $\braket{i}{j}=\delta_{ij}$, we can simplify it to

\begin{align}
F_{n}^{n-1}(\bm{u}, \bm{\theta}) = \sum_{i = 0}^{2^{n-1}-1}\ket{i}\bra{i}\otimes R_{\bm{u}}(\theta_{i}).
\end{align}

§.§ Cost of Multi-Controlled Rotations

Let us now analyze which is the cost of implementing a multi-controlled rotation in terms of single and two-qubit gates.

Let us consider an $n$-qubit system. Then, the multi-controlled rotation $\wedge_{n-1}\left(R_y(\theta)\right)$ can be decomposed employing two one-qubit controlled rotations $\wedge_{1}(R_y(\pm\theta/2))$ and two $\wedge_{n-2}(X)$, as illustrated in Fig. <ref>.

operator/.append style=fill=violet!70!blue!20!

-1 1

R_y(θ)      =

1            1

1            1

⋮     ⋮           ⋮     ⋮

-1 2             -1 2

1           1

R_y(θ/2)            R_y(-θ/2)

Circuit showing the decomposition of $\wedge_{n-1}\left(R_y(\theta)\right)$ employing two controlled rotations $\wedge_{1}(R_y(\pm\theta/2))$ and two $\wedge_{n-2}(X)$ gates.

First, let us see that $R_y(\theta) = XR_y(-\theta/2)XR_y(\theta/2)$:

\begin{align}
&XR_y(-\theta/2)X = \begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}\cdot\begin{pmatrix} \cos(\nicefrac{\theta}{4}) & \sin(\nicefrac{\theta}{4})\\ -\sin\left(\nicefrac{\theta}{4}\right) & \cos(\nicefrac{\theta}{4})\end{pmatrix}\cdot\begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix} \nonumber\\& = \begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}\cdot\begin{pmatrix} \sin(\nicefrac{\theta}{4}) & \cos(\nicefrac{\theta}{4})\\ \cos\left(\nicefrac{\theta}{4}\right) & -\sin(\nicefrac{\theta}{4})\end{pmatrix} = \begin{pmatrix} \cos(\nicefrac{\theta}{4}) & \sin(\nicefrac{\theta}{4})\\ -\sin\left(\nicefrac{\theta}{4}\right) & \cos(\nicefrac{\theta}{4})\end{pmatrix} = R_y(\theta/2)
\end{align}

\begin{equation}
\label{eq:ry}
XR_y(-\theta/2)XR_y(\theta/2) = R_y(\theta/2)R_y(\theta/2) = R_y(\theta).
\end{equation}

The intuition behind this decomposition is that when all the control states are in $\ket{1}$, we obtain the desired rotation, as illustrated in  Eq. (<ref>). Also, when the controlled rotations are activated but the $\wedge_{n-2}(X)$ is not, we obtain the identity, since $R_y(-\theta/2)R_y(\theta/2)=\iden{1}$, as expected. The same occurs when the $\wedge_{n-2}(X)$ gates are activated but the controlled rotations are not. Finally, when none of the controls are triggered, the circuits in both sides are equal to $\iden{n}$.

Now, let us check that the circuits are equivalent. On the left hand side, we have:

\begin{align}
\wedge_{n-1}\left(R_y(\theta)\right) &= \ketone{(n-1)}\braone{(n-1)}\otimes R_y(\theta) \nonumber\\&+ \left(\iden{(n-1)}-\ketone{(n-1)}\braone{(n-1)}\right)\otimes\iden{1}.
\end{align}

On the right-hand side:

\begin{align}
&\wedge_{n-2}(X)\wedge_{1}(R_y(-\theta/2))\wedge_{n-2}(X)\wedge_{1}(R_y(\theta/2)) \nonumber\\&\quad=\left[\ketone{(n-2)}\braone{(n-2)}\otimes\iden{1}\otimes X + \left(\iden{(n-1)}-\ketone{(n-2)}\braone{(n-2)}\otimes\iden{1}\right)\otimes\iden{1}\right]\nonumber\\&\quad\cdot\left[\iden{(n-2)}\otimes\left(\ket{1}\bra{1}\otimes R_y(-\theta/2)+\ket{0}\bra{0}\otimes\iden{1}\right)\right]\nonumber\\
&\quad\cdot\left[\ketone{(n-2)}\braone{(n-2)}\otimes\iden{1}\otimes X + \left(\iden{(n-1)}-\ketone{(n-2)}\braone{(n-2)}\otimes\iden{1}\right)\otimes\iden{1}\right]\nonumber\\&\quad\cdot\left[\iden{(n-2)}\otimes\left(\ket{1}\bra{1}\otimes R_y(\theta/2)+\ket{0}\bra{0}\otimes\iden{1}\right)\right],
\end{align}

where we have ignored the identities on the left hand side of the equation for the sake of simplicity.
If we expand the previous expression, we obtain the following:

\begin{align}
&\left[\ket{1\cdots1}\bra{1\cdots1}\otimes XR_y(-\theta/2) + \ket{1\cdots10}\bra{1\cdots10}\otimes X + F(-\theta/2)\right]\nonumber\\
&\cdot \left[\ket{1\cdots1}\bra{1\cdots1}\otimes XR_y(\theta/2) + \ket{1\cdots10}\bra{1\cdots10}\otimes X + F(\theta/2)\right],
\end{align}

where we have defined the projector

\begin{align}
F(\alpha) &\equiv \left[\left(\iden{(n-1)}-\ket{1\cdots1}\bra{1\cdots1}-\ket{1\cdots10}\bra{1\cdots10}\right)\otimes\iden{1}\right]\nonumber\\
&\cdot\left[\iden{(n-2)}\otimes\left(\ket{1}\bra{1}\otimes R_y(\alpha)+\ket{0}\bra{0}\otimes\iden{1}\right)\right] \nonumber\\
& = \iden{(n-2)}\otimes\ket{1}\bra{1}\otimes R_y(\alpha) + \iden{(n-2)}\otimes\ket{0}\bra{0}\otimes \iden{1} \nonumber\\
&- \ket{1\cdots10}\bra{1\cdots10}\otimes\iden{1} -  \ket{1\cdots1}\bra{1\cdots1}\otimes R_y(\alpha).
\end{align}

It is straightforward to see that both terms $\ket{1\cdots1}\bra{1\cdots1}\otimes XR_y(\alpha)\cdot F(\alpha)$ and $\ket{1\cdots10}\bra{1\cdots10}\otimes X\cdot F(\alpha)$ are equal to $0$, due to the projector $F(\alpha)$. Then, we have

\begin{align}
&\wedge_{n-2}(X)\wedge_{1}(R_y(-\theta/2))\wedge_{n-2}(X)\wedge_{1}(R_y(\theta/2)) \nonumber\\&\quad= \ket{1\cdots1}\bra{1\cdots1}\otimes R_y(\theta)+\ket{1\cdots10}\bra{1\cdots10}\otimes \iden{1} +F(-\theta/2)F(\theta/2),
\end{align}

where we have used that $XR_y(-\theta/2)XR_y(\theta/2) = R_y(\theta)$ and $XX=\iden{1}$. Let us now compute the last term:

\begin{align}
&F(-\theta/2)F(\theta/2) \nonumber\\ &\quad= \iden{(n-2)}\otimes\ket{1}\bra{1}\otimes R_y(-\theta/2)R_y(\theta/2) - \ket{1\cdots1}\bra{1\cdots1}\otimes R_y(-\theta/2)R_y(\theta/2) \nonumber\\ &\quad +\iden{(n-2)}\otimes\ket{1}\bra{1}\otimes\iden{1} - \ket{1\cdots10}\bra{1\cdots10}\otimes \iden{1} - \ket{1\cdots10}\bra{1\cdots10}\otimes \iden{1} \nonumber\\ &\quad +\ket{1\cdots10}\bra{1\cdots10}\otimes \iden{1} - \ket{1\cdots1}\bra{1\cdots1}\otimes R_y(-\theta/2)R_y(\theta/2)\nonumber\\ &\quad +\ket{1\cdots1}\bra{1\cdots1}\otimes R_y(-\theta/2)R_y(\theta/2)\nonumber\\ &\quad = \left(\iden{(n-1)}-\ket{1\cdots1}\bra{1\cdots1}-\ket{1\cdots10}\bra{1\cdots10}\right)\otimes\iden{1},
\end{align}

where we have used that $R_y(-\theta/2)R_y(\theta/2)=\iden{1}$ and

\begin{equation}
\iden{(n-2)}\otimes\ket{1}\bra{1}\otimes\iden{1} + \iden{(n-2)}\otimes\ket{0}\bra{0}\otimes\iden{1}=\iden{(n-1)}\otimes\iden{1}.
\end{equation}

\begin{align}
&\wedge_{n-2}(X)\wedge_{1}(R_y(-\theta/2))\wedge_{n-2}(X)\wedge_{1}(R_y(\theta/2)) \nonumber\\&\quad = \ket{1\cdots1}\bra{1\cdots1}\otimes R_y(\theta)+\ket{1\cdots10}\bra{1\cdots10}\otimes \iden{1} \nonumber\\&\quad +\left(\iden{(n-1)}-\ket{1\cdots1}\bra{1\cdots1}-\ket{1\cdots10}\bra{1\cdots10}\right)\otimes\iden{1}\nonumber\\&\quad = \ket{1\cdots1}\bra{1\cdots1}\otimes R_y(\theta) + \left(\iden{(n-1)}-\ket{1\cdots1}\bra{1\cdots1}\right)\otimes\iden{1}\nonumber\\&\quad= \wedge_{n-1}\left(R_y(\theta)\right),
\end{align}

as we wanted to prove.

The order of complexity of implementing the $\wedge_{n-1}\left(R_y(\theta)\right)$ gate is linear in $n$, for $n>6$. In fact, $\wedge_{n-1}\left(R_y(\theta)\right)$ can be decomposed with $80n-398$ two-qubit gates.

Following Th. <ref>, we know that the multi-controlled rotation $\wedge_{n-1}\left(R_y(\theta)\right)$ can be decomposed into two one-qubit controlled rotations,  $\wedge_{1}\left(R_y(\theta/2)\right)$ and $\wedge_{1}\left(R_y(-\theta/2)\right)$, and two $\wedge_{n-2}\left(X\right)$. Therefore, we need to study the number of two-qubit gates necessary to implement the last two gates $\wedge_{n-2}\left(X\right)$.

Corollary $7.4$ from [51] states that, for $n>6$, $\wedge_{n-2}\left(X\right)$ can be realized with $8(n-5)\wedge_{2}\left(X\right)$. Then, we reduce the problem  to obtain the cost of the Toffoli gate, $\wedge_{2}\left(X\right)$.

Lemma 6.1 from the previous reference [51] states that any one-qubit unitary, $\wedge_{2}\left(U\right)$ can be implemented with 5 two-qubit gates comprising two CNOTs, two $\wedge_{1}\left(V\right)$, and one $\wedge_{1}\left(V^\dagger\right)$, with $V^2=U$. In our case, we have $U=X$. Then, we can choose $V = \frac{1-i}{2}(\iden{1}+iX)$. We provide the exact decomposition of the Toffoli gate using the previous description in Fig. <ref>.

operator/.append style=fill=violet!70!blue!20!

R_y(-π/4)          =

1         1    2

1          1

V       V^†         V

Circuit of the decomposition of the Toffoli gate, $\wedge_{n-2}\left(X\right)$ using five two-qubit controlled operations, where $V^2=X$.

Finally, if $n>6$, then the number of gates required for decomposing $\wedge_{n-1}\left(R_y(\theta)\right)$, considering only two-qubit operators, is

\begin{align}
\#{\text {gates}}\left(\wedge_{n-1}\left(R_y(\theta)\right)\right) = 2 + 2\cdot\left(8(n-5)\cdot5\right) = 80n-398.
\end{align}

§ B  THEOREM PROOFS

§.§ Relation Between the Difference in the Angles and the Target Function

Let $f$ be a continuous function such that $f:[0, 1]\to\mathbb{R}^+$ and consider a block comprising a uniformly controlled rotation of $k$ qubits. Then, the difference between two contiguous angles is bounded by the second derivative of the logarithm of $f$ in the following way:

\begin{equation}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| \leq \frac{\delta_k^2}{4}\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\left.\left(\partial^2_{y}\log f(y)\right)\right\rvert_{y=y'}\right|,
\end{equation}

where $\delta_k=\frac{1}{2^{k-1}}$, $l\in\{1, \dots, 2^{k-1}-1\}$.

The discrete expression of the angles for a block comprised of uniformly controlled rotation with $k$ qubits, first introduced in [28], is given by

\begin{equation}
\theta^{(k-1)}_l = 2\arccos\left(\sqrt{\frac{\int_{l\delta_k}^{(l+1/2)\delta_k}f(x)dx}{\int_{l\delta_k}^{(l+1)\delta_k}f(x)dx}}\right),
\end{equation}

with $\delta_k=\frac{1}{2^{k-1}}$ and $l\in\{1, \dots, 2^{k-1}-1\}$.
We can define the continuous extension of the previous function as

\begin{equation}
\theta^{(k-1)} (y) \vcentcolon= 2\arccos\left(\sqrt{\frac{\int_{y}^{y +\delta_k/2}f(x)dx}{\int_{y}^{y +\delta_k}f(x)dx}}\right),
\label{eq:thetas}
\end{equation}

where the original expression can be recovered replacing $y=l\delta_k$.

Now, given a function $g$ and a displacement $\mu$, we define the numerical derivative of $g$ with respect to $\mu$ as

\begin{equation}
\npartial{y}{\mu}g(y) \vcentcolon= \frac{g(y+\mu)-g(y)}{\mu}.
\label{eq:num-der}
\end{equation}

Notice that for $\mu \rightarrow 0$, $\npartial{y}{\mu}g(y) \rightarrow \partial_yg(y)$, the usual derivative. Then, for two consecutive angles, we have

\begin{equation}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| = \left|\left. \delta_k\left( \npartial{y}{\delta_k}\theta^{(k-1)}(y)\right)\right\rvert_{y=l\delta_k}\right|.
\end{equation}

Now, we require a connection between the difference in the angles and the partial derivative of Eq.(<ref>). The Mean Value Theorem guarantees that the numerical derivative, which corresponds to the slope of the straight line connecting $\theta^{(k-1)}_l$ and $\theta^{(k-1)}_{l+1}$, is constrained by the maximum absolute value of the exact derivative in that interval. Then, we have

\begin{align}\label{eq:der-ineq}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| \leq \delta_k\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\left.\left(\partial_{y}\theta^{(k-1)}(y)\right)\right\rvert_{y=y'}\right|.
\end{align}

Let us now develop the term of the exact derivative of the function $\theta^{(k-1)}(y)$:

\begin{equation}
\partial_{y}\theta^{(k-1)}(y) = \frac{\left(f(y+\delta_k)-f(y)\right)\int_y^{y+\delta_k/2}f(x)dx - \left(f(y+\delta_k/2)-f(y)\right)\int_y^{y+\delta_k}f(x)dx}{\int_y^{y+\delta_k}f(x)dx \sqrt{\int_y^{y+\delta_k/2}f(x)dx \int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}}
\end{equation}

By introducing $\left(f(y+\delta_k/2)-f(y+\delta_k/2)\right)$ in the first term and reorganizing the expression, we have:

\begin{align}
&\partial_{y}\theta^{(k-1)}(y)\nonumber\\&\quad= \frac{\left(f(y+\delta_k)-f(y+\delta_k/2)\right)\int_y^{y+\delta_k/2}f(x)dx - \left(f(y+\delta_k/2)-f(y)\right)\left(\int_y^{y+\delta_k}f(x)dx - \int_y^{y+\delta_k/2}f(x)dx\right)}{\int_y^{y+\delta_k}f(x)dx \sqrt{\int_y^{y+\delta_k/2}f(x)dx \int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}} \nonumber\\
&\quad= \frac{\left(f(y+\delta_k)-f(y+\delta_k/2)\right)\int_y^{y+\delta_k/2}f(x)dx - \left(f(y+\delta_k/2)-f(y)\right)\int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}{\int_y^{y+\delta_k}f(x)dx \sqrt{\int_y^{y+\delta_k/2}f(x)dx \int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}} \nonumber\\
&\quad= \frac{\sqrt{\int_y^{y+\delta_k/2}f(x)dx \int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}}{\int_y^{y+\delta_k}f(x)dx}\left(\frac{f(y+\delta_k)-f(y+\delta_k/2)}{\int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}-\frac{f(y+\delta_k/2)-f(y)}{\int_y^{y+\delta_k/2}f(x)dx}\right),
\end{align}

where the first term is positive. Now, since $\left(\sqrt{a}-\sqrt{b}\right)^2=a+b-2\sqrt{ab}\geq 0$, we have the well-known inequality for the geometric and arithmetic means $\sqrt{ab}\leq\frac{a+b}{2}$. Then, if we select $a=\frac{\int_y^{y+\delta_k/2}f(x)dx}{\int_y^{y+\delta_k}f(x)dx}$ and $b=\frac{\int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}{\int_y^{y+\delta_k}f(x)dx}$, it leads to

\begin{equation}
\frac{\sqrt{\int_y^{y+\delta_k/2}f(x)dx \int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}}{\int_y^{y+\delta_k}f(x)dx} \leq \frac{1}{2}\frac{\int_y^{y+\delta_k/2}f(x)dx + \int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}{\int_y^{y+\delta_k}f(x)dx}=\frac{1}{2}.
\end{equation}

With this result, we get the following inequality for the absolute value of the derivative of $\theta^{(k-1)}(y)$:

\begin{align}
\label{eq:angle-ineq}
\left|\partial_{y}\theta^{(k-1)}(y)\right| \leq \frac{1}{2}\left|\frac{f(y+\delta_k)-f(y+\delta_k/2)}{\int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}-\frac{f(y+\delta_k/2)-f(y)}{\int_y^{y+\delta_k/2}f(x)dx}\right|.
\end{align}

By using again the numerical derivative defined in Eq. (<ref>), we can simplify the previous expression in the following way:

\begin{align}
\frac{f(y+\delta_k)-f(y+\delta_k/2)}{\int_{y+\delta_k/2}^{y+\delta_k}f(x)dx}-\underbrace{\frac{f(y+\delta_k/2)-f(y))}{\int_y^{y+\delta_k/2}f(x)dx}}_{\vcentcolon=h(y)} = h(y+\delta_k/2) - h(y) = \frac{\delta_k}{2}\npartial{y}{\delta_k/2}h(y).
\end{align}

Also, since $f$ must be integrable, its primitive $F$ exists. Then,

\begin{align}
h(y)= \frac{f(y+\delta_k/2)-f(y)}{\int_y^{y+\delta_k/2}f(x)dx}=\frac{f(y+\delta_k/2)-f(y)}{F(y+\delta_k/2)-F(y)} = \frac{\npartial{y}{\delta_k/2}f(y)}{\npartial{y}{\delta_k/2}F(y)}.
\end{align}

Let us now plug it into the inequality

\begin{align}
\left|\partial_{y}\theta^{(k-1)}(y)\right| \leq \frac{1}{2}\left|\frac{\delta_k}{2}\cdot\npartial{y}{\delta_k/2}\left(\frac{\npartial{y}{\delta_k/2}f(y)}{\npartial{y}{\delta_k/2}F(y)}\right)\right| = \frac{\delta_k}{4}\left|\npartial{y}{\delta_k/2}\left(\frac{\npartial{y}{\delta_k/2}f(y)}{\npartial{y}{\delta_k/2}F(y)}\right)\right|.
\end{align}

If we follow the same argument as in Eq. (<ref>), we have that the term of the numerical derivative is upper bounded by the exact derivative, in the absolute value, which corresponds to the second derivative of the function's logarithm. Hence, the following inequality holds

\begin{equation}
\left|\partial_{y}\theta^{(k-1)}(y)\right|_{[l\delta_k, (l+1)\delta_k]} \leq \frac{\delta_k}{4}\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\left.\left(\partial^2_{y}\log f(y)\right)\right\rvert_{y=y'}\right|.
\end{equation}

Finally, by plugging this result in Eq. (<ref>), we obtain

\begin{align}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| \leq\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\frac{\delta_k^2}{4}\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\left.\left(\partial^2_{y}\log f(y)\right)\right\rvert_{y=y'}\right|\right|=\frac{\delta_k^2}{4}\max_{y'\in[l\delta_k, (l+1)\delta_k]}\left|\left.\left(\partial^2_{y}\log f(y)\right)\right\rvert_{y=y'}\right|,
\end{align}

as we wanted to proof.

Let $f$ be a continuous function such that $f:[0, 1]\to\mathbb{R}^+$ and consider a block of uniformly controlled rotations of $k$ qubits. Then, if $\exists \eta \geq 0$ such that $\left|\partial^2_{y}\log f(y)\right| \leq \eta ~\forall y\in[0, 1]$,

\begin{equation}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right|\leq \frac{\delta_k^2}{4}\eta.
\end{equation}

For a non-standardized function $f:[x_{\text{min}}, x_{\text{max}}]\to\mathbb{R}^+$, the result in Corollary <ref> holds with the modification in the bound

\begin{equation}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| \leq \frac{\delta_k^2\cdot\eta}{4L^2},
\end{equation}

where $L=x_{\text{max}}-x_{\text{min}}$ and $\delta_k=\frac{L}{2^{k-1}}$.

The change of variables that map the $x'\in[0, 1]$ with $x\in[x_{\text{max}}-x_{\text{min}}]$ is given by

\begin{equation}
x = x_{\text{min}} + x'L.
\end{equation}

Now, by using the chain rule, we obtain:

\begin{equation}
\frac{\partial f}{\partial y'} = \frac{\partial f}{\partial y}\frac{\partial y'}{\partial y'} =  \frac{\partial f}{\partial y}L \implies \left|\partial^2_{y}\log f(y)\right| \leq \frac{\tilde{\eta}}{L^2}.
\end{equation}

Let $\eta$ be such that $\left|\partial^2_{y}\log f(y)\right| \leq \eta ~\forall y\in[0, 1]$. Then, the difference between any two angles is

\begin{equation}
\left|\theta^{(k-1)}_{l}-\theta^{(k-1)}_{l'}\right|\leq\frac{\delta_k}{4}\eta,
\end{equation}

$\forall l, l'\in\{1, \dots, 2^{k-1}\}$.

The worst scenario is when $l=1$ and $l'=2^{k-1}$. In this case, by using the triangular inequality, we have

\begin{equation}
\left|\theta^{(k-1)}_{1}-\theta^{(k-1)}_{2^{k-1}}\right|\leq\left|\theta^{(k-1)}_{1}-\theta^{(k-1)}_{2}\right|+\cdots+\left|\theta^{(k-1)}_{2^{k-1}-1}-\theta^{(k-1)}_{2^{k-1}}\right|\leq\frac{\delta_k}{4}\eta,
\end{equation}

since there are less than $\delta_k^{-1}$ elements in the sum.

In this situation, we can define the representative angle $\tilde{\theta}^{(k-1)}$ for the clustering process as the one corresponding to the middle part of the interval, satisfying

\begin{equation}\label{eq:representative_ap}
\left|\tilde{\theta}^{(k-1)}-\theta^{(k-1)}_{l}\right|\leq\frac{\delta_k}{8}\eta \quad\forall l\in\{0, 1, \dots, 2^{k-1}-1\}.
\end{equation}

§.§ Error Bound of the Algorithm

§.§.§ Relation Between the Difference in the Angles and the Fidelity

Let us consider a system with $n$ qubits, thus the unitary gate to prepare the quantum state representing the target density function can be written in terms of the blocks as

\begin{equation}
\mathfrak{U}_n = \mathcal{U}_{n-1}\left(\bm{\theta}^{(n-1)}\right)\cdots\mathcal{U}_{0}\left(\bm{\theta}^{(0)}\right),
\end{equation}

where we define

\begin{align}
\mathcal{U}_{k-1}\left(\bm{\theta}^{(k-1)}\right) &\coloneqq F_{k}^{k-1}(\bm{y}, \bm{\theta}^{(k-1)})\otimes\iden{(n-k)}.
\end{align}

Let $\tilde{\mathfrak{U}}_n$ denote the operation $\mathfrak{U}_n$ given a representative with an error $\eta_k$ between the angles for each block, i.e. $|\theta_l^{(k-1)} -\tilde{\theta}^{(k-1)}|\leq \eta_k$ for $l=0, \dots, 2^{k-1}-1$ and $k=1, \dots, n$.

Consider a system of $n$ qubits and an error $\eta_k$ between any angle of the $k$-th block and its representative such that $\eta_k\leq\pi$, with $k=1, \dots, n$. Then, the fidelity between the final states with and without clustering, $F = |\brazero{n}\mathfrak{U}_n^\dagger \tilde{\mathfrak{U}}_n\ketzero{n}|^2$, satisfies

\begin{equation}
F \geq \prod_{k=1}^n\cos^{2}\left(\eta_k/2\right).
\end{equation}

We use induction to prove the inequality. Hence, we start with the elemental case of $n=1$ and later proceed assuming it is satisfied for $n-1$ and check if it holds for $n$:

* $n=1$:

\begin{align}
\bra{0}\mathfrak{U}_1^\dagger \tilde{\mathfrak{U}}_1\ket{0} =& \bra{0}R_y^\dagger(\theta^{(0)}) R_y(\tilde{\theta}^{(0)})\ket{0}\nonumber\\=& \bra{0}R_y^\dagger\left(\theta^{(0)} - \tilde{\theta}^{(0)}\right)\ket{0} = \cos\left(\frac{\theta^{(0)} - \tilde{\theta}^{(0)}}{2}\right)
\end{align}

Notice that the angles, given by

\begin{align}
\theta^{(n-1)}_l = 2\arccos\left(\sqrt{\frac{\int_{x_\text{min}+ l\delta_n}^{x_\text{min} + (l+1/2)\delta_n}f(x)dx}{\int_{x_\text{min}+ l\delta_n}^{x_\text{min} + (l+1)\delta_n}f(x)dx}}\right),
\end{align}

take values between $2\arccos(1)=0$ and $2\arccos(0)=\pi$. Then, since the cosine is a decreasing function in the interval $[0, \pi/2]$ and $\eta\leq\pi$, we have that

\begin{align}
|\theta^{(0)} - \tilde{\theta}^{(0)}|\leq \eta_1 &\implies
\cos\left(\frac{\theta^{(0)} - \tilde{\theta}^{(0)}}{2}\right) \geq \cos\left(\eta_1/2\right) \nonumber\\&\implies F \geq \cos^{2}\left(\eta_1/2\right).
\end{align}

* $n>1$: Assume the condition holds for $n-1$. We aim to recover the expression of $\brazero{n-1}\mathfrak{U}_{n-1}^\dagger \tilde{\mathfrak{U}}_{n-1}\ketzero{n-1}$ so we can use the induction hypothesis. Therefore, as we did for $n=1$, we sandwich the operator

\begin{equation}
\left(F_{n}^{n-1}(\bm{y}, \bm{\theta}^{(n-1)})\right)^\dagger F_{n}^{n-1}(\bm{y}, \bm{\tilde{\theta}}^{(n-1)})
\end{equation}

with the state $\ketzero{1}$ at each side. From now on during this proof, to simplify the notation,  we denote the gate $F_{k}^{k-1}(\bm{y}, \bm{\tilde{\theta}}^{(k-1)})$ as $U_{k-1}(\bm{\theta}^{(k-1)})$. Then, for $n$:

\begin{align}
&\brazero{n}\mathfrak{U}_n^\dagger \tilde{\mathfrak{U}}_n\ketzero{n} \nonumber\\ & =\brazero{n} \mathcal{U}^\dagger_{0}\left(\theta^{(0)}\right)\cdots \mathcal{U}^\dagger_{n-1}\left(\bm{\theta}^{(n-1)}\right)\mathcal{U}_{n-1}\left(\bm{\tilde{\theta}}^{(n-1)}\right)\cdots
\mathcal{U}_{0}\left(\tilde{\theta}^{(0)}\right)\ketzero{n}\nonumber\\
&=\brazero{(n-1)}\otimes\bra{0} \left(U_0^\dagger(\theta^{(0)})\otimes\iden{(n-1)}\right)\cdots
\left(U_{n-2}^\dagger(\bm{\theta}^{(n-2)})\otimes \mathbbm{1}\right)
\cdot\left(U_{n-1}^\dagger(\bm{\theta}^{(n-1)})\right)\nonumber\\&\cdot\left(U_{n-1}(\bm{\tilde{\theta}}^{(n-1)})\right)
\left(U_{n-2}(\bm{\tilde{\theta}}^{(n-2)})\otimes \mathbbm{1}\right)\cdots \left(U_0(\tilde{\theta}^{(0)})\otimes\iden{(n-1)}\right)\ketzero{(n-1)}\otimes\ket{0}\nonumber\\
&=\brazero{(n-1)}\left(U_0^\dagger(\theta^{(0)})\otimes\iden{(n-2)}\right)\cdots \left(U_{n-2}^\dagger(\bm{\theta}^{(n-2)})\right)\nonumber\\&
\cdot\left(\sum_{i_1, \ldots, i_{n-1} = 0}^1 \ket{i_1\ldots i_{n-1}}\bra{i_1\ldots i_{n-1}}\cdot \bra{0}  R_y^\dagger(\theta_{i_1\ldots i_{n-1}}^{(n-1)}-\tilde{\theta}^{(n-1)})\ket{0}\right)\nonumber\\
&\cdot\left(U_{n-2}(\bm{\tilde{\theta}}^{(n-2)})\right)\cdots \left(U_0(\tilde{\theta}^{(0)})\otimes\iden{(n-2)}\right)\ketzero{(n-1)}
\end{align}

Notice that we now we can substitute the rotation terms in the following way:

\begin{align}
\bra{0}  R_y^\dagger(\theta_{i_1\ldots i_{n-1}}^{(n-1)}-\tilde{\theta}^{(n-1)})\ket{0} = \cos\left(\frac{\theta_{i_1\ldots i_{n-1}}^{(n-1)}-\tilde{\theta}^{(n-1)}}{2}\right)\geq\cos\left({\eta_n/2}\right),
\end{align}

since $\theta_{i_1\ldots i_{n-1}}^{(n-1)}, \tilde{\theta}^{(n-1)}, \eta_n\in[0, \pi]~\forall i_1, \dots i_{n-1}\in\{0, 1\}$. However, before introducing the inequality, we need to check that the rest of the terms have the same sign. To do so, we first sandwich the terms corresponding to the block $n-2$ with the state $\ketzero{1}$ as well:

\begin{align}
&\brazero{n}\mathfrak{U}_n^\dagger \tilde{\mathfrak{U}}_n\ketzero{n} \nonumber\\
\left(U_0^\dagger(\theta^{(0)})\otimes\iden{(n-3)}\right)\cdots \left(U_{n-3}^\dagger(\bm{\theta}^{(n-3)})\right)\nonumber\\&
\left(\sum_{i_1, \ldots, i_{n-1} = 0}^1 \ket{i_1\ldots i_{n-2}}\bra{i_1\ldots i_{n-2}}\cdot \bra{0}  R_y^\dagger(\theta_{i_1\ldots i_{n-2}}^{(n-2)})\ket{i_{n-1}}\bra{i_{n-1}}R_y(\tilde{\theta}^{(n-2)})\ket{0} \right.\nonumber\\& \cdot\left.\cos\left(\frac{\theta_{i_1\ldots i_{n-1}}^{(n-1)}-\tilde{\theta}_{i_1\ldots i_{n-1}}^{(n-1)}}{2}\right)\right)
\cdot \left(U_{n-3}(\bm{\tilde{\theta}}^{(n-3)})\right)\cdots
\left(U_0(\tilde{\theta}^{(0)})\otimes\iden{(n-3)}\right)
\ketzero{(n-2)}.
\end{align}

Here, since all the angles are between 0 and $\pi$, we have that

\begin{align}
&\bra{0}R^\dagger_y(\theta_{i_1\ldots i_{n-2}}^{(n-2)})\ket{0}\bra{0}R_y(\tilde{\theta}^{(n-2)})\ket{0}=\cos\frac{\theta_{i_1\ldots i_{n-2}}^{(n-2)}}{2}\cos\frac{\tilde{\theta}^{(n-2)}}{2}\geq 0,\\
&\bra{0}R^\dagger_y(\theta_{i_1\ldots i_{n-2}}^{(n-2)})\ket{1}\bra{1}R_y(\tilde{\theta}^{(n-2)})\ket{0}=\sin\frac{\theta_{i_1\ldots i_{n-2}}^{(n-2)}}{2}\sin\frac{\tilde{\theta}^{(n-2)}}{2}~\geq 0,
\end{align}

$\forall i_1, \dots i_{n-2}\in\{0, 1\}$. Then, by proceeding analogously with the rest of the blocks, we obtain that all terms are positive. Hence, we can apply the inequality:

\begin{align}
&\brazero{n}\mathfrak{U}_n^\dagger \tilde{\mathfrak{U}}_n\ketzero{n} \nonumber\\
&\geq \brazero{(n-1)}\left(U_0^\dagger(\theta^{(0)})\otimes\iden{(n-2)}\right)\cdots\left(U_{n-3}^\dagger(\bm{\theta}^{(n-3)})\otimes\mathbbm{1}\right)\nonumber\\
&\cdot\left(\sum_{i_1, \ldots, i_{n-2} = 0}^1 \ket{i_1\ldots i_{n-2}}\bra{i_1\ldots i_{n-2}}\cdot \otimes \sum_{i_{n-1} = 0}^1 R_y^\dagger(\theta_{i_1\ldots i_{n-2}}^{(n-2)})\ket{i_{n-1}}\bra{i_{n-1}}R_y(\tilde{\theta}^{(n-2)})\right)\nonumber\\
& \cdot\left(U_{n-3}(\bm{\tilde{\theta}}^{(n-3)}\otimes\mathbbm{1})\right)\cdots \left(U_0(\tilde{\theta}^{(0)})\otimes\iden{(n-2)}\right)\ketzero{(n-1)}\cdot\cos\left({\eta_n/2}\right)\nonumber\\
&\cdot\left(\sum_{i_1, \ldots, i_{n-2} = 0}^1 \ket{i_1\ldots i_{n-2}}\bra{i_1\ldots i_{n-2}}\cdot \otimes  R_y^\dagger(\theta_{i_1\ldots i_{n-2}}^{(n-2)})R_y(\tilde{\theta}^{(n-2)})\right)
\nonumber\\
&=\brazero{(n-1)}\mathfrak{U}_{n-1}^\dagger \tilde{\mathfrak{U}}_{n-1}\ketzero{(n-1)}\cdot\cos\left({\eta_n/2}\right).
\end{align}

Finally, using the induction hypothesis, we conclude that

\begin{align}
\brazero{n}\mathfrak{U}_n^\dagger \tilde{\mathfrak{U}}_n\ketzero{n}
&\geq \prod_{k=1}^{n-1}\cos\left({\eta_k/2}\right)\cdot\cos\left({\eta_n/2}\right) = \prod_{k=1}^{n}\cos\left({\eta_k/2}\right) \nonumber\\&\implies F \geq \prod_{k=1}^{n}\cos^{2}\left(\eta_k/2\right).
\end{align}

Notice that one of the conditions of the previous Theorem is that $\eta_k\leq\pi$ $\forall k\in\{1, \dots, n\}$. Recall that, as seen in Eq. (<ref>), $\eta_k\equiv\frac{\delta_k}{8}\eta$. Therefore, we obtain the following condition for the value of $\eta$ for Th. <ref> to be satisfied:

\begin{equation}
\label{eq:eta}
\frac{\delta_k}{8}\eta \leq \pi \implies \eta \leq 2^{k-1}\cdot8\pi.
\end{equation}

In particular, if this holds for the smallest $k$, with value $1$, it will be satisfied for the rest of the blocks. Then, we can write the previous condition simply as $\eta\leq 8\pi$.

§.§.§ Expression of $k_0$

Assume only clusterization for $k > k_0$:

\begin{align}\label{eq:1}
F \geq \prod_{k=k_0+1}^n \cos^2\left(\frac{\eta}{4\cdot2^k}\right) \geq \prod_{k=k_0+1}^n e^{-2\left(\frac{\eta}{4\cdot2^k}\right)^2} = e^{-2\frac{\eta^2}{16}\sum_{k=k_0+1}^n4^{-k}},
\end{align}

since $\cos(x) \geq e^{-x^2}$ for $x\lessapprox\pi/2$. Also, the sum is a geometric series with value

\begin{equation}
\sum_{k=k_0+1}^n4^{-k} = \frac{1}{3}\left(4^{-k_0}-4^{-n}\right).
\end{equation}

\begin{align}
F \geq e^{-\frac{2}{3}\frac{\eta^2}{16}\left(4^{-k_0}-4^{-n}\right)} \coloneqq F_{k_0}.
\end{align}

Let us finally find an expression for $k_0$ given minimum fidelity $F_{k_0}=1-\epsilon$:

\begin{align}
-\frac{2}{3}\frac{\eta^2}{16}\left(4^{-k_0}-4^{-n}\right) = \log{F_{k_0}}\\
\left(4^{-k_0}-4^{-n}\right) = -\frac{3}{2}\frac{16}{\eta^2}\log{F_{k_0}}\\
4^{-k_0} = 4^{-n}-\frac{3}{2}\frac{16}{\eta^2}\log{F_{k_0}}\\
-k_0\log{4} = \log{\left(4^{-n}-\frac{3}{2}\frac{16}{\eta^2}\log{F_{k_0}}\right)}\\
k_0 = -\frac{\log{\left(4^{-n}-\frac{3}{2}\frac{16}{\eta^2}\log{1-\epsilon}\right)}}{\log{4}}.\label{eq:k0s}
\end{align}

Note that in the limit for both $F_0\to 1$ and $\eta\to\infty$, we have that

\begin{equation}
k_0\to -\log(4^{-n})/\log(4)=n,
\end{equation}

as expected.

Also, since $k_0$ must be an integer, we take the ceiling of its value. In addition, we consider the minimum of this parameter to be 2, so the condition of $\eta\leq 8\pi$ is satisfied. Hence, the final expression for $k_0$ is

\begin{equation}
\label{eq:k0_ap2}
k_0 = \max \left \{ \lceil -\frac{1}{2}\log_2 (4^{-n}-\frac{96}{\eta^2}\log (1-\epsilon))\rceil,2 \right \},
\end{equation}

where we have defined $\epsilon \equiv 1-F_0$.

Once we have a fixed number of uniformly controlled rotation blocks, $k_0$, for which there is no clustering, we count the number of necessary gates. The gate $F_{k}^{k-1}(\bm{\hat{y}}, \bm{\theta})$ can be implemented with $2^{k-1}$ CNOTs and $2^{k-1}$ single-qubit gates [42]. Thus, the total number of gates is:

\begin{align}
&\# \text{CNOTs} = \sum_{k=1}^{k_0}2^{k-1}=2^{k_0}-1,\label{eq:cnot-k0}\\
&\# \text{SQGs} \;\;\;= \sum_{k=1}^{k_0}2^{k-1} + \sum_{k=k_0+1}^{n}1 = 2^{k_0}-1 + n-k_0.
\end{align}

§.§ Final Protocol

Let $f\colon [x_\text{min}, x_\text{max}]\to $IR$^{+}$ be a positive function in $L^2([x_\text{min}, x_\text{max}])$. We define the $n$-qubit normalized representative state of $f(x)$ as the $n$-qubit state $|f(x)\rangle_n=\sum_{j=0}^{2^n-1} f(x_\text{min}+j\delta_n)|j\rangle$, with $\delta_n=\frac{x_\text{max}-x_\text{min}}{2^n-1}$ and $\sum_{j=0}^{2^n-1}f^2(x_\text{min}+j\delta_n) = 1$.

Let $f:[0,1]\rightarrow$IR$^+$ be a positive integrable function in $L^2([0,1])$ and $0\leq\eta\leq 8\pi$ a constant such that $\eta=\sup_{x\in[0,1]}\left|\partial^2_x \log f^2(x) \right|$. Then, it is posible to approximate the $n$-qubit representative state of $f(x)$, $|f(x)\rangle_n$, by a quantum state $|\Psi(f)\rangle_n$ such that the fidelity $|\langle \Psi(f)|f(x)\rangle_n|^2 \geq 1-\epsilon$ with at most $2^{k_0(\epsilon)}-1$ two-qubit gates, with

\begin{equation}\label{eq:k0_ap}
k_0(\epsilon) = \max \left \{ \lceil -\frac{1}{2}\log_2 (4^{-n}-\frac{96}{\eta^2}\log (1-\epsilon))\rceil,2 \right \},
\end{equation}

and the circuit to perform it is given in Fig. <ref>.

Quantum circuit performing the protocol presented in this section, based on the Grover-Rudolph method for a system of $n$ qubits.

Using Th. <ref> and Cor. <ref>, we can define the representative angle for each block as in Eq. (<ref>). Then, by applying Th. <ref> with the development between Eqs. (<ref>) and (<ref>), it is clear that the fidelity will be larger than or equal to $1-\epsilon$. Finally, Eq. (<ref>) states that the number of two-qubit gates necessary to realize this protocol is $2^{k_0}-1$, as we wanted to proof.

§ C  THEOREM GENERALIZATION: SINGULAR POINTS

Consider a function $f:[0,1]\rightarrow\mathbb{R}^+$ with

\begin{equation}
\left|\partial^2_{x}\log f(x)\right|\xrightarrow[x \to 0]{} \infty.
\end{equation}

An example is the beta density function, defined as

\begin{equation}
f(x)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)},
\end{equation}

with $\alpha, \beta >0$ and $B(\alpha, \beta)$ the beta function. Then, the second derivative of its logarithm is

\begin{equation}
\partial^2_{x}\log f(x)=\frac{1-\alpha}{x^2}+\frac{1-\beta}{(1-x)^2}.
\end{equation}

Thus, for $\alpha\neq 1$, we have a singularity at $x=0$. Also, we see that if $\beta\neq 1$, there is a singularity at $x=1$.

In this situation, $\left|\partial^2_{x}\log f(x)\right|$ can not be bounded by a finite factor $\eta$ in the whole interval and, hence, Th. <ref> can not be applied. However, under certain circumstances, this issue can be solved.

Recall that, in Th. <ref>, we obtained that the difference between two consecutive angles in a block of $k$ qubits satisfies

\begin{equation}
\left|\theta^{(k-1)}_{l+1}-\theta^{(k-1)}_l\right| \leq \frac{\delta_k^2}{4}\left|\partial^2_{x}\log f(x)\right|,
\end{equation}

with $\delta_k=\frac{1}{2^{k-1}}$, $l\in\{0, 1, \dots, 2^{k-1}-2\}$, and $x=l\delta_k$. Since the singularity is found in $x=0$, there exists a $k_{\text{max}}$ from which the maximum of $\left|\partial^2_{x}\log f(x)\right|$ in $[2^{-k_{\text{max}}+1}, 1]$ is found in $x=2^{-k_{\text{max}}+1}$. Then,

\begin{equation}
\left|\theta^{(k_{\text{max}}-1)}_{l+1}-\theta^{(k_{\text{max}}-1)}_l\right| \leq \frac{1}{4}\frac{\left|\partial^2_{x}\log f(x)|_{x=2^{-k_{\text{max}}+1}}\right|}{2^{2k_{\text{max}}-2}}.
\end{equation}

Next, if for the limit $k\to\infty$ the term $\left|\partial^2_{x}\log f(x)|_{x=2^{-k+1}}\right|\cdot2^{-2k+2}\to 0$, this value is decreasing and there exists a $k^*\geq k_{\text{max}}$ from which an $\eta=\left|\partial^2_{x}\log f(x)|_{x=2^{-k^*+1}}\right|$ can be set so the clustering conditions from Eq. (<ref>) are satisfied. This $k^*$ acts as the $k_0$ of Th. <ref> and represents the last block without clustering. Then, for $k>k^*$, the angles corresponding to the interval $[2^{-k+1}, 1]$ can be clustered, following our protocol. The case for a singularity in $x=1$ is analogous, and in Fig. <ref> it is depicted the resulting gates corresponding to the process of clustering the inner angles for a block of 3 qubits. Given that $R_y(\tilde{\theta})\cdot R_y(-\tilde{\theta})=\mathbbm{1}$, we can complete the identity of the clusterized block by subtracting $\tilde{\theta}$ from the rest of the angles.

[column sep=0.2cm]
1      1     1     1

1      1     1     1

R_y(θ^(2)_00)     R_y(θ̃)     R_y(θ̃)     R_y(θ^(2)_11)

[column sep=0.2cm]
1     1     1[3,steps=4,style=dashed,rounded corners,fill=violet!70!blue!20!, inner xsep=2pt,background]$\mathbbm{1}\otimes R_y(\tilde{\theta})$      1     1     1     1     1

1     1     1      1     1     1     1     1

R_y(θ^(2)_00)     R_y(-θ̃)     R_y(θ̃)     R_y(θ̃)     R_y(θ̃)     R_y(θ̃)     R_y(-θ̃)     R_y(θ^(2)_11)
[column sep=0.2cm]
1           1

1          1

R_y(θ^(2)_00-θ̃)     R_y(θ̃)     R_y(θ^(2)_11-θ̃)
Example of the reduction of gates for the block of three qubits $F_{3}^2(\bm{\hat{y}}, \bm{\theta}^{(2)})$, where it is assumed that both $\theta^{(2)}_{01}$ and $\theta^{(2)}_{10}$ can be approximated by $\tilde{\theta}$.

In this situation, if we consider that singularities exist in $x=\{0, 1\}$, the total number of two-qubit gates required to load $f$ into a quantum state is

\begin{align}
\# \text{TQGs} &= \sum_{k=1}^{k^*}2^{k-1}+2\sum_{k=k^*+1}^{n}(80k-398)\nonumber\\&=2^{k^*}-1+2(n-k^*)(20k^*+20n-179),
\end{align}

where we have used the result of Cor. <ref>. Notice that this result is only valid for $k^*>6$, but since we are interested in the asymptotical behaviour of the protocol, it is not an issue. All in all, we obtain that the complexity of this process is exponentially dependent on $k^*$ with an extra polynomial term.

On the other hand, if a singularity is found in $(0, 1)$, we end up with multiple clusters of angles in each block. The reason behind this is that the clusters are formed with contiguous angles. Therefore, the reduction of gates is not significant, and the protocol can not be performed efficiently. However, if the representatives of the disjointed clusters are equal, then we can create a single cluster, so the reduction is doable. In this sense, we have numerically observed that far from the singularities, all angles converge to a value of $\pi/2$, but we have not found an analytical proof yet.

Let us see an example of a function that meets the previous description and analyze the outcome of the protocol. Consider the function

\begin{equation}\label{eq:poli}
\end{equation}

in $[0, 1]$, where $N$ is the normalization factor, and a system of $n=10$ qubits. The second derivative of its logarithm is

\begin{align}
\log f(x) = x^{\frac{3}{2}} -\log N
&\implies \partial_{x}\log f(x) = \frac{3}{2}\sqrt{x}\\
&\implies \partial^2_{x}\log f(x) = \frac{3}{4\sqrt{x}}.
\end{align}

Hence, $\partial^2_{x}\log f(x)$ has a singularity at $x=0$.
First of all, since $\partial^2_{x}\log f(x)$ is a monotonic decreasing function, its maximum in the interval $[2^{-k+1}, 1]$ is found in $x=2^{-k+1}$, for any $k=1, \dots, n$. Then, we can set $k_{\text{max}}=1$. Next, we need to compute the limit of $\left|\partial^2_{x}\log f(x)|_{x=2^{-k+1}}\right|\cdot2^{-2k+2}$:

\begin{align}
&\left|\partial^2_{x}\log f(x)|_{x=2^{-k+1}}\right|\cdot2^{-2k+2} = \frac{3}{4\sqrt{x}}|_{x=2^{-k+1}}\cdot2^{-2k+2}\nonumber\\
&\quad=\frac{3}{4} 2^{\frac{1}{2}(k-1)}2^{-2k+2}=\frac{3}{4} 2^{-\frac{3}{2}k+\frac{3}{2}}\xrightarrow[k \to \infty]{} 0.
\end{align}

Therefore, the difference in the angles is bounded and we can select a $k^*\geq 1$ for which the conditions of Theorem $\ref{th:k0}$ are satisfied in $[2^{-k^*+1}, 1]$. The first condition we need to check is the inequality given by

\begin{equation}
\frac{\delta_k}{8}\eta \leq \pi,
\end{equation}

with $\eta=\left|\partial^2_{x}\log f(x)|_{x=2^{-k+1}}\right|=\frac{3}{4}2^{\frac{1}{2}(k-1)}$. Then,

\begin{align}
\frac{3}{32}2^{\frac{1}{2}(k-1)}2^{-k+1}=\frac{3}{32}2^{\frac{1}{2}(1-k)}.
\end{align}

Since this term is decreasing, its maximum is found when $k=1$, with a value of $\frac{3}{32}<\pi$. Thus, this condition is met for any $k=1, \dots, n$, so we can select $k^*=1$. Additionally, if we compute the value of $k_0$ with $\epsilon=0.01$ and $\eta=\left|\partial^2_{x}\log f(x)|_{x=2^{-k^*+1}}\right|=0.75$ using Eq. (<ref>), we obtain $k_0=2$. Then, the last block to remain unclustered must be the maximum between $k_0$ and $k^*$, which in this case is 2.

Now, in Fig. <ref>, we have depicted the result of the experiment of the considered function for a system of $n=10$ qubits and the clustering starting with the 3-qubit block, following the protocol described in this section. With a fidelity of $0.99975$, larger than the one required, we have that the final state $\ket{\Psi(f)}_{10}$ successfully captures the features of $f$ with a reduction of the complexity of two-qubit gates from $\mathcal{O}(2^{10})$ to $\mathcal{O}(2^{2})$.

Simulation of $\ket{f(x)}_{10}$ and $\ket{\Psi(f)}_{10}$ for the function defined in Eq. (<ref>), $n=10$, $\epsilon=0.01$ and the clustering process performed for $k=3, \dots, n$.

§ D  VARIATIONAL CIRCUIT

§.§ Ansatz

The circuit proposed as an ansatz, $U(\bm{\theta})$, for our model is based on the Grover-Rudolph circuit [28]. We consider as hyperparameters the number of angles which would remain unclustered in each uniformly controlled block $F_{k}^{k-1}(\bm{\hat{y}}, \bm{\theta})$, corresponding to the positions of the zeros and singular points in the target function to be loaded, see Fig. <ref>. Initially, one/two angle (depending on the slope in the contiguous intervals) per singular point/zero are assumed, see Fig. <ref> . In case more precision is required, we simply add additional angles by considering the nearest angles from both left and right intervals, if possible. In addition, as we considered with our first protocol, we assume the first $k_0$ blocks are not clustered. This parameter also acts as a hyperparameter.In Fig. <ref>, we illustrate how to select the hyperparameters of the variational circuit for a general function containing multiple zeros and singular points.

Let us now count the necessary two-qubit gates to implement such model. We assume that the desired function $f$ has $z$ zeros, $s$ singular points, and the system has $n$ qubits. Then, if we consider $p(k)$ non-clustered angles for each zero and singular point, for $k=k_0+1, \dots, n$, the necessary two-qubit gates are:

\begin{align}
&\# \text{TQGs} = \sum_{k=1}^{k_0}2^{k-1}+(s+z)\sum_{k=k_0+1}^{n}p(k)(80k-398)=2^{k_0}-1+r_{k_0+1}(n, s+z),
\end{align}

where $r_{k_0+1}(n, s+z)$ is a polynomial number in the size of the system $n$. As we can see, the complexity of the algorithm is dominated by the exponential dependency of the hyperparameter $k_0$.

Illustration of proposed variational circuit for a general function. (a) Target function. The proposed landscape generalizes the possible cases of special points that lead to singularities in the second logarithmic derivative. The special points are classified in the legend of the graph, where $z$ denotes the zeros, $s$ the singularities, and the index denotes whether the slope in the contiguous points has the same sign, index$=1$, or not, index$=2$. This is translated to a minimum amount of 7 necessary parameters required to capture the local behaviour of the function in these points, one per index$=1$ and 2 per index$=2$. We add an additional parameter which represents the clusterization. (b) Variational circuit. The blocks which have a number of angles lower or equal than the necessary parameters will remain unclustered. From the first block which has more angles, we proceed to reduce the number of angles to the number of parameters by clustering the ones that do not correspond to the special points. This leads to a significant reduction in the number of gates. (c) Zoom in of the  $F_5^{4}(y, \bm{\theta})$ block clusterization. As we only need $8$ parameters, the $16$ angles of this block are reduced to the half. We keep one local rotation for the clustering parameter and 7 multi-controlled gates corresponding to the intervals where the singularity is located.

§.§ Training Process

We proceed to illustrate the training method for the variational circuit proposed in the Letter. This process consists of iterative steps in which a loss function that measures how far the outcome is from the desired state is recursively minimized to obtain the optimal parameters.

Consider a function $f$ and $n$ qubits. The desired quantum state is

\begin{equation}
\ket{\Psi(f)}_n = \sum_{i=0}^{2^n-1}f_i\ket{i},
\end{equation}

where the $f_i$ terms are discrete approximations to the objective function, as presented in Def. <ref>.
Now, given a set of angles $\bm{\theta}$, our ansatz $U(\bm{\theta})$ returns the state

\begin{equation}
\ket{\Psi\left(\bm{\theta}\right)}_n=U(\bm{\theta})\ket{0},
\end{equation}

where $\ket{0}\equiv\ketzero{n}$, and the components of the obtained state are products of sines and cosines.

Let us now introduce the mean squared error loss function, defined as

\begin{align}
L(n, f, \bm{\theta}) &= \frac{1}{2^n}\sum_{i=0}^{2^n-1}\left(f_i-\Psi_i\left(\bm{\theta}\right)\right)^2=\frac{1}{2^n}\sum_{i=0}^{2^n-1}\left(f_i^2+\Psi_i\left(\bm{\theta}\right)^2-2f_i\Psi_i\left(\bm{\theta}\right)\right),
\end{align}

where $\Psi_i\left(\bm{\theta}\right):=\braket{i}{\Psi({\bm{\theta}})}_n$.

This process aims to find the optimal parameters $\bm{\theta}$ for which $\ket{\Psi\left(\bm{\theta}\right)}_n$ approximates $\ket{\Psi(f)}_n$, which is equivalent to minimizing the loss function. Here, we use the gradient descent method [52] to do so.

Then, given any angle $\theta_{j}^{(k-1)}$, with $k\in\{1, \dots, n\}$ and $j\in\{0, 2^{k-1}-1\}$, its value gets updated after each training step in the following way:

\begin{equation}
\theta_{j}^{(k-1)} = \theta_{j}^{(k-1)} - \gamma\frac{\partial L(n, f, \bm{\theta})}{\partial \theta_{j}^{(k-1)}},
\end{equation}

where $\gamma$ is the learning rate. Let us now compute the expression of the derivative:

\begin{align}
\frac{\partial L(n, f, \bm{\theta})}{\partial \theta_{j}^{(k-1)}}=\frac{1}{2^n}\sum_{i=0}^{2^n-1}2\left(\Psi_i({\bm{\theta}})-f_i\right)\frac{\partial\Psi_i(\bm{\theta})}{\partial \theta_{j}^{(k-1)}}.
\end{align}

Since $\Psi_i(\bm{\theta})$ is a product of sines and cosines, its partial derivative with respect to a given angle is

\begin{align}
\frac{\partial\Psi_i(\bm{\theta})}{\partial \theta_{j}^{(k-1)}}=\begin{cases}
-\frac{1}{2}\frac{\sin\left(\theta_{j}^{(k-1)}/{2}\right)}{\cos\left(\theta_{j}^{(k-1)}/{2}\right)}\Psi_i(\bm{\theta}) ~&\text{if}~j\cdot2^{n-k+1}\leq i< (j+1)\cdot2^{n-k},\\
\frac{1}{2}\frac{\cos\left(\theta_{j}^{(k-1)}/{2}\right)}{\sin\left(\theta_{j}^{(k-1)}/{2}\right)}\Psi_i(\bm{\theta}) ~&\text{if}~(j+1)\cdot2^{n-k}\leq i< (j+1)\cdot2^{n-k+1},\\
0 ~&\text{otherwise}.
\end{cases}
\end{align}

Finally, different stopping criteria exist to put to an end the training process. A simple example is to fix a number of training steps. Another criteria, which is the one we consider for our numerical simulations, is to set a tolerance for the loss function. Therefore, when the difference between the cost function of two consecutive steps is less than the given tolerance, the process is considered satisfactory.

As a summary, in Algorithm 1 we have depicted the pseudocode corresponding to the training process.

§.§ Black–Scholes Distribution

Let us now apply the proposed variational circuit to the Black-Scholes distribution [53], which is given by

\begin{align}
\text{BS}(x)=\begin{cases}
K-e^{-x}/s ~&\text{if}~-\log{(K\cdot s)}\leq x< 0,\\
K-e^{x}/s ~&\text{if}~0<x\leq \log{(K\cdot s)},
\label{eq:bs}
\end{cases}
\end{align}

with $s=K\cdot c$, where $K$ and $c$ are the parameters. The zeros are found in $x=\pm\log{(K\cdot s)}$.

First, in Fig. <ref>, we see the resulting circuit for the block of the uniformly-controlled rotation composed by 4 qubits and $p(k)=1$. As described previously, only the multi-controlled operations corresponding to the intervals containing zeros remain, the rest is clustered.

Variational circuit corresponding to the 4-qubit block and $p(k)=1$ for the Black-Scholes distribution given by Eq.(<ref>), with $K=45$ and $c=3$. The function has two zeros of type 1 ($z1$) located in $x=\pm\log{(K\cdot s)}$, consequently we can chose our hyperparameter $k_0$ to be larger or equal than 2. Therefore, as illustrated, the clustered blocks will present two multi-controlled rotations corresponding to each zero and one single qubit rotation with the global clusterization angle.

We benchmark the performance of our ansatz for loading the Black-Scholes distribution, Eq.(<ref>), into a $5$-qubit system. We choose the parameters of the distribution $K=45$ and $c=3$, and for the training parameters, we consider, $k_0=2$, a learning rate of $\gamma=1.5$ and a tolerance of $10^{-9}$. In Fig. <ref>, we depict the results of the from the numerical simulations of the discretized state $\ket{\Psi\left(\text{BS}(x)\right)}_5$ and the trained ansatz $\ket{\Psi\left(\bm{\theta}\right)}_5$, for different values of $p(k)$ and $n=5$ qubits. We can appreciate that in all cases the circuit is able to capture the main features of the target state, consequently all fidelities are close to 1. We highlight that for $p(k)=3$ the resulting state and the objective one are almost equal, $\epsilon=5\mathrm{e}{-5}$.

We also provide an analysis for larger sytems. In Table <ref>, we depict the numerical results for $\ket{\Psi\left(\text{BS}(x)\right)}_{12}$ and $\ket{\Psi\left(\bm{\theta}\right)}_{12}$ in a system of 12-qubit system, with $k_0=2$, and different values of $p(k)$. We are able to obtain large fidelities with only a few percentage of the original angles, which implies that with only a small portion of the initial parameter space, the circuit can approximately simulate the target state efficiently.

Simulations of $\ket{\Psi\left(\text{BS}(x)\right)}_5$ and $\ket{\Psi\left(\bm{\theta}\right)}_5$ for the Black-Scholes distribution, $n=5$, $k_0=2$, and different values of $p(k)$: (a) $p(k)=1$, (b) $p(k)=2$, and (c) $p(k)=3$. Parameters of the distribution: $K=45$ and $c=3$. Training parameters: $ k_0=2$, learning rate $\gamma=1.5$ and  tolerance of $10^{-9}$. We can appreciate that in all cases the circuit is able to capture the main features of the target state, consequently all fidelities are close to 1.

Numerical data for the simulations of $\ket{\Psi\left(\text{BS}(x)\right)}_{12}$ and $\ket{\Psi\left(\bm{\theta}\right)}_{12}$, $n=12$, and different values of $p(k)$, with $k_0=2$. The value #Angles is the number of independent angles necessary to generate $\ket{\Psi\left(\bm{\theta}\right)}_{12}$, and %Angles is in comparison with the full circuit. We are able to obtain large fidelities with only a few percentage of the original angles, which implies that with only a small portion of the initial parameter space, the circuit can approximately simulate the target state efficiently.

$p(k)$   Fidelity   #Angles   %Angles

1     0.99303      33     0.806

2     0.99838      52     1.269

3     0.99890     70     1.709

$k$     0.99913     142     3.468

Finally, in Fig <ref> and Table <ref> we analyze the performance of initializing the value of the parameters with the angles provided by the Grover-Rudolph method for different sizes of the system. We present the training results of the ansatz initialized with several ramdom angles sets and with the angles provided by the Grover-Rudolph method. We analyze the performance in terms of infidelity and number of steps of the training. We can observe how this initialization reduces drastically the number of steps, as well as enables us to achieve the largest fidelities.

Results of training the variational ansätze for a different number of qubits (n=15,16,17,18) to load the Black-Scholes distribution. We compare the training results (numbers of step and final infidelity) of 40 randomly initialized parameters versus the initialization with Grover-Rudolph method angles. These randomly initialized parameters are generated by drawing samples from a  uniform distribution between the values of 0 and $\pi$. As we can appreciate the initialization with Grover-Rudolph method always achieves the best performance, while most of random initializations get stuck around a fidelity value (vertical lines). Training parameters: $ k_0=2$, $p(k)=1$,  learning rate $\gamma=1.5$ and  tolerance of $10^{-9}$.

Numerical results of training the variational ansätze for a different number of qubits (n=15,16,17,18) to load the Black-Scholes distribution. We compare the training results (numbers of step and final infidelity) of 40 randomly initialized parameters versus the initialization with Grover-Rudolph method angles. These randomly initialized parameters are generated by drawing samples from a uniform distribution between the values of 0 and $\pi$. Training parameters: $ k_0=2$, $p(k)=1$, learning rate $\gamma=1.5$ and  tolerance of $10^{-9}$.

$n$   Fidelity (GR)   Steps (GR)    Avg. Fidelity (Random)    Avg. Steps (Random)     Max Fidelity (Random)    Dif. Fid. GR vs Avg. Random

15     0.99317      13     0.99261     82     0.99295     0.00056

16     0.99316      13     0.99211     59    0.99271    0.00105

17     0.99314     14     0.99052     57    0.99145    0.00262

18     0.99309     13     0.98566    31    0.99199    0.00743