Quantum-inspired algorithms for multivariate analysis: from interpolation to partial differential equations

In this work we study the encoding of smooth, differentiable multivariate functions distributions in quantum registers, using quantum computers or tensor-network representations. We show that a large family of distributions can be encoded as low-entanglement states of the quantum register. These states can be efficiently created in a quantum computer, but they are also efficiently stored, manipulated and probed using Matrix-Product States techniques. Inspired by this idea, we present eight quantum-inspired numerical analysis algorithms, that include Fourier sampling, interpolation, differentiation and integration of partial derivative equations. These algorithms combine classical ideas---finite-differences, spectral methods---with the efficient encoding of quantum registers, and well known algorithms, such as the Quantum Fourier Transform. {When these heuristic methods work}, they provide an exponential speed-up over other classical algorithms, such as Monte Carlo integration, finite-difference and fast Fourier transforms (FFT). But even when they don't, some of these algorithms can be translated back to a quantum computer to implement a similar task.


Introduction
Quantum computers use the exponential capacity of a Hilbert space to process information. A quantum computer with m qubits can store 2 m complex numbers as the components in of the quantum register wavefunction, |ψ = 2 m −1 s=0 ψ s |s . A quantum algorithm creates, manipulates and probes these amplitudes to solve a concrete problem. In this work we discuss algorithms where the amplitudes ψ(s) encode a smooth function defined over some a volume in R N . A common example is storing probability distributions in the quantum register [12] and developing algorithms to extract expected values [18] or conditional probabilities [33]. With this, it becomes possible to perform valuations of complex financial assets [22,26], VaR estimates [7], and other sophisticated interrogations. Using a similar encoding, one may also address radically different problems, such as solving partial differential equations with finite differences [5,6,8].
In this work we discuss how efficient it is to encode discretized functions in a quantum register. We find that for certain distributions-smooth differentiable functions with bounded derivatives or bounded spectrum-, the accuracy of the discretization increases exponentially with the number of qubits, while the bipartite entanglement grows slowly or even remains bounded with the problem size. This implies such distributions may be constructed with quasi-local operations and polynomial resources on quantum computers. But it also opens an exciting possibility: a family of quantum-inspired matrix-product state (MPS) techniques that, under approximations of low entanglement, represent the quantum register efficiently and provide new classical (and quantum) algorithms for interpolating, differentiating, Fourier transforming or solving differential equations of such distributions. These techniques work efficiently because of an implicit renormalization where different qubits work with different length scales, in a way that lends itself to efficient interpolation and compression. This idea, with strong parallelisms to the 2D quantum image processing world [16], gives rise to a performance improvement over earlier techniques based on tensor trains [3,11] or MPS encodings of mode expansions [14].
This paper is structured in three parts. The bulk of the work is preceded by a summary (Section 2) of the main results and heuristic algorithms that are developed in this work. The first section addresses the representation of discretized functions in quantum registers. It presents state-of-the-art techniques (Sect. 3.1) and new discretizations (Sect. 3.2) on an equal footing. We argue that these samplings produce to weakly entangled multi-qubit states (Sect. 3.3). This prediction is confirmed numerically for common distributions in 1, 2 and 3 dimensions, using exact simulations of up to 28 qubits and MPS simulations of up to 36 qubits. The second part of this work introduces the idea of MPS quantum registers, whereby we encode multivariate functions in arrangements of qubits that are represented, manipulated and interrogated using the MPS representation. Section 4 further develops this idea, recalling well known algorithms from the literature and how they specialize for our purposes. With these tools at hand, Section 5 develops new quantum-inspired algorithms for the numerical analysis of multivariate functions. These include the mapping of functions to MPS format (Sect. 5.5), Fourier analysis (Sect. 5.1), interpolation methods (Sect. 5.2), and techniques for approximating derivatives of discretized functions (Sect. 5.3), both through finite-difference (Sect. 5.3.1) and the Quantum Fourier Transform (QFT, Sect. 5.3.2). Sect. 5.4 combines these techniques into higher level algorithms for solving partial derivative equations, demonstrating their performance in the Fokker-Planck model. This work is closed with a discussion of the results, including connections to recent advances in tensor-based numerical analysis, and an outlook of applications.

Summary
Quantum computing has introduced the clever idea of encoding probability distributions in quantum registers. In this encoding, a small number of qubits can store an exponentially large number of function samples. To fix ideas and notation, let us take a function p(x 1 , . . . , x N ) of N variables in a bounded interval. Let us use m qubits to uniformly discretize the domain The non-negative integer s i takes all possible values obtained by grouping m bits s 1 i s 2 i · · · s m i , ordered in decreasing significance. The N integers or N m bits can be associated to different states of a quantum register, enabling two representations of the function p. The first one assumes a non-negative function, |p ∝ s 1 ,...,s N p(s 1 , . . . , s N ) |s 1 ⊗ · · · |s N , s i ∈ {0, 1, . . .
Both representations can be extended to situations where the functions are not normalized, just by keeping track of global prefactors. In this work we argue that both representations are exponentially efficient in many ways. First, the quantum register demands only a logarithmically growing number of qubits N m to store an exponential amount of weights 2 N m in a discretized function. Second, we also need an exponentially small number of qubits m ∼ − log(ε) to reduce the discretization error below a given tolerance ε. Third and finally, we find that for smooth, differentiable functions with bounded derivatives, these states have a small amount of entanglement. Indeed, for many distributions of interest we obtain the scaling O(N ) with the dimension of the problem. We conjecture that this behavior is due to an implicit renormalization that happens in the quantum register, where some bits s 1  2) and for the solution of partial differential equations (Sect. 5.4). Table 1 summarizes the algorithms discussed in this work, paired with alternatives that already exist for quantum computers or in the field of numerical analysis. The table summarizes the costs of those algorithms, expressed in terms of well known quantitiesdiscretization error, number of qubits, estimated entanglement and bond dimension size, etc-. In the case of quantum-inspired numerical analysis, we must emphasize that the performance metrics are heuristic. However, if entanglement remains bounded throughout the simulations, the quantum register method demands a small bond dimension χ, and the algorithms provide an exponential speedup over other classical techniques-from finite differences to the highly performant FFT techniques.

