Implementation of rigorous renormalization group method for ground space and low-energy states of local Hamiltonians
IImplementation of rigorous renormalization group method for ground space andlow-energy states of local Hamiltonians
Brenden Roberts, ∗ Thomas Vidick, and Olexei I. Motrunich
Institute for Quantum Information and Matter,California Institute of Technology, Pasadena, CA 91125 (Dated: February 14, 2018)The success of polynomial-time tensor network methods for computing ground states of certainquantum local Hamiltonians has recently been given a sound theoretical basis by Arad et al. . Theconvergence proof, however, relies on “rigorous renormalization group” (RRG) techniques which dif-fer fundamentally from existing algorithms. We introduce a practical adaptation of the RRG proce-dure which, while no longer theoretically guaranteed to converge, finds MPS ansatz approximationsto the ground spaces and low-lying excited spectra of local Hamiltonians in realistic situations. Incontrast to other schemes, RRG does not utilize variational methods on tensor networks. Rather, itoperates on subsets of the system Hilbert space by constructing approximations to the global groundspace in a tree-like manner. We evaluate the algorithm numerically, finding similar performance toDMRG in the case of a gapped nondegenerate Hamiltonian. Even in challenging situations of criti-cality, or large ground-state degeneracy, or long-range entanglement, RRG remains able to identifycandidate states having large overlap with ground and low-energy eigenstates, outperforming DMRGin some cases. I. INTRODUCTION
Many important techniques for solving lattice modelsin condensed matter physics take the form of tensor net-work algorithms. The seminal such method is White’sdensity matrix renormalization group (DMRG), an op-timization algorithm on the matrix product state (MPS)ansatz of quantum wavefunctions developed as a con-trollable method improving on Wilson’s numerical renor-malization group for impurity systems. DMRG remainsthe most versatile procedure in its class. It has beenheavily used to numerically solve low-dimensional quan-tum models; an early instance was the Haldane phase inthe Heisenberg chain. More recently, a related approachwas used in the classification of all gapped phases in onedimension. Other related techniques include the ten-sor renormalization group and tensor network renormal-ization, which utilize a two-dimensional coarse-grainingprocess to solve quantum systems in one dimension bythe quantum-to-classical correspondence.
Other vari-ational algorithms operate on the multiscale entangle-ment renormalization ansatz (MERA), which efficientlyrepresents states exhibiting logarithmic violation of thearea law by encoding correlations at all scales in an op-timized quantum circuit of logarithmic depth.
In this paper, we present an alternative approach tothe solution of local Hamiltonians in one dimension (1D).The rigorous renormalization group (RRG) is a recentalgorithm developed as part of a proof of the tractabil-ity of computing ground states of gapped local Hamil-tonians in 1D.
The proof employs techniques first in-troduced to establish an improved one-dimensional arealaw. Broadly, the strategy is as follows: partition thesystem into small initial blocks, and, focusing on theHilbert space of the blocks individually, identify sets ofstates that are “extendable” to the rest of the system tocreate a good approximation to the system-wide ground space. This property is termed viability , and formallydefined in Eq. (1). The identification of viable sets is ac-complished with an approximate ground state projector (AGSP), an operator approximately filtering out highlyexcited states on the entire system, whose support is re-stricted to perform this filtering within each block indi-vidually. In this way RRG deviates from a traditionalreal-space blocking scheme, in which each block does nothave access to global information. The next step is tomerge the identified viable sets on adjacent blocks, ob-taining states supported on blocks of larger size. How-ever, this step and the local application of the AGSPresult in an untenable blow-up of the number of states,so a reduction step is performed, returning the numberof states per block (now comprising two blocks of thesmaller size) to a constant value. This procedure is it-erated, merging blocks in a tree-like manner, and at thefull system scale, the identified states are shown to closelyapproximate the low-energy space. In the present work we adapt these techniques to spec-ify a concrete RRG procedure allowing for the explicitcomputation of ground and low-energy states of localHamiltonians. This requires making allowance for com-putational limitations, and generally our modificationsoperate outside of the regime of rigorous guarantee. Still,our algorithm presents a conceptually new approach tothis task. We emphasize that the use of the word “rig-orous” is in reference to the title of Arad et al. , ratherthan in order to establish a contrast with other tensornetwork algorithms.The main conceptual departure of this algorithm fromexisting tensor network methods is that RRG operates onviable sets of states supported on blocks, rather than onvariational states in the full Hilbert space. Two impor-tant features arise from this distinction. First, no localenergy minimization on a particular ansatz state is per-formed. Even though in the RRG procedure described a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b here the basic operations are performed on MPS com-prising an approximate basis of the viable sets, the MPSobjects themselves are incidental, and the concerns aris-ing from the MPS ansatz (e.g., gauge choice, truncation)are external to the fundamental algorithm.Second, the physical degrees of freedom are not coarse-grained. The objective of a coarse-graining strategy is tolimit the dimensionality of the Hilbert space at increasingscale by the introduction of renormalized degrees of free-dom, determined by some local rule, specifying a smallereffective Hilbert space. Instead, RRG achieves this goalby maintaining viable sets of constant dimension at alllevels of the algorithm hierarchy. These processes cannotbe considered equivalent, as the RRG step of applying theAGSP operator changes the relationship between scalesin a complicated way, and does not match the intuition ofan “RG flow” in a small number of parameters. However,this method still allows for fully controllable systematicimprovements in accuracy.The structure of this paper is as follows. We first givea detailed, self-contained description of our algorithm inSec. II. (We refer the reader familiar with the theoret-ical RRG paper to App. A for a precise discussion ofthe differences between the proof and the present work.)We provide an extended discussion of numerical resultsin Sec. III. Given its origins as a highly technical theo-retical algorithm developed in order to obtain provableguarantees, the RRG method performs surprisingly well,often matching the results of standard DMRG implemen-tations and outperforming them in certain difficult casesexhibiting degenerate ground spaces or highly entangledground states. Finally, we conclude and give directionsfor further work in Sec. IV. II. OPERATION OF ALGORITHMA. Overview and notation
The steps of RRG as implemented are listed in Fig. 1for reference and are discussed in detail in subsequentsections. A visual schematic is shown in Fig. 2. Our nota-tion is as follows. Let H = (cid:80) N − i =0 h i be a 2-local Hamil-tonian on a chain of N qubits, with term h i acting onsites i and i + 1. (The generalization to k -local Hamilto-nians and qudits is straightforward.) Denote the Hilbertspace of the system by H , and refer to the low-energyeigenspace of H as T . Let n be a parameter specifying thesize of initial regions of the system, and assume N/n is apower of 2. For each m = 0 , , . . . , log ( N/n ), partitionthe N -site system into contiguous blocks of equal length2 m n . Call these J λm = { λ m n, . . . , ( λ + 1)2 m n − } , for λ = 0 , , . . . , N/ (2 m n ) −
1. The Hilbert space associ-ated with J λm is denoted H λm , and H = (cid:78) λ H λm .Let H λm be the block Hamiltonian on J λm , compris-ing all terms acting only on sites in J λm and excludingboundary terms. Explicitly, H λm = (cid:80) i ∈ J λm ∗ h i , where J λm ∗ = { λ m n, . . . , ( λ + 1)2 m n − } . B. Initialization
The first step is to construct an approximate groundstate projector (AGSP) K , whose action on states in H increases overlap with T , the low-energy subspace of H .Many constructions of AGSP are possible. In the in-terest of efficiency, we use an AGSP obtained as an ap-proximation to a thermal operator at temperature t/k , K ≈ e − kH/t , t, k >
0. Let Q t denote a matrix prod-uct operator (MPO) approximating the thermal opera-tor e − H/t at temperature t ; procedures such as a Trotterdecomposition or cluster expansion can be used to effi-ciently compute this MPO. The AGSP is then obtainedas a power of Q t , contracting the product on the physicalindices k times.Because the AGSP must later be divided into operatorsacting on individual blocks, to compute Q t requires con-traction of the tensor network having terms of the form e − h i /t . After each contraction an SVD is performed be-tween site indices, and the MPO is truncated by eliminat-ing low-weight Schmidt vectors across each bond. Heretruncation is meant in the sense of MPS truncation, rep-resenting the MPO as a state in a higher-dimensionallocal Hilbert space. This amounts to using the Frobe-nius norm to order the terms arising from the SVD, andmay not be an optimal way to approximate operators;we address this issue in more detail later.The second step in the initialization is to identify setsof states V λ ⊂ H λ , for λ = 0 , , . . . , N/n −
1, of con-stant dimension s , where s is a parameter of the algo-rithm which bounds the dimension of the sets manipu-lated throughout. We use the term “viable sets” for the V λ (generally, for V λm ⊂ H λm ) because the intent of the al-gorithm is that each V λm be extendable to include a goodapproximation to the global low-energy eigenspace T .That is, each set V λm is chosen such that if H = H λm ⊗ H λm ,then V λm ⊗ H λm contains a subspace which is a good ap-proximation to T . We identify a set V λm as δ -viable if P T P V λm ⊗ H λm P T ≥ (1 − δ ) P T , (1)where P T is a projector onto a subspace T . More con-cretely, consider the case of a non-degenerate globalground space T = Span {| τ (cid:105)} , | τ (cid:105) ∈ H . The viability of theset V λm is given by δ = 1 − max | x (cid:105)∈ V λm ⊗ H λm |(cid:104) τ | x (cid:105)| , (2)where | x (cid:105) = (cid:80) j a j | v j (cid:105)| v j (cid:105) for a collection of states {| v j (cid:105)} ⊂ V λm along with coefficients a j , and states {| v j (cid:105)} arbitrary in the Hilbert space of the sites in the com-plement. It will be shown in Sec. II C that one neednever explicitly compute the {| v j (cid:105)} . For the case thatdim( T ) > δ is obtained by taking the smallest value The input is a local Hamiltonian H acting on N qubits, specified by an MPO. Let n , s and D be input parameters.1. Initialize:(a) Construct approximate ground state projector (AGSP) K from Hamiltonian H (b) Partition system into contiguous blocks of length n , denoted J λ , λ = 0 , , . . . , N/n −
1. Obtain s -dimensionallow-energy eigenspace V λ of block Hamiltonian H λ for each λ .2. For m = 0 , , . . . , log ( N/n ) −
1, denoting an “RG step” or scale factor:(a) Expand: for λ = 0 , , . . . , N/ (2 m n ) − D operators { A λm,r } r =1 ,...,D from the AGSP K , acting on subsystem J λm . Operate on the viable set,taking V λm → W λm ≡ { A λm,r V λm } r , where dim( W λm ) ≤ sD .ii. Compute the restriction of the block Hamiltonian H λm to W λm ⊂ H λm .(b) Reduce: For λ = 0 , , . . . , N/ (2 m n ) − W λm ⊗ W λ +1 m ⊂ H λ/ m +1 , supported on qubits in J λ/ m +1 = J λm ∪ J λ +1 m .Compute the restriction of H λ/ m +1 to the tensor product set.ii. Obtain s -dimensional low-energy eigenspace of the restriction of H λ/ m +1 to the tensor product space. Use theeigenstates as a basis for viable set V λ/ m +1 in iteration m + 1.3. At m = m ∗ = log( N/n ), the viable set V m ∗ is a candidate for the low-energy space T supported on the full system.Figure 1. Outline of the implemented RRG algorithm. of the maximum in (2), over all | τ (cid:105) ∈ T . The goal of thealgorithm is to construct the viable sets V λm in such away that they are indeed δ -viable for some constant δ less than 1 for all scales m . Note that a small value of δ corresponds to a better approximation, in contrast withmeasures like overlap. We emphasize that the viabilityparameter is not explicitly computed by the algorithm.Instead, it provides a useful metric to evaluate perfor-mance, both in terms of the theoretical results and interms of experimental investigations for cases where wewish to compare with other methods providing estimatesfor the ground space (such as when exact diagonalizationis possible).If n is chosen to be small enough, generic operatorson H λ can be exactly diagonalized. In the initializationstep, the initial viable set V λ is specified to be the spanof the s eigenvectors of H λ of lowest energy, obtained byexact diagonalization. C. Iteration over scale
The algorithm proceeds through a tree-like hierarchy,the levels of which are specified by a scale parameter m = 0 , , . . . , log ( N/n ). At scale m , block J λm consistsof 2 m n sites and the region index λ runs from 0 to N/ (2 m n ) −
1. Note that although the scale of the algo-rithm is increasing, we do not eliminate any of the phys-ical degrees of freedom. At each step we assume that theprevious level has produced a viable set V λm with basis {| v q (cid:105)} q =1 ,...,s represented by MPS, for every λ .The algorithm performs two steps. The first step is the expansion of the viable set, which has the effect of im-proving the viability parameter δ as defined in Eq. (1). This is accomplished using the AGSP constructed in theinitialization step as follows. Let J λm,L denote the qubitsto the left of J λm , and J λm,R those to the right. (Generally J λm has two boundaries with its complement J λm,L ∪ J λm,R .The system-edge cases follow immediately.) Consider theMPO representation of the AGSP K , whose elementarytensors are collections of operators on the local Hilbertspace, as an MPS. The Schmidt decomposition of K across the left boundary, separating J λm,L from J λm ∪ J λm,R ,produces a virtual index of dimension ζ : K = (cid:88) α<ζ σ α L α M α . (3)The L α are the left Schmidt vectors and the M α theright—which are operators on J λm ∪ J λm,R —each with acorresponding Schmidt coefficient σ α . The Schmidt de-composition may then be obtained for each of the M i across the boundary between J λm and J λm,R , producing avirtual index of dimension ξ . That is, M α = (cid:88) β<ξ ν αβ A αβ R αβ (4)Each A αβ is an operator on H λm , with weight γ αβ = σ α ν αβ in the expansion of K . For clarity we makethe algorithm variables explicit: A λm,αβ . Now let D > V λm , act on it with the D operators A λm,r , r = ( α, β ), having highest weight γ r = γ αβ . That is, take V λm → W λm = Span( { A λm,r | v q (cid:105)} r,q ),which we refer to as an expanded viable set with dimen-sion bounded by sD .One expects this operation to produce a set W λm of bet-ter viability than V λm because the A λm,r operators together 𝑊𝑊 𝐻𝐻 𝐻𝐻 𝑉𝑉 𝑊𝑊 𝑊𝑊 𝑉𝑉 𝑉𝑉 𝐴𝐴 𝐴𝐴 ReduceMergeExpand 𝐴𝐴 ReduceMergeExpand 𝐴𝐴 ⟩|𝑣𝑣 𝑎𝑎 ′ , ⟩|𝑣𝑣 𝑏𝑏 ′ , ⟩|𝑣𝑣 𝑎𝑎 ′′ , ⟩|𝑣𝑣 𝑏𝑏 ′′ ⟩|𝑣𝑣 𝑎𝑎 ′ ⟩|𝑣𝑣 𝑐𝑐 ′ , ⟩|𝑣𝑣 𝑎𝑎 ′ ⟩|𝑣𝑣 𝑑𝑑 ′ , ⟩|𝑣𝑣 𝑎𝑎 ′ ⟩|𝑣𝑣 𝑐𝑐 ′′ , ⟩|𝑣𝑣 𝑎𝑎 ′ ⟩|𝑣𝑣 𝑑𝑑 ′′ ,⟩|𝑣𝑣 𝑏𝑏 ′ ⟩|𝑣𝑣 𝑐𝑐 ′ , ⟩|𝑣𝑣 𝑏𝑏 ′ ⟩|𝑣𝑣 𝑑𝑑 ′ , ⟩|𝑣𝑣 𝑏𝑏 ′ ⟩|𝑣𝑣 𝑐𝑐 ′′ , ⟩|𝑣𝑣 𝑏𝑏 ′ ⟩|𝑣𝑣 𝑑𝑑 ′′ ,⟩|𝑣𝑣 𝑎𝑎 ′′ ⟩|𝑣𝑣 𝑐𝑐 ′ , ⟩|𝑣𝑣 𝑎𝑎 ′′ ⟩|𝑣𝑣 𝑑𝑑 ′ , ⟩|𝑣𝑣 𝑎𝑎 ′′ ⟩|𝑣𝑣 𝑐𝑐 ′′ , ⟩|𝑣𝑣 𝑎𝑎 ′′ ⟩|𝑣𝑣 𝑑𝑑 ′′ ,⟩|𝑣𝑣 𝑏𝑏 ′′ ⟩|𝑣𝑣 𝑐𝑐 ′ , ⟩|𝑣𝑣 𝑏𝑏 ′′ ⟩|𝑣𝑣 𝑑𝑑 ′ , ⟩|𝑣𝑣 𝑏𝑏 ′′ ⟩|𝑣𝑣 𝑐𝑐 ′′ , ⟩|𝑣𝑣 𝑏𝑏 ′′ ⟩|𝑣𝑣 𝑑𝑑 ′′⟩|𝑣𝑣 𝑖𝑖 , ��𝑣𝑣 𝑗𝑗 ��𝑣𝑣 𝑖𝑖 ′ , ��𝑣𝑣 𝑗𝑗 ′ , ��𝑣𝑣 𝑖𝑖 ′′ , ��𝑣𝑣 𝑗𝑗 ′′ ��𝑣𝑣 𝑚𝑚 , ��𝑣𝑣 𝑛𝑛 ��𝑣𝑣 𝑎𝑎 , ��𝑣𝑣 𝑏𝑏 𝑉𝑉 𝐴𝐴 𝐴𝐴 ⟩|𝑣𝑣 𝑐𝑐 ′ , ⟩|𝑣𝑣 𝑑𝑑 ′ , ⟩|𝑣𝑣 𝑐𝑐 ′′ , ⟩|𝑣𝑣 𝑑𝑑 ′′ ��𝑣𝑣 𝑐𝑐 , ��𝑣𝑣 𝑑𝑑 ⟩|𝑣𝑣 𝑖𝑖 ′ ⟩|𝑣𝑣 𝑘𝑘 ′ , ⟩|𝑣𝑣 𝑖𝑖 ′ ⟩|𝑣𝑣 𝑙𝑙 ′ , ⟩|𝑣𝑣 𝑖𝑖 ′ ⟩|𝑣𝑣 𝑘𝑘 ′′ , ⟩|𝑣𝑣 𝑖𝑖 ′ ⟩|𝑣𝑣 𝑙𝑙 ′′ ,��𝑣𝑣 𝑗𝑗 ′ ⟩|𝑣𝑣 𝑘𝑘 ′ , ��𝑣𝑣 𝑗𝑗 ′ ⟩|𝑣𝑣 𝑙𝑙 ′ , ��𝑣𝑣 𝑗𝑗 ′ ⟩|𝑣𝑣 𝑘𝑘 ′′ , ��𝑣𝑣 𝑗𝑗 ′ ⟩|𝑣𝑣 𝑙𝑙 ′′ ,⟩|𝑣𝑣 𝑖𝑖 ′′ ⟩|𝑣𝑣 𝑘𝑘 ′ , ⟩|𝑣𝑣 𝑖𝑖 ′′ ⟩|𝑣𝑣 𝑙𝑙 ′ , ⟩|𝑣𝑣 𝑖𝑖 ′′ ⟩|𝑣𝑣 𝑘𝑘 ′′ , ⟩|𝑣𝑣 𝑖𝑖 ′′ ⟩|𝑣𝑣 𝑙𝑙 ′′ ,��𝑣𝑣 𝑗𝑗 ′′ ⟩|𝑣𝑣 𝑘𝑘 ′ , ��𝑣𝑣 𝑗𝑗 ′′ ⟩|𝑣𝑣 𝑙𝑙 ′ , ��𝑣𝑣 𝑗𝑗 ′′ ⟩|𝑣𝑣 𝑘𝑘 ′′ , ��𝑣𝑣 𝑗𝑗 ′′ ⟩|𝑣𝑣 𝑙𝑙 ′′ 𝐻𝐻 𝑉𝑉 𝑊𝑊 𝑊𝑊 𝑉𝑉 𝐴𝐴 𝐴𝐴 𝐴𝐴 𝐴𝐴 ⟩|𝑣𝑣 𝑒𝑒 ′ , ��𝑣𝑣 𝑓𝑓 ′ , ⟩|𝑣𝑣 𝑒𝑒 ′′ , ��𝑣𝑣 𝑓𝑓 ′′ ⟩|𝑣𝑣 𝑒𝑒 ′ ��𝑣𝑣 𝑔𝑔 ′ , ⟩|𝑣𝑣 𝑒𝑒 ′ ⟩|𝑣𝑣 ℎ ′ , ⟩|𝑣𝑣 𝑒𝑒 ′ ��𝑣𝑣 𝑔𝑔 ′′ , ⟩|𝑣𝑣 𝑒𝑒 ′ ⟩|𝑣𝑣 ℎ ′′ ,��𝑣𝑣 𝑓𝑓 ′ ��𝑣𝑣 𝑔𝑔 ′ , ��𝑣𝑣 𝑓𝑓 ′ ⟩|𝑣𝑣 ℎ ′ , ��𝑣𝑣 𝑓𝑓 ′ ��𝑣𝑣 𝑔𝑔 ′′ , ��𝑣𝑣 𝑓𝑓 ′ ⟩|𝑣𝑣 ℎ ′′ ,⟩|𝑣𝑣 𝑒𝑒 ′′ ��𝑣𝑣 𝑔𝑔 ′ , ⟩|𝑣𝑣 𝑒𝑒 ′′ ⟩|𝑣𝑣 ℎ ′ , ⟩|𝑣𝑣 𝑒𝑒 ′′ ��𝑣𝑣 𝑔𝑔 ′′ , ⟩|𝑣𝑣 𝑒𝑒 ′′ ⟩|𝑣𝑣 ℎ ′′ ,��𝑣𝑣 𝑓𝑓 ′′ ��𝑣𝑣 𝑔𝑔 ′ , ��𝑣𝑣 𝑓𝑓 ′′ ⟩|𝑣𝑣 ℎ ′ , ��𝑣𝑣 𝑓𝑓 ′′ ��𝑣𝑣 𝑔𝑔 ′′ , ��𝑣𝑣 𝑓𝑓 ′′ ⟩|𝑣𝑣 ℎ ′′⟩|𝑣𝑣 𝑘𝑘 , ⟩|𝑣𝑣 𝑙𝑙 ��𝑣𝑣 𝑘𝑘 ′ , ⟩|𝑣𝑣 𝑙𝑙 ′ , ��𝑣𝑣 𝑘𝑘 ′′ , ⟩|𝑣𝑣 𝑙𝑙 ′′��𝑣𝑣 𝑒𝑒 , ��𝑣𝑣 𝑓𝑓 𝑉𝑉 𝑊𝑊 𝐴𝐴 𝐴𝐴 ��𝑣𝑣 𝑔𝑔 ′ , ⟩|𝑣𝑣 ℎ ′ , ��𝑣𝑣 𝑔𝑔 ′′ , ⟩|𝑣𝑣 ℎ ′′ ��𝑣𝑣 𝑔𝑔 , ��𝑣𝑣 ℎ Figure 2. Schematic illustration of the RRG algorithm over several length scales m = 0 , ,
2. As shown, ( s, D ) = (2 ,
2) andblock size n = 5. Gray dots represent local Hilbert spaces, and the V , V , etc, are the viable sets over blocks of sites. Thelabeled vectors | v a (cid:105) , etc., are basis states for the viable sets and have no relationship between blocks at a given scale. Theaction of the projector operators A λm,r on the states is represented by primes and double primes (e.g., A , | v a (cid:105) = | v a (cid:48) (cid:105) and A , | v a (cid:105) = | v a (cid:48)(cid:48) (cid:105) ). These generate the expanded viable sets W , and so on. The Merge procedure obtains tensor productstates such as | v a (cid:48) (cid:105)| v c (cid:48) (cid:105) supported on two blocks, and the tensor product set is reduced in dimension via diagonalizing blockHamiltonians like H , producing a viable set supported on two of the previous blocks. are meant to increase overlap with the global low-energyspace T : this is the defining property of the AGSP. Moreprecisely, let {| v j (cid:105)} be a collection of states in V λm suchthat there exists {| v j (cid:105)} ∈ H λm such that for some coeffi-cients a j , the state | x (cid:105) = (cid:80) j a j | v j (cid:105)| v j (cid:105) has good overlapwith T . By construction, K | x (cid:105) has better overlap with T than | x (cid:105) . Using the decomposition of Eqs. (3) and (4), K | x (cid:105) = (cid:88) α,β,j γ λm,αβ a j A λm,αβ | v j (cid:105) ⊗ A λm,αβ | v j (cid:105) , (5)where A λm,αβ = L λm,α R λm,αβ . In this way the viability asdefined in Eq. (2) of the set V λm can be improved whileleaving both the states and the operators supported onthe complement H λm entirely implicit. If all operators A λm,αβ were applied to V λm , the resultingset would contain the collection of states { A λm,αβ | v j (cid:105)} ,which has improved viability. However, instead of ap-plying all A λm,αβ , which would lead to an unmanageableblow-up in the size of the viable set, we introduce an ap-proximation by selecting the D operators A λm,r of high-est weight γ r in order to obtain W λm . There is no formalguarantee that this is the best choice, as the Schmidt de-composition is based on the Frobenius rather than theoperator norm. In practice we found the choice to bequite reasonable: to observe the increase in viability ina nondegenerate gapped model, compare the V and W points in Fig. 3, and in a critical model in Figs. 5, 6.The second step performed at each scale m is that of reduction of the dimension of the expanded viable sets W λm and W λ +1 m to generate V λ/ m +1 . At the cost of a lossof viability, this step restores s -dimensionality, resultingin a viable set suitable to use at the next level. One firstperforms a merge operation on disjoint pairs of blocks( λ, λ +1), with λ = 0 , , . . . , N/ (2 m n ) −
2. Merging refersto computing the tensor product set W λm ⊗ W λ +1 m thathas support on sites J λm ∪ J λ +1 m . One obtains the viableset V λ/ m +1 , a subspace of H λ/ m +1 = H λm ⊗ H λ +1 m , from the s -dimensional low-energy eigenspace of the restriction of H λ/ m +1 to W λm ⊗ W λ +1 m . We note that this step differs fromits counterpart in the theoretical algorithm, which pro-ceeds via random sampling instead of deterministicallyselecting the lowest-energy eigenvectors of H λ/ m +1 , as wedo here. Our choice is based on efficiency considerationsdescribed below; see also Appendix A for further discus-sion. The effect of the operation on the viability of thereduced subspaces can be seen in Figs. 3, 5, and 6.The single viable set V m ∗ generated at m ∗ = log ( N/n )after the reduction step at scale m ∗ −
1, is a constant-dimensional δ -viable subspace with support on the fullsystem. The algorithm returns the s lowest-energy eigen-vectors of the restriction of H to W m ∗ − ⊗ W m ∗ − , whichcomprise a basis for this candidate subspace. D. Scaling and computational considerations
The accuracy with which RRG approximates low-energy eigenstates of H is controlled primarily by twoparameters, s and D . To recapitulate, s bounds the di-mension of the reduced viable sets at each step, and D controls the level of approximation in the application ofthe AGSP via the operators { A λm,r } , r = 1 , . . . , D . Bothparameters are reflected in the bound on the dimension sD of the expanded viable sets W λm .We review the steps in the algorithm and discuss theircomplexity scaling based on these parameters. In addi-tion to s and D , important parameters are the system size N and the bond dimensions χ for MPS and η for MPOthat are manipulated throughout. For physical Hamilto-nians it is reasonable to expect χ and η to be constantin the gapped case, and in gapless systems χ, η ∼ N .See Schollw¨ock for a discussion of the scaling of basicMPS operations. Note that the initial block size n onlyenters this analysis in determining the number of neces-sary layers log( N/n ).The initialization requires obtaining viable sets V λ ofthe Hilbert space H λ on the qubits J λ . For small enoughchoices of n the complexity of this step will be negligi-ble, so we omit it. Similarly, the computation of the fullAGSP K ≈ ( e − H/t ) k can be done efficiently via Trotterdecomposition, and is not an important bottleneck. Inorder to extract the operators { A λm,r } , r = 1 , . . . , D , theAGSP must be obtained as an MPO in canonical form,analogous to that used for MPS. To do so requires a se-quence of Ø( N ) SVD operations, each with cost Ø( η ).For the steps comprising the iterated procedure we give scaling results applicable at the final computational level m = m ∗ −
1. The first step is to apply K to each V λm by means of the Schmidt decomposition of K across theboundary separating J λm from its complement (cid:83) λ (cid:48) (cid:54) = λ J λ (cid:48) m .This yields a set of operators acting on H λm . Applying the D such operators of highest Schmidt weight to a basis ofthe subspace takes V λm → W λm , increasing the dimensionto sD . The total cost of contracting these MPS andMPO is Ø( sD N χ η ).The second step acts on disjoint pairs of neighboringregions, forming the tensor product of expanded viablesets: W λm ⊗ W λ +1 m , with dimension ( sD ) . We com-pute the matrix elements of the restriction of the blockHamiltonian to the tensor product set. The scaling ofthis step is Ø (cid:0) [( sD ) ] N χ (cid:1) . For local Hamiltoniansthe constant can be improved using the decomposition H λ/ m +1 = H λm + H λ +1 m + B λ/ m +1 = H λm + H λ +1 m + (cid:88) p B λm,p ⊗ B λ +1 m,p . (6)The operator B λ/ m +1 contains Ø(1) terms in H actingacross the boundary between J λm and J λ +1 m .Exact diagonalization of the restricted blockHamiltonian in the subspace has complexityØ([( sD ) ] ) = Ø( s D ). After this, the final stepis to explicitly compute the s lowest-energy eigenstates,which has a total cost Ø( s ( sD ) N χ ). These states areused as a basis for the viable set at the next iteration.From this coarse analysis it is clear that the limitingstep with respect to s and D is the diagonalization of therestricted block Hamiltonian. This step is not part ofthe original formulation, which specifies instead that thereduction of viable set dimension take place by randomlyselecting states from the tensor product set. The choiceof our variant is motivated by its effect on the entangle-ment of the intermediate basis states: low-energy excitedstates of a block Hamiltonian may display lower entangle-ment than states chosen randomly. In practice this lowers χ in some systems. It also demonstrates a different pos-sible interpretation of the parameter s , which during theiteration step implicitly defines an energy scale with re-spect to the restricted Hamiltonian. States having blockexcitation energy higher than this scale are inaccessibleto the algorithm for the purposes of the expansion step. III. NUMERICAL RESULTS
We now present results from RRG for some examplemodels with the following goals in mind. We first validatethe algorithm in a simple gapped nondegenerate systemin Sec. III A, demonstrating consistency with DMRG aswell as previous numerical and perturbation theory re-sults. In this case the states obtained by RRG are ofsimilar accuracy to those of DMRG, with run times afactor of 5–10 slower depending on s , D , and N . How-ever, we emphasize that it is not the objective of RRGto obtain a numerically precise ground state; rather, itis to accurately identify states having constant overlapwith the global low-energy subspace. One expects an op-timization algorithm to obtain a more precise state inthe absence of local energy minima or very flat energylandscapes, and for simple models we take the DMRGground state to be exact (in particular, using it to mea-sure viability δ ). The RRG candidate states may laterbe variationally optimized in order to achieve a particularaccuracy, but we do not modify the states here.Our next goal is to demonstrate the practical scalingof the algorithm’s performance and computational costsassociated with the subspace parameters ( s, D ). We usethe familiar case of the Ising model in the transverse fieldin Sec. III B, both away from and at criticality. We findthat for low values of these parameters, often surprisinglygood results can be obtained, with close to unity over-lap between DMRG and RRG ground state candidates.However, neither algorithm scales linearly with systemsize in the critical regime. Here the slowdown of RRG isno longer a simple numerical factor but becomes a signif-icant cost at larger system sizes (beyond a few hundredsites in our implementation) or for larger values of thealgorithm parameters.Finally, we consider somewhat more challenging mod-els demonstrating areas in which RRG may hold an ad-vantage. In Sec. III C we investigate the Bravyi-Gossetmodel, which has Ø( N ) ground state degeneracy, byobtaining a complete basis for the ground space. InSec. III D we consider the XY model with randomly-distributed couplings. The ground state of this model,the random singlet phase, displays long-range entangle-ment in that it supports algebraic decay of correlations.We compare the correlations present in the candidatestates of DMRG and RRG to exact results obtainedby the Jordan-Wigner transformation, finding that RRGmore accurately reproduces observables measured on thestate.All numerical results were obtained using the tensornetwork library ITensor for both the DMRG and RRGcomputations. In all of the following, a Trotter decompo-sition with 60 steps was used to obtain the tensor networkfor Q t ≈ e − H/t , with t = 10, and degree k = 8 used tocompute the AGSP K ≈ ( Q t ) k . Thus the effective tem-perature t/k is of order unity. For reasonable choices ofparameters the accuracy of the approximation Q t is nota limiting factor of the algorithm. Computations wereperformed on standard hardware on a single node of acomputing cluster, with only single threading for the re-ported run times. A single error parameter τ was used tocontrol MPS truncation in ITensor for both DMRG andRRG (usually τ ∼ − –10 − ); in most cases a morelenient value would drastically improve run times withlittle effect on accuracy. DMRG convergence was han-dled using a fixed number of sweeps ≥
20 and relyingon the internal diagonalization routine included in ITen-sor without any modifications specific to the individualsystems. Excited states were found iteratively in DMRG by adding projectors into previously-found states to theHamitonian and using random trial wavefunctions. Of-ten the average viability will be used as a metric; this issimply the average over region label λ of the viability δ of each viable set ( V λm or W λm ) at fixed level m . A. Nonintegrable Ising model
This model refers to a spin-1/2 Hamiltonian H = − J N − (cid:88) i =0 σ zi σ zi +1 − g N − (cid:88) i =0 σ xi − h N − (cid:88) i =0 σ zi . (7)For h (cid:54) = 0 the model is gapped with a nondegenerateground state, and admits no good quantum numbersdue to the longitudinal component of the field. A re-cent numerical study for the parameters ( J, g, h ) =(1 , − . , .
5) found the ground state energy density tobe ε /N ≈ − .
722 and the gap γ = 3 . V W V W V W V W V W V Subspace and scale m − − − − − − − − − A v e r a g e v i ab ili t y δ D = 1 D = 2 D = 3 s = 3 s = 5 Figure 3. (color online) Viability of sets V λm , W λm averagedover λ , for nonintegrable Ising model with N = 256 spins,obtained as the RRG algorithm progresses through the scalehierarchy. Data are shown for parameter values s = 3 , D = 1 , , We run the RRG algorithm for a fixed system size N = 256, initial block size n = 8, and track the aver-age viability δ of the viable sets V m and W m through thesequence of dimensional expansion and reduction at eachscale m (see Fig. 2). Each data point shown in Fig. 3 isthe average over λ at a given length scale m . The param-eters ( s, D ) are varied to demonstrate their influence onthe results. For gapped systems of this size both DMRGand RRG have run times scaling linearly with systemsize, however RRG runs more slowly by a factor of 5–10compared to DMRG. At N = 256, DMRG took 5 min-utes to converge s = 5 states (ground and four excited)and RRG ran in 30 minutes with ( s, D ) = (5 , Index of energy eigenvalue − − − − − − − E n e r gy RRGDMRG − − − − − Figure 4. Energy eigenvalues of the nonintegrable IsingHamiltonian for N = 48 within the subspace obtained byRRG for ( s, D ) = (52 , N = 320 and( s, D ) = (12 , The large improvement in viability from V to W is at-tributable to the AGSP, rather than simple increase indimension. Both dim( V ) = s and dim( W ) = sD areconstant in m and very small compared to the dimen-sions of the block Hilbert spaces. Choosing n vectorswithout bias from an M -dimensional space will producea subspace whose squared overlap with a specific vectoris of order n/M . Since M here is exponentially large, aconstant increase in n would not much affect measured vi-ability. Thus, the AGSP is an effective projector even atlow values of D , which we expect as the model is gapped.A consequence is that the accuracy of RRG for thelargest ( s, D ) is comparable to that of DMRG, but wedo not expect this to be a general feature. Recall fromSec. II B that the criterion the algorithm seeks to main-tain is that the measured viability δ of the V λm (and thusthe average viability) be bounded for all m by some con-stant δ ∗ <
1, rather than approaching unity exponen-tially in m . The viability of the W λm is not necessarilyspecified, but should be sufficiently good for the V m +1 viable sets at the next level to satisfy the bound. For as-sessing the performance of RRG, as in Fig. 3, one seeksthat δ be maintained away from 1 for the V m averages.The final s -dimensional viable sets V m ∗ ( V in Fig. 3)here and in the following examples display much betteraverage viability than that of the previous V m . This isgenerally true: at steps m < m ∗ the viable set is foundby diagonalizing a block Hamiltonian H λm , which omitsterms present in H . The low-energy eigenspace of thisoperator need not be close to T , the global low-energyspace. At m = m ∗ , however, the low-energy eigenspaceof H m ∗ = H coincides with T , resulting in minimal lossof viability from the dimensional reduction.By changing the parameters of RRG, we obtain can-didates for low-energy excited states. The ground state N RRG runtime (s) DMRG runtime (s)32 158 9448 337 13264 866 20896 1871 277128 3912 393Table I. Runtimes of DMRG and RRG for the transverse-field Ising model with
J/g = 0 .
6, using ( s, D ) = (5 , s = 5 states are found by DMRG. of this model is close to a uniform spin-up state, andthe excited band contains a spin-flip excitation. Underopen boundary conditions two nearly degenerate lower-energy states separate from the first band, correspondingto quasiparticles localized at either edge. We obtain thelow-energy spectrum for N = 48 with ( s, D ) = (52 , N = 320 with ( s, D ) = (12 , N both methods find the entire first excited band. Inthe larger system, the localized edge states are more dif-ficult for DMRG, and it does not consistently find theedge states in sequence. The RRG ground state energydensity at N = 320 is ε /N = − .
721 and the gap tothe excited band is γ = 3 . S = 0 .
01 bits,and of the states in the band to be S ≈ .
01 bits, consis-tent with qualitative understanding of these states. ForDMRG and RRG, ground states have bond dimension 4and excited states in the band have bond dimension 31.(The methods do not yield identical bond dimension inall cases.)
B. Transverse-field Ising model
Consider the same Hamiltonian in the regime h = 0;that is, the Ising model in a transverse field. Fig. 5 showsthe result as we approach the critical point J = g fromthe paramagnetic phase for N = 128, measuring aver-age viability throughout the algorithm. One observesa strong deterioration of the measured viability as thegap closes. Approaching the critical point, RRG takesincreasingly more time than DMRG to run: runtimesfor J/g = 0 . s and D at criticality in Fig. 6. The improvement in viabilitywith increasing D is less dramatic than seen in Fig. 3,corresponding to a flatter spectrum of Schmidt valuesacross the cuts between subsystems. Note in this casethat at the critical point, as the algorithm progresses,the average viability of the V m sets visibly approachesunity, in contrast to the gapped case, which appears to V W V W V W V W V Subspace and scale m − − − − − − − − − A v e r a g e v i ab ili t y δ J/g = 1 . J/g = 0 . J/g = 0 . J/g = 0 . Figure 5. Viability of sets V λm , W λm , averaged over λ , for thetransverse-field Ising model both away from and at critical-ity. The number of spins is N = 128. All data points weregenerated using parameter values ( s, D ) = (5 , V W V W V W V W V Subspace and scale m − − − − − − − − A v e r a g e v i ab ili t y δ D = 1 D = 2 D = 4 s = 4 s = 5 Figure 6. (color online) Viability of sets V λm , W λm , aver-aged over λ , for the transverse-field Ising model at criticality,obtained as the RRG algorithm progresses through the scalehierarchy. Data are shown for parameter values s = 4 , D = 1 , , maintain viability bounded away from 1. C. Bravyi-Gosset model
This model was initially introduced as a classificationscheme for frustration-free 2-local Hamiltonians. TheHamiltonian is H = N − (cid:88) i =0 | ψ (cid:105)(cid:104) ψ | i,i +1 , (8) where | ψ (cid:105) is a generic state on two qubits. Up to aglobal phase, such a state can be specified in the form | ψ (cid:105) = R ( θ ) (cid:16) p | (cid:105) + (cid:112) − p | (cid:105) (cid:17) , with R ( θ ) a rota-tion performed on the first qubit. As the spectrum isinvariant under global rotation, the Hamiltonian is fullyspecified by the two parameters θ ∈ [0 , π ), p ∈ [0 , / θ = 0, we may rewrite Eq. (8) in a morefamiliar notation: H = N − (cid:88) i =0 (cid:32) (cid:112) p (1 − p )2 (cid:0) σ xi σ xi +1 − σ yi σ yi +1 (cid:1) + 14 σ zi σ zi +1 (cid:33) + N − (cid:88) i =0 (cid:18) − p σ zi + 18 (cid:19) (9)That is, this model is equivalent to a particular XYZmodel in a fine-tuned field. For any value of p the sys-tem exhibits ( N + 1)-fold ground state degeneracy. Basisstates for the ground space can roughly be thought of ashaving two regions of differing magnetization, with an in-terface which can be located at any site with ground stateenergy ε = 0. (Refer to Bravyi and Gosset for a fulldescription.) Therefore the algorithm choice s ≥ N + 1is sufficient to obtain the full ground space. Index of energy eigenvalue − − − − − − − − − − E n e r gy RRGDMRG
Figure 7. Energy eigenvalues of Bravyi-Gosset model with N = 32 sites within the subspace obtained by RRG for( s, D ) = (36 , The low-energy spectrum obtained by RRG for thismodel at N = 32 is shown in Fig. 7, along with theDMRG results. We use p = 1 /
2; that is, | ψ (cid:105) is a Bellstate. Using ( s, D ) = (36 , τ , the truncation error of the MPS. In contrast,obtaining the full ground space of this model is challeng-ing for DMRG, which becomes hampered by candidatestates of very high entanglement, often requiring a bonddimension an order of magnitude larger than those ofRRG candidate states in order to achieve similar trunca-tion error. These not only are computationally intensiveto optimize, but also present DMRG with difficulty find-ing further excited states, as the modified Hamiltonianincludes nonlocal projectors. Thus, the candidate statesare not accurate eigenstates of the original Hamiltonian.This difficulty is evident in run times as well; to obtainthe data shown took 10 hours for RRG and 40 hours forDMRG. Here we use DMRG without taking into accountthe degenerate ground state manifold, and we considerthese results to be only a point of reference. Use of a spe-cialized approach like multiple targeting could improveaccuracy, or diagonalization of the original Hamiltonianwithin the subspace spanned by the DMRG candidatestates could recover much of the ground space; however,no such specialized approach is needed for RRG. D. Random XY model
The random XY model is an inhomogeneous spin-1/2system with Hamiltonian H = N − (cid:88) i =0 J i ( σ xi σ xi +1 + σ yi σ yi +1 ) , (10)where the position-dependent coupling constants J i aredrawn from a random distribution. If the logarithm ofthe distribution is broad, Dasgupta-Ma real-space renor-malization group analysis identifies the ground state asthe random singlet phase, in which pairs of spins formsinglet states at all length scales. This model istractable by the Jordan-Wigner transformation, whichmaps onto free spinless fermions. We use this system asa benchmark of algorithmic ability to encode long-rangecorrelations in the ground state. | i − j | − − − − | h σ z i σ z j i | ExactRRGDMRG
Figure 8. (color online) Disorder-averaged decay of correla-tions of candidate ground states of the random XY modelfor N = 128, as compared to exact results obtained throughthe Jordan-Wigner transformation. The predicted power-lawbehavior is indicated by the red line. We use the following distribution for the Hamiltonianterms: p ( J i ) = J i − (1 − ) , J i ∈ (0 , We fixΓ = 2, which is sufficiently broad that the ground stateis composed predominantly of localized singlet states on | i − j | − − − − | h σ z i σ z j i | ExactRRGDMRG
Figure 9. A typical “hard” instance contained in the disor-der average above, with energy gap γ ≈ − . This is suf-ficiently large for RRG to track the long-range correlationswith ( s, D ) = (4 , neighboring sites, along with spatially separated corre-lated qubits occurring at all length scales. As a metricwe use the average two-point correlation function (cid:104) σ zi σ zj (cid:105) as a function of separation r = | i − j | in the ground state,which is known to decay algebraically as r − . This quan-tity is compared to exact diagonalization results from theinhomogeneous free fermion description in Fig. 8.These results are intended to present a fair compari-son between DMRG and RRG. Both methods used un-restricted MPS bond dimension to achieve a truncationerror τ ≤ − . Typically the ground state bond dimen-sion is similar for both methods. The RRG parametersare ( s, D ) = (4 , s = 4 states and RRG took 8 hours tocomplete. The average is over 150 disorder realizations.The observed “saturation” of the correlations of Fig. 8to a noise floor arises from the structure of the low-energyexcited states. For a broad initial distribution, the en-ergy gap of a specific disorder Hamiltonian may be verysmall. For any method using MPS, a lower limit on thegap in order to distinguish the ground state (at energy ε ) is γ ∼ τ ε , below which the MPS truncation proce-dure will randomly select a vector from the low-lying sub-space. However, even for realizations with much largergaps a candidate ground state may include substantialcontributions from low-energy excited states. A singletof length l has energy scale ε ∼ e −√ l ; thus, the low-lyingstates involve excitations localized on the long-range en-tangled sites. Choosing a random superposition of theseamounts to white noise at long distances. Instances ofsuch Hamiltonians in the disorder average must necessar-ily eventually overwhelm the decay of correlations; herethe distribution of energy gaps is very broad on a logscale, so these cases are frequent. However in all casesthe RRG candidate state has Ø(1) overlap with the true0 E x a c t RR G
106 111 116 121 126106111116121126 D M R G Figure 10. (color online) Expectation value (cid:104) σ zi σ zj (cid:105) , wheresites i and j are given by the axes for i, j ∈ [106 , |(cid:104) σ zi σ zj (cid:105)| and runs from [ − ,
0] in all plots,with darker color indicating a higher value. The diagonal isomitted. Circles mark particular sites where differences be-tween exact results and candidate states are evident. This dis-order realization is the same “hard” instance shown in Fig. 9. ground state, and typically this overlap is greater than99%.For disorder-averaged correlations at short range upto | i − j | ≈
20, RRG reproduces algebraic decay of cor-relations matching the exact results. In contrast, theDMRG candidate states demonstrate stronger decay ofcorrelations. There is no systematic difference in MPS bond dimension between DMRG and RRG, indicatingthat RRG is not simply using additional resources, butis indeed more sensitive to long-range correlations.Independent of the saturation due to the energy gap,the disorder average comprises both “easy” and “hard”instances. In easy cases both DMRG and RRG matchthe exact results closely at all length scales. In the hardcases both algorithms obtain the correlations only ap-proximately, but DMRG appears to consistently under-estimate correlations. RRG does not demonstrate a ten-dency towards either enhanced or reduced correlations.We provide an example of the spatially averaged correla-tions from a hard disorder realization in Fig. 9. Fig. 10shows an example of measured correlations (cid:104) σ zi σ zj (cid:105) forvarious sites i, j ∈ [106 , (cid:104) σ zi σ zj (cid:105) where ( i, j ) are specified by the axes. Darkersquares indicate a larger magnitude of correlation be-tween these sites. We show the exact results, RRG, andDMRG, and indicate some particular pairs of sites whereeither DMRG (red) or RRG (green) differ visibly fromexact results. These variations in certain entangled sitestend toward reduced correlations in DMRG candidateground states; it is unclear how much additional sweepingis required to compensate. RRG shows similar inaccura-cies, but these are random, due to states missing fromcertain viable sets. Accurate correlations emerge in thedisorder-averaged value, and the performance on indi-vidual disorder realizations can be controllably improvedby tuning the dimension of the viable sets through theparameters s and D . IV. DISCUSSION
DMRG has long been the method of choice for numer-ical calculations involving ground states of 1D systems,and over time both its efficiency and range of applica-bility have gone through multiple improvements and ex-tensions. One of the main findings of our initial numer-ical investigation is that the RRG algorithm, developedfor theoretical purposes, can in fact be made quite effec-tive in practice, to the point of providing a potentiallyviable alternative to DMRG in certain cases of practi-cal interest. We stress that the choices of parametersthat we employ in our numerics are quite far from thetheoretically guaranteed regime. Additionally, many ofthe building blocks required for the proof have been al-tered in our implementation. Therefore the strict guar-antees no longer hold. Regardless, we find that RRG ob-tains ground state candidates having large overlap withthe true ground state in a variety of physically relevantmodels, and surpasses existing techniques in obtaininglow-energy excited states and ground states of particu-lar models demonstrating large degeneracy or long-rangeentanglement.Like another numerical scheme, time-evolving blockdecimation (TEBD), the RRG procedure is a projector1method, relying on operators extracted from the AGSPto guide the choice of states between scales. As a result,given a sufficiently accurate AGSP, RRG will not outputa part of the spectrum strictly excited above the groundspace. This is advantageous relative to variational ansatzmethods which may without warning converge to an ex-cited state rather than the ground state. (For example,if the energy landscape in Hilbert space has local min-ima or is very flat in the low-energy space, as is the casewith the random XY model of Sec. III D.) The downsidesto projector methods are that performance strongly de-pends on the gap and that a random initial state, eventaken from the manifold of low bond dimension MPS,has exponentially small overlap with the ground state.RRG circumvents the latter issue by never choosing atrial wavefunction on the entire system, but rather build-ing global states from wavefunctions supported on blockswhich already have good viability ; thus the projectionstep never has to overcome starting with an exponen-tially small overlap between the initial and the targetstate.At present the run times required by the algorithm re-main a challenge. Thus, the feasibility of RRG as a nu-merical method is essentially determined by the scalingdiscussed previously. This situation does invite futureimprovements. Some are immediate: for example, onemay exploit symmetries of a particular problem (say, re-flection symmetry across the middle of the system) in or-der to reduce duplication of work. Other improvementsto the current implementation are more technical. Forexample, as described here the management of subspacesis clumsy: operations such as addition of MPS necessi-tate keeping careful track of gauge and add computa-tional overhead for what is in principle a simple proce-dure. The use of data structures more appropriate tothese operations could ameliorate scaling problems in allsteps of the algorithm. Indeed, an advantage of RRG is precisely this flexibil-ity, to operate independently of a specific representationof states in Hilbert space. Here we have described anMPS RRG. In order to translate the logic to subspaceswhose basis states are described by MERA—as would benatural for critical phases—one needs only the ability toperform evaluation of observables and addition. The for-mer is a standard contraction which is highly efficient inMERA, and the latter can be seen as a variational processon overlaps, providing a straightforward interpretationas a MERA operation. Systems with periodic boundaryconditions also present an interesting generalization, asuntil the final level the steps of the algorithm are insen- sitive to the system boundaries, provided an appropriateAGSP is given. On a more speculative note, other tensornetwork ansatze may also be amenable: although it isnot known that the RRG algorithm scales efficiently inhigher dimensions, the hierarchical structure does gener-alize in an evident way and it may be the case that thealgorithm gives acceptable results for PEPS representa-tions of some two-dimensional systems.Our numerical results suggest situations in which RRGmay perform well relative to existing algorithms. Thefirst, informed by Sec. III A, is a case in which localizedand delocalized excitations lie close in energy. An opti-mization algorithm operating on local degrees of freedomin a sweeping pattern may exhibit a bias towards delocal-ized excitations, which allow for effective optimization onmany lattice sites. RRG is largely insensitive to such dis-tinctions. The second case is that of Sec. III C, exhibitinghighly degenerate ground states. The full ground spaceis more accurately found in its entirety by RRG thanDMRG. The iterative DMRG procedure of finding statesis susceptible to finding poor or highly entangled candi-dates, which reduce the accuracy of subsequent candi-dates. Such a limitation is not fundamental and couldlikely be eliminated by modification of the procedure;however no such modification is necessary for RRG. Fi-nally, in Sec. III D we observe in the random XY model inthe random singlet phase that long-range correlations areencoded more precisely in the ground state candidate ofRRG than of DMRG, influencing observables computedfor the state.The examples we provide illustrate specific propertiesindicating that a model may be well suited for RRG.However, very little is known about its more general per-formance: other systems with disorder, periodic bound-ary conditions, and higher dimensions all pose interestingchallenges and could constitute exciting new directionswithin this formalism.
ACKNOWLEDGMENTS
We acknowledge useful discussions with M. Fishmanand S. White’s research group, as well as with C. Whiteand C.-J. Lin. The numerical results were computed withthe ITensor library of E. Stoudenmire and S. White.This work was supported by the Institute for QuantumInformation and Matter, an NSF Physics Frontiers Cen-ter, with support of the Gordon and Betty Moore Foun-dation. Additional funding support was provided by theNSF through Grant DMR-1619696. ∗ [email protected] I. Arad, Z. Landau, U. Vazirani, and T. Vidick, Commu-nications in Mathematical Physics , 65 (2017). S. R. White, Phys. Rev. Lett. , 2863 (1992). F. Verstraete, V. Murg, and J. Cirac, Advances in Physics , 143 (2008). U. Schollw¨ock, Annals of Physics , 96 (2011). K. G. Wilson, Rev. Mod. Phys. , 773 (1975). S. R. White, R. M. Noack, and D. J. Scalapino, Phys.Rev. Lett. , 886 (1994). X. Chen, Z.-C. Gu, and X.-G. Wen, Physical Review B , 035107 (2011). M. Levin and C. P. Nave, Phys. Rev. Lett. , 120601(2007). Z.-C. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B ,205116 (2008). G. Evenbly and G. Vidal, Phys. Rev. Lett. , 180405(2015). G. Vidal, Physical Review Letters , 220405 (2007). G. Evenbly and G. Vidal, Physical Review B , 144108(2009). R. N. C. Pfeifer, G. Evenbly, and G. Vidal, Physical Re-view A , 040301 (2009). G. Evenbly and G. Vidal, Phys. Rev. B , 235102 (2010). G. Evenbly and G. Vidal, Phys. Rev. Lett. , 200401(2015). Z. Landau, U. Vazirani, and T. Vidick, Nature Physics , 566 (2015). I. Arad, Z. Landau, and U. Vazirani, Phys. Rev. B ,195145 (2012). M. Suzuki, Progress of theoretical physics , 1454 (1976). S. Bravyi and D. Gosset, Journal of Mathematical Physics , 061902 (2015). E. M. Stoudenmire and S. R. White, “ITensor - IntelligentTensor Library,” http://itensor.org . C.-J. Lin and O. I. Motrunich, Phys. Rev. A , 023621(2017). R. N. Bhatt and P. A. Lee, Phys. Rev. Lett. , 344 (1982). S.-k. Ma, C. Dasgupta, and C.-k. Hu, Phys. Rev. Lett. , 1434 (1979). C. Dasgupta and S.-k. Ma, Phys. Rev. B , 1305 (1980). D. S. Fisher, Phys. Rev. B , 3799 (1994). D. S. Fisher and A. P. Young, Phys. Rev. B , 9131(1998). Since submission of this manuscript, the au-thors have significantly advanced the state ofthe RRG code, which may be found online at https://github.com/brendenroberts/RigorousRG . Appendix A: Differences from Arad et al. In this appendix we give a detailed account of the mainpoints of departure of our numerical procedure from thetheoretically guaranteed algorithm introduced in Arad et al. , giving heuristic justification for our choices. Werefer to the paper for a more thorough introduction tothe main concepts discussed here, such as the notion ofviable set and AGSP.For concreteness we base our comparison on the algo-rithm presented in Arad et al. for the case of a localHamiltonian with degenerate gapped ground space (As-sumption (DG)). The algorithm is stated as Algorithm 1in Arad et al. . It consists of two main steps, Gener-ate and
Merge . The two steps together recursively con-struct a sequence of viable sets V λm for an N -qubit localHamiltonian, where as in the main text m denotes a scaleparameter and λ indexes a subregion.
1. Generate
The goal of the
Generate step is to generate an MPOrepresentation for a suitable AGSP. In Arad et al. afresh AGSP is computed for each scale m and region λ .Given a decomposition H = H L ⊗ H λm ⊗ H R , a global AGSPis defined as K λm = T k ( ˜ H ), where ˜ H is a norm-reducedapproximation of H (which depends on the region decom-position) and T k a suitably scaled Chebyshev polynomialof degree k . The operators A λm,r are then computed froma specific decomposition of K λm across the left and rightboundaries, yielding D terms A λm,r such that the ex-pansion procedure V λm → W λm described in the main textis guaranteed to have a significant improvement on theviability parameter.Here we depart from the theoretical algorithm in twoimportant ways. First we use a simpler construction ofAGSP, which we expect to exhibit similar behavior butis more efficient to compute. Our AGSP takes the formof an approximation K ≈ e − kH/t obtained by Trotter de-composition. (In Arad et al. a similar approach is takento norm-reduce the parts of the Hamiltonian that lie inthe regions L, M and R but are a distance at least (cid:96) > et al. the properties ofthe Chebyshev polynomial are essential to establish thatthe AGSP has sufficiently low bond dimension across theboundaries of region M . Considering only the efficiencyin terms of improvement in viability, however, the use of e − kH/t over the whole chain gives similar guarantees.Using our simpler construction implies a loss of theo-retical control over the bond dimension D of the AGSPoperator across the left and right cuts. This entails asecond main point of departure from the theoretical al-gorithm, as a choice has to be made as to which operators A λm,r to keep. As described in the main text we proceedin a natural way by considering the MPO as an MPSand performing SVD operations to create virtual bondsbetween sites. We then make the choice of keeping op-erators associated with the D highest Schmidt weights.This choice is heuristic: the Schmidt weights control theFrobenius norm of the associated term A λm,r , rather thanthe operator norm of the resulting operator, as would bedesirable. The heuristic nevertheless proved effective: inpractice the magnitude of the Schmidt coefficients oftenfell off quickly, allowing for a relatively aggressive choiceof cutting point.
2. Merge process
The second step in the algorithm is called
Merge . Thegoal of this step is to combine two neighboring viablesets into a single viable set over the union of the tworegions, with similar approximation and size guarantees.The procedure is described as
Merge’ in Arad et al. .Merge’ is provided as input viable sets W λm and W λ +1 m defined over neighboring regions, and returns a viable set3 W λ/ m +1 defined over the union of the two regions. Mergeconsists of three steps: Tensoring , Random Sampling ,and
Error Reduction .1.
Tensoring:
This step is the same as in Arad et al. .2. Random Sampling:
Here as already mentioned inthe main text we depart from Arad et al. in animportant way. In Arad et al. a family of s vec-tors lying in W λm ⊗ W λ +1 m is obtained by randomsampling within the subspace. In practice this pro-cedure is very inefficient: (i) it requires performinghigh-weight (random) linear combinations of MPS,a step that is computationally expensive due tothe MPS renormalization procedure; (ii) the lin-ear combinations formed tend to be arbitrary, andin particular their MPS representations may havehigh MPS bond dimension, as each vector mayinclude an “irrelevant” (with respect to the low-energy eigenspace of the Hamiltonian) componentthat artificially inflates its complexity.Here we replace random sampling by a determin-istic choice of the s lowest-energy eigenvectors ofthe restriction of H to W λm ⊗ W λ +1 m . The idea isthat low-energy eigenstates are likely, due to thelocal structure of the Hamiltonian, to display lessentanglement. Indeed in practice this procedureis much more efficient, and yields MPS with lowerbond dimension, than the random sampling pro-posed in Arad et al. .However, there is a priori no reason for the low-energy eigenstates of the block Hamiltonian to forma viable set for the global low-energy space. A sim-ple heuristic argument can nevertheless be given toargue correctness of our procedure. Recall that theviability criterion Eq. (1) guarantees that the initial tensor product space supports a good approxima-tion to any ground state. Considering the Schmidtdecomposition of this approximation, each of theSchmidt vectors will have a certain energy with re-spect to the block Hamiltonian H λ/ m +1 , which maynot be minimal. The key is thus to argue thatvectors with high energy will not have an impor-tant contribution to the Schmidt decomposition ofthe ground state. In general approximation errorand energy difference can scale with the norm ofthe Hamiltonian, making the argument difficult.However, for the purposes of approximating theground space of a local Hamiltonian two elementsplay in our favor: first, locality of H , and second,the area law. The former allows to show that thelow-energy space of H is well-approximated by anapproximation of H with constant norm, so thatthe error blow-up mentioned above can be con-trolled (see Arad et al. , Proposition 3, for a precisestatement). The latter establishes that the groundstate has low bond dimension, so that few Schmidtvectors need to be considered (see Arad et al. ,Lemma 15, for details on how this can be used).Together these two properties provide a heuristicargument in favor of our modified procedure.3. Error Reduction:
The goal of this step is to im-prove the approximation quality of the viable set.We follow the procedure described in Arad et al. ,except that the operators { A λm,r } are generated dif-ferently, as already described.The final iteration is performed on two viable sets V m ∗ − and V1
The goal of this step is to im-prove the approximation quality of the viable set.We follow the procedure described in Arad et al. ,except that the operators { A λm,r } are generated dif-ferently, as already described.The final iteration is performed on two viable sets V m ∗ − and V1 m ∗ −1