16. Fixed-Point algorithms#

In this section we introduce two algorithms for approximating the ground state of local gapped Hamiltonians using matrix product state techniques. Approximating ground states in a variational manner boils down to minimizing

\[ \min_{\ket{\psi}\in D} \frac{\braket{\psi|H|\psi}}{\braket{\psi|\psi}},\]

over a restricted class of states \(D\). For simplicity, we will assume the Hamiltonian under consideration has an MPO representation of the form

which can encode interactions of arbitrary range as discussed in the previous section. In this formulation, approximating the ground state of \(H\) is equivalent to finding the MPS fixed point the MPO Hamiltonian corresponding to the eigenvalue \(\Lambda\) with the smallest real part,

In the algorithms discussed below we optimize over matrix product states of a fixed finite bond dimension. In the first algorithm known as DMRG (density matrix renormalization group) the states we consider are finite MPS, whereas the second algorithm VUMPS (variational uniform matrix product state algorithm), as the name suggests, optimizes over uniform MPS. Hence, VUMPS enables direct optimization in the thermodynamic limit, without breaking translation invariance.

Our exposition of DMRG closes follows the one in [Bridgeman and Chubb, 2017], and that of VUMPS closely follows the excellent set of lecture notes [Vanderstraeten et al., 2019].

16.1. DMRG#

Starting from a random MPS ansatz, DMRG tries to approximate the ground state by sequentially optimizing over all the MPS tensors one by one and sweeping through the chain, until convergence is reached. Let us discuss this algorithm in a bit more detail step by step.

16.1.1. Algorithm#

Let us consider a random ansatz, by taking random tensors \(\{A_1,A_2,...,A_L\}\), \(L\) being the number of sites. Fixing all tensors but the one at site \(i\), the local tensor \(A_i\) is updated according to

Though seemingly daunting we can turn this problem in a simple eigenvalue problem by making full use of the mixed gauge. By bringing all tensors on the right of \(A_i\) in the right canonical form and those to the left in left canonical form the denominator simply becomes \(\braket{A_i|A_i}\) and the update reduces to

Here the effective Hamiltonian \(\mathcal H_i\), defined as

encodes the effect of the full system Hamiltonian on the current center site \(i\). The variational problem of the local update can then be solved by finding the eigenvector of \(\mathcal{H}_i\) corresponding to the smallest real eigenvalue, and this repeatedly for every site sweeping back and forth through the chain, each time moving the orthogonality center of the MPS. At each update step a large part of the effective Hamiltonian can simply be reused, making the routine very efficient. Notice however that DMRG manifestly breaks translation invariance by updating one tensor at a time. As we will see, VUMPS does not suffer from this artefact.

From this brief explanation it should be clear that DMRG is a surprisingly simple algorithm. Nevertheless DMRG has proven itself time and time again, and is the most successful algorithm for variationally approximating the ground state of local gapped (1+1)d Hamiltonians. DMRG is implemented in MPSKit and can be called by DMRG().

16.1.2. Example#

Let us illustrate the use of DMRG in MPSKit by approximating the ground state of the transverse field Ising model. The Ising model is implemented in MPSKitModels as follows

\[ H = -J\left(\sum_{<i,j>} Z_i Z_j + \sum_i h_x X_i + h_z Z_i\right),\]

where we are free to choose the parameters \(J\), \(h_x\) and \(h_z\), and \(X\) and \(Z\) are the generators of \(\mathfrak{su}(2)\), and thus differ from the usual Pauli matrices by a factor of \(\frac{1}{2}\).

Let us consider 16 lattice sites, bond dimension 12, open boundary conditions and let’s stick to the default critical values of \(h_x=0.5\) and \(h_z=0\). Finding the ground state using DMRG then only takes a handful of iterations!

using TensorKit, MPSKit, MPSKitModels

d = 2 # Physical dimension
L = 16 # Length spin chain
D = 12 # Bond dimension

