Efficient quantum amplitude encoding of polynomial functions

Loading functions into quantum computers represents an essential step in several quantum algorithms, such as quantum partial differential equation solvers. Therefore, the inefficiency of this process leads to a major bottleneck for the application of these algorithms. Here, we present and compare two efficient methods for the amplitude encoding of real polynomial functions on $n$ qubits. This case holds special relevance, as any continuous function on a closed interval can be uniformly approximated with arbitrary precision by a polynomial function. The first approach relies on the matrix product state representation. We study and benchmark the approximations of the target state when the bond dimension is assumed to be small. The second algorithm combines two subroutines. Initially we encode the linear function into the quantum registers with a shallow sequence of multi-controlled gates that loads the linear function's Hadamard-Walsh series, exploring how truncating the Hadamard-Walsh series of the linear function affects the final fidelity. Applying the inverse discrete Hadamard-Walsh transform transforms the series coefficients into an amplitude encoding of the linear function. Then, we use this construction as a building block to achieve a block encoding of the amplitudes corresponding to the linear function on $k_0$ qubits and apply the quantum singular value transformation that implements a polynomial transformation to the block encoding of the amplitudes. This unitary together with the Amplitude Amplification algorithm will enable us to prepare the quantum state that encodes the polynomial function on $k_0$ qubits. Finally we pad $n-k_0$ qubits to generate an approximated encoding of the polynomial on $n$ qubits, analyzing the error depending on $k_0$. In this regard, our methodology proposes a method to improve the state-of-the-art complexity by introducing controllable errors.

There exists a necessity for loading polynomial functions due to the growing number of applications of quantum computing.Specifically, in finance [10][11][12][13][14][15][16][17], the ability to efficiently load first order polynomials, f (j) = aj + b allows for options pricing via quantum amplitude estimation (QAE) without coherent arithmetic [70].In this sense, applying QAE allows us to extract the amplitude j f (j)p(j) = E[f (X)] where X is a random variable with probability mass function p(j).This approach can be generalized to the multidimensional setting where the ability to load multivariate linear functions yields efficient algorithms for pricing basket options [71].Furthermore, quantum circuits for efficiently loading the linear function can be used to construct a block encoding of the identity function thus allowing us to apply the quantum singular value Transformation [72][73][74][75][76][77] (QSVT) in order to obtain any polynomial amplitude encoding [41,42,78,79].This is a powerful tool capable of uniformly approximating the encoding of any continuous real function defined on a closed interval with arbitrary precision [80,81].
In this article we present two methods for implementing the amplitude encoding of real valued polynomials into quantum computers with linear complexity.The first one is based on the matrix product sate (MPS) representation of quantum states and its implementation on a quantum computer.We explore how approximating the MPS affects the achieved fidelity and the resources requirements [82][83][84][85][86][87].On the other hand, the second method consists of two steps, first we propose a novel protocol to efficiently load the linear function on k 0 qubits based of the discrete Hadamard-Walsh transform (DHWT) [88] that we use to achieve a block encoding of the amplitudes with complexity O(k 0 ).Secondly, we use the QSVT to implement a polynomial transformation on the eigenvalues of the block encoding to achieve the desired target state [78,79].Note that as we have implemented a block encoding of the linear function, our method avoids errors due to the polynomial approximation of the arcsin(x) as proposed in previous works [41,42].Furthermore, we reduce the complexity of the state-of-the-art for encoding polynomials functions via QSVT by introducing a controllable error.
The article is structured as follows, first we review the loading of polynomials via MPS.Next, we introduce our new methodology that combines DHWT to load the linear function with the QSVT to implement the polynomial transformation of the amplitudes.Finally we show numerical results and compare our method with previous results in the literature.

Loading of polynomials via MPS
In this section we analyze the methods based on MPSs to encode polynomials into the amplitude of a quantum state according to following definition.
Definition 1 Let P (x) : [a, b] → C, be a polynomial with complex coefficients.We define the nqubit normalized representative state of P (x) as the quantum state |Φ P ⟩ = 1 C P 2 n −1 j=0 P (x j ) |j⟩, with x j = a + j b−a 2 n −1 and C P the normalization factor.
In the particular case of the linear function, i.e.P (x) = x, defined on [0, 2 n − 1], we have the normalization factor is C P = (2 n+1 − 1)2 n /(6(2 n − 1)).For now on for simplicity, we will assume that the polynomials we work with are defined on the interval [0, 1], so that might include rescaling of the domain of the linear function.
The complete description of a quantum state of n linearly connected qubits (sites) can be represented by a tensor, A, with n physical indices.A state of this kind is referred to as a matrix product state (MPS) [82][83][84][85].Each physical index is assigned to a qubit and has a degree of d = 2 i.e., the index is either 0 or 1.For a specific choice of each physical index j i , the tensor's values give rise to a collection of n matrices whose product is equal to the amplitudes of the computational basis state |j n−1 . . .j 0 ⟩.For open boundary conditions, the MPS representation of a quantum state is given by This representation has n − 2 tensors of order 3, denoted as

and 2 external tensors A
[1] j 0 ,α 1 and A [n] j n−1 ,α n−1 of order 2, with j i the physical indexes that range from 0 to d − 1 and α i the virtual indices from 0 to χ i − 1.The virtual dimensions, χ i , connecting each pair of tensors via the virtual indices are referred as the bond dimensions.We define the bond dimension of the entire MPS as χ = max i χ i .
When it comes to representing polynomials as quantum states, Grasedyck [60] proved that for a real valued polynomial of degree d encoded in the amplitudes of a quantum state according to Def. 1, the MPS bond dimension is as much χ = d+1.

