View source for Quantum Simulation/Density-Matrix Renormalization Group

{{lecture}}

==== Variational method ====

While quantum Monte-Carlo methods provide very useful tools to study many-body systems, we have seen that the sign problem constitutes a significant challenge. In addition, the required quantum-classical mapping makes it difficult to implement direct computer algorithms taking only the Hamiltonian as their input. On the other hand, there is a conceptually simple but quite powerful method to approximate the ground state of arbitrary quantum systems which is known as the variational method. In a nutshell, one writes down an ansatz for the ground state as a trial wave function, which contains one or more variational parameters that are chosen such that the minimize the energy of the system.

As a specific example, let us study the Hamiltonian for the Helium atom,
:<math>H = \sum\limits_{i=1,2}\frac{p_i^2}{2} - \sum\limits_{i=1,2} \frac{Z}{r_i} + \frac{1}{|\mathbf{r}_1-\mathbf{r}_2|},</math>
where we have used [[w:Atomic units|atomic units]], <math>Z=2</math> denotes the charge of the nucleus, and the positions <math>\mathbf{r}_i</math> and momenta <math>\mathbf{p}_i</math> satisfy the usual commutation relations <math>[r_{i,\alpha},p_{j,\beta}] = i\delta_{ij}\delta_{\alpha\beta}</math> for the components <math>\alpha,\beta=x,y,z</math>. Physically speaking, it makes sense that the interaction between the electrons shield the Coulomb potential of the nucleus, so we use as an ansatz a product wave function of the form given in position space by <ref name="Schwabl2010">{{Cite book
| publisher = Springer
| isbn = 9783642091070
| last = Schwabl
| first = Franz
| title = Quantum mechanics
| location = Berlin, Germany
| date = 2010
}}</ref>
:<math>\psi(\mathbf{r}_1,\mathbf{r}_2) = \frac{\tilde{Z}^3}{\pi}\exp(-\tilde{Z}[r_1+r_2]),</math>
where the variational parameter <math>\tilde{Z}</math> accounts for the shielding of the core. Computing the expectation value of the Hamiltonian with respect to this trial wave function yields
:<math>\langle H \rangle = \tilde{Z}^2 - 4\tilde{Z} + \frac{5}{8}\tilde{Z}.</math>
This function has a minimum at the effective charge <math>\tilde{Z}=27/16=1.6875</math> and leads to a ground state energy of <math>E_g = -729/256\approx -2.848</math>, which is within 2% accuracy of the experimentally measured value of <math>E_g = -2.903</math> <ref name="Sansonetti2005">{{Cite journal
| doi = doi:10.1063/1.1800011
| issn = 00472689
| volume = 34
| issue = 4
| pages = 1559-2259
| last = Sansonetti
| first = J. E.
| coauthors = W. C. Martin
| title = Handbook of Basic Atomic Spectroscopic Data
| journal = Journal of Physical and Chemical Reference Data
| accessdate = 2013-04-18
| date = 2005-09-28
| url = http://www.nist.gov/pml/data/handbook/index.cfm
}}</ref>.

==== Reduced density matrices ====

