Quantum state preparation using tensor networks

The matrix product state formalism, which emerged from the description of many-body quantum systems [21, 24], can be utilized for efficient quantum state preparation [33–35]. Although such methods are accurate, they require unitary operations acting immediately on the $(\lfloor \log(r) \rfloor + 1)$ qubits, where r is a rank (bond dimension) of the corresponding MPS [24], which is undesirable for modern quantum computers considering the limited connectivity and high gate errors. To reduce the hardware requirements, a matrix product disentangler algorithm was developed in [36], which prepares the state approximately and requires the implementation only of two-qubit gates acting on neighbor qubits. The disadvantage of this algorithm is that updating the gate parameters occurs sequentially layer-by-layer, and not on all gates simultaneously, which causes a higher approximation error [37]. Also, the quantum circuit is fixed and cannot be tuned accordingly to the given quantum hardware topology.

To facilitate the change in the architecture of the quantum circuit and simultaneously optimize all the gates in it, a projected gradient descent algorithm based on the automatic differentiation [38] was introduced in [37]. Firstly, the approximation error is minimized over unconstrained latent gates in a general matrix space, which are automatically differentiable. Secondly, the polar decomposition is applied to each latent gate independently to turn it into the actual (unitary) gate. However, the latter step is not variational, i.e. it does not take into account the gradient direction, and may produce a suboptimal solution.

In this work, we use a more natural update procedure for variational circuits aimed at MPS preparation based on the Riemannian optimization [39], where the gates are updated within the manifold of unitary matrices (see appendix B). This approach was shown to provide a higher speed and accuracy in some regimes [40] (see detailed comparison in supplemental material). Similarly to [41], we prepare vectors obtained by sampling analytical and smooth functions, However, in contrast to [41], we prepare states featuring an arbitrary rank allowing for the encoding of a wider class of functions. We approximate quantum states utilizing variational circuits with the native quantum gate set and [42] artificially restricting the connectivity between qubits to adjust the algorithm for NISQ hardware. We demonstrate that despite the limited connectivity high approximation accuracy is achieved by shallow circuits. Moreover, our method features the logarithmic complexity in the vector size and, in combination with the simple circuit structure, can be easily implemented on modern quantum hardware, to be applied to a variety of problems (see appendix C).

The matrix product state is a compressed representation of a vector that was originally introduced to approximate the ground state of a many-body Hamiltonian in quantum theory. The main characteristic of such a representation is the rank (bond dimension) r, which expresses the correlations and entanglement of the state vector. In fact, using singular value decomposition (SVD) [43] or cross-approximation [44], any 2ⁿ vector can be represented as MPS with the number of elements $2n r^2$. In the worst scenario, the rank r depends exponentially on n. However, considering vectors arising from sampling of analytical functions on an discretized interval, often a small and n-independent rank is sufficient to approximate the function with high accuracy [41]. Moreover, there is an exact MPS-decomposition for trigonometric, exponential, polynomial and rational functions [45]. That is why the idea outlined in this work is not to consider the preparation of arbitrary states, but vectors that can be represented as MPS with a small rank. It is important to note that the set of such states, although it may seem small, nevertheless covers a wide range of the applied problems (see Appendix C).

The key quality of the MPS is weak entanglement. Indeed, considering the entropy of the reduced subsystem $\rho = Tr_2{[\vert \psi\rangle \langle \psi\vert]}$ of an MPS $\vert \psi\rangle$ with rank r, the following bound on entropy $S(\rho)$ holds: $S(\rho) \leqslant \log{r}$ [21]. Since such states feature a small rank, hence weak entanglement, we expect that they can be efficiently approximated by a circuit where only neighbor qubits are coupled and the depth is restricted. Such circuits possess limited entanglement and are similar to circuits specifically generating states belonging to the MPS class [34, 42].

In order to perform the MPS approximation using a quantum circuit, we represent the quantum circuit also as a tensor network allowing us to optimize the gate parameters by processing tensors directly. Leveraging the locality of operations in hardware-efficient circuits, we apply the tensor networks method for circuit simulation and optimization [46, 47]. Such a method allows us to process the circuit with polynomial runtime scaling in the number of qubits.

