Time-evolving block decimation

The time-evolving block decimation (TEBD) algorithm is a numerical scheme used to simulate one-dimensional quantum many-body systems, characterized by at most nearest-neighbour interactions. It is dubbed "time-evolving block decimation" because it dynamically identifies the relevant low-dimensional Hilbert subspaces of an exponentially larger original Hilbert space. The algorithm, based on the matrix product state formalism, is highly efficient when the amount of entanglement in the system is limited, a requirement fulfilled by a large class of quantum many-body systems in one dimension.

Introduction

Time-evolving block decimation is a numerical algorithm that can efficiently simulate the time evolution of one dimensional quantum systems having limited entanglement entropy. Naively, to time evolve a system characterized by a Hamiltonian \hat{H} one would directly exponentiate the Hamiltonian to get the time evolution operator \hat{U}(t)=e^{-i\hat{H}t/\hbar} and apply this to an initial state: \psi(t)=\hat{U}(t)\psi(t=0) However, as the number of degrees of freedom of the system grows it quickly becomes computationally infeasible to perform the associated matrix exponentiation and matrix-vector multiplication. For example, if \psi represents a system of n qubits then the Hilbert space in which \psi resides has dimension 2^n , meaning matrix operations are effectively intractable for all but the smallest values of n . TEBD presents an efficient scheme for performing time evolution by limiting itself to a much smaller subspace of the configuration space. There are several other noteworthy examples of ways to get around this exponential scaling, including quantum Monte Carlo and the density matrix renormalization group. Guifré Vidal proposed the scheme while at the Institute for Quantum Information, Caltech. He asserts that "any quantum computation with pure states can be efficiently simulated with a classical computer provided the amount of entanglement involved is sufficiently restricted". This happens to be the case for a wide suite of Hamiltonians characterized by local interactions, for example, Hubbard-like Hamiltonians. The method exhibits a low-degree polynomial behavior in the increase of computational time with respect to the amount of entanglement present in the system. The algorithm is based on a scheme that exploits the fact that in these one-dimensional systems the eigenvalues of the reduced density matrix on a bipartite split of the system are exponentially decaying, thus allowing one to work in a re-sized space spanned by the eigenvectors corresponding to the selected eigenvalues. The numerical method is efficient in simulating real-time dynamics or calculations of ground states using imaginary-time evolution or isentropic interpolations between a target Hamiltonian and a Hamiltonian with an already-known ground state. The computational time scales linearly with the system size, hence many-particles systems in 1D can be investigated. A useful feature of the TEBD algorithm is that it can be reliably employed for time evolution simulations of time-dependent Hamiltonians, describing systems that can be realized with cold atoms in optical lattices, or in systems far from equilibrium in quantum transport. From this point of view, TEBD had a certain ascendance over DMRG, a very powerful technique, but until recently not very well suited for simulating time-evolutions. With the Matrix Product States formalism being at the mathematical heart of DMRG, the TEBD scheme was adopted by the DMRG community, thus giving birth to the time dependent DMRG, t-DMRG for short. Other groups have developed similar approaches in which quantum information plays a predominant role: for example, in DMRG implementations for periodic boundary conditions, and for studying mixed-state dynamics in one-dimensional quantum lattice systems. Those last approaches actually provide a formalism that is more general than the original TEBD approach, as it also allows to deal with evolutions with matrix product operators; this enables the simulation of nontrivial non-infinitesimal evolutions as opposed to the TEBD case, and is a crucial ingredient to deal with higher-dimensional analogues of matrix product states.

The decomposition of state

Introducing the decomposition of State

