The problem of determining the ground state of a quantum many-body system has wide-ranging applications in physics, chemistry, and materials science. This problem can be viewed as the problem of determining the lowest eigenvalue of a Hermitian operator on a Hilbert space whose dimension grows exponentially with the size of the system or the number of particles, such as electrons in the case of electronic structure. Here we highlight two relevant categories of approaches to taming this curse of dimensionality.
The first category is that of semidefinite relaxations, which rephrase the aforementioned energy minimization problem as an optimization problem in terms of a reduced set of physical observables, almost always a semidefinite program (SDP). In principle these observables satisfy representability constraints, i.e., constraints that ensure that they can be recovered from a bona fide quantum many-body state. However, only a subset of representability constraints can be efficiently enforced, yielding tractable optimization problems that provide lower bounds on the ground-state energy. Such approaches include the 2-RDM theories [28,26,4,27,42,24,1,29,7], as well as methods that may be classified as quantum marginal relaxations such as [23,30,9,16].
Meanwhile, quantum embedding theories take the perspective of dividing a system into local clusters, small enough to be treated with a highly accurate or exact method referred to as the 'high-level' method. Local problems are then stitched together via a reduced set of global quantities or a less-accurate 'low-level' method that operates on the global scale, and the local and global perspectives are constrained to be compatible via some self-consistency condition. Such approaches include dynamical mean-field theory (DMFT) [13,22] and density matrix embedding theory (DMET) [20,21], as well as variants such as the energy-weighted DMET (EwDMET) [10,11] which in a certain sense interpolates between DMFT and DMET [35].
Recently, variational embedding [25] was introduced as a semidefinite relaxation which is also a quantum embedding method. Like other relaxations such as [23,30,9,16], the key optimization variables are quantum marginals for local clusters, but variational embedding additionally includes global constraints which tighten the relaxation and accommodate the treatment of, e.g., long-range interactions.
As an SDP, variational embedding can be solved with black-box methods, as is done in [25], but scalability calls for a solver that is specially adapted to the problem. In this work, we introduce a scalable solver for this SDP which takes advantage of the embedding interpretation of the approach. The aforementioned global constraint is dualized to reduce the problem to a simpler relaxation (similar to those of [23,30,9,16]) in which the global constraints have been exchanged for effective contributions to relevant Hamiltonian operators at the local scale. This problem can then be solved in a fashion in which effective problems for the key variables (the two-cluster marginals) are completely decoupled and can be solved in parallel, with dual variable optimization enforcing the self-consistency of these problems. Furthermore, translation invariance of a lattice system can be used to significantly speed up the running time.
We now briefly mention some further details of our method for solving the simplified problem yielded by partial dualization. Our approach is similar to that of augmented Lagrangian methods (ALM) [17,31,32,37], in that it is derived from an augmented Lagrangian function. However, the primal update for the ALM suggested by our problem is too expensive to solve exactly in each iteration. Inexact ALMs [33] provide one option for circumventing this difficulty, but due to the separable structure of our objective (after partial dualization of the global constraint), it is natural for us to pursue an ADMM-type approach for this subproblem. The alternating direction method of multipliers (ADMM) was introduced by [15,12], and see [40] for recent developments. This approach allows us to solve several decoupled, physically local effective sub-problems completely in parallel, preserving the flavor of a quantum embedding theory.
We note that [41,36] have developed variants of ALM and ADMM, respectively, to solve large-scale SDPs. One key difference of our work is that we dualize a partial set of semidefinite constraints, then we apply ADMM to solve the partially dualized problem, which has convenient revealed structure. On the other hand, [41] and [36] directly work with the full primal and dual problems respectively.
In Section 2 we provide relevant background on the ground-state eigenvalue problem, examples of interest, and the two-marginal relaxation for variational embedding introduced in [25]. (In Appendix A, additional background is provided for the context of fermionic systems.) In Section 3, we describe our optimization approach to this problem, which is an SDP. The section begins with an idealized scheme of projected gradient ascent on the dual variable to the aforementioned global constraint. In order to implement such a scheme, it is necessary to solve an effective problem in terms of the primal variables. In Section 3.1, we introduce an ADMM-type approach to this problem, and in Section 3.2 we integrate this approach with dual ascent to define our practical scheme. In Section 3.3 we explain how translation-invariance can be exploited, and in Section 3.4 we include a detailed discussion of the computational scaling. Finally in Section 4 we present numerical experiments on several model systems of quantum spins and fermions.
In this section we review the formulation of variational embedding for quantum spins, following [25]. In Appendix A, we review the case of fermions (also following [25]), which requires a bit more care but nonetheless yields a semidefinite program of identical form after suitable manipulations.
We consider a model with sites, indexed = 1, … , , each endowed with a classical local state space (which shall be discrete). These in turn yields local quantum state spaces
= ℂ | | . The global quantum state space (i.e., the space of wavefunctions) is then given by ∶= ⨂ =1 ≃ ℂ || , where ∶= ∏ =1 is the global classical state space, so quantum states (wavefunctions) correspond to complex-valued functions on the classical state space. Let (resp. ) denote Hermitian operators → (resp. → ), and let ̂ (resp. ̂ ) denote the corresponding operators → obtained by tensoring (resp.
) by the identity operator on all sites ≠ (resp. ≠ , ). We consider pairwise Hamiltonian ̂ ∶ → of the form
and our interest is in determining the ground-state energy, i.e., the smallest eigenvalue, of ̂ . We denote this eigenvalue by 0 , which is defined variationally by
Now we review some examples of interest. First consider the case of quantum spin-
systems, i.e., the case = {-1, 1}. To construct operators on , one first starts with the Pauli matrices
which (together with the identity 2 ) form a basis for the real vector space of Hermitian operators on ℂ 2 . We let ∕ ∕ denote the operator → obtained by tensoring ∕ ∕ on the -th site with 2 on all other sites. Then in terms of these operators we can define the transverse-field Ising (TFI) Hamiltonian and the anti-ferromagnetic Heisenberg (AFH) Hamiltonian by
2)
ℎ is a scalar parameter and summation over ∼ indicates summation over pairs of indices that are adjacent within some graph defined on the index set {1, … , }, often a rectangular lattice in some dimension. These problems have been considered as prototypical quantum many-body problems, e.g., in [5], as well as models for the study of quantum phase transitions, as in [34].
In [25], the optimization problem (2.1) is reformulated as an optimization over the density operator , which (for nondegenerate ground states) corresponds at optimality to Φ 0 Φ * 0 where Φ 0 is the ground state eigenvector, i.e., the optimizer of (2.1). This optimization is in turn relaxed as a computationally tractable optimization over the quantum two-marginals
< , which are defined as partial traces of , analogous to classical marginals. Indeed, recall [25] that for any subset ⊂ {1, … , }, the partial trace = Tr [ ] may be defined as the unique operator on ⨂ ∈ such that Tr[ ] = Tr[ ̂ ] for all operators on ⨂ ∈ (lifted to operators ̂ on by tensoring with the identity on ). In particular, is an operator on ⊗ , and moreover it is positive semidefinite with unit trace (following from the same properties for ).
Then the two-marginal relaxation of [25] reads in terms of the two marginals (and the analogously-defined one-marginals, which can be obtained from the two-marginals by further partial trace) as the following semidefinite program, whose optimal value we denote by
[{ }, { } ≤ ] ⪰ 0.
(2.8)
Here = [{ }, { } ≤ ] is an operator defined linearly in terms of the one-and two-marginals, subordinate to the specification of an arbitrary collection { , ∶ = 1, … , } of linear operators → at each site = 1, … , . In specific, is specified blockwise, with blocks
The choice of operators only matters up to span{ , ∶ = 1, … , }, and in our numerical experiments we shall consider the complete operator collection spanning all linear maps → . The last constraint (2.8) is called the global semidefinite constraint.
To motivate the relaxation (2.4) further, we examine the problem of finding the lowest energy state of a classical energy function of a pairwise form. (This can be viewed as a special case of the more general quantum ground-state problem by taking the to be diagonal operators.) More concretely, for ∈ , = 1, … , , define an energy function
and observe that the minimizer of can be determined via the linear program arg min
where Π() is the space of probability measures on . Indeed, the optimizer is a -function supported on the minimizer of (provided that it is unique). Exploiting the pairwise structure of , we have
where the two-marginal variables { } < are constrained to be jointly representable, i.e., to be derivable as the two-marginals of a high-dimensional measure . Enforcing this constraint demands exponential complexity, so various convex relaxation approaches have been proposed, where only certain necessary conditions for the { } < are kept; see, for instance, [39] for a review. In particular, the analogous convex relaxation to (2.4) is minimize
where is an all-one vector of appropriate size. Here (2.13) and (2.15) are standard constraints for any discrete probability distributions, and (2.14) are constraints that enforces 'local consistency' [39] among the { } < . The global semidefinite constraint (2.16) is discussed in [19] and [6] in the contexts of multi-marginal optimal transport and energy minimization, respectively.
Given a quantum spin model as above and a decomposition of the sites {1, … , } as a disjoint union
′ to be the classical state space for the ′ -th cluster. One see that any Hamiltonian that is pairwise with respect to sites is pairwise with respect to clusters, so by viewing our clusters as sites and applying the above formalism, we obtain a tighter relaxation [25].
In [25] it was shown that (2.4) admits the minimax formalization (obtained via dualization of the global semidefinite constraint (2.8))
(2.17)
where
5)-(2.7) . (2.18) Here denote the blocks of , and the 'effective' Hamiltonian terms [ ] and [ ] are defined linearly in terms of via [ ] ∶= -∑ , =1 ( ) † , , , [ ] ∶= -∑ =1 ∑ =1 ( ) † , ⊗ , + h.c. , (2.19) where 'h.c.' denotes the Hermitian conjugate term. This partial dual formulation can be obtained from (2.4) by exchanging the global semidefinite constraint (2.8) for an extra term
in the Lagrangian, where ⪰ 0 is a dual variable with respect to which the Lagrangian is to be maximized. Then (2.19) is recovered by breaking this additional term into a blockwise sum, collecting terms, and minimizing over the primal variables, subject to the remaining constraints. Since and are dual variables [25] we have that
where { }, { } < are the minimizers of the infimum in (2.18).
In order to solve the two-marginal relaxation (2.4), our point of departure is the partial dual formulation (2.17). For simplicity we consider the case in which the Hamiltonian and all variables are purely real, and we let Π ⪰0 be defined by
for a symmetric matrix , i.e. the Euclidean projection (in Frobenius norm) of a symmetric matrix onto the set of real symmetric positive semidefinite matrices (equivalent to setting all negative eigenvalues of the argument to zero). As an idealized scheme, we can imagine performing projected gradient ascent on (2.17), which is implemented by Algorithm 1. In practice, the step size > 0 in Algorithm 1 should be chosen sufficiently small to guarantee convergence, but not much smaller.
Algorithm 1 Exact projected gradient ascent.
Require: > 0, ⪰ 0 1: while not converged do 2:
Set { }, { } < to be the minimizer in (2.18), holding fixed 3:
The focus of this section is in the development of algorithm for step 2 for a general Hamiltonian, and we also study the case in the presence of translational invariance. In practice, we will not fully converge a solution to step 2 of Algorithm 1, resulting in an inexact projected gradient ascent scheme. However, in order to motivate our practical scheme, we will first discuss how to solve step 2 exactly for fixed .
We rephrase step 2 as the following optimization problem:
and is fixed for the duration of this subsection. Moreover, for simplicity we have assumed that ∶= | | is constant, and 1 , 2 are defined to be the linear operators Tr {2} and Tr {1} , respectively. Hence 1 , 2 can be realized as sparse matrices of size 2 × 4 We then formulate an equivalent optimization problem via the introduction of dummy variables ̃ and the inclusion of augmented Lagrangian terms in the objective:
Tr
where , > 0 are constant parameters, and Λ , Λ (1) , Λ (2) ∈ End( ⊗ ) are the dual variables for the associated constraints in (3.4) and (3.5). Let [{ }, { , ̃ } < ] denote the objective function (3.3), yielding the Lagrangian Algorithm 2 Augmented Lagrangian method for (3.2).
Require: , > 0, { ′ }, { ′ , Λ , Λ (1) , Λ (2) } < 1: while not converged do 2:
for each pair < do 4:
:
end for 8: end while
In practice, it is difficult to solve step 2 of Algorithm 2 exactly. Therefore, instead of optimizing { }, { , ̃ } < jointly, we consider an ADMM-type [3] substitute, namely Algorithm 3. Notice that in step 2 of Algorithm 3, the are all determined independently Algorithm 3 Pseudo-code for ADMM-type method for (3.2).
Require:
} < 1: while not converged do 2: 1) , Λ (2) } <
for each pair < do 5:
end for 9: end while as the solutions of decoupled optimization problems
After suitable manipulation of the objective (neglecting constant terms), we obtain
which can be exactly optimized via the update
Meanwhile, in step 3, the and the ̃ can all be updated via decoupled optimization problems. In particular, we find that
Finally we turn to the update. Collecting the relevant terms we have that
Observe that the objective may be rewritten as
which we must minimize subject to Tr = 1. This is simply a constrained least squares problem, the solution of which yields the update ← 1 ( -1)
where is a Lagrange multiplier chosen to satisfy the constraint.
Then via (3.7), (3.8), and (3.9), we can rewrite Algorithm 3 concretely as the equivalent Algorithm 4. Observe that all of the for-loops in Algorithm 4 can be run in parallel.
Require: , > 0, { ′ , }, { ′ , ̃ , Λ , Λ (1) , Λ (2) } < 1: while not converged do 2:
for each pair < do 3:
end for 5:
for each pair < do 6:
end for 8:
for each do 9:
end for 13: for each pair < do 14:
:
end for 18: end while
Now we return to the full Algorithm 1 where we optimize via an idealized gradient ascent. Instead of exactly implementing step 2 of Algorithm 1, we replace it with a single iteration of Algorithm 4, yielding our practical approach Algorithm 5 for solving the two-marginal relaxation (2.4).
Note that within Algorithm 5 it is possible to include an additional loop allowing for several repetitions of the block of primal updates (steps 2-18) before the block of dual updates (steps [19][20][21][22][23][24]. If this additional loop were iterated to convergence, we would recover an exact augmented Lagrangian method. However, in practice we find that a single iteration of such a loop is the most effective choice (noting that the cost of the primal updates dominates that of the dual updates), and there is no practical motivation for further complicating Algorithm 5 as written. A practical comparison is provided in Appendix B. Furthermore, inspired by [14], it is natural to consider a dual step size ∈ (1, (1 + √ 5)∕2) appearing as a prefactor in front of the dual variable increments of steps 19-24. This possibility is explored in Appendix C, but due to a lack of demonstrated advantage is left out of our main numerical experiments for simplicity.
We pause now to provide some perspectives on the convergence of Algorithm 5. First, it is worth noting by inspection that our iterative algorithm formally retains the correct fixed point at optimality. However, a proof of convergence is beyond the scope of this paper. Indeed, certain features of the current algorithm differ importantly from the conventional augmented Lagrangian method for SDP. In our outer loop, we perform a dual update via gradient ascent on the dual variable for the global PSD constraint. The inner loop provides a gradient for the outer loop, though it is in practice solved only approximately. For the outer loop, we cannot readily leverage the nice properties of proximal point methods as in the case of ALM [33] to establish convergence. Furthermore, as the inner loop is only roughly solved, the resulting gradient is inexact, which introduces further difficulty.
One of the most expensive step in Algorithm 5 is step 24, where a projection to the positive semidefinite cone is required. We now discuss how translation-invariance of a lattice system can be exploited algorithmically in the solution of the two-marginal relaxation Algorithm 5 Practical ADMM / projected dual gradient ascent method for (2.4).
Require: , , > 0, ⪰ 0, { }, { ̃ , Λ , Λ (1) , Λ (2) } < 1: while not converged do 2:
for each do 3:
end for 5:
end for 8:
for each pair < do 9:
10: end for 11:
for each pair < do 12:
13: end for 14:
for each do 15:
end for 19:
for each pair < do 20:
22:
end for 24:
25: end while (2.4). It is convenient in this section to adopt zero-indexing for the site, i.e., to index the sites by the multi-index = ( 1 , … , ) where is the lattice dimension and = 0, … , -1. Then we assume translation-invariance in that ̂ = ̂ for all and ̂ = ̂ 0, -for all < . The symmetries = and = 0, -are in turn guaranteed to be satisfied by some optimizer of (2.4) [25].
Then we can implement Algorithm 5 (whose iterations preserve this symmetry) without any reference to variables besides and the , . The main challenge is the implementation of step 24 of Algorithm 5. Given and the , , we can only compute the top row ′ , of ′ ∶= + [{ }, { } < ]. However, by translation invariance, the rest of ′ is determined by the property that ′ , = ′ , -. (In the case = 1, ′ is a block-circulant matrix, though a more general term is lacking for the case of arbitrary .)
Via translation invariance, ′ is block-diagonalized by the block-discrete Fourier transform. More precisely, one can write ′ as
, where indicates the appropriate -dimensional discrete Fourier transform matrix, the identity matrix of the appropriate block size, ⊗ the Kronecker product, and ̂ ′ a block-diagonal matrix. Hence to compute the projection Π ⪰0 [ ′ ] we first compute the diagonal block ̂ ′ of ̂ ′ as
Note that the blocks ̂ ′ (concatenated into a block row) can be viewed as the entrywise discrete Fourier transform of the first block row of ′ , hence can be computed simultaneously via FFT. Then project ∶= Π ⪰0 ̂ ′ for all , and set
The final pseudocode for the translation-invariant setting is given in Algorithm 6.
In Algorithm 5, observe that the for-loops run over ( -1) pairs < , and the scaling bottleneck among these loops is the projection ̃ ← Π ⪰0 --1 Λ occurring in step 12. Since this step requires full diagonalization of a matrix of size 2 × 2 , for which the cost is ( 6 ). Meanwhile, suppose for simplicity that = 2 , corresponding to the complete choice of operator collection , ∶ = for each site. Then the size of and is 2 × 2 . Hence step 24, which involves a complete diagonalization of a matrix of this size, costs ( 6 3 ), dominating the ( 6 2 ) cost of the for-loops. If our sites are in fact supersites, each formed from clusters of sites in an underlying spin-
model, then = 2 . Therefore the scaling is (2 6 3 ) per iteration, where is the cluster size and is the number of clusters.
Meanwhile, in the translation-invariant setting of Algorithm 6, the for-loops run only over ( ) sites, so-neglecting the update for -the asymptotic cost per iteration is ( 6 ). Meanwhile, the construction of the ̂ ′ in terms of the 0, can be achieved in time ( 4 log ) via FFT. (Note that we simply treat the lattice dimension as constant.) The cost of each projection of step Algorithm 6 Translation-invariant ADMM / projected dual gradient ascent method for (2.4).
, , , { ̃ , , Λ , , Λ (1) , , Λ (2) , } ≠ 1: while not converged do 2:
3:
5: end for 6:
for each ≠ do 7:
end for 9:
for each ≠ do 10:
end for 12:
, -′ 13:
for each ≠ do 16:
Λ (1) , ← Λ (1) , + ( 1 [ , ] -)
18:
, ← Λ (2) , + ( 2
[ , ] -) 19: end for 20: ′ , ← , + , 21: for each ≠ do 22: ′ , ← , + , , 23: end for 24:
for each do 25:
← Π ⪰0 ̂ ′
end for 28:
for each do
, ←
1 √ ∑ exp 2 ⋅ 30: end for 31: end while 26 is ( 6 ) via full diagonalization, and forming the , in terms of the also costs ( 4 log ) in total via FFT. Hence the cost of updating (i.e., steps 20 through 29) is ( 4 log + 6 ). Hence the total cost per iteration of Algorithm 6, under the assumption that the sites are supersites each composed of spin-1 2 sites, is (2 4 log + 2 6
), where is the cluster size and is the number of clusters. The exponential scaling in the cluster size is unavoidable in our formulation due to the exact treatment of reduced density operators on the clusters (i.e., the cluster marginals). In this work we consider clusters of size no larger than size = 4. Future work will investigate the possibility of treating larger clusters by introducing further relaxation and/or compression of the optimization variables to avoid exponential scaling in the cluster size. As a final comment, observe that every for-loop in Algorithms 5 and 6 can be run fully in parallel.
The numerical experiments were implemented in MATLAB following Algorithm 6. (We shall consider only translation-invariant Hamiltonians.) We present results for the transverse-field Ising (TFI) model (2.2), the anti-ferromagnetic Heisenberg (AFH) model (2.3), the spinless fermion (SF) model (A.1), and the long-range spinless fermion (LRSF) model (A.2). Note that Algorithm 6 can be applied in the fermionic case mutatis mutandi to the problem (A.3). Throughout we fix the value of the algorithmic parameters to be = = 10, = 2 (i.e., we do not tune them specifically to different problems). The dual variables {Λ , Λ (1) , Λ (2) } ≠ are all initialized to be zero, and the primal density operator variables { }, { ̃ } ≠ are all initialized as multiples of the identity with unit trace.
is initialized as the identity. We run Algorithm 6 for 10,000 iterations. (The convergence behavior will be studied in detail below.)
First we consider the TFI model on a periodic 20 × 1 lattice, which is small enough to be solved by exact diagonalization of (2.2). We benchmark the per-site energy error of the two-marginal relaxation with clusters of size 1 × 1, 2 × 1, and 4 × 1. (In these cases the semidefinite matrix variables are each of size 4, 16, and 256, respectively; refer to Section 3.4 for further discussion of scaling.). The results are shown in Fig. 4.1. Observe that the approximations yield lower bounds for the energy as the theory requires, and these lower bounds become tighter as the cluster size is increased. In the same figure we also consider the TFI model on a periodic 4 × 4 lattice. Here we benchmark the energy error of the two-marginal relaxation with clusters of size 1 × 1, 2 × 1, and 2 × 2.
We perform completely analogous experiments for the AFH model with similar conclusions. The results are shown in Tables 1 and 2.
We also benchmark the SF and LRSF models on a 20 × 1 periodic lattice, with results pictured in Fig. 4.2. Note that the fermionic relaxations are exact for = 0, as guaranteed in [25]. Next we consider the TFI model on a periodic 100 × 1 lattice for ℎ = 0.5, 1, 1.5. This problem is too large to solve by exact diagonalization. We report the relaxation energy for several cluster sizes in Table 3. We track convergence behavior of Algorithm 6 on these same problems. We use two different quantities to track convergence. The first is the per-site primal objective (i.e., energy) change between subsequent iterations. The second is the per-cluster feasibility error for the equality constraints, defined as Table 3 Relaxation energy per site for the TFI model with ℎ = 0.5, 1, 1.5 on a 100 × 1 periodic lattice for various cluster sizes.
We plot these quantities as functions of the iteration count in Fig. 4.3. We also recommend that these quantities be used to define a practical stopping criterion for the optimizer, i.e., to stop when both converge within some specified numerical tolerance. In Appendix D we present additional experiments demonstrating the satisfaction of the global semidefinite constraints.
It is possible that tuning the parameters , , to a specific problem and specific choice of clusters could yield smoother convergence profiles. However, even using our fixed choice for all problems, we achieve convergence of the per-site energy within 10 -6 (which is dominated by the relaxation error itself) in a number of iterations that does not seem to grow with the cluster size. Though the asymptotic convergence of ADMM-type methods can be slow, for practical purposes we only need the convergence accuracy to be dominated by the relaxation error of the relaxed problem. If convergence to very high accuracy is required, alternative approaches such as inexact ALM [33] might be more appropriate.
Then we fix clusters of size 2 × 1 and vary the system size of the TFI model with ℎ = 0.5, 1, 1.5, to investigate the effect of system size on convergence. The results are shown in Fig. 4.4. The system size does not appear to have any obvious detrimental effect on the convergence rate.
To conclude, we report the relaxation energies obtained from these last experiments in Table 4. We observe that the relaxation energy approaches a limiting value in this thermodynamic limit, i.e., limit of infinite volume.
Relaxation energy per site for TFI models of several sizes with ℎ = 0.5, 1, 1.