We emphasize that the learning of the circuit parameters for the state preparation does not need to be run on quantum hardware—the parameters can be learned using tensor networks on a classical computer, which are subsequently embedded in a circuit executed on a quantum processor.

Here, we provide a step-by-step description of the protocol, which is depicted in figure 1. In order to prepare a quantum state $\vert \phi_f\rangle$ we utilize a variational circuit, which is represented as an action of a variational unitary $\mathfrak{U}(\mathbf{a})$ over a state $\vert \psi_\mathrm{in}\rangle = \vert 0\rangle_n$. The problem of quantum state preparation is finding a set of parameters a such that the output state $\mathfrak{U}(\mathbf{a})\vert \psi_\mathrm{in}\rangle$ approximates the desired $\vert \phi_f\rangle$:

$$\begin{align} \max\limits_{\mathbf{a}} | \langle \phi_f\vert \mathfrak{U}(\mathbf{a})\vert \psi_\mathrm{in}\rangle |^2. \end{align} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

2.2.1. Approximation of a function by an MPS: steps (a) to (d)

The desired state $\vert \phi_f\rangle$ is produced from an analytical function f(x), which is assumed to have a small MPS rank. One way to encode a function into a quantum state is the amplitude encoding. The function is sampled on a given interval $[x_l, \hspace{3pt} x_r]$ at 2ⁿ equispaced points $x_0,\ldots,x_{2^n-1}$. Then a normalized vector of function values at these points is encoded in the quantum state amplitudes:

$$\begin{equation*} \vert \phi_f\rangle = \sum_{k=0}^{2^n-1} \tilde{f}_k\vert k\rangle_n, \quad \tilde{f}_k = \frac{f(x_k)}{Z}, \quad Z = \sum_{k=0}^{2^n-1} |f(x_k)|^2. \end{equation*}$$

Further, we represent the obtained normalized vector $\vert \phi_f\rangle$ in the MPS format. For this, the vector is reshaped into a tensor of an order n and then sequential SVD applications allow us to represent it as an MPS [21, 48, 49]. However, the disadvantage of this approach is that it requires 2ⁿ function calls, which is impractical in case of large vectors. Alternatively, the tensor-train cross approximation technique [44] allows one to recover an MPS by an adaptive sampling of just $O(nr^2)$ points. In this work, we utilize the cross approximation algorithm.

2.2.2. Processing a circuit as a tensor network: step (e)

In order to prepare states using NISQ devices we consider a unitary $\mathfrak{U}(\mathbf{a})$ consisting of sequentially applied single-qubit variational gates $u_l^k$ mixed with CNOTs, as shown in figure 2.

Zoom In Zoom Out Reset image size

As we mentioned in the previous section, the quantum circuit itself should be represented as a tensor network, with gates acting as tensors. The projection of the desired state on the state generated by this circuit $\langle \phi_f\vert \mathfrak{U}(\mathbf{a})\vert \psi_\mathrm{in}\rangle$ is equivalent to the contraction of two tensor networks: $\vert \phi_f\rangle$ as an MPS and utilized circuit as a tensor network. Weak entanglement in shallow circuits with fixed depth, considered here, leads to the linear scaling of the simulation runtime with respect to the number of qubits. To achieve this, we utilize hyper-optimized tensor network contraction introduced in [46] that allows us to find a near-optimal contraction sequence. We show the numerical runtime of the projection of a random MPS $\vert \phi_f\rangle$ with rank 2 onto a quantum circuit with different number of layers up to 300 qubits in figure 3(a).

Zoom In Zoom Out Reset image size

2.2.3. Optimization of the quantum circuit: step (f)

One of the most commonly used ways to solve the optimization problem equation (1) is to use the gradient-based parameter-shift rule [50]. Usually in that case, variational gates $u_l^k$ depend on the parameters $\theta_l^k$ and have the form of $\exp{(-i\theta_l^k \sigma)}$, with σ being a Pauli operator. However, such a direct optimization of $\theta_l^k$ may suffer from a barren plateau of vanishing gradients due to the nonlinearity of the exponential function.

