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$.
To run this notebook, ALPS, Sympy, Scipy, and SDPA must be installed. A recent version of Ncpol2sdpa is also necessary.
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)
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)
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)
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)
[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.