Lower Bounds for Ground States of Condensed Matter Systems

Standard variational methods tend to obtain upper bounds on the ground state energy of quantum many-body systems. Here we study a complementary method that determines lower bounds on the ground state energy in a systematic fashion, scales polynomially in the system size and gives direct access to correlation functions. This is achieved by relaxing the positivity constraint on the density matrix and replacing it by positivity constraints on moment matrices, thus yielding a semi-definite programme. Further, the number of free parameters in the optimization problem can be reduced dramatically under the assumption of translational invariance. A novel numerical approach, principally a combination of a projected gradient algorithm with Dykstra's algorithm, for solving the optimization problem in a memory-efficient manner is presented and a proof of convergence for this iterative method is given. Numerical experiments that determine lower bounds on the ground state energies for the Ising and Heisenberg Hamiltonians confirm that the approach can be applied to large systems, especially under the assumption of translational invariance.


I. BACKGROUND AND MOTIVATION
The determination of the ground state properties of strongly interacting quantum systems is of considerable interest. As only a very few models are exactly solvable, numerical approximation methods are of significant importance. A key challenge in this respect is the 'curse of dimensionality', namely the exponential growth of the Hilbert space dimension with the number of sites. Ingenious methods such as the density matrix renormalization group (DMRG) approach [1,2] can overcome this problem as they amount to formulating a variational problem over a well-chosen class of quantum states that can be parametrized efficiently, i.e. with a polynomial number of parameters [3]. As the variation is over a subclass of all possible quantum states, variational methods of this type provide upper bounds on the ground state energy. Assessing the quality of the so obtained upper bounds on the ground state energy is then difficult. Hence there is considerable interest in the development of efficient methods for the determination of lower bounds on the ground state energy. One straightforward approach in this direction consists in the determination of Anderson bounds [4], i.e. the decomposition of a Hamiltonian H = k H k and the subsequent determination of the ground state energy of each H k whose sum then yields a lower bound on the ground state energy of H. While useful, this approach does not scale well with the system size as it remains exponential in the size of the support of the H k . The application of DMRG to the H k may overcome the scaling issue but amounts to finding an 'upper bound to a lower bound' and does not resolve the principal issue with DMRG.
As early as 1940, however, it was suggested to avoid the use of wave functions altogether and rather consider reduced density matrices together with conditions that these arise from a valid quantum state [5] -the so-called N-representability problem [6,7]. Varying over the reduced density matrices that are compatible with a valid physical state then allows for the determination of the ground state energy. The exact ground state energy can be obtained if all conditions for the N-representability problem are known. If only partial conditions are specified, then one minimizes over too large a set of reduced density matrices and one obtains a lower bound on the ground state energy. The analytical treatment of this problem is highly challenging and at the time numerical solutions were hard to come by. Further progress was made by Mazziotti [8,9], who pointed out that the problem can be formulated as a semidefinite programme and use can be made of the guaranteed and certified convergence of the primal-dual interior-point algorithm for semi-definite programming [10,11]. A general mathematical framework for this type of optimization problems and applications was given in [12]. Unfortunately, the primal-dual interior-point algorithm and related methods do come with a considerable overhead in computation time and, crucially, memory that is making it hard to scale to very large systems. Lately, new algorithms to address the optimization problems of finding lower bounds on ground state energies have been developed [13,14] and considerable improvements towards the interior-point algorithms have been demonstrated.
In this work, we are exploring further the idea of relaxations of the ground state energy finding problem that allow us to determine lower bounds on the ground state energy in an efficient manner. To this end we formulate the problem, introduce a number of simplifications including the explicit treatment of translation invariance and symmetries and point out additional potential routes for numerical simplifications. In addition, we present a novel algorithm, principally a combination of a projected gradient algorithm with Dykstra's algorithm [15][16][17], for solving the resulting optimization problem in a memory-efficient manner and prove its convergence. Then, in order to demonstrate the capabilities of this approach, we present numerical experiments that determine lower bounds on the ground state energies for two one-dimensional (1D) lattice models, namely the Ising and Heisenberg Hamiltonian, which confirm that the approach can be applied for large systems especially when implementing translational invariance. Note that although we restrict our numerical analysis to 1D lattice models it has been demonstrated that the variational determination of ground state energies considering only constraints on reduced density matrices applies to higher spatial dimensions [18] and to non-periodic systems such as molecules in quantum chemistry [13] as well. We complete this work with a discussion and outlook.