A more general initial circuit can be constructed from arbitrary unitary operators as variational gates. For this purpose, the Riemannian optimization on the Stiefel manifold of $m_1 \times m_2$ isometric matrices $V_{m_1, m_2}$ can be used [51], where $m_1=2,~m_2=2$ in case of single-qubit gates.

In order to compare those methods, we optimize a 5-qubit circuit with 2 layers of the said structure for preparing a random MPS with rank 2 and show how the approximation error decreases with the number of iterations of gradient descent in figure 3(b). The shift rule achieves only an approximation error worse than 10⁻², while the Riemannian optimization provides an error below 10⁻⁴. Since the Riemannian optimization provides a better accuracy compared to the parameter-shift rule, we proceed with using the former.

For convenience, let us reformulate the fidelity maximization problem equation (1) into a minimization problem of infidelity (or approximation error)

$$\begin{align} C(u_l^k) = 1 - | \langle \phi\vert \mathfrak{U}(u_l^k )\vert \psi_\mathrm{in}\rangle |^2 \rightarrow \min, \end{align} \tag{ 2 }$$

for a unitary $2\times2$ matrix $u_l^k$ for each $l = 1 {\ldots} L$ and $k = 1 {\ldots} 2n$. In order to find a local minimum of some differentiable function $C(\mathbf{u})$, where $\mathbf{u}$ is $m_1 \times m_2$ isometric matrix, the gradient descent algorithm can be applied

$$\begin{align} \mathbf{u}_{t+1} = \mathbf{u}_t - \alpha \nabla C(\mathbf{u}_t), \end{align} \tag{ 3 }$$

where $\mathbf{u}_t$ is a matrix $\mathbf{u}$ on step t, and α is a step size. Here, we consider the case of single unitary matrix dependence, however this analysis is valid in case of Cartesian product of such matricies [51].

The optimization problem in equation (2) must be constrained to $\mathbf{u} \in V_{m_1, m_2}$. This prevents us from using a conventional gradient descent algorithm—we cannot guarantee $\mathbf{u}_{t+1} \in V_{m_1, m_2}$ even if $\mathbf{u}_t \in V_{m_1, m_2}$. Addressing this issue, we use the Riemannian gradient and a retraction R on the manifold [51]. Firstly, instead of a standard gradient $\nabla C(\mathbf{u})$ we use the Riemannian gradient $\nabla_R C(\mathbf{u})$ which is a vector in the tangent space to the point $\mathbf{u} \in V_{m_1, m_2}$. The construction of a Riemannian gradient $\nabla_R C(\mathbf{u})$ as an orthogonal projection of the Euclidean gradient $\nabla C(\mathbf{u})$ onto the tangent space is described in the appendix B. Secondly, the next point of approximation $\mathbf{u}_{t+1}$ is determined by the SVD-retraction

$$\begin{align*} \mathbf{u}_{t+1} = R_{\mathbf{u}_t}^\mathrm{SVD} [-\alpha \nabla_R C(\mathbf{u}_t)]. \end{align*}$$

The SVD-retraction is defined as $R_X^\mathrm{SVD}(Y) = U V^{{\,}\dagger}$ after computing the SVD $X + Y = U S V^{{\,}\dagger}$. By definition, this is just a projection of the matrix $(\mathbf{u}_t -\alpha \nabla_R C(\mathbf{u}_t))$ onto the manifold of $m_1 \times m_2$ isometric matrices. Using the retraction instead of a linear step equation (3) allows for the movement along the manifold in a desired direction, satisfying the constraint $\mathbf{u}_{t} \in V_{m_1, m_2}$ automatically for each step t.

Throughout the optimization, we calculate derivatives of the function with respect to gates using the automatic differentiation technique [38]. While for large and complex networks an analytical differentiation becomes impractical, the automatic differentiation computes derivatives of tensor network programs efficiently with machine precision, in contrast to numerical differentiation. The automatic differentiation, being usually applied for computer programs, is also applicable for tensor networks since the contraction of a tensor network can be represented as a computational graph [52]. Its algorithmic complexity is theoretically guaranteed to be not greater than the algorithmic complexity of the original program, which is equivalent to the complexity of a tensor network contraction [52].

