Decomposition of matrix product states into shallow quantum circuits

TNs are low-rank factorizations of large tensors into networks of smaller tensors, which can represent wave functions and quantum circuits on classical hardware [1, 32]. TNs are typically illustrated as undirected graphs, whose nodes represent 'core' tensors containing the parameters of the model, and whose edges indicate contractions between core tensors along particular tensor axes. TN model families are defined by their underlying graph structures, with 1D line graphs giving MPS [33], 2D grid graphs giving projected entangled pair states (PEPS) [34], and tree graphs (or equivalently, acyclic graphs) giving tree tensor networks (TTNs) [26]. The main capacity parameter of TNs is the dimension of the vector spaces being contracted over, known as the bond dimension χ. This parameter controls the expressive power of the model, as well as the computational resources and memory requirements for classical simulations. The bond dimensions needed to encode a target wavefunction can be estimated by the entanglement entropy within the target state [35], with states possessing entanglement entropy S requiring TNs of bond dimensions $\chi \sim e^S$ to exactly represent. MPS implementations typically utilize bond dimensions of at most $\chi \sim 10^3$, although high-performance simulations using MPS are able to reach bond dimensions of $\chi \sim 10^4 - 10^5$ [36].

TN simulations of quantum circuits can be very effective when the graph structure of the TN matches the circuit layout [17, 27], with the singular value decomposition (SVD) enabling efficient approximate simulations by dropping lower magnitude singular values in the Schmidt decomposition of the target state [35]. As a concrete example, TN simulations in [3] were able to achieve higher overall state fidelity when simulating the multi-layer 2D circuits of the so-called quantum supremacy experiments in [37], where the fidelity deteriorated due to quantum hardware errors. On the other hand, computational cost of classical simulations can quickly escalate when the TN graph and quantum circuit layouts differ, such that circuit gates act on qubits that are far apart on the TN graph. In such cases, the bond dimensions necessary for a high-fidelity simulation will generally increase exponentially with the system size, owing to the need to carry entanglement between distant qubits through all intermediate edges.

While simulating quantum circuits on classical hardware via representations as TNs has many fundamental use-cases, our focus here is on the reverse direction, of efficiently representing TNs as the output of quantum circuits. An important tool for carrying out this latter conversion are TN canonical forms, which facilitate the loss-less conversion of all core tensors in a TN into isometries (i.e. inner product preserving transformations between Hilbert spaces) [38]. While weaker canonical forms have been developed for general TNs [39], the ability to efficiently convert core tensors into isometries requires TN architectures to be defined on acyclic graphs, such as MPS and TTNs.

By restricting to bond dimensions which are powers of 2, the isometries arising in a TN in canonical form are equivalent to matrices of shape $2^n\times 2^m$ (assuming $n \unicode{x2A7E} m$ without loss of generality), which can always be expanded to $2^n \times 2^n$ unitary matrices (i.e. n-qubit unitaries). In more detail, if Q is an isometry of shape $2^n\times 2^m$ satisfying $Q^{\dagger}Q = \unicode{x1D7D9}_{2^m \times 2^m}$, the construction of the corresponding n-qubit unitary U involves extending the 2^m columns of Q with $\kappa : = 2^n - 2^m$ basis vectors from the kernel space of $Q^{\dagger}$. These additional basis vectors can be arranged in a $2^n \times \kappa$ matrix X satisfying $Q^{\dagger}X = 0$ and $X^{\dagger} X = \unicode{x1D7D9}_{\kappa \times \kappa}$, which guarantees that the $2^n\times 2^n$ square matrix $U = [Q \quad X]$ satisfies the unitary condition $U^{{\,}\dagger}U = U U^{{\,}\dagger} = \unicode{x1D7D9}_{2^n \times 2^n}$.