Consider a chain of N qubits, described by the function | \Psi \rangle \in H^{{\otimes} N } . The most natural way of describing | \Psi \rangle would be using the local M^N -dimensional basis | i_1,i_2,..,i_{N-1},i_N \rangle : | \Psi \rangle=\sum\limits_{i=1}^{M}c_{i_1i_2..i_N} | {i_1,i_2,..,i_{N-1},i_N} \rangle where M is the on-site dimension. The trick of TEBD is to re-write the coefficients c_{i_1i_2..i_N} : c_{i_1i_2..i_N}=\sum\limits_{\alpha_1,..,\alpha_{N-1}=0}^{\chi}\Gamma^{[1]i_1}_{\alpha_1}\lambda^{[1]}_{\alpha_1}\Gamma^{[2]i_2}_{\alpha_1\alpha_2}\lambda^{[2]}_{\alpha_2}\Gamma^{[3]i_3}_{\alpha_2\alpha_3}\lambda^{[3]}_{\alpha_3}\cdot..\cdot\Gamma^{[{N-1}]i_{N-1}}_{\alpha_{N-2}\alpha_{N-1}}\lambda^{[N-1]}_{\alpha_{N-1}}\Gamma^{[N]i_N}_{\alpha_{N-1}} This form, known as a Matrix product state, simplifies the calculations greatly. To understand why, one can look at the Schmidt decomposition of a state, which uses singular value decomposition to express a state with limited entanglement more simply.

The Schmidt decomposition

Consider the state of a bipartite system \vert \Psi \rangle \in {H_A \otimes H_B} . Every such state |{\Psi}\rangle can be represented in an appropriately chosen basis as: \left\vert \Psi \right\rangle = \sum\limits_{i=1}^{M_{A|B}} a_i \left\vert{\Phi^A_i \Phi^B_i}\right\rangle where |{\Phi^A_i \Phi^B_i}\rangle=| {\Phi^A_i}\rangle\otimes | {\Phi^B_i}\rangle are formed with vectors |{\Phi^A_i}\rangle that make an orthonormal basis in H_A and, correspondingly, vectors |{\Phi^B_i}\rangle , which form an orthonormal basis in {H_B} , with the coefficients a_i being real and positive, \sum\limits_{i=1}^{M_{A|B}}a^2_i = 1 . This is called the Schmidt decomposition (SD) of a state. In general the summation goes up to M_{A|B}=\min(\dim( ),\dim( )) . The Schmidt rank of a bipartite split is given by the number of non-zero Schmidt coefficients. If the Schmidt rank is one, the split is characterized by a product state. The vectors of the SD are determined up to a phase and the eigenvalues and the Schmidt rank are unique. For example, the two-qubit state: | {\Psi}\rangle=\frac{1}{2\sqrt{2}}\left( | {00}\rangle + {\sqrt{3}} | {01}\rangle + {\sqrt{3}} | {10}\rangle + |{11}\rangle\right) has the following SD: \left|{\Psi}\right\rangle = \frac{\sqrt{3}+1}{2\sqrt{2}} \left|{\phi^{A}_1\phi^{B}_1}\right\rangle + \frac{\sqrt{3}-1}{2\sqrt{2}} \left|{\phi^{A}_2\phi^{B}_2}\right\rangle with |{\phi^{A}_1}\rangle=\frac{1}{\sqrt{2}}(|{0_{A}}\rangle+|{1_{A}}\rangle), \ \ |{\phi^{B}_1}\rangle=\frac{1}{\sqrt{2}}(|{0_{B}}\rangle+|{1_{B}}\rangle), \ \ |{\phi^{A}_2}\rangle=\frac{1}{\sqrt{2}}(|{0_{A}}\rangle-|{1_{A}}\rangle), \ \ |{\phi^{B}_2}\rangle=\frac{1}{\sqrt{2}}(|{1_{B}}\rangle-|{0_{B}}\rangle) On the other hand, the state: |{\Phi}\rangle =\frac{1}{\sqrt{3}}|{00}\rangle + \frac{1}{\sqrt{6}}|{01}\rangle- \frac{i}{\sqrt{3}}|{10}\rangle - \frac{i}{\sqrt{6}}|{11}\rangle is a product state: \left|\Phi\right\rangle = \left( \frac{1}{\sqrt{3}} \left|0_A\right\rangle - \frac{i}{\sqrt{3}} \left|1_A\right\rangle \right) \otimes \left( \left|0_B\right\rangle + \frac{1}{\sqrt{2}} \left|1_B\right\rangle \right)