Obtaining the exact MPS
The MPS representation of a quantum state |Ψ⟩ is not unique, as different choices of A [i] j i ,α i−1 α i can yield the same quantum state.We focus on the left canonical form, which implies the following conditions To obtain the MPS representation of a quantum state, the singular value decomposition (SVD) is employed [89].Initially, the quantum state, represented as a tensor A of rank n and dimension 2, is reshaped into a matrix by combining all the indices except one.The SVD is then applied to this matrix, decomposing it into the matrix of left singular vectors, U , the matrix of singular values, Σ, and the matrix of right singular vectors, V † .As we will depict in the following section, it is possible to truncate the smallest singular values by choosing a desired bond dimension χ and keeping only the χ largest values.Then, the matrix of singular values is contracted with the matrix of right singular vectors (left canonical form), and the resulting matrix, ΣV † , is reshaped back into a tensor that now has an extra virtual index.This process, shown in Fig. 1, is iterated for each physical index.Finally, we obtain an MPS that approximates the original quantum state while providing a compact and efficient representation.The computational cost of performing SVD on an m × n matrix is typically O min(mn 2 , m 2 n) and must be taken into account as a preprossessing cost in this methodology.This cost can be reduced for sparse or structured low-rank matrices.In the case of an exact matrix product state (MPS), the bond dimension doubles for each connection between core tensors.This leads to a maximum bond dimension of 2 ⌊n/2⌋ , occurring in the middle of the MPS.As a result, the computational cost of the entire algorithm is primarily determined by the SVD of this central square matrix, i.e.O(2 3n/2 ), which exhibits exponential scaling with the number of qubits.Therefore, this non-negligible classical pre-processing cost must be considered in the overall complexity of the MPS algorithm.
In the particular case of the linear function, the analytical expression of the exact MPS [90], which has bond dimension χ = 2, reads

Approximation of the exact MPS
While, in the worst case scenario, it is possible to exactly represent any quantum state as an MPS by allowing the bond dimensions to grow up to Figure 1: Iterative singular value decomposition (SVD) [89] procedure for obtaining the MPS from a tensor with n = 3 physical indices.The process involves n−1 uses of the SVD.Notably, the matrices containing the singular values are absorbed to the right, resulting in the left canonical form of the MPS.
2 ⌊n/2⌋ , we can potentially achieve an exponential compression by approximating the initial tensor using O(2nχ 2 ) elements, given a fixed bond dimension, χ i = χ ∀ i.However, as the entanglement of the state to be approximated increases, the minimum bond dimension χ required to get a good description of the state using an approximated MPS representation also grows [91,92].Notice that when truncating the maximum bond dimension of the MPS, the complexity of the classical preprocessing becomes O(nχ 3 ).
In order to estimate the error incurred when truncating the bond dimension, it is necessary to consider that during each iteration of the compression protocol we perform a SVD and discard singular values.The error incurred when approximating a matrix by considering its k largest singular values is determined by the Eckart-Young theorem [93], and corresponds to the sum of the omitted singular values.During compression, this rank-k approximation is performed for each core tensor.Thus, the overall error of the MPS approximation can be upper bounded in the Frobenius norm as where Ã denotes the approximated MPS.This equation encompasses the error contributions from the n − 1 singular value matrices Σ i , each characterized by a fixed bond dimension, χ i .In order to keep the approximation error low, it is crucial that the spectrum of each Σ i decays rapidly.
In this regard, the state-of-the-art technique for preparing smooth, differentiable, real-valued functions using matrix product states relies on singular values exhibiting exponential decay as demonstrated in Ref. [61].This allows for good approximations of such a class of functions while maintaining low bond dimension values, thus significantly reducing the required resources to prepare the associated function-encoding quantum state.The method relies on the fact that, for such functions, the entanglement entropy, quantified by the von Neumann entropy, scales logarithmically with the number of qubits, n, and therefore, these functions can be efficiently approximated by an MPS with a low bond dimension, as argued in Ref. [91].Although empirical results with this technique applied to polynomials show good performance with χ = 2, the upper bound of the von Neumann entropy depends on the maximum derivative value of the function within the considered interval.Consequently, as we discuss later in this work, for certain polynomials the truncation to χ = 2 does not yield a satisfactory approximation.
In the particular case of the linear function, the analytical expression from Eq. 2.1 reveals that achieving the exact MPS representation requires a bond dimension of χ = 2.However, one might consider reducing the bond dimension to 1 to investigate how severely accuracy decays.We explore this possibility by comparing the linear function approximated by MPSs of bond dimension χ = 2 (exact MPS) and χ = 1 (product state) in Section 4.

From MPS to circuit
Let us now analyze the resources needed to translate an MPS in either cases, exact or approximated, into a quantum circuit [85,86,94,95].First we will assume that the bond dimensions are powers of two, padding the tensors with zeros if needed.Therefore, we can assume that these tensors are isometries of 2χ i × χ i+1 and thereby can be embedded into a unitary gate acting on n χ i = max{log 2 (2χ i ), log 2 (χ i+1 )} qubits.The factor of two arises from the two possible values that each of the physical indices can take.The arrangement of gates in the MPS conversion process, following a linear topology, results in the formation of a single layer of multi-qubit unitaries, whose sizes n χ i depends on the associated bond dimensions.These unitaries are organized in a staircase topology, commonly referred to as a linear circuit layer [94].Therefore the complexity is equivalent to implementing a cascade of n − 1 multi-qubit unitaries.In the case χ = 2, this complexity scales as O(n) two-qubit unitaries.On top of this, the cost of decomposing these unitaries into two-qubit gates must be taken into account.This cost is considerably larger when decomposing multi-qubit unitaries and might result in an exponential number of two-qubit gates.[96,97].Additionally, one might consider implementing circuits that approximate the exact multi-qubit unitaries [94].Lastly, in Ref. [87] authors proposed a method for loading translation-invariant short-correlation MPS with an error ϵ in depth O(log(N/ϵ)) thereby establishing a class of MPSs that admit efficient circuit implementations.
We can conclude that even though MPSs encoding polynomials admit an analytical construction, we have to consider the cost of obtaining the SVD and the cost of the decomposition of the associated unitaries.These costs can be significantly reduced by truncating down to a bond dimension of χ = 2, although this introduces an uncontrollable source of error [41,61,62,98].Additionally, one can also consider the possibility of approximating the multi-qubit unitaries of a MPS's circuit representation [94], which could result in high fidelity states in some cases although there is not a priori guarantee of success.