II. BASICS
We begin by explaining the general ideas underlying the method for the determination of lower bounds on ground state energies following [8] and then specialize to the case of fermionic systems, for which later sections will provide numerical examples. Analogous ideas can also be implemented for bosonic systems or general spin-systems.

A. Generalities
The determination of the ground state energy E 0 of a many-body system composed of N particles and described by the Hamilton operator H can be formulated as the minimization problem This approach is pushed to its limits rather rapidly as it uses the positivity constraint on ρ, which grows dramatically with the system size, e.g. for particles with spin-1/2 or spinless fermions ρ is a 2 N × 2 N matrix. Hence, this problem cannot be solved efficiently due to the exponentially increasing parameter space. One technique to deal with this difficulty is to parametrize a specific family of candidates |φ ( α) for the ground state by a set of parameters α. One then minimizes the expectation value of the Hamiltonian with respect to this family of states. As this family of states does not encompass the entire state space, this method will find a state |φ 0 which is optimal in {|φ ( α) } α and which yields an upper bound on the true ground state energy. This is the route taken by a variety of variational methods including, for example, DMRG [1][2][3].
To overcome the scaling issues and to provide lower rather than upper bounds on the ground state energy, we need to relax the positivity constraint ρ ≥ 0, replacing it by a weaker set of constraints which scale only polynomially in the system's size [8,12]. This enlarges the class of states over which one varies the energy and hence yields lower bounds on the ground state energy. To this end we use the fact that for all operators O. Demanding only for some operator O does not, in general, imply ρ ≥ 0. For some subset of operators {O k } k , we can then define O = k α k O k and the requirement that O † O ≥ 0 for all choices of α k is then equivalent to the positivity of the matrix X with components X kl = O † k O l , imposing a particular N-representability condition. Hence, replacing in the original optimization problem eq. (1) the condition ρ ≥ 0 with the weaker condition X ≥ 0 and omitting the normalization tr[ρ] = 1, we find that The present formulation still contains the density operator ρ explicitly in the computation of tr Note, furthermore, that the entries of X will generally not be independent but suffer some linear equality constraints, for example due to (anti)-commutation relations between operators. Denoting by C the set of all Hermitian matrices that satisfy X ≥ 0 and any other linear constraints, the energy minimization problem can be written as [8,9] which is evidently a semi-definite programme [11] and efficiently solvable as long as the set {O k } k contains a number of operators that scale polynomially in system size. The decomposition of H into H kl , and hence the choice of the set {O k } k , is not unique and the optimal choice is not obvious. From a practical point of view however, the systematic choice presented in the remainder of the paper leading to the so-called p-positivity conditions [8,9] does perhaps provide the most straightforward formulation.

B. Bosons and Fermions
For concreteness, let us now consider the class of bosonic or fermionic operators defined by where δ k , k , α kl , β kl , γ kl ∈ C are arbitrary complex variables. The N-representability conditions arising from this class of operators are called 1-positivity and 2-positivity conditions, which is a terminology used to reveal the powers of the creation and annihilation operators a † i and a i in O. In general, allowing polynomials of order p in the second-quantized operators leads to restrictions called p-positivity conditions [8,9]. For the particular choice of operators of the order of 2, requiring the positivity of O † O is equivalent to where the N × N matrices comprising the second moments are defined as follows: T k,l = a † k a l , S k,l = a k a † l , U k,l = a k a l (8) and for the third and fourth moments, we find the N 2 × N and N 2 × N 2 matrices A kl,m = a † k a † l a m , M kl,mn = a † k a † l a n a m , B kl,m = a † k a l a m , G kl,mn = a † k a l a n a m , C kl,m = a k a l a m , H kl,mn = a k a l a n a m , In addition to the positivity constraint on the 2N + 3N 2 × 2N + 3N 2 Hermitian matrix X ∈ H, where H denotes the vector space of all Hermitian matrices with appropriate dimension, there are linear constraints relating its entries. In fact, the fermionic anti-commutator relations or the bosonic commutator relations force some entries of X to be equal, to differ by sign or to be a linear function of others. Now, let C be the set of all positive semi-definite Hermitian matrices of the form eq. (7) and satisfying the additional requirement that the entries of X obey the constraints given by the commutation or anti-commutation relations. Further, let E (X) : H → R be a function expressing the expectation value of a Hamiltonian with up to 4-mode terms, i.e. E (X) = H . Then, the original optimization problem eq. (1) can be replaced by where the relaxed optimization problem scales polynomially in the system size N . Note that it is straightforward to extend this formulation to higher moments, e.g. 3-positivity conditions. In the following, we are going to present scenarios and strategies where the number of free variables can be decreased considerably.