Building the decomposition of state

At this point we know enough to try to see how we explicitly build the decomposition (let's call it D). Consider the bipartite splitting [1]:[2..N] . The SD has the coefficients \lambda^{[1]}_{{\alpha}_1} and eigenvectors \left|{\Phi^{[1]}_{\alpha_1}}\right\rangle \left|{ \Phi^{[2..N]}_{\alpha_1}}\right\rangle . By expanding the \left|{\Phi^{[1]}_{\alpha_1}}\right\rangle 's in the local basis, one can write: |{\Psi}\rangle=\sum\limits_{i_1,{\alpha_1=1}}^{M,\chi}\Gamma^{[1]i_1}_{\alpha_1}\lambda^{[1]}_{\alpha_1}|{i_1}\rangle|{\Phi^{[2..N]}_{\alpha_1}}\rangle The process can be decomposed in three steps, iterated for each bond (and, correspondingly, SD) in the chain: Step 1: express the |{\Phi^{[2..N]}_{\alpha_1}}\rangle 's in a local basis for qubit 2: |{\Phi^{[2..N]}_{\alpha_1}}\rangle=\sum_{i_2}|{i_2}\rangle|{\tau^{[3..N]}_{\alpha_1i_2}}\rangle The vectors |{\tau^{[3..N]}_{\alpha_1i_2}}\rangle are not necessarily normalized. Step 2: write each vector |{\tau^{[3..N]}_{\alpha_1i_2}}\rangle in terms of the at most (Vidal's emphasis) \chi Schmidt vectors |{\Phi^{[3..N]}_{\alpha_2}}\rangle and, correspondingly, coefficients \lambda^{[2]}_{{\alpha}_2} : |\tau^{[3..N]}_{\alpha_1i_2}\rangle=\sum_{\alpha_2}\Gamma^{[2]i_2}_{\alpha_1\alpha_2}\lambda^{[2]}_{{\alpha}_2}|{\Phi^{[3..N]}_{\alpha_2}}\rangle Step 3: make the substitutions and obtain: |{\Psi}\rangle=\sum_{i_1,i_2,\alpha_1,\alpha_2}\Gamma^{[1]i_1}_{\alpha_1}\lambda^{[1]}_{\alpha_1}\Gamma^{[2]i_2}_{\alpha_1\alpha_2}\lambda^{[2]}_{{\alpha}_2}|{i_1i_2}\rangle|{\Phi^{[3..N]}_{\alpha_2}}\rangle Repeating the steps 1 to 3, one can construct the whole decomposition of state D. The last \Gamma 's are a special case, like the first ones, expressing the right-hand Schmidt vectors at the (N-1)^{th} bond in terms of the local basis at the N^{th} lattice place. As shown in, it is straightforward to obtain the Schmidt decomposition at k^{th} bond, i.e. [1..k]:[k+1..N] , from D. The Schmidt eigenvalues, are given explicitly in D: |{\Psi}\rangle=\sum_{\alpha_k}\lambda^{[k]}_{{\alpha}_k}|{\Phi^{[1..k]}_{\alpha_k}}\rangle|{\Phi^{[k+1..N]}_{\alpha_k}}\rangle The Schmidt eigenvectors are simply: |{\Phi^{[1..k]}_{\alpha_k}}\rangle=\sum_{\alpha_1,\alpha_2..\alpha_{k-1}}\Gamma^{[1]i_1}_{\alpha_1}\lambda^{[1]}_{\alpha_1}\cdot\cdot\Gamma^{[k]i_k}_{\alpha_{k-1}\alpha_k}|{i_1i_2..i_k}\rangle and |{\Phi^{[k+1..N]}_{\alpha_k}}\rangle=\sum_{\alpha_{k+1},\alpha_{k+2}..\alpha_{N}}\Gamma^{[k+1]i_{k+1}}_{\alpha_k\alpha_{k+1}}\lambda^{[k+1]}_{\alpha_{k+1}}\cdot\cdot\lambda^{N-1}_{\alpha_{N-1}}\Gamma^{[N]i_N}_{\alpha_{N-1}}|{i_{k+1}i_{k+2}..i_N}\rangle