Efficient loading of polynomials via DHWT and QSVT
In this section we present a method for loading polynomials by combining a technique to encode the linear function into the amplitudes of a quantum state via the discrete Hadamard-Walsh transform (DHWT) with the quantum Singular Value Transformation (QSVT) algorithm.The first methodology encodes the linear function with a shallow sequence of multi-controlled gates that loads its Hadamard-Walsh series expansion, followed by the inverse discrete Hadamard-Walsh transform.By truncating the series expansion, we demonstrate that this approach allows for a controllable approximation of the state representing the linear function up to a certain error.We analyze this error in terms of infidelity, ϵ with respect to the exact state and through the deviation δ (2) j in each amplitude from the exact amplitudes.The second step leverages the results of [78,79] to generate a block encoding of the linear function.This involves using the unitary to load the linear function into amplitudes of a k 0 -qubit quantum state.Then, the Quantum Singular Value Transformation (QSVT) is applied to implement a polynomial transformation, P (x), on the amplitudes of the quantum state corresponding to the identity function.This results into a unitary block encoding of the polynomial that applied to an adequate initial state yields an encoding of the desired polynomial.Finally, we pad the remaining n−k 0 qubits, obtaining an approximated encoding.

The Discrete Hadamard-Walsh transform
), ( 7) where n−1 m=0 jmkm is the k-th Walsh function, where we have used the natural order.
We also define the binary norm of an integer as |k| b = n−1 m=0 k m .Note that when |k| b = 1, it is equivalent to say that k is a power of 2. Additionally, when representing a quantum state |j⟩ in terms of a binary notation we will denote the state as |j n−1 . . .j 0 ⟩, taking the order of the tensor product from right to left.
Lemma 1 Let (0, 1 . . ., 2 n − 1) be the discrete sorted input sequence.Then, the coefficients of its Hadamard-Walsh series are given by Note that in general, the sparsity of the state encoding the Hadamard-Walsh series of a polynomial of degree d represented in n qubits, with d ≤ n, is s = d k=0 n k [99].Therefore, one might consider the techniques in Ref. [54] to implement the state with depth O(log(ns)) and O(ns log(s)) ancillary qubits.3, where the angle of every multi-controlled rotation is given by θ k = arcsin 2x

Exact loading of the linear function via DHWT
We denote the unitary corresponding to this circuit as U L,n .
y X p z X n 2 j B y f 8 c w h 8 5 H 9 / 0 J Y s 2 < / l a t e x i t > q 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " y W u N z 9 x x n y V T Y w q u m i g < l a t e x i t s h a 1 _ b a s e 6 4 = " F 1 e 5 a H j R 6 < l a t e x i t s h a 1 _ b a s e 6 4 = " n I q 2 a y 0 r q 2 p L r Y B Truncation up to largest k 0 terms < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 p R N e M r / q r E 1 6 w l A T 1 c 8 5 r 9 d t 6 p X G V P 6 J I j s g x q R K H X J A G u S F N 0 i K c P J B n 8 k b e r c R 6 s l 6 s 1 x 9 r w c p n D s k f W B / f s + i P F g = = < / l a t e x i t >

Ry(✓0)
< l a t e x i t s h a 1 _ b a s e 6 4 = " r z c n J U P Z J 5 3 5 1 / d p G 0 T m r u e e 3 s 9 r R S v 8 q L K K J 9 d I C q y E U X q I 5 u U A M 1 E U U p e k Z v 6 N 1 6 t J 6 s F + v 1 x 1 q w 8 p s y + g P r 4 x u z v p H h < / l a t e x i t > Ry(✓n k0 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " r z c n J U P

r t t l a p X + W H K M I R H M M p e H A B d b i B B j S B g Y R n e I N 3 R z h P z o v z + h M t O P m f Q / g j 5 + M b m r G N 8 w = = < / l a t e x i t > |0i
< l a t e x i t s h a 1 _ b a s e 6 4 = " C H R u X p h 4 c q s w h s x a h u b 0 1 k c h W 5 D 5 1 a 7 I r 7 L N X u A e q q 0 8 q P i a q Y 6 E w J 6 U 3 L I 7   The angle of every multi-controlled rotation is given by θ k = arcsin 2x The first part of the circuit loads the state encoding its Hadamard-Walsh series ΦL n and once it has been loaded, we apply the Hadamard-Walsh transform to achieve |Φ L ⟩ n .If we consider only the first k 0 rotations, then we get the approximated Hadamard-Walsh series Φ(k0:n) L n and the respectively approximated state . With this truncation the angles change its value to θ k0 k = arcsin(2x the new normalization factor and G and the complexity is reduced to O(k 0 ).We denote this last circuit that loads the approximated state as U (k0:n) L,n .See the appendix in Ref. [55] for the details of the decomposition of multicontrolled gates.
Once that we have encoded the state ΦL n that represents the discrete Hadamard-Walsh series of our target state, we can simply uncompute the Hadamard-Walsh transform to obtain |Φ L ⟩ n .In terms of gates on a quantum computer, this operation is a parallel implementation of Hadamard gates on all the qubits.Additionally, we would like to remark that due to normalization, any linear function with an arbitrary slope and 0 offset will lead to the same encoding.However, it is possible to vary the slope by changing the offset of the linear function under the restriction of the quantum state to be normalized.
< l a t e x i t s h a 1 _ b a s e 6 4 = " z 1 d 4 a 3 y N i H i V P T w 4 8 J P W / l K p s m c = " + l q 0 F J 5 8 5 h T 9 w P n 8 A W P K N w Q = = < / l a t e x i t > (ii)  for DHWT method in presence of noise.This is shown in (ii) and (iii), where fidelity and error measured in the L 2 norm for both the ideal and noisy DHWT methods are plotted with respect to k 0 .We have not considered noise for the case MPS χ = 1 as its loading circuit is only a layer of single qubit rotations, and its impact is marginal.