For a general MPS parametrizing an N-qubit wavefunction, each tensor core will have left and right bond dimensions of $\chi_{\text{left}}$ and $\chi_{\text{right}}$, which, by padding with zeros, can be expanded to the nearest powers-of-two, namely $\chi_{\text{left}}^{\prime} = 2^n$ and $\chi_{\text{right}}^{\prime} = 2^m$, where $n = \lceil\log_2(\chi_{\text{left}})\rceil$, $m = \lceil\log_2(\chi_{\text{right}})\rceil$, and $\lceil \cdot \rceil$ denotes the ceiling operation. Placing this MPS in (left) canonical form via iterated QR decompositions (see [33] for details) will then yield a $2^{n+1} \times 2^m$ isometry, which can then be converted into a unitary gate acting on n + 1 qubits. The linear topology of the MPS leads the composition of gates from this conversion process to form one layer of multi-qubit unitaries arranged in a staircase topology, which we refer to as a linear circuit layer. Figure 1 illustrates this exact conversion from MPS to PQC for an MPS with ${\chi_{\text{max}}} = 8$.

Zoom In Zoom Out Reset image size

Despite the simplicity of the conversion process described above, a major deficit of this procedure is the fact that arbitrary n + 1 qubit unitary gates cannot be directly implemented in quantum circuits running on real quantum hardware, which typically only support gates acting on one or two qubits. Overcoming this limitation while still obtaining an efficient and accurate conversion from TNs to PQCs is a primary motivation for our work.

The problem of efficiently representing general MPS as the output of PQCs is crucial for many applications in condensed matter physics and quantum machine learning (QML) which utilize both classical and quantum computational resources. In the former domain, [29, 40–42] explored the use of PQCs as quantum variational ansätze for calculating diverse properties of complex many-body systems, and identified distinct benefits of this method relative to the use of TNs on their own. Within QML, [5] introduced the idea of a joint optimization framework utilizing both TNs and PQCs, which the authors predicted would give improved performance in QML tasks. While this proposal lacked a concrete method for transferring TNs solutions to quantum computers, the PQC pretraining method of [43] gave numerical verification of the benefits of this style of joint training for MPS with ${\chi_{\text{max}}} = 2$ (which exactly decompose into a single linear layer of two-qubit unitaries). The recent work of [10], which utilizes the conversion method described here, illustrates the increased boost in QML task performance that results from the use of larger TN models for initializing PQCs. [44–46] represent other approaches which combine TN and PQC models to solve QML tasks.

At the level of specific algorithms, [28] proposed an analytic technique for decomposing MPS of arbitrary ${\chi_{\text{max}}}$ into multiple linear layers of two-qubit unitaries, where the unitaries for each layer are chosen based on truncations of the MPS to χ = 2 across each bond, determined through iterated SVDs (see also the earlier work of [47]). This has the advantage of requiring relatively small computational resources, but unfortunately leads to relatively small gains with increasing circuit layers, decreasing fidelity past a critical number of layers, and poor performance compared to more sophisticated approaches. In a different direction, [29] proposed decomposing arbitrary MPS into circuits of stacked linear layers via iterative optimization over the two-qubit unitaries forming that circuit. This optimization method is similar to the Evenbly-Vidal algorithm [48], as well as the procedure described in section 3.2, and can support more general circuit architectures and achieve better performance via greater numbers of optimization sweeps. However, without a high-quality initialization for the circuit gates, this method remains vulnerable to barren plateaus [49], and becomes expensive when large numbers of sweeps are needed.

Closer to our work is the recent [30], which utilizes both analytical and optimization-based decomposition methods to efficiently approximate arbitrary MPS using PQCs. This method first initializes the gates of a circuit using the analytical method, before using iterated optimization of the unitary gates to improve fidelity with the target MPS, a combination that was shown to have significant performance benefits over that of [28]. This method is only one of many decomposition protocols we benchmark here (denoted by $D_{all} O_{all}$ in the language of section 3.3), and while we confirm the clear advantages of combining analytical and optimization methods, our empirical findings show that there are yet better ways of utilizing these two styles of decompositions to approximate arbitrary MPS using quantum circuits.

[28] presents a technique to sequentially decompose an MPS into K layers of two-qubit unitaries with a linear next-neighbor topology. The approach is aimed at maximizing the fidelity