As the name implies, the Density-Matrix Renormalization Group (DMRG) method does not operate on pure quantum states, but on density matrices, which were originally discussed in the context of open quantum systems. Generally, we can think of a system decribed by a density matrix <math>\rho</math> as a statistical mixture of pure quantum states <math>|\psi_i\rangle</math>, i.e.,
:<math> \rho = \sum\limits_i p_i |\psi_i\rangle\langle\psi_i|,</math>
where <math>p_i</math> can be interpreted as the probability to find the system in the pure quantum state <math>|\psi_i\rangle</math>. Density matrices have unit trace and are positive-semidefinite, and their time evolution is governed by the von Neumann equation,
:<math>\frac{\mathrm{d}}{\mathrm{dt}}\rho = -\frac{i}{\hbar} [H, \rho].</math>
Density matrices play an important role when looking at a smaller subsystem embedded in an environment. Consider a bipartite system composed of the subsystems <math>A</math> and <math>B</math>, with pure states being represented by
:<math>|\psi\rangle = \sum\limits_{ij} c_{ij} |i\rangle_A|j\rangle_B \equiv \sum\limits_{ij} c_{ij}|ij\rangle.</math>
Then, we can define the reduced density matrix of <math>A</math> as the partial trace over <math>B</math>, given by
:<math>\rho_A = \mathrm{Tr}_B\{\rho\} = \sum\limits_{i,i'}\sum\limits_j \langle ij|\rho|i'j\rangle |i\rangle\langle i'|.</math>

==== The DMRG algorithm ====

The key element of DMRG is to think of the many-body system of interest being composed of a system of size <math>l</math>, attached to an environment of the size. Suppose we have already found a set of <math>D</math> states <math>\{|\phi_S\rangle\}</math>, which allows us to give a good approximation to both the Hamiltonian and its ground state for this particular system size. Then, we can increase the system size to <math>l+1</math> by the following procedure <ref name="Schollwock2005">{{Cite journal
| doi = 10.1103/RevModPhys.77.259
| volume = 77
| issue = 1
| pages = 259-315
| last = Schollwöck
| first = U.
| title = The density-matrix renormalization group
| journal = Reviews of Modern Physics
| accessdate = 2013-04-24
| date = 2005-04-26
| url = http://link.aps.org/doi/10.1103/RevModPhys.77.259
}}</ref>:
# Form a new system <math>S'</math> of size <math>l+1</math> by combining the <math>D</math> states describing the system with the full Hilbert space describing the additional site. For spin <math>1/2</math>, we have a new set of states according to <math>\{|\phi^S_l\uparrow\rangle,|\phi^S_l\downarrow\rangle\}</math>.
# Build a superblock of size <math>2l+2</math> by combining the enlarged system with the enlarged environment. The Hilbert space of this dimension will be <math>4D^2</math> in the case of spins. If the Hamiltonian is reflection symmetric (i.e, left and right are identical), the basis states for the system and the environment will be the same.
# Find the ground state <math>|\psi\rangle</math> of the Hamiltonian of this superblock by [[Quantum Simulation/Exact diagonalization|exact diagonalization]].
# Compute the reduced density matrix for the enlarged system
#:<math>\rho_S' = \mathrm{Tr}_E'\{|\psi\rangle\langle\psi|\}</math>.
# Diagonalize the reduced density matrix and keep only the eigenvectors corresponding to the <math>D</math> largest eigenvalues. This set of states <math>\{|\psi_{l+1}^S\rangle\}</math> will serve as the input for the next iteration step.
# Express the Hamiltonian for system size <math>l+1</math> in this new basis.
# Continue this iteration until sufficient convergence of observables (e.g, of the ground state energy) is obtained.

This procedure is commonly referred to as "infinite-system DMRG" <ref name="Schollwock2005"/>. In practice, one can often achieve significantly better results by accounting for finite size effects. This can be done in a relatively straightforward way by keeping track of the system and the environment separately. Once the infinite-system algorithm has reached the desired size, one can successively increase the system at the expense of the environment. When the environment reaches a minimum size, the procedure is reversed and the environment is increased at the expense of the system. Multiple sweeps of this kind can be performed until the ground state energy has converged to an even better value. During each step, we can also use the estimate for the ground state obtained during the previous step as the initial state of our exact diagonalization procedure, which will result in a drastically improved performance <ref name="Schollwock2005" />. In general, the error of the DMRG procedure will be given by the truncation error <math>\varepsilon</math>, which is the weight of all the states that have been dropped during the DMRG step. Typical values are <math>\varepsilon \sim 10^{-10}</math> for <math>D\sim 100</math>.

As in the case of exact diagonalization, DMRG is not restricted to ground states, but due to the constraints imposed by the convergence of the exact diagonalization step, it works best for low-energy eigenvalues. It is also possible to compute the time evolution of a many-body system using DMRG, for instance by combining it with the Runge-Kutta scheme for the numerical integration of ordinary differential equations <ref name="Feiguin2005">{{Cite journal
| doi = 10.1103/PhysRevB.72.020404
| volume = 72
| issue = 2
| pages = 020404
| last = Feiguin
| first = Adrian E.
| coauthors = Steven R. White
| title = Time-step targeting methods for real-time dynamics using the density matrix renormalization group
| journal = Physical Review B
| accessdate = 2013-04-24
| date = 2005-07-18
| url = http://link.aps.org/doi/10.1103/PhysRevB.72.020404
}}</ref>. This is done by using the DMRG procedure to converge towards four different target states,
:<math>\begin{align}
|k_1\rangle &= \tau H(t)|\psi(t)\rangle\\
|k_2\rangle &= \tau H(t+\tau/2)\left[|\psi(t)\rangle + \frac{1}{2}|k_1\rangle\right]\\
|k_3\rangle &= \tau H(t+\tau/2)\left[|\psi(t)\rangle + \frac{1}{2}|k_2\rangle\right]\\
|k_3\rangle &= \tau H(t+\tau)\left[|\psi(t)\rangle + |k_3\rangle\right]\end{align}.</math>
The state at the time <math>t+\tau</math> then follows from <ref name="Press1992">{{Cite book
| publisher = Cambridge Univ. Pr.
| isbn = 9780521431088
| last1 = Press
| first1 = William H
| last2 = Teukolsky
| first2 = Saul A
| last3 = Vetterling
| first3 = William T
| last4 = Flannery
| first4 = Brian P 
| title = Numerical recipes in C Book, C version / William H. Press.
| location = Cambridge
| date = 1992
}}</ref>
:<math>|\psi(t+\tau)\rangle = \frac{1}{6}[|k_1\rangle + 2|k_2\rangle + 2|k_3\rangle + |k_4\rangle] + O(\tau^5).</math>

==== Matrix product states ====

To understand the success of DMRG, it is instructive to look at the set of states that can be created using this method. It turns out that these are given by so-called "matrix product states" of the form
:<math>|\psi\rangle = \sum\limits_{\{i_j\}} \mathrm{Tr}\left\{\prod\limits_{k=1}^N A_k^{i_k}\right\}|i_1,i_2,\ldots\rangle,</math>
where <math>A_k^{i_k}</math> is a <math>D\times D</math> matrix associated with every site <math>k</math> and the corresponding basis vector <math>|i_k\rangle</math>. For translational invariant systems, the matrix elements will be identical for all sites, but each matrix will still act on its own (auxiliary) Hilbert space, i.e., we can express each <math>A_j</math> as
:<math>A_j = \sum\limits_{\alpha\beta}(A)_{\alpha\beta}|\alpha\rangle_j\langle\beta|_j,</math>
for some orthonormal basis <math>\{|\alpha\rangle\}</math>. During each DMRG step, we perform basis transforms of the form
:<math>|m_l\rangle = \sum_{m_{l-1},\sigma_l} = (A_l)_{m_l,m_{l-1}}|m_{l-1}\rangle|\sigma_l\rangle,</math>
which tell us how to go from the block of size <math>l-1</math> to size <math>l</math> by including an additional site, with <math>\sigma_l = \uparrow,\downarrow</math> for spin <math>1/2</math> <ref name="Schollwock2005" />. Applying this construction recursively, we see that the basis states obtained within DMRG indeed belong to the class of matrix product states. The variational character of the DMRG algorithm becomes apparent when we consider the state of a superblock,
:<math>|\psi\rangle = \sum\limits_{mn}^D \sum\limits_i \psi_{m\sigma_{l/2}\sigma_{l/2+1}n} (A_{l/2-1}\cdots A_M)_{m,(\sigma_1\cdots \sigma_M)} (A_{l/2+2}\cdots A_{N-M})_{n,(\sigma_{N-M}\cdots\sigma_N)}|i\rangle,</math>
where <math>M</math> forms the largest block that can be solved exactly, i.e, <math>\dim(M) = D</math>. Then, we can interpret the coefficients <math>\psi_{m\sigma_{l/2}\sigma_{l/2+1}n}</math>as the variational parameters according to which the energy is being minimized. Further information about the relationship between DMRG and matrix product states can be found in a more recent review article <ref name="Schollwock2011">{{Cite journal
| doi = 10.1016/j.aop.2010.09.012
| issn = 0003-4916
| volume = 326
| issue = 1
| pages = 96-192
| last = Schollwöck
| first = Ulrich
| title = The density-matrix renormalization group in the age of matrix product states
| journal = Annals of Physics
| accessdate = 2013-04-25
| date = 2011-01
| url = http://www.sciencedirect.com/science/article/pii/S0003491610001752
}}</ref>.

==== Entanglement and the area law ====

Having seen that DMRG implements a variational method based upon matrix product states, it is important to pose the question how generic these states will be in terms of ground states of many-body Hamiltonians. A key concept in answering this question is entanglement, i.e., nonclassical correlations between different parts of the system. Consider the two-spin state
:<math>|\psi\rangle = \frac{1}{\sqrt{2}}(|\uparrow\uparrow\rangle+|\downarrow\downarrow\rangle,</math>
which can be represented by the density matrix
:<math>\rho = \left(\begin{array}{cccc}1/2 & 0 & 0 & 1/2\\0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\1/2 & 0 & 0 & 1/2\end{array}\right).</math>
As can be seen from the diagonal elements, this state possesses classical correlations, but its off-diagonal contributions are what is crucial here. It is instructive to look at the reduced density matrix of a single spin,
:<math>\rho_1 = \mathrm{Tr}_2\{\rho\} = \left(\begin{array}{cc}1/2 & 0\\0 & 1/2\end{array}\right).</math>
This state is proportional to the identity, so it is invariant under all unitary transformations, i.e., <math>U\rho U^\dagger=\rho</math>. Consequently, despite the two-spin state being pure, its local density matrix of a single spin is completely random and does not carry any information at all!

We can quantify this behavior using the von Neumann entropy, defined as
:<math>S = -\mathrm{Tr}\{\rho\log\rho\},</math>
which vanishes for pure states and reaches its maximum <math>\log d</math>, with <math>d</math> being the Hilbert space dimension, for the "maximally mixed state",
:<math>\rho = \left(\begin{array}{ccc}1/d \\ & 1/d \\ & & \ddots\end{array}\right),</math>
which for a single spin is precisely the state we have obtained above. The von Neumann entropy is invariant under unitary transformations and therefore does not change under Hamiltonian dynamics. Additionally, it is subadditive with respect to its subsystems, i.e.,
:<math>S(\rho_{AB}) \leq S(\rho_A) + S(\rho_B),</math>
where <math>\rho_A</math> and <math>\rho_B</math> are the density matrices of the individual subsystems <math>A</math> and <math>B</math> <ref name="Araki1970">{{Cite journal
| issn = 0010-3616, 1432-0916
| volume = 18
| issue = 2
| pages = 160-170
| last = Araki
| first = Huzihiro
| coauthors = Elliott H. Lieb
| title = Entropy inequalities
| journal = Communications in Mathematical Physics (1965-1997)
| accessdate = 2013-05-02
| date = 1970
| url = http://projecteuclid.org/euclid.cmp/1103842506
}}</ref>. Consequently, by looking only at parts of the full system, we can only obtain partial information, as seen in the example above. Remarkably, for pure states of the full system, the entropies <math>S(\rho_A)</math> and <math>S(\rho_B)</math> are identical. This can be seen by looking at the Schmidt decomposition <ref name="Nielsen2000">{{Cite book
| publisher = Cambridge University Press
| isbn = 9781107002173
| last1 = Nielsen
| first1 = Michael A
| last2 = Chuang
| first2 = Isaac L 
| title = Quantum computation and quantum information
| location = Cambridge, UK
| date = 2010
}}</ref> of the wave function of the full system,
:<math>|\psi\rangle = \sum\limits_i c_i|\phi_i\rangle|\chi_i\rangle,</math>
which results in the reduced matrices
:<math>\rho_A = \sum\limits_i |c_i|^2 |\phi_i\rangle\langle\phi_i|,\;\;\;\;\rho_B = \sum\limits_i |c_i|^2 |\chi_i\rangle\langle\chi_i|.</math>
As the coefficients <math>|c_i|^2</math> are identical for both subsystems, the reduced density matrices have the same nonzero eigenvalues and therefore the same entropy. Therefore, we can define an "entanglement entropy" for any pure state <math>|\psi\rangle</math> simply as the entropy of the reduced density matrix <math>S(\rho_A)</math>. In general, the entanglement entropy will depend on the choice of the subsystems, however, we are mostly interested in the partitioning that maximizes the entanglement entropy. For translationally invariant systems, this is achieved simply by cutting the system in half.

Coming back to the question on the usefulness of DMRG, we take a look at the entanglement entropy of matrix product states. As the states are encoded in terms of a <math>D\times D</math> matrix, the entanglement entropy of matrix product states satisfies the relation
:<math>S(\rho_A) = 2\log D,</math>
which does not change with the size of the subsystem. So how does this figure compare to typical ground states of many-body systems? For gapped one-dimensional systems with local interactions, it can be shown that the entanglement entropy of the ground state is constant, which can be understood as a consequence of a finite speed of information transfer <ref name="Eisert2010">{{Cite journal
| doi = 10.1103/RevModPhys.82.277
| volume = 82
| issue = 1
| pages = 277-306
| last = Eisert
| first = J.
| coauthors = M. Cramer, M. B. Plenio
| title = Colloquium: Area laws for the entanglement entropy
| journal = Reviews of Modern Physics
| accessdate = 2013-05-02
| date = 2010-02-04
| url = http://link.aps.org/doi/10.1103/RevModPhys.82.277
}}</ref>. Consequently, for this large class of systems, DMRG can be expected to produce accurate results. Furthermore, for critical (gapless) models or models with strong disorder, one typically finds a logarithmic divergence of the entanglement entropy with system size, which still allows for an efficient simulation with DMRG <ref name="Eisert2010" />. On the other hand, in higher-dimensional systems, the entanglement entropy typically satisfies a relation of the form
:<math>S(\rho_A) \lesssim \mathcal{A}(A)</math>,
where <math>\mathcal{A}(A)</math> is the surface area of the subsystem <ref name="Eisert2010" />. Thus, the constant scaling of the entanglement entropy in one-dimensional systems can be seen as the surface area being constant in this case. The consequences of this "area law" of the entanglement entropy for DMRG are disastrous: we cannot expect DMRG to be a useful method for anything beyond one-dimensional systems. While there have been various approaches to extend DMRG-style methods to higher dimensions, their success has been fairly limited.

==== References ====

{{Reflist}}