Infinitesize density matrix renormalization group with parallel Hida’s algorithm
Abstract
We report a parallel algorithm for the infinitesize density matrix renormalization group (iDMRG) method applicable for onedimensional (1D) quantum systems with longperiod structures, which combines the Hida’s iDMRG applied to 1D random spin systems with the McCulloch’s wave function prediction. This algorithm can be regarded as an infinitesize version of the realspace parallel DMRG proposed by Stoudenmire and White. We describe details of our combined method and show benchmark calculations for physical quantities of a ground state of the spin1/2 Heisenberg model on the YC8 kagome cylinder in the thermodynamic limit. We find that our iDMRG can be parallelized efficiently for both sharedmemory and distributed memory systems and properly evaluate various observables including the total energy per site, bond strength on nearest neighbor spins, spinspin correlation function along the infinite YC8 cylinder and its correlation length.
pacs:
05.10.Cc, 02.70.c, 05.30.d, 71.27.+aI Introduction
The density matrix renormalization group White () (DMRG) is a powerful tool to analyze ground states of onedimensional (1D) quantum systems. After the success of DMRG, a lot of algorithms for improving DMRG have been proposed Schollwöck (2005, 2011) to attack the analysis of nearly twodimensional (2D) quantum systems, where the essential parts of the algorithms are categorized into two kinds, the finitesize algorithm and the infinitesize algorithm.
The finitesize DMRG (fDMRG) is a variational method and optimizes finitelycorrelated states Fannes et al. (1992) or matrixproduct states Östlund and Rommer (1995); Rommer and Östlund (1997) (MPS) on a finitesize lattice. It is known that numerical calculations in fDMRG can be accelerated by use of the structure of MPS White (1996). Many researchers sophisticated this acceleration method and proposed the singlesite DMRG White (2005) to avoid converging into a metastable state even if the 1D system has longperiod structures. This improvement gives us a way to analyze ground states of nearly twodimensional (2D) quantum systems by use of DMRG Stoudenmire and White (2012). Even now the acceleration method is improved, such as the strictly singlesite DMRG Hubig et al. (2015).
On the other hand, the infinitesize DMRG (iDMRG) grows the lattice in each iteration by a unit cell, which contains one or more sites, and produces a uniform MPS with a nontrivial unit cell at the fixed point Östlund and Rommer (1995); Rommer and Östlund (1997). We can also accelerate the numerical calculation in iDMRG by using the structure of MPS, such as the product wave function renormalization group Nishino and Okunishi (1995); Ueda et al. (2006, 2008) and a wave function prediction proposed by McCulloch McCulloch (); Ueda et al. (2010).
Unfortunately, applying iDMRG to 1D quantum systems with a longperiod unit cell, such as multileg ladder/cylinder systems, is practically not so easy. Of course, in principle, we can apply the iDMRG to such systems, but the number of orthonormal bases for representing local states on a unit cell increases exponentially. This is too inconvenient for iDMRG in terms of numerical costs, because the calculation cost of iDMRG is proportional to the cube of the number of bases White (). To avoid this exponential increase, a combination of iDMRG and fDMRG is often used. In such combined method, we have to pay attention to an additional numerical convergence with respect to fDMRG in each iteration process of iDMRG. It is desirable that both of the exponential increase of numerical costs and the additional convergence process are avoided simultaneously in an iDMRG algorithm for multileg ladder/cylinder systems.
In this article, we report an extension of iDMRG with a wave function prediction which can be practically applied to quantum systems with longperiod structures. The numerical cost of this iDMRG is proportional to the number of sites contained in a unit cell, and the algorithm is efficiently parallelized for both sharedmemory and distributed memory systems. In addition, even if the system is a nearly 2D system, this iDMRG can provide a proper ground state without using fDMRG. Therefore, any additional convergence procedures are not included in our parallel iDMRG. We demonstrate above numerical performance when the method applied to the spin1/2 Heisenberg model on the YC8 kagome cylinder Yan et al. (2011), which can be mapped to 1D quantum systems with a site period. The essential algorithm of our parallel iDMRG is Hida’s iDMRG for 1D random systems Hida (1996), which can be regarded as a prototype algorithm of the realspace parallel DMRG Stoudenmire and White (2013).
The organization of this article is as follows. In the next section, we recall the Hida’s iDMRG algorithm in terms of the matrix product formalism. In Section III, we introduce the parallel iDMRG algorithm for quantum 1D systems with longperiod structures. We demonstrate the numerical efficiency of our algorithm in Section IV, where we discuss fidelity error McCulloch () of the wave function predictions, a subtraction procedure, and an extrapolation procedure for estimating the ground state energy of the YC8 kagome cylinder in the thermodynamic limit. Also we discuss bond strength on nearest neighbor spins, a spinspin correlation function and its correlation length. We show all of our results are consistent with reported observables given by largescale fDMRG calculations Yan et al. (2011); Depenbrock et al. (2012); Kolley et al. (2015). We summarize conclusion and remarks of this article in the last section, where we mention the relation between Hida’s iDMRG and the realspace parallel DMRG Stoudenmire and White (2013).
Ii Hida’s iDMRG
In this section, we review the algorithm of Hida’s iDMRG Hida (1996) in terms of a matrix product formalism. Target systems of this algorithm are site systems represented by positiondependent Hamiltonians, where is an even value. Using a formalism of matrix product operator (MPO) McCulloch (), we can express a position dependent Hamiltonian as
(1) 
where is a lowertriangular matrix defined for the outer product of local states and is depends on the position . Left and right boundary vectors, and , are identical with the last row and the first column of the matrix , respectively. Hereafter, unless otherwise noted, we abbreviate the subscripts on for simplifying notations. Figure 1 shows the graphical representation of , , and , where vertical and horizontal lines coming from rounded squares represent physical and auxiliary variables, respectively. An integer means the number of rows (columns) of (). The number of lines without indexes of states represents the rank of tensor, which means that is a matrix, and are row and column vectors, respectively.
The Hida’s iDMRG proceeds as follows:

Give in Eq. (1) and prepare pairs of blocks, and , where . Then, set the iteration number .

Expand each pair of blocks as follows
(2) (3) where the range of integers and are and . Parameters and are the degree of freedom of local state and the maximum number of renormalized sates for each block, respectively. Graphical representations of Eq. (2) and (3) are shown in Fig. 2(a), where we take configuration sum for those links connecting the neighboring diagrams.

Solve an eigenvalue problem of each super block Hamiltonian White () as by the Lanczos or JacobiDavidson iteration method, where the groundstate energy and a corresponding eigenvector are represented by and , respectively. An initial vector is required for starting the iteration method, which is often given randomly. Graphical representations of in Fig. 2(b) show that , and can be reshaped as rankeight, zero and four tensors, respectively.

Apply the singular value decomposition (SVD) to as
(4) where and are unitary matrices. The diagonal matrix contains the singular values, which are normalized as due to . A graphical representation of this decomposition is shown in Fig. 2(c).

If , finish the calculations because the vector represents a variational/exact wave function of the original Hamiltonian in Eq. (1).

Apply blockspin transformations to each expanded block as follows
(5) (6) where the range of is . The graphical representation of these transformations are depicted in Fig. 2(d). Truncations of the degree of freedom of blocks can be introduced in this step. The calculation of DMRG becomes monotonically precise with increasing .

Set and go to step 2.
Through the processes stated above, we obtain a variational/exact ground state, , of the original Hamiltonian using the matrix product formalism as follows
(7)  
where and are matrices defined for the local state . In the same equation, and are a column and a row vectors defined for and , respectively. A graphical representation of Eq. (7) is shown in Fig. 2(e), and we describe the state given by Eq. (7) as MPS.
To see the whole picture of the Hida’s algorithm, we show a schematic procedure of this algorithm for in Fig. 3.
Keys of this algorithm are preparing and growing up to pairs of blocks, and , for providing proper environments in each positiondependent DMRG calculations. Because of this careful treatment, this infinitesize algorithm can be effectively applied to analyze a ground state of random quantum 1D systems without using fDMRG Hida (1996, 1997a, 1997b, 1999).
In addition, the expanding and blockspin transformation process of each pair of blocks can be parallelized easily. When we have a computing machine including nodes, the th node may store and update only a pair of blocks and . A onetoone communication between “Node ” and “Node ” is required for in Eq. (3) and its cost per node does not change with respect to . This property is suitable for the message passing interface parallel (MPI) programming, and the parallel number of nodes becomes in the Hida’s iDMRG.
Iii Parallel DMRG algorithm for systems with site period structure
Here, we discuss a combination of the Hida’s iDMRG Hida (1996) with McCulloch’s wave function prediction McCulloch (); Ueda et al. (2010). A target Hamiltonian containing sites can be represented by an MPO as
(8) 
The parallel iDMRG calculation proceeds as follows:

Initialization steps are identical with the step 1 of the Hida’s algorithm in Section II.

Expand each pair of blocks as
(9) (10) where and . The range of is always

An initial vector for iteration methods is a random vector if and . When , is given by the wave function prediction methods McCulloch (); Stoudenmire and White (2013) as follows
(11) A graphical representation of is shown in Fig. 4.

Steps for solving eigenvalue problems and applying SVDs are identical with the step 3 and 4 of the Hida’s algorithm in Section II.

If , estimate the ground state energy per site as to subtract boundary effects, where set . This subtraction method for estimating the energy per site is identical with that in Ref. [Stoudenmire and White, 2012]. Then, if converges with respect to , finish the iDMRG calculation.

Steps for blockspin transformations and setting are identical with the step 6 and 7 of the Hida’s algorithm in Section II.
As shown in Fig. 5, this parallel iDMRG for systems with site period can be implemented by introducing slight modifications to the Hida’s iDMRG for site stems.
After the calculations, we obtain an MPS for the ground state of the original Hamiltonian using the wave function prediction method iteratively, namely
(12)  
where and . As an example, in Fig. 6, we show a graphical representation of a unit of the uniform part of matrix product structures in Eq. (12) for .
The wave function prediction methods in the step 3 is used to accelerate the iteration methods for eigenvalue problems. A degree of acceleration of solving the problems by the wave function prediction can be discussed using a fidelity error
(13) 
where we assume that the predicted initial vector is unnormalized. When this fidelity error becomes zero, the prediction method produces an eigenvector of the super block Hamiltonian . However, since the Schmidt rank of initial state in Eq.(11) is up to , the fidelity error must be not less than the best fidelity error given by
(14) 
where is an approximated eigenvector defined by with a number of renormalized states . This behavior can be confirmed in benchmark calculations for the kagome cylinder in the next section [see the caption of Fig. 9]. Due to this prediction methods, the numerical error of the eigenvalue with respect to the Lanczos or JacobiDavidson iterations becomes less than up to 30 iterations in the benchmark calculations.
Iv Benchmark calculations
To confirm numerical performance of our parallel iDMRG, we estimate the basic physical quantities, namely, bond strength on nearest neighbor spins, spinspin correlation functions, and correlation lengths, of the spin1/2 Heisenberg model on the YC8 kagome cylinder with the infinite length. The Hamiltonian is given as
(15) 
where the sum runs over nearest neighbor sites. The shape of YC8 cylinder is shown in the inset of Fig. 7. The Hamiltonian of the cylinder can be represented by an MPO with as in Eq. (8). The path of matrix product of our MPO and MPS are denoted by bold red lines on the YC8. The blue broken line separates the MPO into each period of the MPO. The ground state energy of this system is nicely estimated by the fDMRG with keeping states up to .Yan et al. (2011); Depenbrock et al. (2012) Also the spinspin correlation lengths of the YC8 with finite lengths are estimate properly by the fDMRG with keeping states up to .Kolley et al. (2015) We will show that the parallel iDMRG can estimate consistent physical quantities with only up to .
iv.1 Parallel performance
First, we check a parallel performance of the combined method of the Hida’s iDMRG with the McCulloch’s wave function prediction for the ground state of the spin1/2 Heisenberg model on the YC8 kagome cylinder as shown in Fig. 7. Each pair of numbers located nearby each plots, respectively, means the number of nodes and threads per node in the hybrid MPI/OpenMP parallel calculations. Calculation time is defined by the mean time of our parallel iDMRG keeping states per iteration, where the iteration number is up to . In this paper, we do not introduce block diagonalizations with respect to typical quantum numbers, the total spin and its component. Of course, our parallel iDMRG is compatible with the use of abelian and nonabelian symmetries McCulloch and Gulácsi (2002). The calculation time can be fitted by a linear function , where , and are the reciprocal of the parallel cores, the parallel and serial processed parts of our calculation, respectively. A parallel efficiency is defined by , and we obtain in the calculations. This means the parallelization works well up to several hundred cores in this system.
iv.2 Ground state energy per site
Taking to the double limit, namely, the number of iteration of calculations and the number of kept states of iDMRG, we can address the true physical quantity of the cylinder in the thermodynamic limit. In this subsection, we show how to take the double limit properly when we estimate the ground state energy per site of the cylinder.
To estimate the ground state energy, first, we focus on the convergence of the energy per site with respect to under a fixed . As shown in the step 6 of the section III, we use a subtraction method for the convergence acceleration to obtain the energy per site in the limit . If this treatment becomes a proper method for the convergence acceleration, can rapidly converge to the energy per site of . Figure 8 shows the energies per site of finite systems and the subtracted energies of the YC8 kagome cylinder as a function of the reciprocal of the cylinder length with . The energies per site have an almost linear dependence with respect to depicted as the black dashed line. We find the energy per site agrees with the converged values of subtracted energies up to , where the error comes from a common cancellation of significant digits due to the subtraction analysis. Using this subtraction method, thus, we can avoid a careful extrapolation of the energy per site with respect to the length of the cylinder.
Secondly, to obtain the true ground state energy, we have to extrapolate to the limit . In this algorithm, we define a truncation error , where
(16) 
and extrapolate to the limit . This truncation error is reduced by increasing and must be zero in the limit . The reasons for taking and maximum of are as follows:

Reflecting the site period structure of the system, the value has same periodicity with respect to as shown in Fig. 9. We assume that the biggest truncation error mainly determines a quality of MPS.
If this definition of truncation error is appropriate, we can expect that the expectation values of operators are well fitted with an quadratic form of in the region of , as discussed in Ref. [White and Chernyshev, 2007].
Using each pair of energy per site and truncation error up to and , we extrapolate to the limit . The process of the extrapolations with respect to the truncation error is shown in Fig. 10. The linear fit and quadratic fit are denoted by the black line and the dashed red carve, and give extrapolated values and , respectively, where errors are due to the standard deviation of the fits. Both values are consistent with the reported value with up to Depenbrock et al. (2012). Thus, the combination of Hida’s infinitesize DMRG and McCulloch’s wave function prediction effectively searches the ground state energy per site of the kagome cylinder in the thermodynamic limit without fDMRG procedures.
iv.3 bond strength on nearest neighbor spins
In order to better understand the convergence behavior of the parallel iDMRG, next, we discuss bond strength of the nearest neighbor spins on the YC8 kagome cylinder as one of typical local observables. Our parallel iDMRG can predict and assume an infinitesize MPS (iMPS) with 12site period structures like Eq. (12) as a variational wave function. Therefore, the observables are equivalent to each other in our MPS, where a numbering of sites for YC8 cylinder is shown in Fig. 11 (a). Also, if the numerical calculation is executed exactly, four observables are the same, reflecting a translational symmetry along the circumferential direction of the cylinder. However, in our parallel iDMRG, these equalities are not the same because of the finite effect. Figure 11 (b) shows the difference of the bond strength from a mean values of , namely
(17)  
(18) 
where is the arithmetic average of . The difference approaches zero in the limit . Actually, we confirm the extrapolated values by the quadratic fits for data of are less than . This behavior is consistent with the fact that the translational symmetry along the circumferential direction must be recovered to the limit . Thus, it is allowed to focus only on if we discuss the values to the limit .
Due to the twotype of translational symmetries discussed above, we estimate only a set of bond strength to discuss the bond strength of nearest neighbor spins on the YC8 cylinder. Figure 11 (c) shows the bond strength of the set versus the truncation error . Broken curves are obtained by the quadratic fitting for data of . Extrapolated values to the limit can be grouped into two, namely and . The configuration of strength of nearestneighbor spins corresponding with this result are shown in 11(a). It is found that the bond intensity along the circumferential direction are slightly weaker than that along the diagonal bonds of the cylinder. A similar configuration result is reported in the XC8 kagome cylinder Yan et al. (2011).
iv.4 Spinspin correlations and correlation length
As the final demonstration, we estimate a spinspin correlation function along the YC8 cylinder and determine its correlation length. We set , and sweep along the axis of the cylinder [see Fig. 11(a)]. Every extrapolated value is given by the quadratic fitting as we discussed above and fitting error of are smaller than the symbol sizes shown in Fig. 12. The absolute values of the correlation function decays exponentially with respect to the distance between site and . We fit the data for with an exponential function and obtain the correlation length , where the error is the standard deviation of the fit. The spinspin correlation length of the YC8 cylinder with finite length up to is already evaluated by the nonAbelian DMRG Kolley et al. (2015) (inset of Fig. 12). We confirm that the correlation length of obtained by our parallel iDMRG agrees with the extrapolated value from the data of and with the linear fit.
V Conclusion and Remarks
We have investigated a parallel iDMRG method applied to onedimensional quantum systems with a longperiod structure. This parallel iDMRG is based on the Hida’s iDMRG Hida (1996) for 1D random quantum systems and the McCulloch’s wave function prediction McCulloch (). Numerical efficiency of our parallel iDMRG is demonstrated in a spin1/2 Heisenberg kagome cylinder. We have succeeded in obtaining proper observables, including the ground state energy per site, bond strength on nearest neighbor spins, spinspin correlation functions and its correlation length, without using fDMRG.
Finlay, several remarks are in order. First, the Hida’s iDMRG is intimately related to the realspace parallel DMRG Stoudenmire and White (2013). Figure 13 shows the whole picture of the realspace parallel DMRG, starting from the Hida’s iDMRG, where the diagrams, including overbraces and arrows, have the same meaning as shown in Fig. 3 and Fig. 5. The green shaded region is identical with the Hida’s iDMRG. In the procedures of Figure 13, the initial MPS is no longer needed to starting the parallel DMRG calculations.
Second, the physical back ground of the wave function prediction in iDMRG is wellunderstood from the view point of the 2D classical vertex model Ueda et al. (2010). Applying the quantumclassical correspondence discussed in Ref. [Ueda et al., 2010], we can easily find that our parallel iDMRG algorithm is also applicable to analyze 2D classical vertex models with arbitrary periodic structures along only the horizontal (vertical) direction.
Third, our parallel iDMRG is compatible with other parallel algorithms, such as parallelizing over different terms in the Hamiltonian Chan (2004) and block diagonalization of the matrix with respect to the quantum number Hager et al. (2004); Kurashige and Yanai (2009).
We expect that our parallel iDMRG and its extension are widely used to various other quantum systems.
Acknowledgements.
The author thanks S. Onoda, T. Nishino, and S. Yunoki for fruitful discussions. The work was partially supported by GrantsinAid for Scientific Research under Grant No. 25800221 from Japan Society for the Promotion of Science. The computation was performed in the computing facilities at the HOKUSAIGreat Wave system of RIKEN.References
 S. R. White, Phys. Rev. Lett. 69 (1992) 2863; Phys. Rev. B 48 (1993) 10345 .
 U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
 U. Schollwöck, Annals of Physics 326, 96 (2011).
 M. Fannes, B. Nachtergaele, and R. F. Werner, Comm. Math. Phys. 144, 443 (1992).
 S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995).
 S. Rommer and S. Östlund, Phys. Rev. B 55, 2164 (1997).
 S. R. White, Phys. Rev. Lett. 77, 3633 (1996).
 S. R. White, Phys. Rev. B 72, 180403 (2005).
 E. Stoudenmire and S. R. White, Annu. Rev. Condens. Matter Phys. 3, 111 (2012).
 C. Hubig, I. P. McCulloch, U. Schollwöck, and F. A. Wolf, Phys. Rev. B 91, 155115 (2015).
 T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 64, 4084 (1995).
 K. Ueda, T. Nishino, K. Okunishi, Y. Hieida, R. Derian, and A. Gendiar, J. Phys. Soc. Jpn. 75, 014003 (2006).
 H. Ueda, T. Nishino, and K. Kusakabe, J. Phys. Soc. Jpn. 77, 114002 (2008).
 I. P. McCulloch, arXiv:0804.2509 .
 H. Ueda, A. Gendiar, and T. Nishino, J. Phys. Soc. Jpn. 79, 044001 (2010).
 S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
 K. Hida, J. Phys. Soc. Jpn. 65, 895 (1996).
 E. M. Stoudenmire and S. R. White, Phys. Rev. B 87, 155137 (2013).
 S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012).
 F. Kolley, S. Depenbrock, I. P. McCulloch, U. Schollwöck, and V. Alba, Phys. Rev. B 91, 104418 (2015).
 K. Hida, J. Phys. Soc. Jpn. 66, 3237 (1997a).
 K. Hida, J. Phys. Soc. Jpn. 66, 330 (1997b).
 K. Hida, Phys. Rev. Lett. 83, 3297 (1999).
 I. P. McCulloch and M. Gulácsi, EPL (Europhysics Letters) 57, 852 (2002).
 S. R. White and A. L. Chernyshev, Phys. Rev. Lett. 99, 127004 (2007).
 G. K.L. Chan, J. Chem. Phys. 120 (2004).
 G. Hager, E. Jeckelmann, H. Fehske, and G. Wellein, J. Comput. Phys. 194, 795 (2004).
 Y. Kurashige and T. Yanai, J. Chem. Phys. 130, 234114 (2009).