GR discretization
One of the earliest works suggesting the encoding of functions in quantum registers is the unpublished manuscript by Grover and Rudolph [12]. This designed a unitary operator U p that encodes a probability distribution p(x) in an empty quantum register with m qubits |p (m) The original construct assumes a random variable x in a bounded interval [a, b] subdivided into 2 m smaller intervals, labeled by the quantum register states |s = |s 1 s 2 · · · s m . A practical application of this encoded state would be the computation of expected values for any observable or function f (x). This requires engineering an observableÔ Typically, we approximateÔ |s s| , and apply a uniform discretization (1), to have an integration error that decays exponentially with register size, ε int ∼ O(δ m ).
As found by A. Montanaro [18], using amplitude estimation with U p and the operator O f , one may estimatef with a precision that scales better than Monte Carlo algorithms. If the cost of implementing U p is T GR , and we aim for a sampling precision ε sample , the asymptotic time cost of the ideal amplitude estimation algorithm is roughly This represents a favorable scaling when compared with traditional Monte Carlo, where the sampling uncertainty goes as O(ε −2 sample ), but only if the cost of encoding the probability state T GR remains small or weakly dependent on the integration error ε int . The unitary by Grover and Rudolph's U p is a recursive construct that adds one more qubit of precision in each step. As sketched in Fig. 1a, the procedure reads where we identify |s |0 =: |2s and |s |1 =: |2s + 1 . The rotation angle 0 ≤ θ i ≤ π/2 is obtained from two identities Unfortunately, this algorithm requires an exponentially large number of angles and involves a highly non-local unitary with a potentially bad decomposition. The GR algorithm must be therefore considered more a proof of existence, than a practical recipe that can be used when bounding the resources of Monte Carlo analysis. We will study the GR construct and other discretizations, demonstrating that there are efficient alternatives to equation (7). Our analysis centers on the complexity of the sampled states. Using the bipartite entanglement as quantifier, we will show that the cost of adding one more qubit of resolution decays exponentially. This will help us understand that there one quasi-local unitary procedure that builds the GR state with a cost that is polynomial in the number of qubits, T GR ∼ O(−m log(ε)). We will confirm numerically this result using various well-known probability distributions.

Other discretizations
The integral representation by Grover and Rudolph reproduces exactly the probability that is contained inside each interval, but it requires computing 2 m integrals. In practice, this is unnecessary because the estimation of expected values already introduces a discretization error (5) ε sample ∼ O(δ m ) in the operator definition. It is not difficult to find simpler representations that have the same or better scaling. The obvious one is the uniform sampling of the probability distribution The standard error bound for this Riemann-type state is ε sample ≤ max d dx (f p) δ m , which depends on the derivatives of both the sampled observables and probability.
The first order scaling is good enough for the simulations that we will show below, because the interval size decreases exponentially with the number of qubits m. We will therefore stick to the GR states or to Eq. (9), unless otherwise noted. However, if we need to save some qubits, we can try variations, such as a probability state that implements the trapezoidal or the Simpson rule The discretization error of the Simpson state decreases faster with the interval size, ε sample ∼ O(δ 2 m ). This scaling means that the qubits required to achieve a given precision m ∼ log 2 (ε sample ) can be half those required by the uniform ansatz. These savings may be interesting in resource-starved architectures, such as NISQ computers, and also in the algorithms to be considered later in Sect. 4.

Existence of efficient constructs
We can show that adding k new qubits to a GR probability state with m qubits demands a vanishingly small amount of entanglement, that decreases exponentially with 2 m . This small amount of entanglement suggests that the GR unitary U (m) p can be replaced with simpler decomposition in terms of quasi-local gates. The reasoning proceeds as follows. In App. A.1 we study the reduced density matrix ρ (m,k) of the k new qubits in a GR state with m + k bits of resolution. Assuming that the probability distribution is smooth, differentiable and has an upper bound on its derivative, we obtain upper bounds for the entropy of the extra qubit. For one added qubit, k = 1 our bound reads If k is larger, we have roughly Note that this argument easily extends to the Riemann-type (9) or Simpson sampling (10), because the differences between them decays exponentially with the number of qubits. The fact that adding every new bit requires a small amount of entanglement implies that the probability state |p (m) has an efficient matrix-product state (MPS) representation. A matrix product state is a decomposition of a wavefunction as a contraction of matrices that are labeled by the physical indices of a composite quantum system Each of the tensors is in general different, and has three indices There is one physical index s of dimension 2, and two bond dimensions of sizes χ i and χ i+1 . In our particular case, the fact that S[ρ (i,m−i) ] decays exponentially with i means that the bond dimensions also decrease extremely fast, making the representation (14) efficient.
It is known that MPS's admit an efficient, sequential construct [24,25]. This algorithm consists of m steps. On each step, an ancilla with log 2 χ i qubits is correlated with a fresh new qubit using a unitary operation of dimension (d max{χ i , χ i+1 }) 2 . This unitary operation is potentially smaller and more efficient than the GR unitary, because the whole process has bounded time-cost T GR ∼ O(mχ 2 ).