Approximated loading of linear polynomials via DHWT
In the preceding section, we have conclusively demonstrated the efficient loading of linear functions.Now, the inevitable query arises: can the Hadamard-Walsh series be deliberately truncated while maintaining control over the process?In other words, is it possible to strike a balance between the number of terms truncated and the resulting error, thereby achieving a quantum state approximation that accurately encodes our intended target?
We now proceed to illustrate how the loading of the linear function can be approximated by truncating the Hadamard-Walsh series.Let us assume that we have the DHWT for n qubits, then the non zero coefficients are with |x (1) n |.We keep x 0 and the largest k 0 values of the coefficients with |k| b = 1, i.e. ⃗ h Next, we construct the circuit to generate the state encoding these renormalized coefficients of the truncated series.Then, the fidelity of the resulting state, |Φ More details about how to derive this expression are given in the Appendix A.
Note that while the structure of the circuit is the same, see Fig. 3, the angles have changed its value to θ k 0 k = arcsin(2x the new normalization factor and G Now, assuming an infidelity ϵ = 1 − F , it is possible to obtain the expression which establishes the trade off between infidelity and truncation.In the asymptotic limit n → ∞, we obtain Additionally to this analysis, we also study the deviation of the amplitudes of the approximated state, Φ with respect to the ideal state, |Φ L ⟩ n .
We now focus on describing the state Φ that comes from the truncation of the series, denoted as ⃗ h We compare the k 0 qubit state obtained from ⃗ h ) with the state produced by the series corresponding on doing the DHWT to the linear function encoded on k 0 qubits given by ⃗ h k 0 = (x k 0 x (0) From this comparison we can obtain that

Therefore we can write
with β = 2 (n−k 0 )/2−1 (2 n−k 0 − 1) the offset that will induce the maximum deviation in the amplitude between both the ideal and approximated quantum state.Thus, we can conclude that when truncating the sate, we are loading the k 0 -qubit state Φ the corresponding normalization factor.When adding the remaining n − k 0 qubits, this will lead to an state with a degeneracy 2 n−k 0 on every of this amplitudes (graphically a step wise function) Finally, merging both sums we finally get From these expression we can define the deviation of the amplitude of |j⟩ as In the asymptotic limit the upper bound of δ (2) j given by β/C 1/2 , (17) which converges to 0 when k 0 tends to infinity as O( 12 k 0 ).We depict the complexities of loading the exact and approximated versions of linear function either via MPS or DWHT in Tab. 1.
Alternative to this set up, one can also consider a MPS based methodology to load the linear function, which in the exact case scales linearly on the number of qubits, see Section 2.3.Thus, from now on, we will assume that we have access to a state preparation oracle of the exact linear function on k 0 qubits denoted as U L,k 0 , for 1 ≤ k 0 ≤ n, along with its adjoint and their controlled variants.This allows us to load the linear function either in an exact manner with a circuit depth that scales O(k 0 ).Note that in the worst case scenario, the controlled version can be achieved by controlling every gate of the oracle, which keeps the complexity invariant.

Polynomial transformation of amplitudes via quantum singular value transformation
Once we have introduced the quantum circuit that loads the linear function into k 0 qubits via the unitary operator U L,k 0 with a complexity of O(k 0 ), our method uses the quantum singular value transformation (QSVT) [72][73][74][75][76][77] to achieve the polynomial transformation of the amplitudes.One also could consider the unitary to load the MPS of the linear function with χ = 1 on k 0 qubits, although the final results are considerable worst as we show later.
In this work we follow the procedure detailed in Ref. [78] and in its subsequent exponential improvement in Ref. [79], which provides a generalization that presents a diagonal block encoding and introduces the importance sampling.The remarkable insight of these works lies in the authors' demonstration of how it is possible to use the QSVT polynomial transformation to explicitly construct quantum circuits that apply the polynomial transformation to the amplitudes of any quantum state of interest, provided that the qubits method Exact Aprx  encoding unitary is known.
The steps needed to achieve the state amplitude encoding of the polynomial are: the block encoding of the amplitudes, the polynomial transformation of the block encoding and the quantum amplitude amplification algorithm.

Block encoding
The first step is the block encoding of the amplitudes of the linear function given the unitary U L,k 0 that prepares the state that represents the linear function.
According to Theorem 2 in Ref. [79], given a k 0 -qubit state |ψ⟩ = queries to a controlled version of U .Note that if the amplitudes are real the number of ancillary qubits required for the block encoding is k 0 + 2.
Here we use this result to create from U L,k 0 the unitary U A L,k 0 which is a (1, k 0 + 2, 0) block encoding of the Hermitian operator An alternative (1, 1, 0)-block encoding U B resulting by following the idea of Ref. [41] could be implemented by applying the unitary dilation technique [101] This operation would require an efficient simulation of the Hamiltonian H = arccos(B) [102].Discussions of more possible block encoding for the linear function are shown in Appendix C.

Polynomial transformation of the block encoding
Once that we have obtained the amplitude block encoding, we show how to implement polynomial transformations of complex amplitudes via the Quantum Singular Value Transformation (QSVT) [41,42,75,78,79].The construction and efficiency of explicit quantum circuits for applying the polynomial transformation to amplitudes of a k 0 -qubit quantum state rely on the complexity of implementing the block encoding of the amplitudes, denoted as U L,k 0 .It is crucial that the complexity of U L,k 0 remains of the order O(k 0 ) to prevent it from dominating the overall algorithmic cost.
By applying Lemma 2 with a polynomial P (x) to U L,k 0 one obtains a block encoding U P,k 0 of P (A L,k 0 ).Additionally, note that due to the normalization, the values of the amplitudes of the state encoding the linear function have been renormalized to the interval [0, (2 . Therefore, to standardize the application of polynomials to the interval [0,1], we must compose the polynomial with the rescaling transformation, which keeps the degree of the polynomial invariant.For simplicity's sake, we will ignore this transformation from now on.

Amplitude amplification
We now apply the unitary U P (A) to the initial state |+⟩ ⊗k 0 |0⟩ ⊗k 0 +5 .This results into a quantum state The success probability of measuring given by 0.0625F 2 P , with F P being the filling ratio defined as and Therefore the quantum amplitude amplification algorithm would require O(4/F) rounds of the oracle that prepares the state to boost the success probability to 1 [41].
Note that the complexity of this protocol depends on the polynomial encoding that we aim to achieve.In this sense, some particular transformations will head to an efficient circuit, while others will introduce an exponential consumption of resources, depending on the filling ratio F.
Once that we have obtained the state |Φ P ⟩ k 0 we can tensor it with the state |+⟩ n−k 0 to obtain a quantum state Φ (k 0 :n) P n =|Φ P ⟩ k 0 ⊗|+⟩ n−k 0 that approximately encodes the polynomial P on n qubits.
Before concluding this section, we would like to remark that for the particular case the polynomial transformation satisfies P (0) = 0, one can use the importance sampling technique recently presented in Ref. [79] (Theorem 3) to achieve an exponential enhancement in the overall complexity do to the improvement in the step that encodes the transformed amplitudes into the final state.

Amplitude deviation in the polynomial encoding
Finally, once the polynomial transformation has been achieved, we can analyze how the error coming from the approximation of loading the exact k 0 qubit polynomial into n qubits.In order to set a proper comparison with Ref. [41,42], we compute the error in the same way, (20) for our particular case we assume ||P exact (x)|| ∞ = M and therefore If we now assume with D = max k {|c k |}.In the asymptotic limit this upper bound decays as O( 1 2 k 0 ).Note that one might think that it can be much cheaper to prepare an appropriate product state and then try to achieve an approximation to the desired polynomial encoding by applying a polynomial function to the product state via QSVT.However, since the complexity of the block encoding is not dominated by the oracle that loads the linear function, a reduction in complexity would not be achieved.Additionally, due to the approximation of MPS,this leads to a non-uniform grid discretization would result in a higher overall error when transforming the amplitudes.

Numerical Results
Having established the theoretical framework, in this section we present numerical simulations

Method
. Methods are arranged in ascending order based on fidelity of the polynomial encoding (fifth column).We also show the fidelity of loading the linear functions for those 2 steps methods (second column) and the filling ratio F (third column) for the different k 0 values.
comparing the methods outlined in this paper.The primary objective is to provide empirical support for our analytical findings.

Linear function
We begin our analysis by evaluating the performance of the exact and approximated loading of the linear function and its its robustness in terms of fidelity against noise.
The analysis of loading the linear function is depicted in Fig. 4, where we have considered n = 6 qubits.In Fig. 4 (i) we have depicted a comparison of the ideal case χ = 1 for the MPS technique, the noisy and ideal cases with k 0 = 3 for the DHWT, and the exact encoding.In Fig. 4 (ii) and (iii), we observe a trade-off between the truncation error and the experimental error for different values of k 0 when using the DHWT technique.The two figures of merit considered are the fidelity and the L 2 norm, which is defined as , where ⃗ v corresponds to the difference between the approximation and the ideal quantum state.We observe that the highest fidelity and the smallest L 2 norm error are achieved at truncation levels k 0 = 2 and k 0 = 3.In the comparison, we also include the ideal MPS with χ = 1.

Example of polynomial function
We present a second analysis where we consider the encoding of a polynomial function.In Fig. 5, we present the loading of the polynomial function with C p the normalization factor, by using n = 6 qubits and the two methods studied in this paper.Our approach involves two distinct loading strategies: first, loading the linear function and then applying the polynomial transformation with the QSVT method, ensuring no induced errors (Figs. 5 (i) and (ii)), and second, using the matrix product state (MPS) representation of the polynomial itself (Fig. 5 (iii)).For additional examples see Appendix D.
In order to load the linear function, we utilize either the DHWT method introduced in this work for different truncation values k 0 , or the MPS approach with χ = 1.In Tab. 6 we show the L 2 norm and fidelity of the final state with respect to the exact state for each methodology.We observe that for k 0 ≥ 5, the value of the fidelity < l a t e x i t s h a 1 _ b a s e 6 4 = " z 1 d 4 a 3 y N i H i V P T w 4 8 J P W / l K p s m c = " + l q 0 F J 5 8 5 h T 9 w P n 8 A W P K N w Q = = < / l a t e x i t > (ii) < l a t e x i t s h a 1 _ b a s e 6 4 = " R Y 0 I 1 n l q Q p f H l v s e I K r L 4 3 x a S ) for n = 6 qubits using different methods.For this particular example the filling ratio for k 0 = 6 is F = 0.6184.The worst and best implementation methods in which the linear function is loaded first, followed by the application of QSVT are displayed in (i) and (iii), respectively.In (ii), results for the loading using MPS are shown, achieving perfect loading with χ = 4, although theory predicts χ = 5 as the upper bound.The L 2 norm and fidelities for each case are displayed in Tab. 6.
achieved by the combined protocol is better than the direct polynomial MPS technique with χ = 2.
Additionally, the fidelities resulting from combining the approximate loading of the linear function with the QSVT to perform the polynomial transformation are presented.Remarkably, the loading of the linear function via a χ = 1 MPS achieves a fidelity of F = 0.9885.This suggests that the quantum state of the linear function is nearly a product state.To gain deeper insights, we performed an analysis of the single-qubit rotations that generate this approximate state.We fitted the angles of the rotations to an analytical expression and leveraging this fitting information, we trained a variational circuit aimed at preparing a product state that optimizes the fidelity with respect to the exact linear function.It is crucial to emphasize that the effectiveness of the fitting is not guaranteed across different number of qubits, hence we undertook this approach as a validation method.For additional details, please refer to the Appendix B.

Comparison with other methods
The first statement we would like to highlight is that we have put large part of our efforts in achieving an approximate loading of functions as simple as the linear function, by introducing a controllable error that reduces the depth of an already efficient circuit.As far as our knowledge is, there is no previous result in the literature that can do the same task with a comparable performance.That being said, we proceed to compare our method with similar results.
The QSVT enables the application of polynomial transformations to the amplitudes of a quantum state.However, when utilizing this technique to load a polynomial function encoded in the quantum state's amplitudes, the efficient loading of the linear function is crucial.Previously in literature, authors in Ref. [41,42] explored the possibility of applying the QSVT to block encoding of the sine function.Our method posses the same advantages of using QSVT that they mention: 'we avoid discretizing the values the function can take, providing instead a continuous approximation to the function.Our method is straightforward and versatile, as the same circuit template can be used for a wide range of functions.In contrast, our method avoids the error that propagates into the final loaded state due to the polynomial approximation of the arcsin function by efficiently loading the block encoding of linear function, and the subsequent polynomial transformations applied to load the desired polynomials into the amplitudes.In this context, our approach focuses on implementing the block encoding of the linear function rather than the sinusoidal function.To achieve this, we have proposed a method based on the Hadamard-Walsh Transform.The cost of this replacement in the block encoding can me mainly expressed in terms of adding n + 1 extra ancillary qubits to the circuit, see Tab. 4 for comparison.From this table we can observe that or the encoding of a polynomial without error, our algorithm scales with the same complexity as [41,42]   ror resulting from the arcsin approximation.In order to explore the polynomial degree overhead needed in the methodology presented in Ref [41] to encode a polynomial into amplitudes of a quantum state we present Tab. 5. Additionally, our methodology allows the introduction of a controllable error that reduces the complexity of the exact case for the encoding of polynomial functions.In the general case of loading an arbitrary function the trade off between the additional resources of our protocol versus the approximation of arcsin should be taken into account.For instance, [41,42] are more efficient for loading trigonometric polynomials.
that for such functions, favorable outcomes can be obtained by employing a fixed bond dimension of two, attributed to the logarithmic scaling of entanglement entropy in these cases [62].In this paper, we have explored the resource requirements of this approach and compare it to our linear function+QSVT approach, specially in the case of loading of the linear function, for which we have analytical expressions for fidelity and error propagation.
Considering a different perspective, it was recently proposed the use of the Hadamard-Walsh series [58] for amplitude encoding.This proposal leverages the fact that for functions whose derivative is bounded, the error resulting from truncating their discrete Hadamard-Walsh Series is exponentially suppressed with the index of the truncation [103].The authors use Hamiltonian simulation techniques presented in [99] to achieve a Hadamard-Walsh approximated simulation of the unitary U = e −i f ϵ 0 with error ϵ 1 , where f = x f (x) |x⟩ ⟨x| is the operator corresponding to the target function to be encoded.Finally they use an ancillary qubit to generate the operator −i(I − e −i f ϵ 0 ) acting on the state x |x⟩ |1⟩, which approximates the target state to first order of ϵ 0 and introduces a protocol conditioned to the probability of measuring the ancillary qubit in the state |1⟩.Therefore the total protocol introduces two sources of error, ϵ 1 corresponding to the truncation of the Hadamard-Walsh series and ϵ 0 corresponding to the Taylor expansion.This technique introduces an error for loading functions even for those cases in which ϵ 1 = 0, as in the linear case.By contrast our subroutine that uses the DHWT leverages the sparsity of the series corresponding to polynomial functions to efficiently generate quantum circuits that encode the states directly into the amplitudes.Additionally, this state loading algorithm based on the Hadamard-Walsh transform is not as efficient as our methodology for the application of loading the linear function, which is an important factor when loading the linear function for the derivative pricing use case.
Finally, we remark that our method is not limited by the condition of a bounded derivative of the target function, as is the case for the MPS [61] or the Hamiltonian simulation based on DHWT [58].

Conclusions
In this article, we have considered the problem of loading real polynomials into a quantum computer, with a particular focus on the encoding of the linear function.We have presented and reviewed two methodologies based on different approaches.The first method is based on matrix product states, and even though it offers competitive results for some particular cases [61], it is difficult to precisely control the error the method incurs.The second algorithm introduced in this paper utilizes the discrete Hadamard-Walsh transform (DHWT) to achieve the block encoding of the amplitudes which is fed into the quantum singular value transformation (QSVT) in order to apply a polynomial transformation to the amplitudes.This technique is able to exactly encode polynomials with same complexity as previous works [41,42] that incurred into an approximation error.Furthermore we have been able to reduce the complexity of the protocol via introducing a controllable error.
Using this technique based on the DHWT, the coefficients of the Hadamard-Walsh series of a given function are loaded into a quantum state and the inverse discrete Hadamard-Walsh transform (DHWT) is applied to achieve the amplitude encoding of the target function.This idea constitutes a novel and promising approach for functions whose Hadamard-Walsh series can be efficiently encoded, as for instance when it sparse [54] or efficiently truncated [99].
We would like to remark that even though our work has been focused on encoding real polyno-mials, it could be easily extended to load polynomials valued on complex values, i.e.P : C → C , multivariate polynomials or even non-linear functions approximated with polynomials [41,78,79].Future work will consider the feasibility of using our DHWT-based method to load highly discontinuous square wave-like functions.

Appendix D Additional numerical examples
Here we provide two more numerical experiments for n = 6 qubits for two different polynomials.

Figure 2 :
Figure 2: Discrete Hadamard-Walsh Transform for n = 3 qubits in the natural order representation.

Definition 2
The discreteHadamard-Walsh transform (DHWT) is a linear, orthogonal and symmetric operation that transforms discrete signals or sequence of sorted data into a new representation given by the Hadamard-Walsh Series

<
l a t e x i t s h a 1 _ b a s e 6 4 = " i b I e 7 H D 6 A l z V w s m t E d H J w / S Y 5 6 I = " > A A A B 5 H i c b Z D N T g I x F I X v 4 B / i H + r S T S M x c U V m D F G X R D c u M c p P A h P S K X e g o Z 0 Z 2 4 4 J m f A G u j L q z i f y B X w b C 8 5 C w b P 6 e s 9 p c s 8 N E s G 1 c d 0 v p 7 C y u r a + U d w s b W 3 v 7 O 6 V 9 w 9 a O k 4 V w y a L R a w 6 A d U o e I R N w 4 3 A T q K Q y k B g O x h f z / z 2 I y r N 4 + j e T B L 0 J R 1 G P O S M G j u 6 e + i 7 / X L F r b 2 m X m l c 5 o c o w h E c w y l 4 c A 4 N u I Y m t I D B E J 7 h D d 6 d 0 H l y X p z X n 2 j B y f 8 c w h 8 5 H 9 / 3 I Y s 4 < / l a t e x i t > q 2 w b C 8 5 C w b P 6 e s 9 p c s 8 N E s G 1 c d 0 v p 7 C y u r a + U d w s b W 3 v 7 O 6 V 9 w 9 a O k 4 V w y a L R a w 6 A d U o e I R N w 4 3 A T w b C 8 5 C w b P 6 e s 9 p c s 8 N E s G 1 c d 0 v p 7 C y u r a + U d w s b W 3 v 7 O 6 V 9 w 9 a O k 4 V w y a L R a w 6 A d U o e I R N w 4 3 A T 8 5 h D 9 g f X x D S / 6 k Q M = < / l a t e x i t > Ry(✓n 1) < l a t e x i t s h a 1 _ b a s e 6 4 = " J s m z 2 9 u Y S W H 2 j e 6 8 G i D O l 2 D d 2 / Q = " > A A A B 9 H i c b V D L T s J A F J 3 i C / F V d O l m I j H B h a Q l v p Z E N y 7 R y C O B p p k O t z B h + s j M F N M 0 / I m u j L r z S / w B / 8 Y B u 1 D w r M 6 9 5 9 z c e 4 8 X c y a 7 B c b 6 s w t L y y u p a c b 2 0 s b m 1 v W P v 7 r V 0 l C j K m j Q S k e o E R D P B J W s C B 8 E 6 s W I k D A R r B + P r q d 5 + Y E r z S N 5 D G j M v J E P J B 5 w S M C v f L t / 5 a b U H I w b E z + T x 2 H c m R 7 5 d c W r O D H i R u D m p o B w N 3 / 7 s 9 S O a h E w C F U T r r u a t e x i t > H < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 p R N e M r / q r E 1 6 w l A T 1 6 i 17 H A f Y U = " > A A A B 6 n i c b Z D N T g I x F I X v 4 B / i H + r S T S M x c U V m D F G X R D c u M Z E f A x P S K X e g o Z 2 Z t B 0 T M v I S u j L q z s f x B Xw b C 8 5 C w b P 6 e s 9 p c s 8 N E s G 1 c d 0 v p 7 C y u r a + U d w s b W 3 v 7 O 6 V 9 w 9 a O k 4 V w y a L R a w 6 A d U o e I R N w 4 3 A T r 7 x 7 t l M 6 O h 4 V M U P W y D r Z J B 7 Z J 0 f k l F R J j X B y T x 7 J C 3 l 1 7 p w H 5 8 l 5 / r K O O a O b V f I D z v s n Z V m e a w = = < / l a t e x i t > Loads | ˜ (k0:n) L i n < l a t e x i t s h a 1 _ b a s e 6 4 = " r y 8 i e 9 M h 5 0 l 9 a O S c 1 o 6 u T 4 u l i 8 m R e T I L t k j h 8 Q h Z 6 R M r k i V 1 A g n j + S Z v J F 3 6 8 F 6 s l 6 s 1 x / r n D W 5 2 S F / Y H 1 8 A 2 s D m r 4 = < / l a t e x i t > Loads | ˜ L i n

Figure 3 :
Figure 3: Circuit implementation for loading the quantum state |Φ L ⟩ n , denoted as U L,n with complexity O(n) multi-controlled gates.The angle of every multi-controlled rotation is given by θ k = arcsin 2x

25 < l a t e x i t s h a 1 _ b a s e 6 4 =
" w g l r n M O / 0 e C 8 y X 1 4 t 0 w p V J h y O n I= " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o M Q L 2 F X 4 u M Y 8 O I x o n l A s o T Z S W 8 y Z H Z 2 m Z k V Q s g n e P G g i F e / y J t / 4 y T Z g y Y W N B R V 3 X R 3 B Y n g 2 r j u t 5 N b W 9 / Y 3 M p v F 3 Z 2 9 / Y P i o d H T R 2 n i m G D x S J W 7 Y B q F F x i w 3 A j s J 0 o p F E g s B W M b m d + 6 w m V 5 r F 8 N O M E / Y g O J A 8 5 o 8 Z K D 2 V + 3 i u W 3 I o 7 B 1 k l X k Z K k K H e K3 5 1 + z F L I 5 S G C a p 1 x 3 M T 4 0 + o M p w J n B a 6 q c a E s h E d Y M d S S S P U / m R + 6 p S c W a V P w l j Z k o b M 1 d 8 T E x p p P Y 4 C 2 x l R M 9 T L 3 k z

45 <
l a t e x i t s h a 1 _ b a s e 6 4 = " R Y 0I 1 n l q Q p f H l v s e I K r L 4 3 x a S x A = " > A A A B 7 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S L U S 0 m k f h w L X j x W M G 2 h D W W z n b R L N 5 u w u x F K 6 W / w 4 k E Rr / 4 g b / 4 b t 2 0 O 2 v p g 4 P H e D D P z w l R w b V z 3 2 1 l b 3 9 j c 2 i 7 s F H f 3 9 g 8 O S 0 f H T Z 1 k i q H P E p G o d k g 1 C i 7 R N 9 w I b K c K a R w K b I W j u 5 n f e k K l e S I f z T j F I K Y D y S P O

Figure 4 :
Figure 4: Comparison of different methods for loading the linear function on n = 6 qubits.(i)Illustrates the resulting state by using different encoding protocols, DHWT with k 0 = 3 and MPS χ = 1, in noisy an ideal scenarios.The choice of k 0 = 3 is due to the fact this case has a high fidelity and the minimum error measured in the L 2 norm for DHWT method in presence of noise.This is shown in (ii) and (iii), where fidelity and error measured in the L 2 norm for both the ideal and noisy DHWT methods are plotted with respect to k 0 .We have not considered noise for the case MPS χ = 1 as its loading circuit is only a layer of single qubit rotations, and its impact is marginal.

25 <
l a t e x i t s h a 1 _ b a s e 6 4 = " w g l r n M O / 0 e C 8 y X 1 4 t 0 w pV J h y O n I = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o M Q L 2 F X 4 u M Y 8 O I x o n l A s o T Z S W 8 y Z H Z 2 m Z k V Q s g n e P G g i F e / y J t / 4 y T Z g y Y W N B R V 3 X R 3 B Y n g 2 r j u t 5 N b W 9 / Y 3 M p v F 3 Z 2 9 / Y P i o d H T R 2 n i m G D x S J W 7 Y B q F F x i w 3 A j s J 0 o p F E g s B W M b m d + 6 w m V 5 r F 8 N O M E / Y g O J A 8 5 o 8 Z K D 2 V + 3 i u W 3 I o 7 B 1 k l X k Z K k K H e K3 5 1 + z F L I 5 S G C a p 1 x 3 M T 4 0 + o M p w J n B a 6 q c a E s h E d Y M d S S S P U / m R + 6 p S c W a V P w l j Z k o b M 1 d 8 T E x p p P Y 4 C 2 x l R M 9 T L 3 k z 8 z + u k J r z x J 1 w m

Figure 5 :
Figure 5: Different methods to load the polynomial function P(x) = 1 Cp (x − 1/(2 n − 1))(x − 20/(2 n − 1))(x − 50/(2 n − 1))(x − 60/(2 n − 1)) for n = 6 qubits using different methods.For this particular example the filling ratio for k 0 = 6 is F = 0.6184.The worst and best implementation methods in which the linear function is loaded first, followed by the application of QSVT are displayed in (i) and (iii), respectively.In (ii), results for the loading using MPS are shown, achieving perfect loading with χ = 4, although theory predicts χ = 5 as the upper bound.The L 2 norm and fidelities for each case are displayed in Tab. 6.

Table 1 :
Complexities for loading the linear function into a different number of qubits k 0 and n, with k 0 < n, either in an exact or approximated form.k0(ϵ) corresponds to Eq. (12).

Table 2 :
Noise parameters description and their value.We have estimated the numerical values from the calibration data provided for the IBM device 'ibm_jakarta'.

Table 3 :
Error in L 2 Norm (fourth column) and Fidelities (fifth column) for Different Loading Methods of P with the er-

Table 4 :
[41,42]son between the methodology proposed in Refs.[41,42]and the current work, with||P exact (x)|| ∞ = M , F the filling ratio, D = max k {|c k |} and δ ∞ = ||P (j)/||P (j)|| ∞ − P (j)/|| P (j)|| ∞ || ∞ , where P = d k=0 c k x kis the exact polynomial and P its approximation, and B ≥ k |c k |.We neglect amplitude amplification and precomputational errors and complexities.In both cases where QSVT is used we have assumed no error in the classical preporcesing of the angles needed to implement the polynomial transformation.

Table 6 :
Error in L 2 Norm and Fidelities for Different Loading Methods of two different polynomials P 1 (x) = 1 Cp