Real-space parallel density matrix renormalization group

We demonstrate how to parallelize the density matrix renormalization group (DMRG) algorithm in real space through a straightforward modification of serial DMRG. This makes it possible to apply at least an order of magnitude more computational power to challenging simulations, greatly accelerating investigations of two-dimensional systems and large parameter spaces. We discuss details of the algorithm and present benchmark results including a study of valence-bond-solid order within the square-lattice $Q_2$ model and Néel order within the triangular lattice Heisenberg model. The parallel DMRG algorithm also motivates an alternative canonical form for matrix product states.

Figure 1

Timing of parallel DMRG ground state calculations for the spin $1/2$ Heisenberg model on the $24\times 8$ square lattice (cylindrical boundary conditions[14]). A speedup of $S$ indicates that the calculation using $n$ nodes was $S$ times faster than the same calculation using only one node. Each calculation consisted of 10 sweeps and reached a relative energy accuracy of ${10}^{-5}$ by keeping $m=2000$ states in the final two sweeps. The speedup is so close to ideal for the two-node case because the reflection symmetry permits optimal load balancing and minimizes waiting time at the communication step.

Figure 2

The parallel DMRG algorithm for the simplest case of dividing a single calculation over two machines. Circles represent lattice sites (or matrix product state site tensors) and the diamonds each represent the same matrix $V$ on the shared bond as in Eq. (4). In step (a), the local information needed to sweep left is copied to the left machine and similarly for the right. (b) Each machine then performs a DMRG sweep in parallel over its half of the system. In step (c) the wave functions are merged together using the prescription preceding Eq. (5). Before repeating the algorithm, (d) the merged state is optimized using Lanczos or Davidson steps on the shared bond.

Figure 3

Prediction fidelity $1-\langle {\Psi }_n^{\prime }|{\Psi }_n\rangle$ at the center (shared) bond of a DMRG calculation parallelized across two nodes. $|{\Psi }_n^{\prime }\rangle$ is the initial merged wave function following sweep $n$ as in Eq. (5). $|{\Psi }_n\rangle$ is the wave function resulting from a fully converging Davidson calculation initialized with $|{\Psi }_n^{\prime }\rangle$. The systems are, from top to bottom, the 200 site $S=1/2$ Heisenberg chain without warmup sweeps (blue dashed line and open circles); $S=1/2$ chain with warmup sweeps (blue solid line and solid circles); $S=1$ Heisenberg chain without (black dashed line and open squares) and with warmup sweeps (black solid line and solid squares). The warmup consisted of five serial DMRG sweeps keeping $m=50$ states. For the parallel sweeps, $m$ was increased after every other sweep to a maximum of 600.

Figure 4

Sweeping pattern for one full sweep of the parallel DMRG algorithm split over four computational nodes. First, (a) the nodes are positioned in a spatially staggered pattern and sweep to the other end of their block. (b) When the nodes reach the end of their block they wait for their neighboring node to arrive, then communicate. (c) Finally the nodes sweep back to their starting positions and (d) communicate with their other neighbor.

Figure 5

Columnar VBS order parameter $\langle {\text{D}}_x\rangle$ on open cylinders of the square-lattice $Q_2$ model [Eq. (6)] with strong-pinning boundary conditions.[14] Each order parameter value was estimated by increasing the states kept $m$ after every other sweep and extrapolating in the truncation error. Error bars are smaller than symbol sizes except for the $L=10$ system.

Figure 6

Staggered magnetization $M$${}_z$ of the antiferromagnetic $S=1/2$ Heisenberg model on open cylinders of the triangular lattice with strong-pinning boundary conditions.[4] The aspect ratio of each cylinder is defined as $\alpha =L_x/L_y$. Each magnetization value was computed by gradually increasing the number of states kept and extrapolating in the truncation error. Error bars are smaller than symbol sizes except for the $L_y=9$ systems.

Figure 7

Matrix product state in (a) the left-canonical gauge where the OC is the first site. Contracting (b) the first two site tensors, computing a singular-value decomposition, and then multiplying the singular-value matrix (shaded diamond) into the second site tensor transforms the MPS into (c) a mixed-canonical gauge with the second site as the OC.[30]

Figure 8

Matrix product state in (a) the left-canonical gauge with the first site an OC. In panel (b) combine the first two site tensors and compute a singular-value decomposition but now with two copies of the singular-value matrix $\Lambda$ (shaded diamonds) and its inverse $V$ (white circle). Multiplying each singular-value matrix $\Lambda$ into the neighboring site tensor transforms the MPS into (c) a new gauge having two OCs.

Figure 9

Inverse canonical matrix product state gauge. Each site tensor $\psi$ is an OC. This gauge is similar to the canonical gauge ($\Gamma$-$\Lambda$ form)[30, 36] where each bond tensor $\Lambda$ is an OC and is a diagonal matrix containing the Schmidt decomposition weights at that bond. Here each bond tensor $V={\Lambda }^{-1}$ is a diagonal matrix containing the inverse Schmidt weights. One can directly map between the two gauges via $V_j={\Lambda }_j^{-1}$ and ${\psi }_j={\Lambda }_{j-1}{\Gamma }_j{\Lambda }_j$.