The Density Matrix Renormalization Group (DMRG) method

The block decimation idea

Let us try a simple idea, using the recursion scheme described above. At every step in the recursion, we add one spin on the right, and our basis dimension grows by a factor $2$. At some point during this recursion, the matrix will be too large to deal with. So let us fix a maximum number of states that we want to keep, $m$. At certain point during the process, the basis dimension will become larger than $m$. It is here that we start applying the truncation rule: diagonalize the Hamiltonian matrix exactly, and keep only the $m$ states with lowest eigenvalues(see Fig. [fig:nrg]).

As the system grows, the basis of the left block changes as we rotate to the new basis of eigenstates of the Hamiltonian. This is done by using a unitary transformation $U$. This matrix $U$ is nothing else but the matrix with the eigenstates ordered in columns. Therefore, adding a spin to the block now involves two steps: (i) we need to build the ‘tilde’ operators as before, and (ii) rotate the Hamiltonian matrix and the tilde operators to the new basis.

Let us assume that our old block before adding a site has a basis $\{|\alpha_{i-1}\rangle\}$, of dimension $D_{i-1}$, and the site has a basis $\{|s_{i}\rangle\}$ of dimension $d$. The new block basis $\{|\alpha_{i-1},s_{i}\rangle\}$ has dimension $d\times D_{i-1}$, such that we can easily diagonalize it to obtain all the eigenvalues and corresponding eigenvectors $\{|\alpha_{i+1}\rangle\}$. We build the matrix $U$ as the $D_{i-1} \times D_{i}$ unitary matrix with the $D_{i}=m$ eigenvectors with largest eigenvalues in the columns: $$U_{\alpha_{i-1}s_i,\alpha_i} = \langle \alpha_{i-1}s_i|\alpha_i\rangle.$$ Before the rotation, the operators had matrix elements: $$\tilde{O}_{\alpha_{i-1}s_i,\alpha'_{i-1}s'_i} = \langle \alpha_{i-1}s_i |\hat{O}|\alpha'_{i-1}s'_i\rangle.$$ We can now rotate all the tilde operators to the new basis as: $$\begin{aligned} \tilde{O}_{\alpha_i,\alpha'_i} & = & \langle \alpha_i |\hat{O}|\alpha'_i\rangle = \sum_{\alpha_{i-1},s_i} \sum_{\alpha'_{i-1}s'_i} \langle \alpha_i|\alpha_{i-1}s\rangle\langle \alpha_{i-1} s_i|\hat{O}|\alpha'_{i-1}s'_i\rangle\langle \alpha'_{i-1}s'|\alpha_i\rangle \nonumber \\ & = & \sum_{\alpha_{i-1},s_i} \sum_{\alpha'_{i-1}s'_i} (U^\dagger)_{\alpha_i,\alpha_{i-1}s_i} \tilde{O}_{\alpha_i \alpha'_i} U_{\alpha'_{i-1}s'_i,\alpha'_i}\end{aligned}$$ where the new matrices will have dimensions $m \times m$. we can now use these matrices to continue to block-growing process by adding another site. This can be repeated until the energy per site converges, or until we reach a desired system size.

It may seem that this new basis would be a natural choice is we assume that the physics of the problem is described by different manifolds with different energy scales. If we keep the lowest energy states and we get rid of the high energy states we can expect to get the low energy physics right. This in fact is the case in problems such as the Kondo and Anderson impurity problems@hewson. However, in strongly correlated, many-body problems such as the Heisenberg chain, this scheme performs poorly.

The density matrix truncation: the kernel of the DMRG

The problem was solved by Steve White by using what he called the ‘density matrix truncation’. He realized (without knowing at the time) that instead of getting rid of high energy states, one has to redistribute the ‘entanglement’ and minimize the loss of information. However, the way he formulated the problem did not incorporate the idea of entanglement, a concept that entered the picture much later after quantum information ideas were used to understand why and when the DMRG actually works. Before introducing these ideas, we shall describe the original formulation of the density matrix truncation@dmrg.

In order to introduce this new concept, we are going to use a new scheme: We are going to use two blocks instead of one, a left block, and a right block, as shown in Fig.[fig:two_blocks]. We are going to grow both blocks simultaneously using the same procedure outlined previously: at every step we add one site at the right of the left block, and one site to the left of the right block. The ground state can then be written as: $$|\Psi\rangle = \sum_{i,j} \Psi_{ij} |i\rangle |j\rangle,$$ where the sum runs over all the states of the left block $|i\rangle$ and right block $|j\rangle$, with the corresponding coefficients $\Psi_{ij}$.

Now the idea is as follows: once we reach the desired basis dimension $m$, we shall rotate the left block to a new basis $|i\rangle \rightarrow |\alpha\rangle$. We want to pick these states $|\alpha\rangle$ in such a way that when we truncate the basis, the “distance” between the original ground state $|\Psi\rangle$, and the new, truncated, variational approximation $|\tilde{\Psi}\rangle$, is minimized: $$S=\left||\Psi\rangle - |\tilde{\Psi}\rangle \right|^2,$$ where $$|\tilde{\Psi} \rangle = \sum_{\alpha=1}^m\sum_{j} \Psi_{\alpha j} |\alpha\rangle|j\rangle.$$

We are going to anticipate the solution: pick the basis $|\alpha\rangle$ given by the $m$ eigenvectors of the reduced density matrix of the left block with the $m$ largest eigenvalues. In order to justify this result, we first need to introduce some important concepts. The reduced density matrix

Imagine that we have a bipartite system, composed by subsystem $A$ and subsystem $B$, as shown in Fig.[fig:partition]. The Hilbert space of the system $A+B$ will be given by the tensor product of the Hilbert spaces of the two subsystems: $H_{A+B}=H_A \otimes H_B$, and will have dimension $D_{A+B} = D_A\times D_B$. Assume that the state of our system is described by a normalized wave-function $|\Psi\rangle$ that has support on $H_{A+B}$. We define the reduced density matrix of subsystem $A$ as $$\hat{\rho}_A = \mathrm{Tr}_B |\Psi\rangle\langle\Psi|.$$ Its corresponding matrix form is $${\rho_A}_{ii'} = \langle i|\hat{\rho}_A|i'\rangle = \sum_j \langle ij|\Psi\rangle\langle \Psi|i'j\rangle = \sum_{j} \Psi_{ij}\Psi^*_{i'j}$$ This operator describes the density matrix of a mixed state, in which the system $A$ is in contact with a bath or environment $B$. This is the price we have to pay for our ignorance of subsystem $B$.

The reduced density matrix has some nice properties:

• It is Hermitian (or symmetric in case of real matrices). This means that its eigenvalues are real.

• Its eigenvalues are non-negative.

• The trace equals to unity: $\mathrm{Tr \rho_A} = 1$.

• Its eigenvectors $|\alpha\rangle$ and eigenvalues $\omega_\alpha$ form an orthonormal basis.

This means that the reduced density matrix can be re-defined in the new eigenvector basis as: $$\hat{\rho}_A = \sum_\alpha \omega_\alpha |\alpha\rangle\langle\alpha|;$$ with $\omega_\alpha \ge 0$ and $\sum_\alpha {\omega_\alpha} = 1$.

These same considerations are valid for the block $B$.

The Singular Value decomposition (SVD)

Consider an arbitrary matrix $\Psi$ of dimensions $D_A \times D_B$. One can prove that $\Psi$ can de factorized as $$\Psi = UDV^{\dagger},$$ where $U$ is a ($D_A \times D_B$) unitary matrix, $V$ is a $(D_B \times D_B)$ unitary matrix, and $D$ is a $(D_B \times D_B$ diagonal matrix with real non-negative numbers along the diagonal, and zeroes elsewhere. Since $U$ and $V$ are unitary, they satisfy: $$\begin{aligned} UU^\dagger = {1}; \nonumber \\ VV^\dagger = {1}. \nonumber \end{aligned}$$ Their columns are orthonormal vectors, so $U$ and $V$ can be regarded as rotation matrices. The diagonal matrix elements $\lambda_\alpha$ of $D$ are known as the “singular values” of $\Psi$.

The Schmidt decomposition

Let us apply the SVD to our quantum wave-function $|\Psi\rangle$ ([psi]), and for illustration, let us assume that $D_B \le D_A$. The coefficients $\Psi_{ij}$ define a matrix $\Psi$. After a SVD, they can be re-written as: $$\Psi_{ij} = \sum_\alpha^{D_B} U_{i\alpha}\lambda_\alpha(V^{\dagger})_{\alpha j} = \sum_\alpha^{D_B} U_{i\alpha}\lambda_\alpha V^{*}_{\alpha j}.$$ The wave-function can now be expressed as: $$\begin{aligned} |\Psi\rangle & = & \sum_{i}^{D_A} \sum_{j}^{D_B} \sum_\alpha^{D_B} U_{i\alpha}\lambda_\alpha V^*_{\alpha j}|i\rangle|j\rangle \nonumber \\ & = & \sum_\alpha^{D_B} \left(\sum_i^{D_A} U_{i\alpha}|i\rangle\right)\lambda_\alpha\left(\sum_j^{D_B} V^*_{\alpha j}|j\rangle\right) \nonumber \\ & = & \sum_\alpha^{D_B} \lambda_\alpha |\alpha\rangle_A|\alpha\rangle_B, \nonumber \end{aligned}$$ where we have defined the states $|\alpha\rangle_A = \sum_i U_{i\alpha}|i\rangle$ and $|\alpha\rangle_B = \sum_j V^*_{\alpha j}|j\rangle$. Due to the properties of $U$ and $V$, these states define a new orthogonal basis. This final expression is known as the “Schmidt decomposition” of the state $\Psi$, and the bases ${|\alpha\rangle}$ as the “Schmidt bases”..

In general, we have that the state $\Psi$ can be written in the new basis as: $$|\Psi\rangle = \sum_{\alpha}^r \lambda_\alpha |\alpha\rangle_A |\alpha\rangle_B; \,\,\, r=\mathrm{min}(D_A,D_B).$$

In the Schmidt basis, the reduced density matrices for the subsystems $A$ and $B$ are $$\rho_A = \mathrm{Tr}|\Psi\rangle\langle\Psi| = \sum_\alpha \lambda_\alpha^2|\alpha\rangle_{A}\ _A\langle \alpha|,$$ and $$\rho_B = \sum_\alpha \lambda_\alpha^2|\alpha\rangle_{B}\ _B\langle \alpha|$$

At this point, we realize some interesting observations:

• The eigenvalues of the reduced density matrices are $\omega_\alpha = \lambda_\alpha^2$, the square of the singular values.

• The two reduced density matrices share the spectrum.

• The Schmidt bases are the eigenvectors of the reduced density matrices.

Optimizing the truncated wave-function

We here go back to the original problem of optimizing the wave-function in a reduced basis. In order to solve it, we are going to reformulate the question as: Given a matrix $\Psi$, what is the optimal matrix $\tilde{\Psi}$ with fixed rank $m$ that minimizes the Frobenius distance between the two matrices? It turns out, this is a well known problem called the “low ranking approximation”.

If we order the eigenvalues of the reduced density matrix in decreasing order $\omega_1,\omega_2,...\omega_m,...\omega_r$, it is straightforward to see that the Frobenius distance between the two matrices is given by $$S=\left|\Psi - \tilde{\Psi}\right|^2 = \sum_{m+1}^r \omega_i$$ This proves that the optimal basis is given by the eigenvectors of the reduced density matrix with the $m$ largest eigenvalues.

Infinite-size DMRG

The above considerations allow us now to introduce the DMRG algorithm in a very natural way. We are going to present it in the traditional formulation, starting with the infinite-size algorithm, followed by the finite-size scheme.

The main idea behind the infinite-size algorithm consists in growing the left and right blocks by adding one site at a time. As we add sites, the basis of the blocks will grow, until we reach the desired maximum number of states $m$. At this point we need to start applying the density matrix truncation on both blocks. This process is repeated until we reach a desired system-size, or the error in the energy is below a pre-defined tolerance.

[fig:infinite]

The algorithm illustrated in Fig.[fig:infinite] could be outlined as below:

• Build all the operator matrices for a single-site Hamiltonian, and the operators involved in the interactions between the site an the rest of the system.

• Start growing the blocks by adding single-sites, as outlined in the exact diagonalization section. We assume that the Hilbert space for the single site has dimension $d$. When the size of the bocks become larger than $d \times m$, we start applying the density matrix truncation as follows:

  1. Using a suitable library routine (Lanczos,Davidson), diagonalize the full Hamiltonian (sometimes called super-Hamiltonian) of the two blocks combined (sometimes refereed to as superblock), to obtain the ground state $$|\Psi\rangle = \sum_{ij}\Psi_{ij}|i\rangle|j\rangle$$.

  2. Calculate the reduced density matrix of the left block, and right blocks. When the system is symmetric under reflections, we only need one of them.

  3. For each of the blocks, diagonalize the density matrix to obtain the full spectrum and eigenvectors.

  4. Truncate the basis by keeping only the $m$ eigenvectors with the largest eigenvalues.

  5. Rotate the Hamiltonian and the operators involved in the interactions between blocks to the new basis.

  6. Add a new site to the left and right blocks, to build new blocks of dimension $d \times m$, and reiterate the diagonalization and truncation steps. Stop when we reach the desired system-size, or the error in the energy is below a pre-defined tolerance.

In the early days of DMRG it was assumed that this scheme would lead to a good approximation of the system properties in the thermodynamic limit. Today we know that he best way to reach the thermodynamic limit is by using the finite-size algorithm on systems of fixed length, and doing a careful finite-size analysis of the results.

Let us now explain some of these steps in more detail.

Adding a single site to the block

[fig:add_site]

Same as we did in the exact diagonalization section, we can add sites to the blocks by performing tensor products of the “tilde” operators on the block, and single-site operators.

Assume that we are in the $i^\mathrm{th}$ iteration of the algorithm, with our left and right blocks having length $i$. Let us label our $D_L$ basis states for the left block $\{|\alpha_{i}\rangle \}$, and our $d$ basis states for the single site that comes to the right $\{|s_{i+1}\rangle \}$ (See Fig.[fig:add_site]). When we add the site to the block, we obtain a new basis for the new combined block as $|\alpha_{i+1}\rangle = |\alpha_i\rangle \otimes |s_{i+1}\rangle$.

Let us assume for illustration purposes that we are dealing once more with the Heisenberg chain. All these ideas can be easily generalized to arbitrary models. Same as we did in the exact diagonalization section, we obtain the new Hamiltonian matrix for the combined block as: $$H_{L,i+1} = H_{L,i} \otimes {1}_2 + \tilde{S}^z_{L,i} \otimes S^z + \frac{1}{2} \left( \tilde{S}^+_{L,i} \otimes S^- + \tilde{S}^-_{L,i} \otimes S^+ \right).$$ In this expression, the “tilde” operators are in the $|\alpha_i\rangle$ basis, while the others are defined in the single-site basis.

A similar expression applies to the right block, which is obtained from the single site at position $i+2$, with basis $\{|s_{i+2}\rangle\}$ and dimension $d$, and the right block with basis $\{|\beta_{i+3}\rangle\}$ and dimension $D_R$: $$H_{R,i+2} = {1}_2 \otimes H_{R,i+3} + S^z \otimes \tilde{S}^z_{R,i+3} + \frac{1}{2} \left( {S^+} \otimes \tilde{S}^-_{R,i+3} + {S^-} \otimes \tilde{S}^+_{R,i+3} \right).$$

Building the super-Hamiltonian

We now need to combine the left and right blocks to form the super-Hamiltonian: $$\hat{H}=\hat{H}_{L,i+1}+\hat{H}_{R,i+2} + S^z_{i+1} S^z_{i+2} + \frac{1}{2}\left( S^+_{i+1}S^-_{i+2} + S^-_{i+1}S^+_{i+2} \right)$$ where $\hat{H}_{L(R)}$ where obtained above, and only involve terms in the left (right) block. The single sites at positions $i+1$ and $i+2$ were absorbed by the left and right blocks, respectively, so in order to build the interactions, we have to rotate the corresponding operators to the new basis of the blocks. This is again done in the same spirit of the “tilde” transformation: $$\begin{aligned} H & = & H_{L,i+1} \otimes {1}_{D_R\times2} + {1}_{D_L \times 2} \otimes H_{R,i+2} \nonumber \\ & + & {1}_{D_L} \otimes S^z \otimes S^z \otimes {1}_{D_R} \nonumber \\ & + & \frac{1}{2} {1}_{D_L} \otimes S^+ \otimes S^- \otimes {1}_{D_R} \nonumber \\ & + & \frac{1}{2} {1}_{D_L} \otimes S^- \otimes S^+ \otimes {1}_{D_R} \nonumber\end{aligned}$$ or: $$\begin{aligned} H & = & H_{L,i+1} \otimes {1}_{D_R\times2} + {1}_{D_L \times 2} \otimes H_{R,i+2} \nonumber \\ & + & \tilde{S}^z_{L,i+1} \otimes \tilde{S}^z_{R,i+2} \nonumber \\ & + & \frac{1}{2} \tilde{S}^+_{L,i+1} \otimes \tilde{S}^-_{R,i+2} \nonumber \\ & + & \frac{1}{2} \tilde{S}^-_{L,i+1} \otimes \tilde{S}^+_{R,i+2} \nonumber \end{aligned}$$

Density matrix truncation and the rotation to the new basis

The truncation process is similar to the one use in numerical renormalization group, but instead of using the matrix of eigenvectors of the Hamiltonian, we use the eigenvectors $\{|\alpha\rangle\},\{|\beta\rangle\}$ of the left and right reduced density matrix. Therefore, the new basis states for the left and right block are related to the states in the previous step as: $$\begin{aligned} |\alpha_{i+1}\rangle & = & \sum_{s_{i+1},\alpha_i}\langle \alpha_i s_{i+1}|\alpha_{i+1} \rangle |\alpha_i s_{i+1} \rangle = \sum_{s_{i+1},\alpha_i} (U_L)_{\alpha_i s_{i+1},\alpha_{i+1}} | \alpha_i s_{i+1} \rangle \nonumber \\ |\beta_{i+2}\rangle & = & \sum_{s_{i+2},\beta_{i+3}}\langle s_{i+2} \beta_{i+3}|\beta_{i+2} \rangle |s_{i+2} \beta_{i+3} \rangle = \sum_{s_{i+2},\beta_{i+3}} (U_R)_{s_{i+2} \beta_{i+3},\beta_{i+2}} | s_{i+2} \beta_{i+3} \rangle\end{aligned}$$ where $$(U_L)_{\alpha_i s_{i+1},\alpha_{i+1}} = \langle \alpha_i s_{i+1}|\alpha_{i+1} \rangle$$ and $$(U_R)_{s_{i+2} \beta_{i+3},\beta_{i+2}} = \langle s_{i+2} \beta_{i+3}|\beta_{i+2} \rangle.$$ If we keep only $m$ states, the matrices $U_{L(R)}$ will have dimensions $D_{L(R)} d\times m$. If the basis had already been truncated in the previous step, then $D_L(R)=m$.

We can now use these transformations to obtain the matrix elements for all the operators in the new truncated basis. For instance, an operator acting on a site inside the left block will be transformed as: $$\begin{aligned} (\tilde{O}_{L,i+1})_{\alpha_{i+1},\alpha'_{i+1}} & = & \langle \alpha_{i+1} |\hat{O}|\alpha'_{i+1}\rangle \nonumber \\ & = & \sum_{\alpha_i,s_{i+1}} \sum_{\alpha'_i,s'_{i+1}} \langle \alpha_{i+1}|\alpha_i s_{i+1}\rangle\langle \alpha_i s_{i+1}|\hat{O}|\alpha'_i s'_{i+1}\rangle\langle \alpha'_i s'_{i+1}|\alpha'_{i+1}\rangle \nonumber \\ & = & \sum_{\alpha_i s_{i+1}}\sum_{\alpha'_i s'_{i+1}} (U_L^\dagger)_{\alpha_{i+1},\alpha_is_{i+1}} (\tilde{O_{L,i}})_{\alpha_i s_{i+1},\alpha'_i s'_{i+1}} (U_L)_{\alpha'_is'_{i+1},\alpha'_{i+1}},\end{aligned}$$ and a similar expression can be obtained for operators in the right block.

The Finite-size DMRG

[fig:finite]

As we mentioned before, the proper way to reach the thermodynamic limit with DMRG is by studying finite systems and performing a finite-size analysis. In order to study finite system, a generalization of the above ideas needs to be applied. The finite-size DMRG (illustrated in Fig.[fig:finite]) can be summarized as follows:

• Run the infinite-size algorithm until the desired system size is reached. During this process, store all the left and right blocks, with their corresponding operators and basis transformations. This step is typically referred to as the ``warmup''.

• Once the desired system size is reached we start performing ``DMRG sweeps'', from right-to-left, and left-to-right to optimize the bases and improve accuracy. A left-to-right sweep is described as follows:

  1. Add a site to the left block using the same idea of the infinite-size DMRG. Since the total size of the system needs to be kept fixed, we need to ``shrink'' the right block. This is done by using the smaller right block from the infinite-size step, or from the previous right-to-left sweep.

  2. Using a suitable library routine (Lanczos,Davidson), diagonalize the super Hamiltonian (sometimes called super-Hamiltonian) of the two blocks combined, same as for the infinite-size DMRG.

  3. Calculate the reduced density matrix of the left block.

  4. Diagonalize the density matrix to obtain the full spectrum and eigenvectors.

  5. Truncate the basis by keeping only the $m$ eigenvectors with the largest eigenvalues.

  6. Rotate the Hamiltonian and the operators of the left block involved in the interactions between blocks to the new basis.

  7. Iterate until reaching the far right end of the system, with a right block containing a single site. This completes the left-to-right sweep.

• Perform a right-to-left sweep, by growing the right block one site at a time, and using the left block from the previous left-to-right sweep.

• Re-iterate the sweeping. Stop when the change in the energy is below a pre-defined tolerance. One typically stops at a point when both blocks have the same size, the ``symmetric configuration''.

This sweeping process works in a similar fashion as a self-consistent loop, where we iteratively improve the solution. In fact, the DMRG can be formulated as a variational method, in which the variational parameters are continuously improved to minimize an energy functional. Intuitively a way to see it is by imagining a “demon” probing the environment around the block for the optimal states to improve the basis to represent the ground-state. These states are “absorbed” inside the block by the density matrix and its eigenvectors.

As described above, the shrinking block is replaced by the block from the previous sweep in the opposite direction. This means that all the information about the block and its operators needs to be stored, either in memory, or dumped on disk or stored in memory.