Rationale

Now, looking at D, instead of M^N initial terms, there are {\chi}^2{\cdot}M(N-2) + 2{\chi}M + (N-1)\chi . Apparently this is just a fancy way of rewriting the coefficients c_{i_1i_2..i_N} , but in fact there is more to it than that. Assuming that N is even, the Schmidt rank \chi for a bipartite cut in the middle of the chain can have a maximal value of M^{N/2} ; in this case we end up with at least M^{N+1}{\cdot}(N-2) coefficients, considering only the {\chi}^2 ones, slightly more than the initial M^N ! The truth is that the decomposition D is useful when dealing with systems that exhibit a low degree of entanglement, which fortunately is the case with many 1D systems, where the Schmidt coefficients of the ground state decay in an exponential manner with \alpha : \lambda^{[l]}_{{\alpha}_l}{\sim}e^{-K\alpha_l},\ K>0. Therefore, it is possible to take into account only some of the Schmidt coefficients (namely the largest ones), dropping the others and consequently normalizing again the state: |{\Psi}\rangle=\frac{1}{\sqrt{\sum\limits_{{\alpha_l}=1}^{{\chi}_c} ^2}}\cdot\sum\limits_{{{\alpha}_l}=1}^{{\chi}_c}\lambda^{[l]}_{{\alpha}_l}|{\Phi^{[1..l]}_{\alpha_l}}\rangle|{ \Phi^{[l+1..N]}_{\alpha_l}}\rangle, where \chi_c is the number of kept Schmidt coefficients. Let's get away from this abstract picture and refresh ourselves with a concrete example, to emphasize the advantage of making this decomposition. Consider for instance the case of 50 fermions in a ferromagnetic chain, for the sake of simplicity. A dimension of 12, let's say, for the \chi_c would be a reasonable choice, keeping the discarded eigenvalues at 0.0001 % of the total, as shown by numerical studies, meaning roughly 2^{14} coefficients, as compared to the originally 2^{50} ones. Even if the Schmidt eigenvalues don't have this exponential decay, but they show an algebraic decrease, we can still use D to describe our state \psi . The number of coefficients to account for a faithful description of \psi may be sensibly larger, but still within reach of eventual numerical simulations.

The update of the decomposition

One can proceed now to investigate the behaviour of the decomposition D when acted upon with one-qubit gates (OQG) and two-qubit gates (TQG) acting on neighbouring qubits. Instead of updating all the M^N coefficients c_{i_1i_2..i_N} , we will restrict ourselves to a number of operations that increase in \chi as a polynomial of low degree, thus saving computational time.

One-qubit gates acting on qubit k