H = transverse_field_ising()

algorithm = DMRG(); # Summon DMRG
Ψ = FiniteMPS(L, ℂ^d, ℂ^D) # Random MPS ansatz with bond dimension D
Ψ₀,_ = find_groundstate(Ψ, H, algorithm);

┌ Info: DMRG iteration:
│   iter = 1
│   ϵ = 0.0005443461129231371
│   λ = -20.016387900453072 - 1.1529224368640848e-16im
└   Δt = 0.656620179
┌ Info: DMRG iteration:
│   iter = 2
│   ϵ = 2.1371523099817977e-7
│   λ = -20.016387900459865 + 1.140294036246174e-15im
└   Δt = 0.088715297

┌ Info: DMRG iteration:
│   iter = 3
│   ϵ = 6.595416787971344e-8
│   λ = -20.016387900460188 + 2.5945145280068877e-16im
└   Δt = 0.024979348
┌ Info: DMRG iteration:
│   iter = 4
│   ϵ = 2.669445530707093e-8
│   λ = -20.016387900460284 + 7.319051500185613e-17im
└   Δt = 0.026395541
┌ Info: DMRG iteration:
│   iter = 5
│   ϵ = 1.064795309649792e-8
│   λ = -20.016387900460295 + 7.659573627428371e-16im
└   Δt = 0.016926126
┌ Info: DMRG iteration:
│   iter = 6
│   ϵ = 4.236406898772371e-9
│   λ = -20.01638790046029 - 3.8200293026848695e-16im
└   Δt = 0.01972514
┌ Info: DMRG iteration:
│   iter = 7
│   ϵ = 1.6857942437367281e-9
│   λ = -20.016387900460288 + 1.3407727342998587e-15im
└   Δt = 0.013797992
┌ Info: DMRG iteration:
│   iter = 8
│   ϵ = 6.712202527793425e-10
│   λ = -20.016387900460302 - 6.13475308479391e-17im
└   Δt = 0.017379941

┌ Info: DMRG iteration:
│   iter = 9
│   ϵ = 2.674137185440222e-10
│   λ = -20.016387900460277 - 2.9640928893728194e-16im
└   Δt = 0.011330565
┌ Info: DMRG iteration:
│   iter = 10
│   ϵ = 1.0659748891009541e-10
│   λ = -20.016387900460288 - 5.55726973232157e-16im
└   Δt = 0.010242441
┌ Info: DMRG iteration:
│   iter = 11
│   ϵ = 4.251328009840777e-11
│   λ = -20.016387900460288 + 2.839604757125549e-16im
└   Δt = 0.017049266
┌ Info: DMRG iteration:
│   iter = 12
│   ϵ = 1.6962775156682162e-11
│   λ = -20.016387900460288 - 1.5176907562937642e-16im
└   Δt = 0.009095099
┌ Info: DMRG iteration:
│   iter = 13
│   ϵ = 6.769276367248231e-12
│   λ = -20.016387900460273 - 4.340862311951199e-16im
└   Δt = 0.008646554
┌ Info: DMRG iteration:
│   iter = 14
│   ϵ = 2.703510855181987e-12
│   λ = -20.016387900460263 - 3.773077811684178e-16im
└   Δt = 0.012864275
┌ Info: DMRG iteration:
│   iter = 15
│   ϵ = 1.0792915205545097e-12
│   λ = -20.016387900460266 - 8.46443216414052e-16im
└   Δt = 0.007770156
┌ Info: DMRG iteration:
│   iter = 16
│   ϵ = 4.315208327470975e-13
│   λ = -20.01638790046027 + 9.480960494303697e-16im
└   Δt = 0.007647437

┌ Info: DMRG summary:
│   ϵ = 2.0e-12
│   λ = -20.01638790046027 + 9.480960494303697e-16im
└   Δt = 1.203674018

16.2. VUMPS#