$$\begin{align} \begin{aligned} f\left(\left\{\text{L}\left[U\right]^{\left(k\right)}\right\}_{k = 1}^K\right) & = \left|\langle 0^{\otimes N} |\prod_{k = 1}^K \textrm{L}\left[U\right]^{\left(k\right)\dagger}| \psi_{{\chi_{\text{max}}}} \rangle\right|\\ & = \left|\langle 0^{\otimes N} |\text{L}\left[U\right]^{\left(K\right)\dagger}\dots\textrm{L}\left[U\right]^{\left(1\right)\dagger}| \psi_{{\chi_{\text{max}}}} \rangle\right|\\ & = \left|\langle \psi_{QC}^{\left(K\right)} | \psi_{{\chi_{\text{max}}}} \rangle\right| \end{aligned} \end{align} \tag{ 1 }$$

between the original MPS state $| \psi_{{\chi_{\text{max}}}} \rangle$, and a quantum state

$$\begin{align} \begin{aligned} | \psi_{QC}^{\left(K\right)} \rangle & = \prod_{k = K}^1 \textrm{L}\left[U\right]^{\left(k\right)}| 0^{\otimes N} \rangle\\ & = \textrm{L}\left[U\right]^{\left(1\right)}\dots\textrm{L}\left[U\right]^{\left(K\right)}| 0^{\otimes N} \rangle \end{aligned} \end{align} \tag{ 2 }$$

which is meant to approximate the MPS. The resulting quantum circuit is $\prod_{k = K}^1 \textrm{L}[U]^{(k)}$, and it consists only of two-qubit U(4) unitaries. Equivalently, one could view the inverse layers $\prod_{k = 1}^K \textrm{L}[U]^{(k)\dagger}$ as the circuit that disentangles the MPS. To avoid ambiguity, we will at times refer to $| \psi_{{\chi_{\text{max}}}} \rangle$ as the 'target state' when it is being viewed as a target for approximation, and as the 'partially-decomposed state' when it is being disentangled by circuit layers.

Algorithm 1 describes the analytic decomposition technique, which is also illustrated in figure 2(top).

Zoom In Zoom Out Reset image size

Algorithm 1. Analytic Decomposition.
Input: MPS $\psi_{{\chi_{\text{max}}}}$, Maximum layers K, Target fidelity $\hat{f}$
Output: Quantum Circuit layers $\prod_{k = K}^1 \textrm{L}[U]^{(k)}$
1: $k \gets 0$
2: $| \psi^{(0)} \rangle \gets | \psi_{{\chi_{\text{max}}}} \rangle$
3: while k < K or $|\langle 0^{\otimes N} | \psi^{(k)} \rangle| \lt \hat{f}$ do
4:  Truncate $\psi^{(k)}$ to $\psi^{(k)}_{\chi = 2}$ via SVD
5:  Convert $\psi^{(k)}_{\chi = 2}$ to $\text{L}[U]^{(k+1)}$
6:  $|\psi^{(k+1)}\rangle \gets \textrm{L}[U]^{(k+1)\dagger} |\psi^{(k)}\rangle$
7:  $k \gets k+1$
8: end while

In this algorithm, one truncates a copy of the original MPS $\psi_{{\chi_{\text{max}}}}$ to χ = 2 via SVD, converts the truncated MPS $\psi_{\chi = 2}$ to one layer L$[U]$ of two-qubit unitaries, and applies the inverse of that layer L$[U]^{\dagger}$ to the original MPS. The resulting partially-disentangled MPS is expected to have less entanglement and the bond dimensions may decrease depending on the predefined singular value thresholds. The conversion of the χ = 2 MPS to a linear layer of two-qubit unitaries is explained in section 2.2 and [28, 41, 43]. This process can be repeated iteratively to create several layers $\prod_{k = K}^1 \textrm{L}[U]^{(k)}$ which are indexed from K to 1 owing to the reversal of the disentangling layers to form a circuit for approximating the target MPS. In particular, the Kth layer of the disentangling circuit ends up being used as the (adjoint of the) first layer of the quantum circuit for approximating $| \psi_{{\chi_{\text{max}}}} \rangle$.