The OQGs are affecting only the qubit they are acting upon, the update of the state |{\psi}\rangle after a unitary operator at qubit k does not modify the Schmidt eigenvalues or vectors on the left, consequently the \Gamma^{[k-1]} 's, or on the right, hence the \Gamma^{[k+1]} 's. The only \Gamma 's that will be updated are the \Gamma^{[k]} 's (requiring only at most (M^2\cdot\chi^2) operations), as \Gamma^{'[k]i_k}_{\alpha_{k-1}\alpha_k}=\sum_{j}U^{i_k}_{j_k}\Gamma^{[k]j_k}_{\alpha_{k-1}\alpha_k}.

Two-qubit gates acting on qubits k, k+1

The changes required to update the \Gamma 's and the \lambda 's, following a unitary operation V on qubits k, k+1, concern only \Gamma^{[k]} , and \Gamma^{[k+1]} . They consist of a number of ({M\cdot\chi}^3) basic operations. Following Vidal's original approach, |{\psi}\rangle can be regarded as belonging to only four subsystems: {{H}=J H_C{\otimes}H_D{\otimes}K}.\, The subspace J is spanned by the eigenvectors of the reduced density matrix \rho^{J} = Tr_{CDK}|\psi\rangle\langle\psi| : \rho^{[1..{k-1}]}=\sum_{\alpha}{(\lambda^{[k-1]}_{\alpha})}^2|{\Phi^{[1..{k-1}]}_{\alpha}}\rangle\langle{\Phi^{[1..{k-1}]}_{\alpha}}|=\sum_{\alpha}{(\lambda^{[k-1]}_{\alpha})^2}|{\alpha}\rangle\langle{\alpha}|. In a similar way, the subspace K is spanned by the eigenvectors of the reduced density matrix: \rho^{[{k+2}..{N}]}=\sum_{\gamma}{(\lambda^{[k+1]}_{\gamma})^2}|{\Phi^{[{k+2}..N]}_{\gamma}}\rangle\langle{\Phi^{[{k+2}..N]}_{\gamma}}|=\sum_{\gamma}{(\lambda^{[k+1]}_{\gamma})^2}|{\gamma}\rangle\langle{\gamma}|. The subspaces H_C and H_D belong to the qubits k and k + 1. Using this basis and the decomposition D, |{\psi}\rangle can be written as: |{\psi}\rangle=\sum\limits_{\alpha,\beta,\gamma=1}^{\chi}\sum\limits_{i,j=1}^{M}\lambda^{[C-1]}_{\alpha}\Gamma^{[C]i}_{\alpha\beta}\lambda^{[C]}_{\beta}\Gamma^{[D]j}_{\beta\gamma}\lambda^{[D]}_{\gamma}|{{\alpha}ij{\gamma}}\rangle Using the same reasoning as for the OQG, the applying the TQG V to qubits k, k + 1 one needs only to update \Gamma^{[C]} , \lambda and \Gamma^{[D]}. We can write |{\psi'}\rangle=V|{\psi}\rangle as: |{\psi'}\rangle=\sum\limits_{\alpha,\gamma=1}^{\chi}\sum\limits_{i,j=1}^{M}\lambda_{\alpha}\Theta^{ij}_{\alpha\gamma}\lambda_{\gamma}|{{\alpha}ij\gamma}\rangle where \Theta^{ij}_{\alpha\gamma}=\sum\limits_{\beta=1}^{\chi}\sum\limits_{m,n=1}^{M}V^{ij}_{mn}\Gamma^{[C]m}_{\alpha\beta}\lambda_{\beta}\Gamma^{[D]n}_{\beta\gamma}. To find out the new decomposition, the new \lambda 's at the bond k and their corresponding Schmidt eigenvectors must be computed and expressed in terms of the 's of the decomposition D. The reduced density matrix \rho^{'[DK]} is therefore diagonalized: \rho^{'[DK]}=Tr_{JC}|{\psi'}\rangle\langle{\psi'}|=\sum_{j,j',\gamma,\gamma'}\rho^{jj'}_{\gamma\gamma'}|{j\gamma}\rangle\langle{j'\gamma'}|. The square roots of its eigenvalues are the new \lambda 's. Expressing the eigenvectors of the diagonalized matrix in the basis: \ , or, alternatively, simulate an isentropic evolution using a time-dependent Hamiltonian, which interpolates between the Hamiltonian H_1 , which has the product state |{\psi_P}\rangle as its ground state, and the Hamiltonian \tilde{H} ; the evolution must be done slow enough, such that the system is always in the ground state or, at least, very close to it. | 4 = finally, make the time-evolution of the state |{\psi_0}\rangle towards |{\psi_{T}}\rangle using the Hamiltonian H_n : |{\psi_ }\rangle=e^{-iH_nT}|{\psi_0}\rangle }}

Error sources

The errors in the simulation are resulting from the Suzuki–Trotter approximation and the involved truncation of the Hilbert space.