C. Consequences of Superselection Rules and Particle Number Conservation for Fermions
For the remainder we will consider fermionic systems. In that case, superselection rules require that the Hamiltonian commutes with the parity operator given by . As a consequence, only Hamiltonians comprising terms with an even number of creation and annihilation operators are physically allowed. This in turn ensures that only expectation values of products of an even number of creation and annihilation operators are non-vanishing. Hence, we find that we can restrict X to take the form and to satisfy the constraints dictated by the fermionic anti-commutator relations. See appendix B for a list of constraints on the entries for a formulation for up to fourth moments in the fermionic case.
If the Hamiltonian conserves the particle number, a further simplification ensues. In this case, each term in the Hamiltonian will consist of an equal number of annihilation and creation operators. Hence, the Hamiltonian is invariant under the transformation a k → a k e iφ and a † k → a † k e −iφ and we find that expectation values vanish unless they concern operators with the same number of annihilation and creation operators so that X can be assumed to take the form Note that the elements of the matrices T, S are coefficients of the one-electron reduced density matrix (1-RDM), while the matrices M, R, Q are representations of the two-electron reduced density matrix (2-RDM) [9]. Restricting matrices of the form eq. (14) to be positive semi-definite is equivalent with the positivity of the individual matrices T, S and M, R, Q. This constraint on the matrices containing the second moments defines the 1-positivity conditions, while the positivity of the matrices containing the fourth moments generates the 2-positivity conditions [8,9]. For concreteness we list the matrices comprising the constraints on a system composed of two spinless fermions and the additional assumption of particle number conservation in appendix B.

D. Exploiting Translational Invariance
So far we considered optimization problems for up to 2-positivity conditions which reduced the number of free variables to scale as N 4 . Apart from superselection rules and the particle number conservation, no specific assumption on the structure of the underlying Hamiltonian has been exploited such that in principle this procedure is applicable to a wide range of different systems, e.g. quantum chemistry calculations [13] or higher spacial dimensions [18]. Since our intention is to apply this method to large lattice models, we now discuss a further simplification of the system to decrease the number of parameters, namely the translational invariance. Translation invariance with periodic boundary conditions of the physical system reduces the number of free variables considerably. In fact for the second moments, i.e. for the matrices T, S and U , this can be exploited immediately as these matrices then obey T i,j = t i−j , etc. Consequently, for translational invariant systems with periodic boundary conditions, the matrices comprising the second moments become circulant matrices [19] and can therefore be diagonalized by a discrete Fourier transform for k, l = 1, . . . , N . Reducing the number of variables for the fourth moments is not so transparent as for the second moments, but manageable. To do so, note that the ordering of the entries of the matrices M, G, H, R, I and Q is not fixed but can be chosen arbitrarily. Remember that the indices of these matrices are of the form [(k, l) , (m, n)] for k, l, m, n = 1, . . . , N . Rearranging the two-digit indices labeling the rows and columns in the following way: the matrices comprising the fourth moments decompose in a block structure where each block can be diagonalized by the unitary matrix V k,l given above. For example, we find for the matrix M where U = diag (V, V, . . . , V ) with V given by eq. (15) and with diagonal matrices D i,j , i, j = 1, . . . , N . Hence, for translational invariant systems the matrices contained in set C will have the following form: where the matrices denoting the second moments in the upper left block are diagonal matrices and the matrices considering the fourth moments have a block diagonal structure as described above. These considerations reduce the complexity of the constraints considerably and hence make the numerical study feasible for large systems.

E. Varying System Size for Higher Moments
As the treatment of the fourth-order correlations is computationally costly, it may in some cases be advantageous to consider fourth-order constraints on a smaller lattice, while the computationally cheaper second-order constraints are treated on the full lattice. In particular, while T, U and S are N × N matrices, M, G, H, R, I and Q are N 2 × N 2 matrices leading to the unfavourable scaling. To overcome this issue, it might therefore be beneficial to restrict higher moments to smaller subsystems; that is, for the second moments one would consider system size N 2 , while for fourth moments one would consider N 4 where N 2 N 4 . Reducing the support of the matrices M, G, H, R, I and Q, on the one hand, will lead to disproportionate savings in memory and computation time and, on the other hand, will certainly lead to a decrease of the lower bounds.