The analytical decomposition algorithm is relatively inexpensive to perform, and the decomposition can be improved sequentially up until a desired fidelity, with every additional layer likely giving better decomposition performance. However, it is crucial to note that for a fixed number of layers K, it is not possible to iteratively increase the quality of the decomposition. This leads to our observation that K sometimes needs to be orders of magnitudes larger for high-fidelity decomposition than predicted by the lower bound $K_\mathrm{min} \sim \log_2({\chi_{\textrm{max}}})$ [50]. For this reason, this analytical decomposition algorithm is better utilized in conjunction with an effective constrained optimization algorithm for the resulting quantum circuit unitaries.

We now describe the constrained optimization algorithm for the circuit unitaries used throughout this work. It was demonstrated in [51] for learning quantum circuits to prepare quantum states. Given an initial quantum circuit $\prod_{m = 1}^MU_m$ consisting of M unitaries U_m , the goal of the optimization algorithm is to increase the fidelity

$$\begin{align} \begin{aligned} f\left(\left\{U_m\right\}_m\right) & = \left|\langle 0^{\otimes N}|\prod_{m = M}^1 U^{{\,}\dagger}_m|\psi_{{\chi_{\text{max}}}}\rangle\right|\\ & = \left|\langle 0^{\otimes N} |U^{{\,}\dagger}_{1}\dots{U}^{\dagger}_{M}| \psi_{{\chi_{\text{max}}}} \rangle\right|\\ & = \left|\langle \psi_{QC}^{\left(M\right)} | \psi_{{\chi_{\text{max}}}} \rangle\right| \end{aligned} \end{align} \tag{ 3 }$$

between the target MPS state $| \psi_{{\chi_{\text{max}}}} \rangle$, and a quantum state

$$\begin{align} \begin{aligned} | \psi_{QC}^{\left(M\right)} \rangle & = \prod_{m = 1}^M U_{m}| 0^{\otimes N} \rangle\\ & = U_{M}\dots{U}_{1}| 0^{\otimes N} \rangle \end{aligned} \end{align} \tag{ 4 }$$

which is meant to approximate the MPS. In our case, the unitaries U_m act only on two qubits, i.e. the unitaries are U(4) matrices, and the circuit layout is a linear next-neighbor topology for efficient classical simulation using MPS, as well as for compatibility with algorithm 1 in section 3.1.

Algorithm 2 describes the optimization algorithm used in this work, which is also illustrated in figure 2(bottom).

Algorithm 2. Decomposition by Optimization.
Input: Initial quantum circuit $\prod_{m = 1}^M U_{m, 0}$, number of sweeps T, target fidelity $\hat{f}$, learning rate $r \in [0, 1]$
Output: Optimized quantum circuit $\prod_{m = 1}^M U_m$
1: $t \gets 1$
2: While t < T or $\left | \langle 0^{\otimes N}| \prod_{m = M}^1 U^{{\,}\dagger}_m |\psi_{\chi_{max}}\rangle \right | \lt \hat{f}$ do
3:   for m in 1 ... M do
4:    Calculate $\hat{\mathcal{F}}_m$ in equation(5)
5:    SVD $\hat{\mathcal{F}} = \mathcal{U}\mathcal{S}\mathcal{V^{\dagger}}$
6:    $U_{\text{new},m} \gets \mathcal{U}\mathcal{V^{\dagger}}$
7:    Apply learning rate in equation (6)
$U_{m} \gets U_{m}(U^{{\,}\dagger}_{m} U_{\text{new},m})^r$
8:   end for
9:   $t \gets t+1$
10: end While

The optimization algorithm iterates over all unitaries in the circuit and updates them one by one. The number of sweeps can either be predefined through a maximum number T, or until a target fidelity $\hat{f}$ is reached with $f\big(\{U_m\}_m\big) \gt \hat{f}$ (compare equation (3)). The update for one unitary U_m can be calculated by first calculating the so-called environment tensor for U_m :

$$\begin{align} \hat{\mathcal{F}}_m = \textrm{Tr}_{\bar{U}_m}\left[\prod_{i = M}^{m+1} U_{i}| \psi_{{\chi_{\text{max}}}} \rangle\langle 0^{\otimes N} |\prod_{j = 1}^{m-1} U^{{\,}\dagger}_{j} \right]. \end{align} \tag{ 5 }$$