Numerical study of 1D distributions
Let us put these ideas to the test using three paradigmatic distributions We will also consider an unconventional function that also has a finite bandwidth, but which is not log-convex σ We have chosen all these distributions because we can compute the functions p (m) (x) exactly for all sampling sizes, constructing the GR states for up to m = 14 qubits 1 , analyzing their structure and complexity. In Fig. 2a we plot the entanglement required to enlarge the quantum register by one bit, from m to m + 1. This is the entropy S[ρ (m,1) ] in the notation above. As shown in Fig.  2b, the bipartition never exceeds one e-bit of entanglement, and exhibits an exponential decay at large sizes that goes as 2 −γm with γ between 1.73 and 1.84, depending on the simulation and probability distribution. The behavior is therefore more favorable than the bound from Eq. (12), which overestimates the entanglement. Fig. 2c shows that just like the entropy increase of adding one more qubit is small, the cumulative entropy obtained by studying all bipartitions of m + k qubits also remains small, and with a similar tendency. Moreover, for a detailed enough sampling with m = 14, there is no difference between the GR states and the simpler discretizations from Sect. 3.2, shown here with dash-dot lines.
To test whether these favorable dependencies are artifacts of our choice of distributions, we have varied the parameters of the distributions and also tested situations where one of them (18) acquires more features. Fig. 1d shows the maximum entanglement entropy over all bipartitions as we change the parameter σ. The first three probability distributions have a bounded entanglement below an e-bit. The non-log-convex distribution p nc (x) behaves slightly different: increasing σ leads to the appearance of more peaks, that are more difficult to describe. This causes a steady increase in entropy, but this is slow enough that still facilitates an efficient construct.

Extensions to more variables
Challenging applications to fields such as Physics, fluid dynamics or finance [20] require the study of states with many more, p(x 1 , x 2 . . .). To better understand the scaling of entanglement and the complexity of the state, we have studied the discretization (9) of two-and three-dimensional Gaussian distributions with covariance matrix Σ and zero mean, using the same number of qubits in all dimensions. As in the 1D problem, we will treat the quantum register as a one-dimensional arrangement of qubits, studying the entanglement over all 1D bipartitions. Naturally, the complexity of this discretization will depend on how we arrange the qubits. The simple straightforward order (A) distributes the qubits sequentially, first by coordinate, then by significance. In this order, Gaussian states with a diagonal covariance matrix Σ = diag{σ 1 , σ 2 , . . .} become products states of one-dimensional distributions, such as those studied in Sects. 3.3 and 3.4. We also introduce the order (B), where qubits are first sorted by significance and only then by coordinate. This order is inspired by the renormalization group, and deeply similar to the multi-scale representation of 2D quantum image encodings [16]. For a distribution with two random variables, x 1 and x 2 , using three qubits per variable, the two orders read Order (A): |s 1 s 2 → |s 1 1 |s 2 1 |s 3 1 |s 1 2 |s 2 2 |s 3 2 , and Order (B): |s 1 s 2 → |s 1 1 |s 1 2 |s 2 1 |s 2 2 |s 3 1 |s 3 2 .
x 1 / max Figure 3: (a) Two-mode Gaussian probability distribution with variances (σ min , σ max ), rotated an angle θ. We work with a discretization interval [−7σ max , 7σ max ] ⊗2 , compressing σ min down to 0.1, which is a 20dB squeezing of the variance. (b) Maximum entanglement of all bipartitions, for different variances, orientations and orders. In solid lines, we plot the trivial sampling (9) with N × m = 28 qubits in total. We compare those plots with the same simulation using an MPS at a higher sampling (circles). We probe two angles θ = 0 and π/4 in Eq. (21), and two orders (20). (c) Entanglement entropies for all bipartitions of N n = 24 qubits into (k, N m − k), for the highly squeezed state σ min = 0.1σ max with θ = π/4. (d) Similar plot but for the MPS algorithm using a N m = 36 qubits.
Let us begin the discussion using a general two-dimensional covariance matrix One quadrature is squeezed by a factor σ min /σ max , while rotating the frame of reference an angle θ, as shown in Fig. 3a. We sample this probability, reconstructing the exact wavefunction (9) with N m = 28 qubits. As shown in Table 2, this size approaches the limits of a decent computer, using 2 gigabytes of data in real double precision form. In contrast, the same distributions using MPS, the techniques from Sect. 5.5 and up to N m = 36 qubits consume less than 1 megabyte. The output of these simulations is shown in Fig. 2. When θ = 0, the wavefunction |p As shown in Fig. 2b (blue solid), the maximum entanglement over all one-dimensional bipartitions is less than one e-bit, consistent with Sect. 3.4. To grow the entanglement we must combine squeezing and rotation, recreating a two-mode squeezed state. From the theory of Gaussian states, the entanglement should be maximal for θ = π/4 and it should diverge with the squeezing. Our simulations confirm this prediction for the (A) order.  Table 2: Summary of resources to describe numerically the Gaussian probability distributions with covariance matrices given by (21) and (22), using m qubits per dimension, either in an exact formstoring all values p (m) (x s )-or in the compact MPS representation (14). Fig. 3b (orange, solid) shows that the maximum bipartite entanglement between our qubit variables grows as (σ max /σ min )1/4. This entanglement is also spread all along the chain of qubits, as seen in Fig. 3c. We can even improve on these results. As shown in Fig. 3b, if we adopt the renormalization order (B), all states can be described with N = 2 e-bits of entanglement. Moreover, the entanglement distribution concentrates around the most significant qubits [cf. Fig. 3c (dashed, green)], producing significantly smaller tensors. Despite the growth of entanglement, the quantum states that we create always admit a compact MPS representation that beats the classical approach of storing the full wavefunction. As shown in Table 2, a highly-correlated discretization in the (A) order with 14 and 18 bits per coordinate, requires 1Mb of floating point real numbers in MPS form. The same states stored using the (B) order, take, in the worst scenario σ min = 0.1σ max with θ = π/4, just about 80 kilobytes of information. All this is to be compared to the 2 Gb and 524 Tb of data required to write down the wavefunctions of 28 and 36 qubits.
We have performed the same study using three-dimensional Gaussian states. For concreteness, we have focused on a three-mode squeezed state that starts with a diagonal matrix diag(σ min , σ max , σ min ) and performs two identical rotations around the X and Z = /4% axes, with angles θ x = θ y = θ We show results for a large, exact simulation with N m = 21 qubits, which amounts to 256 3 points and 128 Mb of data, together with an MPS that is directly built with 33 qubits, a sampling of 2048 3 points. As in the two-dimensional case, the combination of squeezing and rotation deforms the state, which as shown in Fig. 4b becomes a titled rugby ball. Once more, the unrotated state remains weakly entangled. It is a product of three onedimensional probabilities, and the entanglement never exceeds one 1 e-bit [cf. Fig. 4a (blue, solid)]. The squeezing and rotation leads to a divergence of the maximum entanglement, but this divergence is once more cured by the (B) order. This importance-based structure brings down the entanglement to about 3 e-bits and a reduction of 54,000 in the information required. We conjecture that for these Gaussian states-and other smooth functions-the (B) order consumes at most N e-bits, giving a scaling of resources O(2 2N × N ) in both the time and memory costs of reproducing the probability distribution.