III. THE PROJECTED GRADIENT ALGORITHM
So far we have described how to relax the original optimization problem eq. (1) and discussed how to include symmetries to reduce the number of free variables. Now we are going to describe an iterative algorithm for addressing this minimization problem. Recall the minimization problem we have to face, that is, where again the function E (X) : H → R expresses H and the set C consists of all positive semi-definite Hermitian matrices X fulfilling the constraints dictated by the (anti)-commutator relations and potential other symmetries of the system. Note that the function E (X) is indeed linear in the elements of matrix X and hence could be written as where G is a Hermitian matrix and c ∈ R. Furthermore, the set C is the intersection of the convex cone consisting of positive semi-definite matrices and the affine set defined by the linear constraints on the entries. Hence C is convex, too [11]. It is well known that this optimization problem is in fact a semi-definite programme and as such can be solved efficiently. Standard primal-dual interior point methods for semi-definite programming come with a considerable overhead in memory and time which makes their application to large systems challenging. Recently, two algorithms addressing the problem of ground state energy estimation were developed which scale better than the conventional methods [13,14]. Here, we continue the research for efficient algorithms and formulate an alternative numerical method whose efficacy we demonstrate in later sections. To address the issues described above, we suggest to solve the minimization problem by a projected gradient algorithm [20]. Since this algorithm relies heavily on projections onto convex sets, we need to define the latter. Let H be a vector space and A ∈ H be a convex set. The projection of X ∈ H onto A is defined by where ||·|| F : H → R denotes the Frobenius or Hilbert-Schmidt norm. As the name suggests, the projection determines the nearest point in the convex set A with respect to a given norm. With this, we claim that the solution of problem eq. (20) could be found by iterating 1. Illustration of the gradient projection algorithm (23). The convex sets A and B have a nonempty intersection, which is shown shaded. The straight parallel lines clarify the level sets of the function E (X) such that E * < E < E , where E * is the minimum value the function attains on the intersection. The projections are represented by arrows to demonstrate their directions. Note that these projections are computed using Dykstra's algorithm, which itself uses the projections onto A and B. The algorithm is initialized with the matrix X0.
where the algorithm is initialized by a Hermitian matrix X 0 ∈ H and α > 0 is an arbitrary real and positive number. Furthermore, P C : H → C denotes the projection onto set C, i.e. the intersection of the positive semi-definite Hermitian matrices and the set given by the linear constraints. The scheme of this algorithm is visualized in figure 1. We find that the following lemma holds.
Lemma 1 Let E : H → R be a linear function of the form E (X) = G, X + c where G ∈ H n×n and c ∈ R is a constant. Assume that the function E (X) is lower bounded on the closed convex set C, i.e. there is anX ∈ C such that E X ≤ E (X) for all X ∈ C. Let E k = E (X k ) be the sequence determined by the iterations of the algorithm eq. (23).
A detailed proof of this lemma is presented in appendix A. So far we have shifted the problem of minimizing the function E (X) with respect to the convex set C to the problem of computing the projection onto C. As pointed out above, this set is the intersection of two convex sets. This complicates the computation of the projection because in principle one has to solve the optimization problem eq. (22) in each step of the algorithm. The most obvious method at hand to overcome this difficulty is to use Dykstra's algorithm [17,21]. Dykstra's algorithm is an elaborate technique to compute the projection of X onto the intersection of several convex sets by means of modified iterative projections onto the individual sets. To apply this procedure to the computation onto C we need to determine the projections onto the set of positive semi-definite Hermitian matrices and onto the affine set given by the constraints imposed by the canonical (anti)-commutation relations and other symmetries. The former can be implemented straightforwardly. For the latter we were able to determine formulae for the submatrices in eq. (13) and thus gain a method to compute P C : H → C efficiently. Note that since C is convex, P C (X) is unique for all X ∈ H [17]. Further, the proposed algorithm finds a lower bound on the ground state energy by only considering boundary points of the feasible set and hence is slightly related to the algorithm proposed in [14].
There are two parameters defining the accuracy of the estimated lower bound. In principle the algorithm (23) converges for perfect projections P C (X) to a matrix X ∈ C minimizing the function E (X) in the limit k → ∞, where k is the iteration number. As for all iterative algorithms we need to introduce stopping criteria, on the one hand, for the projected gradient algorithm and, on the other hand, for Dykstra's algorithm. For the latter we use the robust stopping criterion introduced in [21] and determined by a small number τ Dykstra ∈ R. For the former one can exploit the fact that the dual problem of the optimization problem eq. (20), which we identify as the primal problem, can be solved using the projected gradient algorithm, too. Note that the primal problem can be written as where we identify the matrices G and X with the vectors g and x such that E (X) − c = G, X = g † x. The Hermitian matrices F 0 , {F i } i comprise the linear constraints dictated by the fermionic anti-commutator relations such that the matrix F 0 + i x i F i is an element of the affine set for all x i ∈ R. The dual problem of eq. (24) is then for the dual variable Z ∈ H [10] maximize − F 0 , Z subject to g i = F i , Z Z ≥ 0. (25) Note that the constraints of the dual problem force the solution Z * to lie both in the positive semi-definite cone and in the affine set given by the linear equations. The intersection of these two convex sets is again convex and the dual problem becomes where D denotes the intersection of the positive semi-definite cone with the set given by the linear constraints g i = F i , Z . Finally, lemma 1 shows that the projected gradient algorithm can be applied to solve the dual problem as well.
The benefit of the formulation as a primal-dual problem is that the solution of the dual problem is a lower bound on the solution of the primal problem. Hence, solving the primal and dual problems simultaneously using the projected gradient algorithm provides us with an ingenious certificate for the accuracy of the solution. Consequently, one way to introduce a stopping criterion for the projected gradient algorithm is to initiate a small quantity τ primal−dual and stop the computations when where H (Z) = − F 0 , Z + c is basically the dual function up to a constant. Note that eq. (27) is equivalent to the duality gap of the primal and dual problem and that for strictly feasible primal problems, i.e. for optimization problems of the form eq. (24) for which there exists an x such that F 0 + i x i F i > 0, the duality gap is zero if evaluated at optimal points [10]. For the optimization problem eq. (20), one can easily find a point in the feasible set which is strictly feasible. Consequently, we have E (X * ) − H (Z * ) = 0, where X * and Z * are solutions of the primal and dual problems, respectively, and hence eq. (27) suits perfectly as a stopping criterion. Due to the numerical cost to solve the primal and dual problems simultaneously, we decide on using a different stopping criterion for the projected gradient algorithm. We stop when the improvement per iteration step, i.e.
falls below a predetermined small number τ E ∈ R. We implemented both routines for systems considering only second moments and numerical experiments suggest that for a sufficiently small number τ E the latter stopping criterion yields results that are as good as in the primal-dual formalism. Another degree of freedom in the formulation of the algorithm is the step length α. Numerical examples propose that for the first step in the direction of the negative gradient α should be chosen to be large. For the remaining steps the computational cost could be reduced by choosing a small but appropriate value of α [22].