Errors coming from the Suzuki–Trotter expansion

In the case of a Trotter approximation of {\it{p^{th}}} order, the error is of order {\delta}^{p+1} . Taking into account n = \frac{T}{\delta} steps, the error after the time T is: \epsilon=\frac{T}{\delta}\delta^{p+1}=T\delta^p The unapproximated state |{\tilde{\psi}_{Tr}}\rangle is: |{\tilde{\psi}_{Tr}}\rangle = \sqrt{1-{\epsilon}^2}|{\psi_{Tr}}\rangle + {\epsilon}|{\psi^{\bot}_{Tr}}\rangle where |{\psi_{Tr}}\rangle is the state kept after the Trotter expansion and |{\psi^{\bot}_{Tr}}\rangle accounts for the part that is neglected when doing the expansion. The total error scales with time T as: \epsilon(T) = 1 -|\langle{\tilde{\psi_{Tr}}}|{\psi_ }\rangle|^2 = 1 - 1 + \epsilon^2 = \epsilon^2 The Trotter error is independent of the dimension of the chain.

Errors coming from the truncation of the Hilbert space

Considering the errors arising from the truncation of the Hilbert space comprised in the decomposition D, they are twofold. First, as we have seen above, the smallest contributions to the Schmidt spectrum are left away, the state being faithfully represented up to: \epsilon({ }) = 1 - \prod\limits_{n=1}^{N-1}(1-\epsilon_n) where \epsilon_n = \sum\limits_{\alpha=\chi_c}^{\chi}(\lambda^{[n]}_{\alpha})^2 is the sum of all the discarded eigenvalues of the reduced density matrix, at the bond {\it{n}} . The state |{\psi}\rangle is, at a given bond {\it{n}} , described by the Schmidt decomposition: |{\psi}\rangle = \sqrt{1-\epsilon_n}|{\psi_{D}}\rangle + \sqrt{\epsilon_n}|{\psi^{\bot}_{D}}\rangle where |{\psi_{D}}\rangle = \frac{1}{\sqrt{1-\epsilon_n}}\sum\limits_{{{\alpha}_n}=1}^{{\chi}_c}\lambda^{[n]}_{{\alpha}_n}|{\Phi^{[1..n]}_{\alpha_n}}\rangle|{ \Phi^{[n+1..N]}_{\alpha_n}}\rangle is the state kept after the truncation and |{\psi^{\bot}_{D}}\rangle = \frac{1}{\sqrt{\epsilon_n}}\sum\limits_{{{\alpha}_n}={\chi}_c}^ \lambda^{[n]}_{{\alpha}_n}|{\Phi^{[1..n]}_{\alpha_n}}\rangle|{ \Phi^{[n+1..N]}_{\alpha_n}}\rangle is the state formed by the eigenfunctions corresponding to the smallest, irrelevant Schmidt coefficients, which are neglected. Now, \langle\psi^{\bot}_{D}|\psi_{D}\rangle=0 because they are spanned by vectors corresponding to orthogonal spaces. Using the same argument as for the Trotter expansion, the error after the truncation is: \epsilon_n = 1 - |\langle{\psi}|\psi_{D}\rangle|^2 = \sum\limits_{\alpha=\chi_c}^{\chi}(\lambda^{[n]}_{\alpha})^2 After moving to the next bond, the state is, similarly: |{\psi_{D}}\rangle = \sqrt{1-\epsilon_{n+1}}|{{{\psi}'}_{D}}\rangle + \sqrt{\epsilon_{n+1}}|{{\psi'}^{\bot}_{D}}\rangle The error, after the second truncation, is: \epsilon = 1 - |\langle{\psi}|\psi'_{D}\rangle|^2 = 1 - (1-\epsilon_{n+1})|\langle{\psi}|\psi_{D}\rangle|^2 = 1 - (1-\epsilon_{n+1})(1-\epsilon_{n}) and so on, as we move from bond to bond. The second error source enfolded in the decomposition D is more subtle and requires a little bit of calculation. As we calculated before, the normalization constant after making the truncation at bond l ([1..l]:[l+1..N]) is: R = {\sum\limits_{{\alpha_l}=1}^{{\chi}_c} ^2} = {1 - \epsilon_l} Now let us go to the bond {\it{l}}-1 and calculate the norm of the right-hand Schmidt vectors \|{\Phi^{[l-1..N]}_{\alpha_{l-1}}}\| ; taking into account the full Schmidt dimension, the norm is: n_1 = 1 = \sum\limits_{\alpha_l=1}^{\chi_c}(c_{\alpha_{l-1}\alpha_{l}})^2(\lambda^{[l]}_{\alpha_l})^2 + \sum\limits_{\alpha_l=\chi_c}^{\chi}(c_{\alpha_{l-1}\alpha_{l}})^2(\lambda^{[l]}_{\alpha_l})^2 = S_1 + S_2, where (c_{\alpha_{l-1}\alpha_{l}})^2 = \sum\limits_{i_l=1}^{d}(\Gamma^{[l]i_l}_{\alpha_{l-1}\alpha_{l}})^{*}\Gamma^{[l]i_l}_{\alpha_{l-1}\alpha_{l}} . Taking into account the truncated space, the norm is: n_{2}=\sum\limits_{\alpha_l=1}^{\chi_c} (c_{\alpha_ \alpha_{l}})^2\cdot({\lambda'}^{[l]}_{\alpha_l})^2=\sum\limits_{\alpha_l=1}^{\chi_c}(c_{\alpha_ \alpha_{l}})^2\frac{(\lambda^{[l]}_{\alpha_l})^2}{R} = \frac{S_1}{R} Taking the difference, \epsilon = n_2 - n_1 = n_2 - 1 , we get: \epsilon = \frac{S_1}{R} - 1 \leq \frac{1-R}{R} = \frac{\epsilon_l}{1-\epsilon_l} {\to}0\ \ as\ \ Hence, when constructing the reduced density matrix, the trace of the matrix is multiplied by the factor: |\langle{\psi_{D}}|\psi_{D}\rangle|^2 = 1 - \frac{\epsilon_l}{1-\epsilon_l} = \frac{1-2\epsilon_l}{1-\epsilon_l}