Here, the expression $\text{Tr}_{\bar{U}_m}[\dots]$ denotes the partial trace over the qubit indices that do not interact with U_m . The environment tensor $\hat{\mathcal{F}}_m$ is represented as a $4\times 4$ matrix, and can in practice be calculated by 'removing' U_m from the circuit and contracting the remaining tensor network (see bottom of figure 2) while retaining the MPS structure throughout. Interestingly, $\hat{\mathcal{F}}_m$ is the operator which is locally optimal for the circuit to increase the fidelity with the MPS if it were replacing U_m [51], but this operator is not generally unitary. In order to keep all quantum circuit operations unitary, we can compute the SVD $\hat{\mathcal{F}}_m = \mathcal{U}\mathcal{S}\mathcal{V}^{\dagger}$, then use the fact that the product $\mathcal{U}\mathcal{V}^{\dagger}$ is the unitary matrix that minimizes the Frobenius distance with $\hat{\mathcal{F}}_m$ [52], which in turn maximizes the fidelity of the PQC with the MPS among all unitaries. This leads $U_{\text{new}} = \mathcal{U}\mathcal{V}^{\dagger}$ to be the locally optimal unitary for increasing the fidelity with the target MPS. Since this is a rather strong local update, we introduce a learning rate $r\in[0, 1]$, which influences the unitary update rule via

$$\begin{align} U_{m}^{^{\prime}} = U_{m}\left( U^{{\,}\dagger}_{m} U_{\text{new}} \right)^r. \end{align} \tag{ 6 }$$

We find the optimization process to give excellent performance using a value of r = 0.6. We note that the fractional matrix exponential in $\left( U^{{\,}\dagger}_{m} U_{\text{new}} \right)^r$ can be implemented by simply applying the exponent to the diagonal elements in the eigendecomposition of $U^{{\,}\dagger}_{m} U_{\text{new}}$. Finally, $U_{m}^{^{\prime}}$ is the unitary that replaces U_m , which gives improved fidelity with the MPS.

We note that it is vital to the computational feasibility of this algorithm to perform optimization sweeps consisting of consecutive values for the index m, because in that case the states used to compute $\hat{\mathcal{F}}_m$, i.e. $\prod_{i = M}^{m+1} U_{i}| \psi_{\chi_{max}} \rangle$ and $\langle 0^{\otimes N} |\prod_{j = 1}^{m-1} U^{{\,}\dagger}_{j}$, do not require full re-calculation. With a proper choice of the optimization order, most of these states can be cached and reused. For example, after $U_{m}^{^{\prime}}$ has been obtained, we apply $U_{m}^{^{\prime}\dagger}$ to $\langle 0^{\otimes N} |\prod_{j = 1}^{m-1}$ and $U^{{\,}\dagger}_{m+1}$ to $\prod_{i = M}^{m+1} U_{i}| \psi_{\chi_\mathrm {max}} \rangle$, which forwards the two sides of the TN to step m + 1. While one can analogously perform backward sweeps, we restrict ourselves to forward sweeps in this work.

The issue with decomposing an MPS with this optimization algorithm alone (what we call the brute-force method), is that when the initial circuit unitaries are random, the application of the circuit to the fully-optimized MPS state $| \psi_{{\chi_{\text{max}}}} \rangle$ can potentially increase the maximum bond dimension ${\chi_{\text{max}}}$ required to approximate the resulting quantum state without significant truncation of singular values. Therefore, we now combine algorithms 1 and 2. Let algorithm 1 provide the initial circuit guesses that tend to decrease the bond dimension of the MPS, and let algorithm 2 then optimize the unitaries to achieve the best approximation of the target MPS for a fixed number of unitaries and layers.

In sections 3.1 and 3.2 we have introduced the fundamental algorithmic building blocks for decomposing MPS into quantum circuits consisting of only two-qubit unitaries. All steps of both algorithms are to be performed on classical hardware. Both approaches on their own have severe limitations, but complement each other very nicely. To briefly summarize, while relatively inexpensive to perform, the main disadvantage of the analytic decomposition in section 3.1. is that the number K of resulting quantum circuit layers can be very large to achieve a desired fidelity, and there is no way to increase fidelity for a fixed number of layers. Disadvantages of the brute-force decomposition algorithm in section 3.2 include that random initial guesses for the unitaries may not only require many sweeps T to optimize, but the random circuit can cause the bond dimension ${\chi_{\text{max}}}$ of the partially-disentangled MPS to increase. However, optimization is a very valuable approach, because one spends classical resources to increase the TN decomposition quality without increasing circuit depth and potentially hampering the performance of the quantum circuits on noisy quantum hardware. The performance of optimization-based decomposition methods can furthermore be improved by using high-quality initializations of the gate parameters being optimized [53].

