Simulating Noisy Quantum Circuits with Matrix Product Density Operators

Simulating quantum circuits with classical computers requires resources growing exponentially in terms of system size. Real quantum computer with noise, however, may be simulated polynomially with various methods considering different noise models. In this work, we simulate random quantum circuits in 1D with Matrix Product Density Operators (MPDO), for different noise models such as dephasing, depolarizing, and amplitude damping. We show that the method based on Matrix Product States (MPS) fails to approximate the noisy output quantum states for any of the noise models considered, while the MPDO method approximates them well. Compared with the method of Matrix Product Operators (MPO), the MPDO method reflects a clear physical picture of noise (with inner indices taking care of the noise simulation) and quantum entanglement (with bond indices taking care of two-qubit gate simulation). Consequently, in case of weak system noise, the resource cost of MPDO will be significantly less than that of the MPO due to a relatively small inner dimension needed for the simulation. In case of strong system noise, a relatively small bond dimension may be sufficient to simulate the noisy circuits, indicating a regime that the noise is large enough for an `easy' classical simulation. Moreover, we propose a more effective tensor updates scheme with optimal truncations for both the inner and the bond dimensions, performed after each layer of the circuit, which enjoys a canonical form of the MPDO for improving simulation accuracy. With truncated inner dimension to a maximum value $\kappa$ and bond dimension to a maximum value $\chi$, the cost of our simulation scales as $\sim ND\kappa^3\chi^3$, for an $N$-qubit circuit with depth $D$.


I. INTRODUCTION
Quantum computer has the potential to outperform the best possible classical computers in many tasks such as factoring large numbers. It relies on the fact that wavefunctions represent amplitudes that grow exponentially in terms of the system size [1]. At the heart is quantum coherence, which is fragile and easily destroyed by noise. In principle, this drawback may be overcome by the techniques of quantum error correction and fault-tolerance, which however require hundreds of thousands of qubits to perform computing tasks of practical relevance [2,3]. For near-term hardware systems, precision needs to improve while systems size grows, to be able to perform reasonable computing tasks before decoherence, whose performance may be measured by the so-called quantum volume [4]. It is claimed that systems with quantum volume as large as 32 have been achieved [5], and systems of quantum volume 64 is on the way [6].
Real world quantum computers battling noise recently achieved the so-called quantum supremacy, at Google, for implementing random quantum circuits in a 53-qubit system and a circuit depth of 25 [7]. It is well-known that noiseless random circuits are hard to simulate on classical computers [8][9][10][11][12][13], and simulations of (noiseless) random quantum circuits on supercomputers have been implemented for more than 40 qubits [14][15][16][17]. It still remains unclear, however, whether there are classical algorithms running on available supercomputers that may be able to simulate the behavior of the Google system, due to physical noise that results in low fidelity compared to noiseless systems. For instance, a method based on second storage has been proposed, which may simulate the system in a few days [18].
It is recently proposed in [19] that a method based on Matrix Product States (MPS) (or Tensor Networks in higher spatial dimensions) can approximate the behavior of real quantum systems. The MPS has been a powerful method that faithfully represents ground states of local Hamiltonians [20][21][22]. The Singular Value Decomposition (SVD) method for truncating the bond dimension for MPS has been shown great success for finding ground states of 1D local Hamiltonians, both gapped and gapless. It is unclear, however, what is the error model the MPS method represents, for simulating circuit output distribution of real quantum computers.
Since MPS cannot represent mixed states of quantum systems, a natural idea is instead to use Matrix Product Operators (MPO) [23][24][25]. The MPO method has been used to simulate quantum circuits of Shor and Grover algorithm with noise [26]. Very recently, the MPO method has also been used to simulate 1D random circuits [27]. For simulating two-qubit gates, the MPO tensors with a 'canonical' update bring a factor of DN 2 in simulating an N-qubit random circuit with depth D in terms of complexity.
In this work, we simulate noisy 1D random quantum circuits with Matrix Product Density Operators (MPDO), based on the MPDO construction proposed in [25]. Recently, it is also shown that MPDO can describe many-body thermal states efficiently [28]. Compared with the MPO method, the inner indices in the MPDO method capture the classical information of the noise simulation, which also reduces the computational and memory complexity under the condition of weak noise. The MPDO model consists of two parts that are conjugated to each other. By so, the simulation of the two-qubit gates can be done in a similar way as the MPS simulation, which is taken care of by the bond indices. In case of weak system noise, a small inner dimension may be sufficient for the simulation, so the resource cost of MPDO could be significantly less than that of the MPO. In case of strong system noise, a relatively small bond dimension may be sufficient to simulate the noisy circuits, indicating a regime that the noise is large enough for an 'easy' classical simulation.
Moreover, we propose a more effective canonical tensor update scheme, performed after each layer of the circuit, which would truncate the inner dimension to some maximum value κ and the bond dimension to some maximum value χ with a canonicalization of the MPDO for improving simulation accuracy. The complexity of this scheme only proportional to DN for an N-qubit circuit with depth D. The cost of our entire simulation scales as ∼ NDκ 3 χ 3 .
We apply our method to simulate the random quantum circuit with different noise models, including the dephasing noise, the depolarizing noise, and the amplitude damping noise. We demonstrate that MPDO approximates the noisy output quantum states well, while the method based on Matrix Product States (MPS) fails to approximate the noisy output quantum states for any of the noise models considered. This indicates that the bond dimension truncation method of the MPS simulation might not represent any local noise model in real physical systems. With a further look into the deviation from the Porter-Thomas distribution for the ideal random circuit case, relatively small bond dimension for the MPDO method already grasp some 'qualitative behavior' of the noisy output distribution.
We organize our paper as follows. In Sec. II, we discuss the error models we use for our circuit simulation. In Sec. III, we discuss the MPDO method for simulating noisy quantum circuits and its complexity. In Sec. IV, we present our results based on the MPDO method, and compare with an exact noise simulation based on density matrices, and the MPS method based on bond dimension truncation, for different error models. In Sec. V, we study the effect of truncation on the bond and inner dimensions. Summary of the results and discussions on future directions will be given in Sec. VI.

II. NOISE MODELS
Physical noise E for quantum systems with a quantum state ρ are generally characterized by the Kraus representation, as given by where E k s are the Kraus operators and ∑ k E † k E k = I [1]. For the quantum circuit of N qubits, with single-qubit and two-qubit quantum gates, normally the fidelity of single-qubit gates are much higher than that of the two-qubit gates. We will then assume that all single-qubit gates are ideal, and model the noise of a two-qubit gate U by where E k s acting on the two qubits that the two-qubit U is acting on. Alternatively, for the noisy channel by E , denote the channel corresponding to the two-qubit gate U by U , we can hence re-write Eq. (2.2) as In this work, we consider the following noise models.
• Dephasing noise The dephasing noise on a single qubit can be modeled by where Z is the Pauli operator.
For a noisy two-qubit gate U, we model the dephasing noise by • Depolarizing noise We model the depolarizing noise by where M = 2 N , and N is the number of qubits in the system. The output state of the circuit under this depolarizing noise, in fact, will also be given in the form of depolarizing noise, i.e.
for some α ∈ [0, 1], where |ψ is the noiseless output state of the quantum circuit.
• Amplitude damping noise The amplitude damping noise on a single qubit can be modeled by , acting on one of the qubit U is acting on.
For a noisy two-qubit gate U, we model the amplitude damping noise by

III. MODELING NOISE SIMULATION BY MATRIX PRODUCT DENSITY OPERATORS (MPDO)
By applying the Matrix Product Operators(MPO) instead of the Matrix Product States(MPS), one can extended the noise model to represent mixed quantum states, which allows us to introduce typical random noise in a more direct and efficient way. As Fig. 1 shown, we use the MPO to represent the density matrix ρ of N qubits where ρ denotes the density matrix of a mixed quantum state. s k and s k are known as the "physical" indices, which expand into the 2 N × 2 N density matrix ρ. L k and R k denote the left and right "bond" indices, which carry some information of entanglement between qubits. However, there are at least two reasons that directly using the tensor M to form an MPO representation of the density matrix may not be the best choice. First, since the density matrix is a Hermitian matrix, we must guarantee that M = M † , which will make at least half of the M parameters invalid, resulting in additional computational overhead. Second, considering the density matrix ρ = ∑ k λ k |φ k φ k |, when k is not too many, the density matrix can be regarded as the sum of the direct products of several state vectors. In this case, it is very inefficient to use tensor networks to model the direct product of state vectors instead of the sum of several states itself.
Therefore, a proper way to ensure the hermiticity of density matrix ρ is the Matrix Product Density Operators(MPDO) [29], which design the M By so, the MPDO will automatically satisfy the hermiticity of a density matrix, although Ts are general 4th order tensors. In this way, the parameters of the model can be fully used on the one hand, and the inner indices between the T [b] and T [r] have emerged on the other hand. Those inner indices can be used to describe the classic entropy of the system. At the same time, when the classic entropy is small, the model's computational cost will become significantly smaller. The new tensors would be obtained by the singular value decomposition (SVD). The SVD would increase the bond dimension between those two tensors, while we could find the global optimal truncation of tensors based on the singular value of SVD. (c) Applying noise models on the corresponding tensor is equivalent to applying each noise operator on tensor separately and then direct sum them on the inner indices. The increased dimensions of inner indices can also be truncated using SVD.
In others words, the MPDO is composed of a general Matrix Product Operators(MPO) and its conjugation. The indices L k and R k of M are correspond to the direct product of indices l k ⊗ l k and r k ⊗ r k respectively. In practice, the dimension of l k and r k bond indices D k are restricted to a certain maximum value χ.
The two conjugate MPOs were connected by the "inner" indices a k , which could carry the classical information of the system. This 'inner' dimensions will increase by adding statistical noise. The dimension of the inner indices d k would be no more than 2D k D k+1 . While in some cases, we still truncate the d k to a smaller number κ < 2D k D k+1 . The memory cost of the MPDO would only proportional to 2Nκχ 2 . Note that if we use the M matrix directly, the memory cost will be 4Nχ 4 , which means that if κ reaches the upper bound, those two models cost almost the same memory. While if the system noise is small, the dimensions of the inner indices will be fixed at a smaller value, and the model will consume significantly fewer resources than using M directly.
To be more intuitive, here we give two particular examples of the MPDO. The first example is the density matrix of a pure quantum state ρ pure = |s s|, where s k ∈ {0, 1}. The corresponding T (k) is {T s k ,a k l k ,r k : T s k ,0 0,0 = 1, others : 0}. For the density matrix of a maximum mixed state , others : 0}. Applying gates and noise on the MPDO is straightforward. Considering the symmetry of the density matrix, one only needs to apply the gate to half of the density matrix, which is the MPO, the rest automatically becomes the conjugation of updated tensors. As Fig. 2 (a) shown, the one-qubit gate wouldn't change any topology or dimensions of the MPDO, we could do it exactly with complexity ∼ κχ 2 .
A two-qubit gate would be represented as a 4th order tensor U. As shown in Fig. 2 (b), we first contract the gate with the corresponding qubits, which form a 6th order tensor W. It can be separated into two new MPO tensors by the Singular Value Decomposition (SVD). In general cases, applying a two-qubit gate would increase the dimension of bond indices between the two qubits to min(2D k d k , 2D k+2 d k+2 ). If we truncate it to χ, the approximation error is where λ i are the singular values in descending order. More details about the strategy of this truncation are placed at the end of this section. The computation cost is ∼ κ 2 χ 3 for contraction, and ∼ κ 3 χ 3 for SVD.
Compared with the MPS method, the MPDO has a clear efficiency advantage in adding noise. Unlike the state representation, the density matrix can directly express noise as the Kraus representation in (2.1), which avoids repeated Monte Carlo sampling of different noise. Let's take the amplitude damping noise as an example. What we need to do is to apply the A 0 and A 1 in (2.8) to ρ respectively and then add them directly. Same to the gate operators, we only need to apply noise to half of ρ to save the computational cost. Note that the summation of noise and the contraction of the conjugate tensors are not interchangeable. To avoid unnecessary cross-terms, as shown in Fig. 2 (c), here we need direct sum the different noise parts on the inner indices a k . which writes, thus, The computational cost of applying noise is ∼ mκχ 2 , m is the number of terms of the noise model. There is an interesting fact behind the structure of the MPDO. Note that the two-qubit gate introduces entanglement entropy into the system, and it only increases the dimension of the bond indices; meanwhile the single-qubit noise introduces classical statistical entropy into the system, and only increases the dimension of the inner indices. Therefore, if there is no truncation, in the final structure of the MPDO, the dimensions of the bond indices and inner indices will be related to the quantum entanglement entropy and classical statistical entropy of the system, respectively. The exact relationship between bond/inner dimensions and entanglement/classical entropy remains an open question.
In many cases, especially those with significant noise, a small dimension of bond indices is enough to ensure high fidelity. So in order to remove unnecessary parameters, we may truncate the bond indices and inner indices, simultaneously. More specifically, we will first do a local optimal approximation on inner indices by SVD, which separate tensor T into where S µ is the singular value of T, U and V are unitary matrix. Then we only keep κ largest S µ and corresponding orthogonal vectors in U and V. The approximate M can be write as We repeat this process on each inner indices. Then, before we truncate bond indices, we first apply QR decomposition from the left of MPDO to the right to form a canonical form of the MPO, which writes After applying this on all qubits, all tensors except the rightmost one are left canonicalized. In other words, they fulfill ∑ l k ,s k ,a k T s k ,a k l k ,r k T s k ,a k l k ,r k = I r k ,r k . (3.8) Where I denotes the identity matrix. This canonicalization will ensure the SVD truncation on rightmost tensor is globally optimal. We then use the SVD to truncate each bond indices from right to left.
∑ l k+1 T s k ,a k l k ,l k+1 T s k+1 ,a k+1 l k+1 ,r k+1 Note that each time the SVD changes the right tensor from left canonicalized to right canonicalized, which results in the following SVD truncations are all globally optimal. Therefore, the most economical way is first to complete a layer of two-qubit gates and noise (see Fig. 3 for an illustration of a layer in the dotted line circuit), then to perform a canonicalization from left to right and the following SVD decomposition from right to left. Compared with canonicalizing on each qubit independently, the order of this scheme reduces the complexity of canonical truncation from a factor of N 2 to a factor of N without loss of the accuracy. The total cost of simulating an N-qubit circuit with depth D is ∼ DNκ 3 χ 3 . Note that while the canonical form of the MPO can not be used to find the global optimal truncation of inner indices. we could still use the full update method similar to that used in the higher dimensional tensor network to find its global optimal truncation, in case there are some people be willing to tolerate excessive calculation costs.

IV. COMPARISON MPDO SIMULATION WITH DIFFERENT MODELS
We consider 1D random circuits as illustrated in Fig. 3. The colored boxes represent different kinds of single-qubit gates, and the lines and the boxes connecting two qubits are either CNOT or Control-Z gates with equal probability. The single-qubit gates are sampled randomly from the set of universal single-qubit quantum gates. It is known that such kinds of (pseudo-)random circuit with big enough depth D could yield an approximately Haar-distributed unitary, and generate entanglement efficiently [30][31][32]. This kind of (pseudo-)random quantum circuits has been discussed extensively for demonstrating quantum supremacy [7,13,33,34]. The gray dotted line outlines the structure of one layer. The circuit is interleaved by layers of single-qubit gates (colored squares) and two-qubit gates (the dots connected to a square). The single-qubit gates are randomly generated from the rotation operator e ia(σ x sin θ cos φ+σ y sin θ sin φ+σ z cos θ) with three uniformly distributed parameters a, θ, φ ∈ U(0, 2π). The two-qubit gates are either control-NOT or control-Z with equal probability. The noise would apply to each two-qubit gate.

We perform numerical experiments based on
• A circuit simulator for noiseless simulation.
• A density matrix simulator for exact noise simulation.
• An MPS simulator based on SVD truncation.
• An MPDO simulator based on the method discussed in Sec. III.
We consider different error models, including dephasing, depolarizing, and amplitude damping, as discussed in Sec. II.

A. Comparison based on fidelity
We start to consider random circuits with 10 qubits and depth 24. For the output density matrices in the following cases, we define the following notations.
• The output density matrix of the circuit simulator is denoted by ρ 0 , which is a pure state representing the exact output the noiseless random circuit.
• The output density matrix of the density matrix simulator is denoted by ρ e , which gives the exact result of simulating the circuit with a given noise model.
• The output density matrix of the MPDO simulator is denoted by ρ 1 , which is subject to different truncation up to some maximum bond dimension χ and maximum inner dimension κ.
• The output density matrix of the MPS simulator is denoted by ρ 2 , which is a pure state subject to different truncation up to some maximum bond dimension χ.
We use the fidelity function between two quantum states as given by For the MPS method, we define a parameter r as which corresponds to a certain bond dimension truncation χ. For each of the error model, we set r = F (ρ 0 , ρ e ), (4.3) and find the corresponding error rate . We summarize the values of χ and s in Table I.  We then use the MPDO simulator to simulate the noisy random quantum circuit for different error models with error rate as given in Table I. We set cut for max bond dimension χ = 32, and max inner dimension κ = 48.
For each r and each error model, we calculate the fidelity of F (ρ e , ρ 2 ), which demonstrates how the MPS method approximates the exact result of the noisy output density matrix, given the same fidelity of ρ e and ρ 2 with the exact noiseless output density matrix. We also calculate the fidelity of F (ρ e , ρ 1 ), which demonstrates how the MPDO method approximates the exact result of the noisy output density matrix. Our results are shown in Fig. 4.
From Fig. 4, it is clearly shown that the MPS simulator fails to simulate any of the noise model considered, even if it gives the fidelity with the exact noiseless state with ρ e . This indicates that the MPS truncation approximate is not simulating physical noise in real systems. The MPDO method approximates ρ e well, which can, in fact, simulate any physical noise as given by the MPDO model construction. We have used the maximum bond dimension (32) in the simulation but truncates the inner dimension to 48.

B. Deviation from the Porter-Thomas distribution
In this section, we consider random circuits with 15 qubits and depth 24. We focus on analyzing how the Porter-Thomas distribution changes due to the effect of noise, and compare the different method of simulation.
For a density matrix ρ, consider a random variable for x i is the i-th bit-string from {0, 1} N , thus |x i is one of the computational basis.
If there is no noise, the output pure state |ψ is resulted from the random circuit U: Then for a pure state, the probability of getting a certain base |x i is For a random circuit with sufficiently large depth, the distribution of {p = p i (x)} is known to follow the Porter-Thomas distribution, with expectation 1/M [33]. For depolarizing noise, the output state will also be given in the form of depolarizing noise, i.e.
where |ψ is the noiseless output state of the quantum circuit. In this case, the distribution of {p = p i } is given by The detailed derivation can be seen in appendix A. For random with 15 qubits and depth 24, we calculate the cumulated p distribution for different noise models with: • The Porter-Thomas distribution (corresponding to a noiseless circuit simulator) • Distribution given by an exact noise simulator.
• Distribution given by an MPS simulator based on SVD truncation.
• Distribution given by an MPDO simulator based on the method discussed in Sec. III.
Our results are summarized in Fig. 5, which clearly shows that the MPDO distribution can approach its corresponding exact noisy distribution with the increase of the bond dimension χ and the inner dimension κ. While in the case of amplitude damping noise, the MPS method gives a relatively good approximation of the accumulated distribution when the error rate is large (this makes sense since in that case, the system is approaching some product pure state), per the discussion in Sec. IV A, it does not approximate the actual output distribution of ρ e . Notice that a relatively small bond and inner dimensions χ = κ = 32 for the MPDO simulation already grasp some 'qualitative behavior' of the cumulated p distribution.

V. TRUNCATION OF BOND AND INNER DIMENSIONS
In this section, we study the effect of truncation in our MPDO simulation for bond and inner dimensions. We consider random circuits with 10 qubits and depth 24.
We first study the effect of truncating the bond dimension. Notice bond dimension, in fact, puts an upper bound for the inner dimension. That is, when the bond dimension is truncated to a maximum value of χ, the inner dimension is upper bounded by 2χ 2 . For each of the error models with different error rates, we choose to truncate the bond dimension to a maximum value χ and also truncate the inner dimension to a maximum value κ = 2χ. and computes F (ρ e , ρ 1 ), to see how well the MPDO method approximates the exact noisy result ρ e . Our results are shown in Fig. 6.
As shown in Fig. 6, when the gate error is larger than some threshold (approximately 0.01 for all the error models), smaller bound dimension suffices to simulate the noisy circuits, indicating it is the regime that the noise is large enough for the quantum circuit to be 'easy' for classical simulation.
We further study the effect of truncating the inner dimensions. For each of the error models with different error rates, we choose to truncate the bond dimension to a maximum value χ = 32 and to truncate the inner dimension to a maximum value of κ. We compute F (ρ e , ρ 1 ), to see how well the MPDO method approximates the exact noisy result ρ e . Our results are shown in Fig. 7.
As shown in Fig. 7, in case of weak system noise, smaller inner dimension suffices to simulate the noisy circuits. In this case, the memory cost of MPDO is significantly less than the MPO method by directly using the M tensor.

VI. DISCUSSION
In this work, we have developed a method to simulate noisy quantum circuits based on MPDO. We show that our method approximates the noisy output states well. While the method based on MPS bond dimension truncation, failed to approximate the noisy output states for any of the noise model considered, indicating that the MPS method might not represent any local noise model in real physical systems.
Our MPDO method exhibits the following advantages.
• It reflects a clear physical picture, with inner indices taking care of noise simulation, and bond indices taking care of two-qubit gate simulation.
• Both bond and inner dimensions can be truncated using the SVD method, adaptive to the need of different situations of the noise simulation.
• In case of strong system noise, small bond dimensions are sufficient to simulate the noisy circuits.
• In case of weak system noise, the memory cost of MPDO is significantly less than the MPO method.
• With an effective tensor update scheme that truncates the inner dimension up to a maximum value κ and bond dimension up to some maximum value χ, performed after each layer of the circuit, the cost of our simulation scales as ∼ NDκ 3 χ 3 , for simulating an N-qubit circuit with depth D.
It remains an interesting open question to understand further the relationship between bond/inner dimensions and entanglement/classical entropy. It is also highly desired to generalize our method to simulate noisy circuits in two spatial dimensions so that we can more directly compare it to experimental data from e.g. Google's experiments [7].