IV. APPLICATIONS
In the remainder of this work, we provide a number of numerical examples that demonstrate that the proposed method can be applied to rather large systems and is capable of extracting accurate information concerning the ground state energy as well as the correlation functions. The purpose of what follows is not to improve lower bounds on unknown ground state energies, but rather to analyse the behaviour of the proposed algorithm for finite but large lattice models.

A. Second Moments and the Ising Model
We begin our exposition of the numerical application of the algorithm with a simple, exactly solvable model for which the relaxation of the energy minimization problem has been analysed in great detail exploiting different numerical methods in [26]. The Hamiltonian we are considering as a first application of the projected gradient algorithm is given by where the operators describe spinless fermions. The ground state of this model coincides with that of the Ising model in a transverse magnetic field with periodic boundary conditions when j c < 0 and N is even. This can be seen from the connection of eq. (30) with eq. (29) via a Jordan-Wigner transformation [23,24]. To analyse a translational invariant system with periodic boundary conditions, we simply change the signs of the second and the fourth terms in the Hamiltonian eq. (29). Now we require only constraints that are due to second moments, i.e. we require that and that the entries fulfil the constrains dictated by the anti-commutator relations. Solving the relaxed minimization problem eq. (20) by applying the projected gradient algorithm and comparison of the so-obtained lower bound to the exact ground state energy [23,24] yields the results presented in table I. Note that the ground state energy for the system given by the Hamiltonian eq. (29) is computed by the formula (32) with parameter l k = 2k + 1, while one chooses l k = 2k for the energy of a translational invariant system with periodic boundary conditions.  Hamiltonian with parameters jc = −1 and h = 1/2 for N modes. The Hamiltonian is identical to eq. (29) but the signs of the second and fourth terms are changed to make the Hamiltonian translational invariant. The relative deviation ∆E is of the order of the parameters of the stopping criteria τ Dykstra = τE = 10 −8 , which limits the accuracy of this method. We initialized the algorithm with α = N for the first step and α N for the remaining steps until convergence.
In fact, it can be shown directly that the ground state energy of the Ising model is exactly reproduced with the present method as the Ising model can be diagonalized by an orthogonal mode transformation mapping second moments to second moments, this making the transformed optimization problem straightforward to solve [25]. Thus the so obtained lower bounds are, in this case, superior to other approaches such as Anderson lower bounds. In particular, solving the optimization problem for the non-translational invariant Hamiltonian eq. (29) for even N and negative j c yields the exact ground state energies for the translational invariant Ising Hamiltonian eq. (30). Even for this Hamiltonian one is able to solve the optimization problem for large particle numbers since one is only considering the 1-positivity conditions which scale linearly in the system size.