As mentioned above, VUMPS optimizes uniform MPS directly in the thermodynamic limit. Since the total energy becomes unbounded in this limit, our objective should be to rather minimize the energy density. When working in the mixed gauge, this minimization problem can be represented diagrammatically as

where we have introduced the left- and right fixed points \(F_L\) and \(F_R\) defined as

which obey the normalization condition

The VUMPS algorithm offers the advantage of global optimalization by design, since the algorithm, contrary to DMRG, does not rely on individual updates of local tensors.

Given a Hamiltonian of the form mentioned above and an intial random uniform MPS defined by \(\{A_L, A_R,C\}\), VUMPS approximates the ground state by finding an approximate solution to the fixed-point equations

\[\begin{split} A_C' = H_{A_C}(A_C), \\ C' = H_C(C), \\ A_C = A_LC = CA_R.\end{split}\]

A detailed derivation that these equations characterize the variational minimum in the manifold of uniform MPS is beyond the scope of these notes, but see [Vanderstraeten et al., 2019].

In these equations the effective Hamiltonians \(H_{A_C}\) and \(H_{C}\) acting on \(A_C\) and \(C\) respectively are given by is given by

The last equation then simply states that \(C\) intertwines the left - and right-orthonormal form of the tensor \(A\).

16.2.1. Algorithm#

Let us now explain step-by-step how VUMPS finds an approximate solution to the fixed-point equations in an iterative way.

We initialize the algorithm with the random guess \(\{A_L, A_R,C\}\), and chose a tolerance \(\eta\).
We first solve the first two eigenvalue equations

\[\begin{split} A_C = H_{A_C}(A_C), \\ C = H_C(C),\end{split}\]

using for example an Arnoldi algorithm with the previous approximations of \(A_C\) and \(C\) as initial guess. This yields two tensors \(\tilde A_C\) and \(\tilde C\).
From \(\tilde A_C\) and \(\tilde C\) we compute \(\tilde A_L\) and \(\tilde A_R\) that minimize following two-norms

\[\begin{split} \epsilon_L = \min_{A_L^\dagger A_L=1} ||\tilde A_C-\tilde A_L\tilde C||_2, \\ \epsilon_R = \min_{A_R A_R^\dagger=1} ||\tilde A_C-\tilde C\tilde A_R||_2,\end{split}\]

and thus approximately solve the last equation. Note that the minimum is taken over respectively left - and right isometric matrices. We comment below on the analytic soltuion of these equations and how this analytic solution can be approximated efficiently.
Update \(A_L\leftarrow\tilde A_L\), \(A_R\leftarrow\tilde A_R\) and \(C\leftarrow\tilde C\).
Evaluate \(\epsilon=\max(\epsilon_L,\epsilon_R)\) and repeat until \(\epsilon\) is below the tolerance \(\eta\).

Let us finally comment on solving the minimization problem to approximate \(\tilde A_{L/R}\).

A beautiful result in linear algebra states that the minimum is exactly given by \(\tilde A_L=U_lV_l^\dagger\) where \(U_l\) and \(V_l\) are the isometries arising from the singular value decomposition of \(\tilde A_C\tilde C^\dagger=U_l\Sigma_lV_l^\dagger\), and similarly \(\tilde A_R=U_rV_r^\dagger\), where \(\tilde C^\dagger\tilde A_C=U_r\Sigma_rV_r^\dagger\). Even though this approach will work well for the first iteration steps, this might not be the best solution close to convergence. When approaching the exact solution \(A^s_C=A^s_LC=CA^s_R\) the singular values in \(\Sigma_{l/r}\) become really small so that in finite precision arithmetic the singular vector in the isometries \(U_{l/r}\) and \(V_{l/r}\) are poor approximations of the exact singular vectors. A robust and close to optimal solution turns out to be

\[ \tilde A_L = U^l_{A_C}(U^l_C)^\dagger,\qquad \tilde A_R = (U^r_C)^\dagger U^r_{A_C},\]