Fortunately, both algorithms can be combined into several different decomposition protocols, towards achieving the best trade-off between the required classical computing resources and the resulting quantum circuit depth. In fact, one way of understanding their interplay is that the analytical decomposition step serves as a high-quality initialization for each layer in the subsequent constrained optimization. As quantum hardware improves, one can more reliably implement deeper circuits on quantum hardware, which could in turn alleviate some of the burden of the MPS decomposition from classical hardware.

In figure 3, we illustrate the combinations of algorithms 1 and 2 that are studied in this work. The letter D stands for the analytical decomposition algorithm, i.e. algorithm 1 in section 3.1. We make use of the possibility to either decompose one layer (indicated by i), or decompose to the full circuit depth K (indicated by all). The letter I represents an alternative approach to D, where linear layers are added to the quantum circuit in the form of identity unitaries. The letter O stands for the application of the optimization algorithm, i.e. algorithm 2 in section 3.2. We either optimize the last layer that was created (indicated by i), or all layers in the circuit so far (indicated by all). The 'Iter' notation implies that there is an iterative procedure between creating new layers (either with D or I), and optimizing the newest layer or all layers so far (indicated by i or all, respectively).

Zoom In Zoom Out Reset image size

Our first protocol is D_all , which is equivalent to applying only algorithm 1. The second protocol is Iter$[D_i O_i]$, and it represents a minimal extension of D_all where only the unitaries in the layer $\textrm{L}[U]^{(k+1)}, k\lt K$ are optimized before the disentangling step $| \psi^{k+1} \rangle \gets \textrm{L}[U]^{(k+1)\dagger}| \psi^{(k)} \rangle$. This aims to ensure that $|\langle \psi^{(k)}_{\chi = 2} | \psi^{(k)} \rangle|$ is maximal in every iteration before disentangling, which the SVD only ensures approximately. The third protocol is O_all , which is equivalent to applying only algorithm 2. The circuit unitaries are initialized as $4 \times 4$ matrices with random entries from a normal distribution with zero mean. These are then converted into unitaries by selecting the unitary component Q in the QR-decomposition of the random matrix. The fourth protocol, $D_{all}O_{all}$, constitutes the first practical combination of both analytical decomposition and optimization. In this protocol, the K circuit layers are first created using algorithm 1, and then all unitaries are optimized in sweeps in further increase fidelity. The analytically decomposed layers provide a good starting point for the optimization algorithm as compared to random unitaries, and, crucially, the layers tend to already have a disentangling effect on the MPS, which reduces the bond dimension of the MPS. The fifth protocol, Iter$[I_i O_{all}]$, is a sequential protocol where new layers are iteratively added as identity operations, and and all existing layers are optimized using algorithm 2. The sixth and final protocol studied in this work is Iter$[D_i O_{all}]$, where new layers are added via one analytical decomposition step, and all existing layers are optimized using algorithm 2. The difference between the protocols Iter$[D_i O_i]$ and Iter$[D_i O_{all}]$ is that, while both maximize $|\langle \psi^{(k)} | \psi_{{\chi_{\text{max}}}} \rangle|$ at every decomposition step k, Iter$[D_i O_{i}]$ keeps layers $1\dots k-1$ fixed, however, these layers are also optimized in Iter$[D_i O_{all}]$.

The decomposition protocols studied in this work are constructed such that the protocol Iter$[D_i O_{all}]$ represents the most sophisticated combination of algorithms 1 and 2. The other protocols can be seen as reference benchmarks with different components of Iter$[D_i O_{all}]$ being stripped away—either the analytical decomposition, the optimization, or the sequential growing scheme of the circuit. As such, we compare this protocol to previous work, where D_all was proposed in [28], O_all is similar to the algorithm proposed in [29], and $D_{all} O_{all}$ is similar to the method implemented in [30].