B. Fourth Moments and the Heisenberg model
While the Ising model is quasi-free it is of considerable interest to study interacting models, such as the Heisenberg Hamiltonian, that contain, for example, fourth-order fermionic operators (see [27] for the study of another interacting fermionic system where a different numerical approach is employed). In the following, we consider the performance of the proposed algorithm in this situation. The Heisenberg Hamiltonian can be transformed to spinless fermions by the Jordan-Wigner transformation. In this picture, the Hamiltonian reads Now we are in a position to compare the exact ground state energies of the Heisenberg model with the lower bounds obtained by applying the suggested algorithm. Indeed, for the Heisenberg Hamiltonian with parameter J = 1/2 the lower bounds on the ground state energy lie very close to the actual energies.   The range of applicability of the projected gradient algorithm can be extended considerably by analysing a system that is translationally invariant. To this end, we make the Hamiltonian in eq. (34) translationally invariant with periodic boundary conditions. The results for the ground state energy of this Hamiltonian are presented in the lefthand side of figure 2. As a reference, we depict the exact ground state energies for this model for up to N = 22. Additionally, the right-hand side of figure 2 presents the nearest-neighbour and next-nearest-neighbour correlations determined by the gradient projection algorithm.
To illustrate the rate of convergence of the algorithm for fourth moments, we present the improvement per iteration of the function E (X k ) = G, X k + c for the translational invariant Jordan-Wigner transformed Heisenberg Hamiltonian describing N = 50 sites in figure 3. It is noticeable that the algorithm converges to a good estimation in only a few steps and that the most time-consuming factor is given by the projection onto the set C.
Note that in all calculations we initialized the algorithm with the zero matrix since we wanted to test the algorithm without imposing prior knowledge. Generally, for translational invariant systems it is advantageous to complete the circulant submatrices obtained from a run with site number N 1 to circulant matrices describing N 2 sites and use these matrices as an initial matrix for N 2 > N 1 .
Furthermore, since the matrices comprising the fourth moments in the translational invariant case decompose into blocks of diagonal matrices, it is possible to parallelize the projection onto the positive semi-definite matrices. In fact, this step in the computation is the most expensive since one has to diagonalize the matrices in each step. Parallelizing this computation might speed up the calculations. So far, no effort has been made to do so since our intention is rather to demonstrate that this technique yields good results. for α = x, y, z, which is possible since the system is rotationally invariant. For the next-nearest-neighbour correlations, we proceeded the same way. . Energies E k for each iteration step for the Hamiltonian eq. (34) extended to a translational invariant system for N = 50 particles illustrating the good behaviour of the gradient projection algorithm with respect to the iteration step k. Note that the data is plotted double logarithmically.

V. CONCLUSIONS
In this work, we have presented a method for the determination of lower bounds on the ground state energy of quantum many-body systems. To this end, we followed a large body of work to relax the general ground state energy finding problem to obtain a semi-definite programme that involves a number of variables that scale polynomially with the system size. To overcome unfavourable scaling in time and memory of standard primal-dual interior point solvers, we combined a projected gradient algorithm with Dykstra's algorithm to obtain access to larger system sizes and proved its convergence. For quasi-free systems the method yields the exact ground state energies, while for general interacting models it provides very good lower bounds and direct access to correlation functions between sites. Different properties of the underlying system such as higher spatial dimensions [18], bosonic degrees of freedom [28], symmetries [29] and translational invariance can be included straightforwardly and the relative weight of higher order correlation functions can be adjusted to optimize for computational performance.
We hope that further investigations will help to optimize this approach and thus elucidate its power and limitations.