The total truncation error

The total truncation error, considering both sources, is upper bounded by: \epsilon({ }) = 1 - \prod\limits_{n=1}^{N-1}(1-\epsilon_n) \prod\limits_{n=1}^{N-1}\frac{1-2\epsilon_n}{1-\epsilon_n} = 1 - \prod\limits_{n=1}^{N-1}(1-2\epsilon_n) When using the Trotter expansion, we do not move from bond to bond, but between bonds of same parity; moreover, for the ST2, we make a sweep of the even ones and two for the odd. But nevertheless, the calculation presented above still holds. The error is evaluated by successively multiplying with the normalization constant, each time we build the reduced density matrix and select its relevant eigenvalues.

"Adaptive" Schmidt dimension

One thing that can save a lot of computational time without loss of accuracy is to use a different Schmidt dimension for each bond instead of a fixed one for all bonds, keeping only the necessary amount of relevant coefficients, as usual. For example, taking the first bond, in the case of qubits, the Schmidt dimension is just two. Hence, at the first bond, instead of futilely diagonalizing, let us say, 10 by 10 or 20 by 20 matrices, we can just restrict ourselves to ordinary 2 by 2 ones, thus making the algorithm generally faster. What we can do instead is set a threshold for the eigenvalues of the SD, keeping only those that are above the threshold. TEBD also offers the possibility of straightforward parallelization due to the factorization of the exponential time-evolution operator using the Suzuki–Trotter expansion. A parallel-TEBD has the same mathematics as its non-parallelized counterpart, the only difference is in the numerical implementation.

References

Text derived from https://en.wikipedia.org/w/index.php?oldid=1335311278 under CC BY-SA 4.0; modifications have been made.