Ultimately, optimization produces a set of unitary matrices that minimize $C(\textbf{u})$. Such matrices define single-qubit gates in a quantum circuit for the approximation of the desired state.

In order to demonstrate the capabilities of our method, we consider the frequently-used analytical functions from different classes and prepare them in a many-qubit quantum state. In this work, we consider the polynomial x⁵, the power function $\sqrt{x}$, the trigonometric function $\sin{x}$, and the Gaussian probability density function. Polynomials and trigonometric functions are often encountered when considering partial differential equations [28] and two-electron integral calculation [30]. Distribution functions are needed for sampling tasks [31] and the calculation of expectation values [6], which is necessary, e.g. in finance [12]. We sample the function over N points on an equidistant grid in $[0,1]$ interval and form the function values in a N-dimensional vector, which we approximate using N amplitudes of $n = \log{N}$-qubit state.

In figure 4(a), we show the runtime of the whole state preparation routine as a function of the number of qubits. Here, we set the fidelity to $99\%$ and the number of layers to 3, which is sufficient to approximate the considered vectors having low MPS ranks (also depicted in figure 4(a)). For all studied functions, our quantum state preparation protocol features a polylogarithmic scaling in the function discretization accuracy: the runtime scales as $O(\mathrm{poly}(\log{N}))$. The number of layers can be increased which leads to an improvement in the approximation error—the infidelity decrease with the number of layers increase is presented in figure 4(b) for a 20-qubit circuit. The runtime breakdown is shown in figure 4(c). As an example, the approximation of a Gaussian density function sampled over 1024 points is depicted in figure 4(d), where a 10-qubit circuit with 3 layers was used to achieve a $99.94\%$ fidelity.

Zoom In Zoom Out Reset image size

The most significant advantage of our method is the ability to achieve a low approximation error using a shallow circuit. While deeper circuits with more parameters improve the accuracy, as shown in figure 4(b), a low error of 0.01%–0.1% can be achieved using only 5 layers. Such an error is of the same order of magnitude as an error of state-of-the-art single two-qubit gates [53].

Here, we analyze the complexity of the whole algorithm.

In order to investigate how the rank of MPS affects the approximation of our algorithm, we consider a $\sin(x^2)$ function specified on an interval $[0, b]$. As b increases, the function features faster oscillations and becomes more complex, as shown in figure 5(a), which leads to an increase of MPS ranks. In figure 5(b), we plot the approximation error for various intervals and number of layers. As expected, larger rank requires more layers and the error is larger for a given number of layers for higher ranks. As our goal is to approximate a function, we can fix the number of layers that provides certain accuracy. Generally, it is challenging to specify the exact number of layers for a given approximation as it strongly depends on the approximated function. Practically, for the functions investigated above shallow circuits are sufficient to provide decent accuracy.

Zoom In Zoom Out Reset image size

With fixed number of layers $L$, the remaining parameters are the number of qubits n and the MPS rank r, which we assume to be restricted, i.e. $r \ll 2^{n/2}$. Each step of the algorithm, then, features the following complexity.

• (i)  
  Representation of a function as an MPS is done in $O(n r^2)$ time using the cross-approximation algorithm [44]. Moreover, since this procedure is required to be done only once in the whole algorithm (in contrast to optimization), the runtime contribution of this step can often be neglected.
• (ii)  
  Circuit simulation depends linearly on the number of qubits n with fixed number of layers according to figure 3(a). At the same time, the contraction time almost does not depend on the rank r, because all the complexity mainly contributed by the contraction of the quantum circuit, since the dimensions arising from the contraction of the quantum circuit are much larger than reasonable values of ranks (supplemental material, section 3). The treewidth [54] of the used circuit is constant with a fixed number of layers (supplemental material, section 3), thus the complexity of this step is reduced to be proportional to the total number of gates, i.e. $O(n L)$.
• (iii)  
  Riemannian optimization features the complexity of $O(n L)$, where the most computationally demanding part is the SVD-retraction, which has to be done $n(2L +1)$ times for $2 \times 2$ matricies. For the fixed number of layers the runtime scaling is O(n).

