Comparing the ground state energies obtained by density matrix renormalization group, exact diagonalization, and an SDP hierarchy

We would like to compare the ground state energy of the following spinless fermionic system [1]:

$H_{\mathrm{free}}=\sum_{<rs>}\left[c_{r}^{\dagger} c_{s}+c_{s}^{\dagger} c_{r}-\gamma(c_{r}^{\dagger} c_{s}^{\dagger}+c_{s}c_{r} )\right]-2\lambda\sum_{r}c_{r}^{\dagger}c_{r},$

where $<rs>$ goes through nearest neighbour pairs in a two-dimensional lattice. The fermionic operators are subject to the following constraints:

$\{c_{r}, c_{s}^{\dagger}\}=\delta_{rs}I_{r}$

$\{c_r^\dagger, c_s^\dagger\}=0,$

$\{c_{r}, c_{s}\}=0.$

Our primary goal is to benchmark the SDP hierarchy of Reference [2]. The baseline methods are density matrix renormalization group (DMRG) and exact diagonalization (ED), both of which are included in Algorithms and Libraries for Physics Simulations (ALPS, [3]). The range of predefined Hamiltonians is limited, so we simplify the equation by setting $\gamma=0$.

Prerequisites

To run this notebook, ALPS, Sympy, Scipy, and SDPA must be installed. A recent version of Ncpol2sdpa is also necessary.

Calculating the ground state energy with DMRG and ED

DMRG and ED are included in ALPS. To start the calculations, we need to import the Python interface:


In [1]:
import pyalps

For now, we are only interested in relatively small systems, we will try lattice sizes between $2\times 2$ and $5\times 5$. With this, we set the parameters for DMRG and ED:


In [2]:
lattice_range = [2, 3, 4, 5]
parms = [{ 
  'LATTICE'        : "open square lattice",  # Set up the lattice
  'MODEL'          : "spinless fermions",    # Select the model 
  'L'              : L,                      # Lattice dimension
  't'              : -1 ,                    # This and the following
  'mu'             : 2,                      # are parameters to the
  'U'              : 0 ,                     # Hamiltonian.
  'V'              : 0,
  'Nmax'           : 2 ,                     # These parameters are
  'SWEEPS'         : 20,                      # specific to the DMRG
  'MAXSTATES'      : 300,                    # solver.
  'NUMBER_EIGENVALUES' : 1,          
  'MEASURE_ENERGY' : 1
} for L in lattice_range ]

We will need a helper function to extract the ground state energy from the solutions:


In [3]:
def extract_ground_state_energies(data):
    E0 = []
    for Lsets in data:
        allE = []
        for q in pyalps.flatten(Lsets):
            allE.append(q.y[0])
        E0.append(allE[0])
    return sorted(E0, reverse=True)

We invoke the solvers and extract the ground state energies from the solutions. First we use exact diagonalization, which, unfortunately does not scale beyond a lattice size of $4\times 4$.


In [4]:
prefix_sparse = 'comparison_sparse'
input_file_sparse = pyalps.writeInputFiles(prefix_sparse, parms[:-1])

res = pyalps.runApplication('sparsediag', input_file_sparse)
sparsediag_data = pyalps.loadEigenstateMeasurements(
                     pyalps.getResultFiles(prefix=prefix_sparse)) 

sparsediag_ground_state_energy = extract_ground_state_energies(sparsediag_data)
sparsediag_ground_state_energy.append(0)


sparsediag comparison_sparse.in.xml

DMRG scales to all the lattice sizes we want:


In [5]:
prefix_dmrg = 'comparison_dmrg'
input_file_dmrg = pyalps.writeInputFiles(prefix_dmrg, parms)
res = pyalps.runApplication('dmrg',input_file_dmrg)
dmrg_data = pyalps.loadEigenstateMeasurements(
                  pyalps.getResultFiles(prefix=prefix_dmrg)) 
dmrg_ground_state_energy = extract_ground_state_energies(dmrg_data)


dmrg comparison_dmrg.in.xml

Calculating the ground state energy with SDP

The ground state energy problem can be rephrased as a polynomial optimiziation problem of noncommuting variables. We use Ncpol2sdpa to translate this optimization problem to a sparse SDP relaxation [4]. The relaxation is solved with SDPA, a high-performance SDP solver that deals with sparse problems efficiently [5]. First we need to import a few more functions:


In [6]:
from sympy.physics.quantum.dagger import Dagger
from ncpol2sdpa import SdpRelaxation, generate_operators, \
                       fermionic_constraints, get_neighbors

We set the additional parameters for this formulation, including the order of the relaxation:


In [7]:
level = 1
gam, lam = 0, 1

Then we iterate over the lattice range, defining a new Hamiltonian and new constraints in each step:


In [8]:
sdp_ground_state_energy = []
for lattice_dimension in lattice_range:
    n_vars = lattice_dimension * lattice_dimension
    C = generate_operators('C%s' % (lattice_dimension), n_vars)
    
    hamiltonian = 0
    for r in range(n_vars):
        hamiltonian -= 2*lam*Dagger(C[r])*C[r]
        for s in get_neighbors(r, lattice_dimension):
            hamiltonian += Dagger(C[r])*C[s] + Dagger(C[s])*C[r]
            hamiltonian -= gam*(Dagger(C[r])*Dagger(C[s]) + C[s]*C[r])
    
    substitutions = fermionic_constraints(C)
        
    sdpRelaxation = SdpRelaxation(C)
    sdpRelaxation.get_relaxation(level, objective=hamiltonian, substitutions=substitutions)
    sdpRelaxation.solve()
    sdp_ground_state_energy.append(sdpRelaxation.primal)

Comparison

The level-one relaxation matches the ground state energy given by DMRG and ED.


In [9]:
data = [dmrg_ground_state_energy,\
        sparsediag_ground_state_energy,\
        sdp_ground_state_energy]
labels = ["DMRG", "ED", "SDP"]
print ("{:>4} {:>9} {:>10} {:>10} {:>10}").format("", *lattice_range)
for label, row in zip(labels, data):
    print ("{:>4} {:>7.6f} {:>7.6f} {:>7.6f} {:>7.6f}").format(label, *row)


             2          3          4          5
DMRG -8.000000 -18.828427 -33.708204 -52.928203
  ED -8.000000 -18.828427 -33.708204 0.000000
 SDP -8.000000 -18.828427 -33.708204 -52.928207

References

[1] Corboz, P.; Evenbly, G.; Verstraete, F. & Vidal, G. Simulation of interacting fermions with entanglement renormalization. Physics Review A, 2010, 81, pp. 010303.

[2] Pironio, S.; Navascués, M. & Acín, A. Convergent relaxations of polynomial optimization problems with noncommuting variables. SIAM Journal on Optimization, 2010, 20, pp. 2157-2180.

[3] Bauer, B.; Carr, L.; Evertz, H.; Feiguin, A.; Freire, J.; Fuchs, S.; Gamper, L.; Gukelberger, J.; Gull, E.; Guertler, S. & others. The ALPS project release 2.0: Open source software for strongly correlated systems. Journal of Statistical Mechanics: Theory and Experiment, IOP Publishing, 2011, 2011, P05001.

[4] Wittek, P. Ncpol2sdpa -- Sparse Semidefinite Programming Relaxations for Polynomial Optimization Problems of Noncommuting Variables. arXiv:1308.6029, 2013.

[5] Yamashita, M.; Fujisawa, K. & Kojima, M. Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0). Optimization Methods and Software, 2003, 18, 491-505.