where the \(U\)’s are the unitaries appearing in the polar decomposition of

\[\begin{split} \tilde A_C = U^l_{A_C}P^l_{A_C},\qquad \tilde C = U^l_CP^l_C,\\ \tilde A^r_C = P^r_{A_C}U^r_{A_C},\qquad \tilde C = P^r_CU^r_C.\end{split}\]

16.2.2. Example#

Let us demonstrate the algorithm using MPSKit by estimating the ground state energy density of the spin 1 XXX model. The VUMPS algorithm is called in the same way as we called DMRG. We initialize a random initial MPS with bond dimension 12 and physical dimension 3 (because the spin 1 representation of SU(2) is \(2\cdot1+1=3\)-dimensional). Obviously we don’t have to specify a system size because we work directly in the thermodynamic limit.

H = heisenberg_XYZ()

Ψ = InfiniteMPS(ℂ^3, ℂ^D)
algorithm = VUMPS()
Ψ₀, envs = find_groundstate(Ψ, H, algorithm);

┌ Info: VUMPS iteration:
│   iter = 1
│   ϵ = 0.5610100646483375
│   λ = -0.0865104407127983 - 1.2027136760735278e-17im
└   Δt = 0.084018603
┌ Info: VUMPS iteration:
│   iter = 2
│   ϵ = 0.4310799261082684
│   λ = -0.7618901087124069 - 6.792980042905938e-17im
└   Δt = 0.01130146
┌ Info: VUMPS iteration:
│   iter = 3
│   ϵ = 0.14979807151611316
│   λ = -1.3466352449583159 + 4.608048670280587e-19im
└   Δt = 0.034563615
┌ Info: VUMPS iteration:
│   iter = 4
│   ϵ = 0.015646810678441914
│   λ = -1.4006901307171526 + 1.3434458815938065e-16im
└   Δt = 0.008281546
┌ Info: VUMPS iteration:
│   iter = 5
│   ϵ = 0.005604286150776656
│   λ = -1.4013108403346757 + 9.647759176561666e-17im
└   Δt = 0.010333781
┌ Info: VUMPS iteration:
│   iter = 6
│   ϵ = 0.0023036744034752757
│   λ = -1.4013693555925446 + 4.854413600486681e-17im
└   Δt = 0.010374276
┌ Info: VUMPS iteration:
│   iter = 7
│   ϵ = 0.0008337434819378151
│   λ = -1.4013791484785758 - 1.2122138124051279e-16im
└   Δt = 0.010182469
┌ Info: VUMPS iteration:
│   iter = 8
│   ϵ = 0.0002936596335017501
│   λ = -1.4013804527428924 + 9.092701101546693e-17im
└   Δt = 0.009179364