MPS quantum registers
In the previous section we have seen that it is possible to encode single and multimode probability distributions in quantum registers, that these states are typically weakly entangled and admit an efficient MPS representation. We now review the representation of wavefunctions and operators, and the basic ingredients in the MPS toolbox-tensor contraction and reordering, tensor renormalization, time evolution, etc-that will be used in Sect. 5 to implement actual algorithms. More precisely, Sect. 4.1 introduces two representations of multivariate functions: one following the precepts from Sect. 3.2, and another one that improves the computation of expected values (Sect. 4.3) and equation solving. We will also discuss the algorithm of MPS simplification (Sect. 4.4), which is an essential tool to implement all numerical analysis approximation schemes.

Quadratic and linear representation
Our overarching goal for the rest of this work is to encode multivariate functions using a virtual quantum register, and to store this register efficiently as an MPS. If the function we want to encode is non-negative-such as a probability distribution-, we can follow Sect. 3, identify √ p with a quantum register wavefunction (2) and use the MPS representation (14) to compress it. In this representation, observables f (x) become Matrix-Product operators (MPO's)Ô We call this the quadratic MPS representation because the mean value of an observablē f ∼ p|Ô f |p is a quadratic function of any of the tensors in the MPS state [cf. Fig. 5].
In the alternative encoding from equation (3), both observables f (x) and probability distributions p(x) become unnormalized vectors in a Hilbert space In this representation, each probability state is a linear form that maps observables to expected values, or vice versā We will call this strategy the linear MPS representation becausef is a linear function with respect to any of the tensors in |p or |f .

Memory cost
Let us denote χ f , χ p or simply χ the largest bond dimensions to encode those functions as MPS. Assuming an N -dimensional volume, discretized with m qubits per dimension, the space required by an MPS and an MPO scales as O(N mχ 2 ). Since χ dominates this scaling, we need to understand how this dimension behaves in typical problems.
The answer to this question is connected to the bipartite entanglement that is stored in the MPS representations. Since MPS's are obtained through a recursive Schmidt decomposition [30], the entanglement entropy over any sequential bipartition of the state is bounded by S ≤ log 2 (χ). Conversely, if the maximum entanglement is S, we will expect that the bond dimension scales as χ ∼ 2 S . Thus, for the smooth distributions from Sect. 3.4, where S ∼ N, we expect a scaling of resources of the form This represents an exponential saving over the space O(2 N m ) ∼ O(ε −1 int ) required to store the full wavefunction in general classical algorithms.
This exponential improvement is similar in origin to the one in quantum computers, as it also exploits the rapid growth of the Hilbert space with the number of qubits. However, in the MPS quantum register the gain is heuristic: it only appears for distributions with nice sampling properties. There are infinitely many problems where this improvement vanishes, due to the growth of S and the exponential blowup of the bond dimensions ξ. However, it seems that there are still many problems of interest where the MPS quantum register approach is useful, as we see below.

Integrals and expected values
In the linear MPS representation, the scalar product f |p is an approximation to the integral between both functions f (x)p(x)dx. As sketched in Fig. 6a, there is an optimal contraction scheme that starts from one boundary-the left-most site in Fig. 6a-, and sequentially contracts with one tensor from |p and one from |f . The optimal sequence of contractions in Fig. 6b demands O(dχ 3 ) operations in 2N m steps, with an estimated cost O(2dN mχ 3 ). We can perform a similar analysis for the quadratic MPS representation, where mean valuesf = p|Ô f |p involve contracting the state |p twice with the MPOÔ f that represents the observable. In this case the optimal procedure is slightly more involved, with a higher asymptotic cost O (N m × 4χ 4 ).
In order to compare these algorithms with other classical methods, we have to introduce the accuracy of the estimate. In the MPS representation the only error we make is the discretization error of discretizing N variables with m qubits, which scales as ε int ∼ N 2 −m . Thus, the MPS approximation to the integral demands a time T constr + O(−N log 2 (ε int /N ) × 2χ 3 ), where T constr is the time to build the MPS and the rest is a logarithmically growing cost associated to the contraction. We can compare this with Monte Carlo sampling, a good and general method for integration. The errors in this technique arise from the statistical uncertainty ε sample ∼ 1/ √ M . This error decays slowly with the number of iterations M, giving a time cost O(1/ε 2 ). The MPS therefore has the potential of providing an exponential speedup, given that (i) the bond dimension χ remains small and (ii) the cost of constructing the MPS states is also bounded.

Approximating states
In working with the quantum register, we will frequently need to apply operators that distort the MPS representation, increasing the size of the tensors. This is corrected by a process known as MPS simplification, which seeks the closest matrix-product state with the smallest bond dimensions, within a prescribed time and error tolerance. The simplification is an optimization typically defined with respect to the norm-2 distance between states Here p is the state we wish to approximate and φ is the new MPS. The distance d(φ, p) = p|p + φ|φ − 2Re φ|p is a quadratic form with respect to the tensors in φ, which is optimized iteratively, sweeping across the MPS [10,29] in a two-site DMRG-like process. The cost of this optimization has two parts. The estimation of the linear form φ|p involves a contractions like the ones shown in Fig. 6b, involving O(dχ 3 ) operations per site. On top of this, we apply a two-site simplification algorithm that optimizes pairs of tensors simultaneously, dynamically adapting the bond dimension. This has an extra cost O(4(dχ) 3 ) due to the singular value decomposition. Thus, assuming that we need T sweeps for convergence, the simplification time cost grows as O(T sweeps 4d 3 χ 3 ). In practical examples we find that T sweeps is very small: one or two sweeps reach the numerical precision of the computer, giving efficient results.
Finally, note that we can use the algorithm to construct a new state φ that approximates a linear combination of MPS This has a linear increase in the cost, O(kT sweeps 4d 3 χ 3 ), that has been exploited in other algorithms such as time evolution [10].

Function multiplication
In many algorithms below we will need to construct a state |f p that approximates the product of two sampled functions f (x)p(x). This operation can be implemented efficiently in at least four cases. Given is an arbitrary complex constant c and a discretization of N variables with m qubits per dimension, there exist MPO's for f (x) = cx, cx 2 and exp(cx) using bond dimensions 2, N m and 1, respectively. As illustration, let us discuss the implementation of cx. This operator is an MPO (23) with a bond dimension of size ξ = 2 that keeps track of whether any operator has been applied. The n − th tensor reads More generally, we can write the exponential of any QUBO formula as a product of MPO's with bond dimension 2 The MPO's inside the product are constructed with simple tensors. In particular, for the k−th step, the tensors read B s n sn a,b = δ s n sn      δ a0 δ b0 , n < k, δ bs 1 δ a0 exp(Q 11 s n ), n = k, δ ab exp(2aQ an s n ), k < n. (32) Using this idea, we can decompose any exponential or Gaussian, such as exp(cx 2 + dx) or exp( nm c nm x n x m ), using N m layers of MPO's. This technique is discussed in Sect. 5.5 when we construct MPS.
For other formulas, we find other alternative approaches. One is to find an MPO representation of the operatorÔ f . This is expected to be a simple tensor contraction for smooth functions, but the cost of computing it in general may be significant, unless we combine it with strategies such as interpolation (cf. Sect. 5.2.2). Another alternative is to approximate f (x) in Taylor series and use multiplications by constants and by x, combined with simplification stages, to approximate the product. A final approach is to use algorithms from quantum computing to implement the equivalent transformations [4,32]. The interesting thing in this approach is that, since we are working with artificial registers, we can forego many of the ancillas, imposing success in all operations.

Multivariate analysis on the MPS quantum register
After this brief review of MPS quantum register, we will now introduce specific, higher level algorithms for interpolation, Fourier analysis, differentiation, solving PDEs and creating MPS for a sampled function. Unless otherwise stated, all algorithms are implemented using the linear MPS encoding from equations (3) and (25). Whenever possible we provide also a comparison with state-of-the-art classical algorithms, with regard to efficiency, accuracy and numerical stability.

Fourier transforms and frequency analysis
Assume a function p(x s ) uniformly sampled over 2 m points in an interval x s ∈ [a, b] of size L = |b − a|. We can define 2 m distinct Fourier modesψ r (x s ) ∝ exp(iq r x s ) for q r = 2πr/L and r = 0, 1, . . . 2 m . We can also define a Fourier transform which returns the expansion of the original function in the Fourier basis.
The Quantum Fourier Transform (QFT) is the extension of the classical operator F to the quantum realm. The QFT operatorF takes a quantum state with wavefunction p(s) and converts it into its Fourier transform The QFT is defined by the action on the quantum register stateŝ  This operation can be implemented using m+1 layers of quantum gates, as F = U flip F m · · · F 1 . The last operation, U flip just reverses the order of the qubits in the quantum register. The i-th layer F i contains a Hadamard gate H i on the i-th qubit, followed by m − i controlled rotations with respect to the following m − i qubits Each of the F i is an MPO with bond dimension χ = 2 and tensors B s n ,sn Applying the QFT on an MPS requires contracting the N m operators F i and simplifying the resulting states. As in other algorithms [cf. Table 1], the simplification dominates the asymptotic cost and the performance is strongly dependent on the amount of entanglement in the transformed states. For smooth functions we expect that the entanglement will be bounded, both because of the estimates above, and because smooth functions will also tend to be concentrated in Fourier space. In that case we expect an overall scaling of time as O(m 2 N 2 2 3N ), which is exponentially faster than the classical FFT. It is interesting to note that other authors have considered using the QFT in a classical context [19,27], but they never observed a real speed-up-more like a 20-fold slow down-because of working with the complete wavefunction and not with the MPS quantum register.
As example, take the Gaussian probability distribution from Fig. 7a. Its QFT is a highly concentrated state, another Gaussian in momentum space shown in Fig. 7b. We know that the entanglement entropy of the transformed state is upper-bounded by 1 e-bit (see Sect. 3.4). This is what we see not only for the final state-orange line in Fig. 7c-, but also when we analyze all stages of the transform F 1 ψ, F 2 F 1 ψ, etc [cf. Fig. 7d.] Let us inspect more carefully the Fourier transformed wavefunction Fig. 7b. Note how the wavefunction concentrates on both sides of the interval, s 0 and s ∼ 2 m − 1. This is caused by the mapping from the non-negative quantum register states s, to the actual momenta, and which has the representation k = It is useful to implement a two's complement operation that flips all qubits conditioned on the state of the sign qubit, U 2c = |s 1 , s 2 , . . . , s m → |s 1 , s 1 ⊕ s 2 , . . . , s 1 ⊕ s m . When we apply this operator to the Fourier transformed state, U 2c Fψ, we find that the amount of entanglement surprisingly drops down, almost close to zero for all bipartitions-see Figs. 7c-d. This hints at the fact that the signed quantum register is a much better variable for describing this (and probably other) symmetric probability distributions.

Interpolation
Interpolating means approximating values of a discretized function on points that were not initially considered. We discuss two techniques for interpolating from a quantum register with m qubits to a new register and discretization with m + k qubits. The first method uses finite differences to extrapolate new points as linear combinations of previous values. The second method is a spectral technique based on Fourier transformations that, as we will show, can be exponentially more accurate for finite bandwidth functions.

Taylor expansions
Let us consider a scenario in which we have discretized the interval [a, b] uniformly with 2 m points, and we want to add 2 m extra points, moving to m + 1 qubits. Let us call x s the original variables and x r the new sampling, which satisfies We assume that the values at x s determine those at x 2s and we only need to extrapolate the values at the odd sites x 2s+1 . In order to do so, we can assume that our function is analytic and admits a Taylor expansion to some finite order, which gives the following approximation x s [ p](s) We can translate equations (41) and (42) into an algorithm that extends a sampling with m qubits into another one with m + 1. In a way that resembles very much the Grover and Rudolph protocol, but which is definitely not unitary, we add one least significant qubit at the end The ladder operatorsŜ ± increase or decrease the quantum register by one, displacing the function we encodedŜ Instead of using reversible operations (Toffoli, CNOT), we implement the MPO as a smaller, irreversible and classical circuit with one carry bit that propagates through the bond-dimension. This requires a single tensor Ĉ The tensor C s ,s a,b is nonzero only for s = s ⊕ b, and a = s ∧ b. The MPOŜ − is obtained by simply exchanging the indices aŝ and the linear combination algorithm (43) can be implemented using the simplification techniques from Sect. 4.4.

Fourier interpolation
The problem with linear interpolation is that the accuracy of the approximation is constrained by the initial sampling, O(δ 2 m ). This seems to contradict the Nyquist-Shannon theorem, according to which a function with a bounded spectrum only needs to be sampled with a frequency of 2 × ν max for a perfect interpolation. Take for instance a Gaussian probability distribution (15) with width σ, which has been sampled in the interval [−8σ, 8σ]. Its Fourier transform is a normal distribution with center ν = 0 and width 1/σ. We can say that the information beyond ν = 4/σ, is exponentially suppressed. According to Nyquist's theorem, we can exactly reconstruct a Gaussian function by sampling it with period σ/8. For the conditions above, that means 8/σ × 16σ ∼ 128 points stored in 7 qubits. However, if we attempt linear interpolation, we typically will make a bounded error that is fixed by the initial sampling, O(δx ∼ 1/16).
A well known solution is to do the interpolation in frequency space. When we perform a discrete Fourier transform, we are decomposing the sampled function p(x s ) as a sum of discrete Fourier modes ψ(ks ) ∝ exp(iks x s ). We can use this to reconstruct a continuous approximation to the original function p(x) in what is known as Fourier interpolation This continuous approximation can then be resampled with as fine a grid as needed.
We can implement the approximation and re-sampling very efficiently. If we wish to enlarge the number of qubits from m to m + k, we compute We start with a Fourier transform and two's complement over m qubits. We then insert k qubits in positions 2 to k + 1. These are bits that encode very high frequencies and which are populated with zeros, as we do not need any finer details in the sampled function. We finally take the enlarged register and invert both the two's complement and the Fourier transform. As illustrated in Fig. 8 this is a powerful technique that can reconstruct a Gaussian using 2 10 = 1024 points out of a discretized Gaussian with 2 5 = 128 points, with negligible error.

Differentiation
In numerical analysis, there are two main ways to estimate the derivative of a discretized function. The first method is called finite differences, because it relies on linear combinations of the function p(x s ) and its displacements p(x s ± nδx). The second type of methods is called a spectral method, because it works with the Fourier expansion from (47). Both methods have simple translations to the language of MPO's and MPS's.

Finite differences
We will work out this technique by example. Our starting point are two standard finitedifference approximations to the spatial derivatives, with different degrees of approximation The step δx will be the discretization (b − a)/2 m of the uniform sampling (1). Small changes x s ± δx map to increments and decrements of the variable s and using the ladder operators (44). Our finite-difference formulas become These formulas can be implemented as MPO's of bond dimension 3, which can be efficiently contracted and simplified using the algorithm from Sect. 4.4. Note that hhigher order approximations are also possible using the same operator or powers of it. Roughly, the bond dimension of the MPO that implements a finite difference formula grows linearly with the order of the approximation, just like the number of non-diagonals in its matrix representation. This is a moderate cost that makes differentiation an approachable routine in higher level algorithms.

Fourier approximations to derivatives
Since Fourier interpolation works so well, we can use it to approximate the action derivatives at arbitrarily high orders. The action of a general differential operator G(∂ x ) on the interpolated state (47) is very simple If we resample this function with m bits, we obtain a representation of the differential operator in terms of the momentum operator (38)

Solving partial differential equations
One of the main applications of all these techniques is the study of how multivariate functions evolve in time, when subject to one of many partial differential equations. We focus our discussion on the Fokker-Planck equation in one dimension using uniform drift and diffusion, µ(x, t) = µ and D(x, t) = D. This equation governs the evolution of probability distributions for random variables undergoing a Wiener process in the Îto representation-a recurrent problem in quantum optics and finance, for instance. Moreover, this is already a challenging toy model from numerical analysis that is subject to numerical instabilities and demands state-of-the-art integration techniques. We will provide two techniques to solve this equation using the MPS quantum register approach, matching the two methods to work with differential operators from Sects. 5.3.1 and 5.3.2.

Finite differences
Let us assume that we have the MPS representation of the initial value p(x, 0) and we need to estimate the evolution of this probabilty distribution at later times, p(x, t). Following Sect. 5.3.1, we write a finite-difference approximation to the Fokker-Planck model, replacing derivative operators with ladder operators. This transforms the Fokker-Planck equation into a first order differential equation generated by the linear operatorĜ, whose action we approximate with an MPO. For concreteness, we will use the combined first-second order approximation from (51) We also need to build an implicit integration method that works around the fact that equation (54) is not unitary and has the potential to develop exponentially growing numerical instabilities. We have chosen a second order implicit method, which translates into our integration recipe SinceĜ has a simple representation in terms of MPO's, we build our algorithm around the repetition of two elementary steps. First, compute the MPS for the product |φ 1 (t) = (1 + 1 2 δtĜ) |p(x, t) , using standard simplification techniques to get the simplest and best approximation. Second, estimate the MPS |φ 2 (t) that best approximates the equation In practice, we implement this step using a conjugate gradient method, but one could also write quadratic optimization techniques that minimize the distance between both states, as in DMRG's correction vector techniques [15,21,23]. As illustration, Fig. 9a shows limitations. The operator δtĜ ∝ 2 2m diverges as we add more and more qubits. Our simulation is therefore very limited both in time and in the spatial discretization. This intrinsic instability is only partially cured with the implicit methods, but not always and also not for very long simulations. If we wish to have higher precision and numerical stability, we need to develop slightly better techniques, described in the following section.

Spectral split-step method
The spectral methods, and in particular the Fourier transform and the split-step method techniques [31], have been traditionally used in many nonlinear Optics and quantum mechanical problems, due to their efficiency, stability and accuracy. The method is optimally designed to solve equations of the form where G(∂ x ) is a function of the differential operator ∂ x and the coordinates are defined over a regular interval. It works by moving to Fourier space, where the generator of the evolution is a function of the momentum k Herep is the Fourier transform of the original function,p =Fp, over the real line. In our discrete scenario with uniformly sampled, regular intervals, the spectral method can be implemented using the Quantum Fourier Transform (34) and the interpolation techniques from Sect. 5.3.2. The approximate solution to (60) is expressed as We have implemented the recipe (62) using MPS and MPO's. We representF,F −1 and exp(−µk + Dk 2 ) using 3(m + 1) MPO's of bond dimension 2. Provided that states remain weakly entangled, we can exactly solve the time evolution, without any numerical instabilities or truncation error, for any time and coefficients. As illustration, Fig. 10 reproduces once more the solution of the Fokker-Planck equation for a state that is initially Gaussian with variance σ(0) = 1 and centerx(0) = 0, living in an interval [−10σ, 10σ] discretized with m = 14 bits (16384 points in spacee). We can afford larger number of qubits and longer times than in the finite difference method, because the algorithm is orders of magnitude more efficient and very stable numerically. Notice also how the algorithm implements periodic boundaries by default; this is a feature that is very useful to avoid boundary reflections and simulate in small intervals the dynamics of propagating fields.
This recipe can be extended to problems that include dependencies on both the spatial derivatives and the spatial coordinates. We explicitely refer to equations of the form or higher-dimensional equivalents. The solution to this problem is no longer exact, but relies on a Trotter-Suzuki expansion of the generator This technique is known as a split-step method because it combines evolution steps in real space, with stages that are implemented in Fourier space [31]. All ingredients in this formula (64) are well known from earlier pages.

Construction of the MPS
The algorithms from previous sections assume that we already have an MPS wavefunction to begin with. We will now discuss different ways in which to obtain such initial conditions.

Exact discretization
The trivial way to construct an MPS representation is to begin with the discretized function p (m) (s 1 , s 2 . . . s m ) and perform a sequential Schmidt decomposition. This has an exponential cost O(2 3m/2 ) and fails when the number of qubits goes above 30-something, due to memory and time constraints, as we have seen in Sect. 3.5.

Explicit formulas
As discussed in Sect. 4.5, we could write the whole MPS using the exponential of a QUBO formula. In practice, this only works well when the state Σ is not very squeezed 2 . More generally, we need to reconstruct the whole probability state a progressive refinement of the uniform distribution with each step implemented by an MPO that is a better behaved Gaussian function. We have used this technique to reconstruct the MPS representations of the 2D (21) and 3D (22) Gaussian states. Figs. 3b and 4a show the maximum entanglement entropy over all bipartitions of the resulting MPS state, while Table 2 discusses the amount of memory and the size of the tensors. (67) is nothing but a discrete approximation to the imaginary time evolution of the unnormalized state |p u

Imaginary time evolution Equation
with uniform initial condition |p u (β = 0) = |e . This equation can be solved for any partition function distribution (65) generated by an operatorĤ with an efficient MPO representation. This approach is reminiscent of how thermal state density matrix have been simulated using the MPS formalism [28], but it is much simpler, since we do not introduce any auxiliary purification degrees of freedom.

Machine learning The challenge of writing an MPS representation of a function p(x)
is comparable to the challenge of constructing that function in a quantum computer: there are efficient protocols, but we do not know them a priori. The quantum computing field has developed different strategies to address this problem. A promising one is to approximate the state |p using a parameterized quantum circuit that is trained using the techniques of machine learning, such as generative adversarial networks [34]. These techniques can be extended to our domain. Instead of using a quantum-classical approach with a circuit that generates the probabilities and a neural network that discriminates, we can use the MPS as generator and apply similar training techniques.

Porting back to the quantum computers
This work has introduced various classical MPS algorithms that use the tools of quantum computing-quantum registers, function encodings, quantum gates and algorithms-to solve efficiently various numerical analysis tasks. These techniques, and the whole line of research, can feed back to the world of quantum computing. For instance, the Fourier interpolation algorithm from Sect. 5.2.2 has an immediate translation to a real quantum computer.
Other algorithms require some fine tuning. In linear interpolation and finite differences, the theŜ ± operator was designed to use irreversible arithmetic. However, we know that a similar operator can be implemented or approximated using ancillas and reversible arithmetic-half and full adders-, or other ideas from quantum simulation [5].
Something similar happens in the case of PDE solvers. Our discussion has focused on the non-unitary evolution induced by Fokker-Planck techniques, studying the implementation of the non-unitary operator exp[G(k)]. Algorithms (62) and (64) can be trivially generalized to solve Schrödinger equations such as The generator of this equation is anti-Hermitian and both exp[δt G(ik)] and exp[−iV (x)δt/ ] can be implemented as a unitary gate in the quantum register. As in the MPS case, the scaling of the algorithm is problem dependent. The exponential exp[iαk 2 ] can be implemented with O(N 2 m 2 ) steps, but the exponential of V (x) may have a more complicated scaling, strongly dependent on function arithmetic. However, we expect that the smoothness of usual potentials will also lead to simple approximations with quasi-local gates.

Discussion and outlook
This work has presented many numerical algorithms for constructing, manipulating and interrogating multivariate functions in MPS-encoded quantum registers. We have shown that, heuristically, the renormalization provided by the quantum register representation is key to the creation of states with low entanglement, capable of encoding smooth, differentiable functions. The use of tensor network states is a modern development in numerical analysis [3,11]. In the field of low-rank tensor approximations, a multivariate function ψ(x 1 , x 2 . . . x N ) is approximated by a contraction of tensors, labeled either by the continuous variables x i , or discretized versions of them x s i . In other works, the spatial degrees of freedom are replaced by labels in some local mode expansion [14]. Nevertheless, all these approaches preserve a notion of local degrees of freedom that translates into the tensor structure. In this work we are instead using a quantum register discretization, where each local coordinate is exploded into m qubits of information, each of them probing the function at a different length scale. As we have seen in Sect. 3.5, these qubits need not be kept together, and there may be more efficient renormalization schemes when we group common length scales that are more strongly correlated. This implicit renormalization scheme is a discovery of the quantum computing community that has not been sufficiently exploited so far and which may empower many recent developments in the field, specially in the field of quantum finance and probability distribution analysis [20,22,33].
The algorithms in this paper bend themselves to a broad family of problems which have been considered in the quantum computing world. We have illustrated the solution of time-dependent partial differential equations, but the same techniques can be extended to stationary problems. This way, the MPS quantum register becomes a natural tool to solve the Poisson equation [5], the wave equation [1], the fluid equations [27], or even the Schrödinger equation itself. This implies not just abstract, fundamental studies in Physics, but practical applications in fields such as aerodynamics or finance. We expect new applications of quantum-inspired finance that reach beyond the state-of-the-art [20], providing new schemes for evaluating financial products [22], performing risk analysis [33] and even more sophisticated time-dependent simulations and tracking.
Let us also remark that the algorithms developed in this work are of a heuristic nature. All methods and techniques in Sections 4 and 5 assume states with low entanglement. This approximation is bound to break at some point, either because of functions with broad spectra and complex structure, or because of increased dimensionality. Quantum computers become a valuable tool that still has an exponential advantage over classical algorithms, and which may profit from the ideas and developments associated to MPS quantum register techniques, as discussed before [cf. Sect. 5.6].