Figure 4(c) shows our experimental investigation of the complexity and confirms the scaling is indeed O(n). It is clear that the most time-consuming part is the circuit simulation since it is has to be done in all iterations and requires an order of magnitude more operations than Riemannian optimization step. Thus, for further improvement of the algorithm, first of all, one should focus on this particular step.

The problem, which remains to be considered, is that hardware-efficient quantum circuits suffer from the barren plateaus problem [32]. Namely, the number of circuit parameters to be optimized grows linearly with the number of qubits, but the gradients decrease exponentially. If random values of the parameters in variational gates are chosen as an initial guess for the optimization, the initial fidelity can be less than the machine precision (${\sim}1 / 2^{n})$.

However, as was proposed in [32], selecting the appropriate starting optimization gates may help to avoid this issue. Combining the said idea with our MPS-based method, we modify the state preparation routine in the 'cut once, measure twice' way, shown in figure 6:

• (i)  
  We cut the MPS into two parts and truncate the rank at the border of the partition to 1.
• (ii)  
  For each of the resulting MPS vectors, we apply the state preparation algorithm for circuits with the same number of layers as for the initial MPS, but with fewer qubits. As a result, gates for the preparation of the 'small' MPS vectors are obtained.
• (iii)  
  We use thus obtained gates as an initial guess for the optimization of the original circuit for preparation of the initial MPS.

Zoom In Zoom Out Reset image size

If a single cut is not enough, we cut each MPS again into halves and apply this algorithm recursively.

In the numerical experiments, we have found that a single run of this procedure is sufficient to optimize the 100-qubit circuit with 3 layers for the preparation of the Gaussian density function $\mathcal{N}(0.5,0.01)$ on the interval $[0,1]$ with a $99.6\%$ approximation fidelity. At the same time, the state preparation algorithm implemented on the full 100-qubit circuit without any cuts prepares the same state with the fidelity close to zero.

The most straightforward way to apply quantum state preparation is the sampling task [31]. Our algorithm allows for the preparation of probability distribution functions, and then measurement of qubit's states provides us samples. For instance, the function of interest can be the Gaussian distribution, that has restricted MPS ranks in case of good structure of the covariance matrix (small ranks of off-diagonal blocks of the matrix) [74].

Linear systems of equations usually arise considering partial differential equations after their discretization. For instance, simple one-dimensional Poisson equation

$$\begin{align*} -{u^{^{\prime\prime}}}(x) = f(x) \quad x \in [0,1] \,\,\, u(0) = u(1), \end{align*}$$

gives the system of linear equations that can be discretized on the uniform grid

$$\begin{align*} Ax = b, \end{align*}$$

where $b_i = f(x_i)$ is the function value at the corresponding point.

Most of quantum algorithms for linear systems assumes that it is possible to efficiently prepare any desired state b [1, 3, 5]. However, there is still no general and established framework for this task and we expect that our approach facilitates the vectors encoding.

Calculation of the expectation value is one of the main routine in the computational finance, and many quantum algorithms are devoted to this problem [7, 12]. At the same time, the quantum state preparation is one of the main subroutine in these approaches, so our method is applied here as well.

For example, in the risk analysis, the Value at Risk needs to be calculated [6]. That is to find the smallest X such that:

$$\begin{align*} P[x \leqslant X] \geqslant (1 - \alpha) \end{align*}$$

for fixed α. This can be reformulated to calculation of the expectation value of the Heaviside step function $\theta(x - X)$ with various values of X. Suppose we have distribution function p(x) defined on the interval (a, b), then:

$$\begin{align*} P[x \leqslant X] = \int_a^X p(x) \mathrm{d}x = \int_a^b p(x) \theta(x - X) \mathrm{d}x= \mathbf{E} (\theta(x - X)) .\end{align*}$$

Thus, the task is reduced to finding the expectation value. If one wants to calculate it on a quantum computer one can use the amplitude estimation algorithm [75] which provides quadratically better error scaling than classical Monte Carlo. But first of all, one needs to prepare the distribution function p(x) discretized on the interval (a, b):

$$\begin{align*} \vert p\rangle = \sum_i \sqrt{p(x_i)} \vert i\rangle. \end{align*}$$