┌ Info: VUMPS iteration:
│   iter = 9
│   ϵ = 0.00010509743552438332
│   λ = -1.4013806183201412 + 1.0246997139160656e-16im
└   Δt = 0.010018064
┌ Info: VUMPS iteration:
│   iter = 10
│   ϵ = 3.7494334247535774e-5
│   λ = -1.40138064015559 + 3.939797737721078e-18im
└   Δt = 0.009206416
┌ Info: VUMPS iteration:
│   iter = 11
│   ϵ = 1.367659602945338e-5
│   λ = -1.4013806430554205 + 7.858504907951533e-17im
└   Δt = 0.009297735
┌ Info: VUMPS iteration:
│   iter = 12
│   ϵ = 4.96721271892144e-6
│   λ = -1.401380643454414 - 9.294966409855741e-17im
└   Δt = 0.020732121
┌ Info: VUMPS iteration:
│   iter = 13
│   ϵ = 1.8420587042124429e-6
│   λ = -1.401380643509248 - 5.460206790974683e-17im
└   Δt = 0.006727026
┌ Info: VUMPS iteration:
│   iter = 14
│   ϵ = 6.785086789276065e-7
│   λ = -1.4013806435169813 - 6.438660477620306e-17im
└   Δt = 0.00660538
┌ Info: VUMPS iteration:
│   iter = 15
│   ϵ = 2.547473039126722e-7
│   λ = -1.4013806435180676 - 5.113601617677797e-18im
└   Δt = 0.006534318
┌ Info: VUMPS iteration:
│   iter = 16
│   ϵ = 9.485713314017868e-8
│   λ = -1.4013806435182241 - 6.761588005036087e-17im
└   Δt = 0.00655174
┌ Info: VUMPS iteration:
│   iter = 17
│   ϵ = 3.5938667493256724e-8
│   λ = -1.4013806435182465 - 1.3114288608388948e-16im
└   Δt = 0.006598416
┌ Info: VUMPS iteration:
│   iter = 18
│   ϵ = 1.3496801361477126e-8
│   λ = -1.401380643518248 + 6.122891167156048e-18im
└   Δt = 0.006548384
┌ Info: VUMPS iteration:
│   iter = 19
│   ϵ = 5.148257018855053e-9
│   λ = -1.4013806435182479 + 4.55259286979783e-17im
└   Δt = 0.006686911
┌ Info: VUMPS iteration:
│   iter = 20
│   ϵ = 1.9469060433781618e-9
│   λ = -1.4013806435182494 - 4.860066918499831e-18im
└   Δt = 0.006603276

┌ Info: VUMPS iteration:
│   iter = 21
│   ϵ = 7.464848416391231e-10
│   λ = -1.401380643518249 + 8.888490607361698e-17im
└   Δt = 0.007566857
┌ Info: VUMPS iteration:
│   iter = 22
│   ϵ = 2.839357240304302e-10
│   λ = -1.4013806435182492 + 2.907032789921015e-17im
└   Δt = 0.006776227
┌ Info: VUMPS iteration:
│   iter = 23
│   ϵ = 1.0930867955227511e-10
│   λ = -1.4013806435182508 - 4.557461146198056e-17im
└   Δt = 0.01718122
┌ Info: VUMPS iteration:
│   iter = 24
│   ϵ = 4.178118511402224e-11
│   λ = -1.4013806435182505 + 2.76571897107553e-17im
└   Δt = 0.00674009
┌ Info: VUMPS iteration:
│   iter = 25
│   ϵ = 1.6136658374819503e-11
│   λ = -1.4013806435182483 - 1.2684218222261536e-16im
└   Δt = 0.006618524
┌ Info: VUMPS iteration:
│   iter = 26
│   ϵ = 6.193587315436871e-12
│   λ = -1.4013806435182492 + 1.2530217727534505e-16im
└   Δt = 0.006731504
┌ Info: VUMPS iteration:
│   iter = 27
│   ϵ = 2.3984731713030978e-12
│   λ = -1.4013806435182485 + 1.8700481222026167e-17im
└   Δt = 0.006587777
┌ Info: VUMPS iteration:
│   iter = 28
│   ϵ = 9.235210216959459e-13
│   λ = -1.4013806435182485 + 4.674011821362125e-17im
└   Δt = 0.005651065
┌ Info: VUMPS summary:
│   ϵ = 9.235210216959459e-13
│   λ = -1.4013806435182485 + 4.674011821362125e-17im
└   Δt = 1.15866895

It takes about 30 iterations and a second or two to reach convergence. Let us gauge how well the ground state energy density was approximated by calling

expectation_value(Ψ₀, H)

1-element PeriodicArray{ComplexF64, 1}:
 -1.4013806435182485 - 1.6964694180142878e-17im

The value we obtain here is to be compared with the quasi-exact value -1.401 484 038 971 2(2) obtained in [Haegeman et al., 2011]. As you can see, even with such a small bond dimension we can easily approximate the ground state energy up to 3 decimals.