VI. ACKNOWLEDGEMENTS
We acknowledge support from the Bundesministerium für Bildung und Forschung (FK 01BQ1012), the EU Integrating Project Q-ESSENCE and the Alexander von Humboldt Foundation. Interesting discussions with K. Audenaert, J. Cai, M. Cramer, D.A. Mazziotti and M. Navascues are acknowledged.

Appendix A: Proof of Convergence
The analysis of lemma 1 requires some basic properties of the projection operators defined in eq. (22). The following results can be found, for example, in [17], theorems 4.1 and 5.5.
Lemma 2 Let P C : H → C be the projection onto the convex set C.
1. For all Z ∈ C and all X ∈ H we have P C (X) − X, Z − P C (X) ≥ 0.

For all
The following proof is motivated by [20]. Proof of lemma 1 Before we start with the main proof, we need to establish the relation for α ≥ 0 and where the sequence {X k } k is defined via i.e. the projected gradient algorithm. To verify relation (A1), note that with the first inequality of lemma 2, we find and relation (A1) is established. Note that with this expression the following holds for all α > 0: In a first step, we show that the sequence E (X k ) is monotonically decreasing. We find that where E (X) = G, X + c, which holds for all α ≥ 0. Hence the sequence E (X k ) is monotonically decreasing. Since we assume that the function E (X) is lower bounded on C, we then conclude that this sequence converges for all α ≥ 0, i.e. E (X k ) → E * for k → ∞. It remains to be shown that E * is indeed the minimum of the function in C.
In a second step, we show that Since E (X) is linear in X, we have and hence Since the sequence E (X k ) converges, i.e. E (X k ) → E * , we find that lim k→∞ G, X k − X k+1 = 0. (A10) Inequality (A1) then yields and hence Thus we have In a third step, we show that the elements X k are bounded. LetX ∈ C be such that E X ≤ E (X) for all X ∈ C. Then since G, P C (X k − αG) −X = G, X k+1 −X = E (X k+1 ) − E X . Note that with lemma 2 and the fact that E X ≤ E (X) for all X ∈ C, one can show that the last two terms are greater than or equal to zero. Therefore and consequently which leads to We conclude that the elements X k are bounded. As a consequence we know that, since C is closed, there is at least one converging subsequence X n k → X * for k → ∞, where X * is an accumulation point of the sequence {X k } k . In a last step, we show that for any accumulation point X * it holds that E (X * ) ≤ E (X) for all X ∈ C. With lemma 2, we find for all Z ∈ C X k+1 − X k + αG, Z − X k+1 ≥ 0, ⇔ α G, X k+1 − Z ≤ X k+1 − X k , Z − X k+1 .
(A18) Further, we have where we exploited the Cauchy-Schwarz inequality. Finally, we find that holds for all Z ∈ C. Now choose a converging subsequence X n k such that X n k → X * . We already know that G, X n k − X n k +1 → 0 due to equation (A10). Further, equation (A13) claims that ||X n k +1 − X n k || F → 0 for k → ∞. Additionally, since the sequence {X k } k is bounded, we know that forX ∈ C, where E X ≤ E (X) for all X ∈ C, where we exploited that the sequence ||X − X n k || F is upper bounded by ||X − X 0 || F , see equation (A17) , and where C (Z) is a constant which may depend on arbitrary Z ∈ C. Hence we find for all Z ∈ C and all accumulation points X * lim k→∞ G, X n k − Z ≤ 0, ⇒ G, X * − Z ≤ 0, ⇒ G, X * ≤ G, Z , ⇒ G, X * + c ≤ G, Z + c, ⇒ E (X * ) ≤ E (Z) .

(A22)
Since we know that the sequence E (X k ) converges, we can choose the subsequence defined by X n k → X * to see that E (X n k ) → E * = E (X * ). With the argumentation above, we know that E * ≤ E (Z) for all Z ∈ C and hence the algorithm (A2) finally ends up at a point minimizing the function E (X) in the convex set C.