It is crucial to the decomposition protocols presented in this work that all necessary steps can feasibly be performed on classical hardware. A rigorous derivation of the computational resources required for a practical decomposition of an MPS of interest is not within the scope of this work, but we emphasize that the use of MPS as intermediate states during the decomposition procedure allows the classical runtime to remain polynomial with respect to the circuit depth and qubit count. Therefore, in this section, we merely present a general overview of the expected complexities of algorithms 1 and 2, as well as the protocols presented in section 3.3.

For an MPS with a maximal bond dimension ${\chi_{\text{max}}}$, the memory requirements of the MPS scale like $\sim N\chi^2_{\text{max}}$ [32]. However, including the cost of bringing the MPS into canonical form, the computational complexity $\mathcal{C}_{\text{MPS}}$ is

$$\begin{align} \mathcal{C}_{\text{MPS}} \sim N\chi^3_{\text{max}}. \end{align} \tag{ 7 }$$

Now approaching our MPS decomposition protocols, the protocol $D_{all} O_{all}$, for example, requires K steps of analytical decomposition (see section 3.1), whose complexity $\mathcal{C}_{\text{Decomp}}$ scales like

$$\begin{align} \begin{aligned} \mathcal{C}_{\text{Decomp}} &\sim K \cdot \mathcal{C}_{\text{MPS}}\\ &\sim K N \chi^3_{\text{max}}. \end{aligned} \end{align} \tag{ 8 }$$

We note that each linear layer has N − 1 two-qubit unitaries, but the computational cost $\mathcal{C}_{\text{MPS}}$ already includes performing SVD on all bonds of the MPS. Then performing the T optimization sweeps over all K layers requires results in the complexity $\mathcal{C}_{\text{Opt}}$ of the optimization

$$\begin{align} \begin{aligned} \mathcal{C}_{\text{Opt}} &\sim K T \cdot \mathcal{C}_{\text{MPS}}\\ &\sim N K T \chi^3_{\text{max}}. \end{aligned} \end{align} \tag{ 9 }$$

The sequential protocols Iter$[O_i D_{i}]$, Iter$[I_i O_{all}]$ and Iter$[D_i O_{all}]$ apply T optimization sweeps not only on K layers, but on k < K layers for all previous iterations k. However, these previous layers have fewer unitaries. The total number of layers that are optimized in a protocol that concludes at K layers is $\frac{K(K+1)}{2}$. Thus, the sequential protocols (indicated by 'Iter') receive an additional scaling factor of K:

$$\begin{align} \begin{aligned} \mathcal{C}_{\text{Opt}}^{\text{Iter}} \sim N K^2 T \chi^3_{\text{max}}. \end{aligned} \end{align} \tag{ 10 }$$

The question now becomes, how many layers K are needed for a high-fidelity approximation of the MPS? A convenient property of our main proposed protocol Iter$[D_i O_{all}]$ is that the circuit can grow sequentially, and one does not have to guess the circuit depth that can feasibly be optimized on classical hardware. However, a lower-bound of $K_{\text{min}} \sim \log_2({\chi_{\text{max}}})$ is presented in [50]. We note that this bound provides a best-case circuit depth to exactly decompose a generic MPS, but exact decomposition is likely not the goal of a practical decomposition protocol. Therefore, we merely understand the scaling to provide a general sense of how deep quantum circuits consisting of two-qubit unitaries need to be to represent an MPS with maximum bond dimension ${\chi_{\text{max}}}$. In this case, we arrive at scalings of

$$\begin{align} \begin{aligned} \mathcal{C}_{\text{Decomp}} & \sim N \chi^3_{\text{max}} \log_2\left({\chi_{\text{max}}}\right)\\ \mathcal{C}_{\text{Opt}} & \sim N T \chi^3_{\text{max}} \log_2\left({\chi_{\text{max}}}\right) \\ \mathcal{C}_{\text{Opt}}^{\text{Iter}} & \sim N T \chi^3_{\text{max}} \log_2\left({\chi_{\text{max}}}\right)^2. \end{aligned} \end{align} \tag{ 11 }$$