Mitigating the Sign Problem Through Basis Rotations

TL;DR

This paper introduces a method to optimize the basis in quantum Monte Carlo simulations to mitigate the Fermion sign problem, significantly improving simulation efficiency for large two-dimensional Hubbard systems.

Contribution

The authors demonstrate how to use sign-free quantum Monte Carlo to optimize basis rotations, reducing the sign problem and accelerating simulations of large 2D Hubbard models.

Findings

Optimized basis rotations improve the average sign in Hubbard models.

Sign problem reduction leads to exponential speedup in simulations.

Method applicable across various U and density parameters.

Abstract

Quantum Monte Carlo simulations of quantum many body systems are plagued by the Fermion sign problem. The computational complexity of simulating Fermions scales exponentially in the projection time $\beta$ and system size. The sign problem is basis dependent and an improved basis, for fixed errors, lead to exponentially quicker simulations. We show how to use sign-free quantum Monte Carlo simulations to optimize over the choice of basis on large two-dimensional systems. We numerically illustrate these techniques decreasing the `badness' of the sign problem by optimizing over single-particle basis rotations on one and two-dimensional Hubbard systems. We find a generic rotation which improves the average sign of the Hubbard model for a wide range of $U$ and densities for $L \times 4$ systems. In one example improvement, the average sign (and hence simulation cost at fixed accuracy) for the $16\times 4$ Hubbard model at $U/t=4$ and $n=0.75$ increases by $\exp\left[8.64(6)\beta\right]$. For typical projection times of $\beta\gtrapprox 100$, this accelerates such simulation by many orders of magnitude.

Equations

$$⟨s⟩\equiv \frac{1}{N}\sum _i{s}_i$$

$$|H(R){|}_{ij}$$

$$=H(R{)}_{ij}\text{ otherwise }$$

$$H=-t\sum _{\sigma ⟨i,j⟩}{c}_{j\sigma }^{†}{c}_{i\sigma }+h.c.+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow }$$

$$\overline{E(R;\beta )}=\frac{⟨{\Psi }_{\text{init}}|\overline{H(R)}\exp [-\beta \overline{H(R)}]|{\Psi }_{\text{init}}⟩}{⟨{\Psi }_{\text{init}}|\exp [-\beta \overline{H(R)}]|{\Psi }_{\text{init}}⟩}$$

$$v_{i+1}=v_i+\alpha \gamma \frac{\left|\partial \overline{E}/\partial v_i\right|}{\partial \overline{E}/\partial v_i},$$

$$\overline{E_V(R)}\equiv ⟨{\Psi }_{\text{u(R)}}|\overline{H(R)}|{\Psi }_{\text{u(R)}}⟩$$

$$\begin{pmatrix}b_{(n,4)}\\ b_{(n,3)}\\ b_{(n,2)}\\ b_{(n,1)}\end{pmatrix}^{(\dagger)}=\dfrac{1}{\sqrt{2}}\begin{pmatrix}0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{5.50005pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{5.50005pt}{0.0pt}{$-$}}1\\ 1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{5.50005pt}{0.0pt}{$-$}}1&0\\ 0&1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{5.50005pt}{0.0pt}{$-$}}1\\ 1&0&1&0\end{pmatrix}\begin{pmatrix}c_{(n,4)}\\ c_{(n,3)}\\ c_{(n,2)}\\ c_{(n,1)}\end{pmatrix}^{(\dagger)}$$

$$\left(\begin{smallmatrix}b_{(2,4)}\\ \vdots\\ b_{(2,1)}\\ b_{(1,4)}\\ \vdots\\ b_{(1,1)}\end{smallmatrix}\right)^{(\dagger)}=\dfrac{1}{2}\left(\begin{smallmatrix}0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&1&0&1\\ 1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&1&0\\ 0&1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&1\\ 1&0&1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0\\ 0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1\\ 1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0\\ 0&1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1&0&1&0&\leavevmode\hphantom{0}\hbox to0.0pt{\hss\leavevmode\resizebox{3.85002pt}{0.0pt}{$-$}}1\\ 1&0&1&0&1&0&1&0\end{smallmatrix}\right)\left(\begin{smallmatrix}c_{(2,4)}\\ \vdots\\ c_{(2,1)}\\ c_{(1,4)}\\ \vdots\\ c_{(1,1)}\end{smallmatrix}\right)^{(\dagger)}.$$

$$U=\left(\begin{array}{ccc}-\frac{1}{\sqrt{2}}& 0& \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}& 0& \frac{1}{\sqrt{2}}\\ 0& 1& 0\end{array}\right),$$

$$U=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}-1& 0& 0& 1\\ 0& 1& 1& 0\\ 1& 0& 0& 1\\ 0& -1& 1& 0\end{array}\right).$$

$$U=\left(\begin{array}{cccccccc}-0.314786& 0.125218& -0.326724& 0.0175546& 0.465212& 0.366267& 0.415591& 0.504861\\ -0.0370841& -0.251754& -0.321497& -0.511366& 0.00856541& -0.497735& -0.322363& 0.467608\\ 0.390442& -0.48436& -0.326193& 0.0390196& -0.529225& 0.362701& 0.258697& 0.162699\\ -0.527768& -0.447829& -0.355797& 0.509121& 0.0331845& -0.0526765& -0.284827& -0.223855\\ 0.432137& -0.141862& -0.390688& -0.0378131& 0.500662& -0.290055& 0.312168& -0.454778\\ -0.110209& 0.243491& -0.398565& -0.493372& -0.047652& 0.494441& -0.344209& -0.40134\\ -0.258416& 0.472052& -0.373214& 0.0537856& -0.497251& -0.393736& 0.394274& -0.102314\\ 0.450426& 0.428415& -0.325755& 0.479303& 0.0461818& 0.0122356& -0.451591& 0.267417\end{array}\right).$$

$$b_{(i,j)}=\sum _k{U}_{jk}\frac{1}{\sqrt{2}}\left(c_{(i+1,k)}-s_{(i-1)}{c}_{(i-1,k)}\right)$$

$$U_{jk}=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}0& -1& 0& -1\\ 1& 0& -1& 0\\ 0& 1& 0& -1\\ 1& 0& 1& 0\end{array}\right)$$

Full text

The power of quantum computers and the difficulty of simulating quantum many-body ground states stem from similar sources: non-trivial entanglement and paths of interfering signs. Low entanglement states can be solved using the density matrix renormalization group (DMRG) Schollwöck (2011); White (1992). Paths of interfering signs are the source of the fermion sign problem in quantum Monte Carlo (QMC) simulations. Systems, such as bosonic Hamiltonians, which have only positive paths and hence no sign problem, can be simulated efficiently. For some classes of systems with ‘naive’ sign-problems, approaches have been found which also allow for efficient simulations. One common technique involves finding a sign-problem free basis in which to perform the simulation Nakamura (1998); Okunishi and Harada (2014); Honecker et al. (2016); Alet et al. (2016); Wessel et al. (2018). Examples where this has also been successful include Hubbard models at half-filling Wang et al. (2015); Li and Yao (2019) and employing a Majorana representation Li et al. (2015). Note that while there is always a sign-problem free basis (i.e. the eigenstate basis), it is often as (or more) difficult to find this basis as doing the simulation. Even in cases where the sign problem can’t be removed, the choice of basis can affect the ‘badness’ of the sign problem Kolodrubetz et al. (2013). Interestingly, there has not yet been a significant amount of work in automatically searching for a good basis to help mitigate or remove the sign problem. It is this problem that we approach in this Letter.

Our approach is as follows: given a quantification of the sign problem, search over a class of unitary rotations of the Hamiltonian minimizing the ‘badness’ of the sign problem. The key question, then, is the development of an algorithm to efficiently accomplish this. As a proof of principle we will consider hole-doped Hubbard models finding a rotation that decreases the ‘badness’ of the sign problem for $L\times 4$ lattices over a wide range of $U/t$.

I Quantification of Sign Problem

Let $H$ be the Hamiltonian represented in a standard basis (i.e. real space) and $H(R)\equiv RHR^{\dagger}$ be the Hamiltonian rotated by the unitary rotation $R$. Our choice for quantifying the sign problem of $H(R)$ has two goals: (1) The objective function should be proportional to the cost of QMC on $H(R)$. (2) We can measure (and hence optimize) this objective function efficiently on large bulk two-dimensional systems.

The natural objective function to maximize is the relative variance of the average sign

[TABLE]

where $s_{i}$ is the sign for each Monte Carlo sample $i$ and $N$ is the total number of Monte Carlo samples. While the efficiency of different observables vary, their values $\langle O\rangle=\sum_{i}s_{i}O_{i}/\langle s\rangle$ all involve the average sign. Here $O_{i}$ is the observable measured on Monte Carlo sample $i$. Given the variance of $\langle s\rangle$ is $O(1)$, minimizing the relative variance comes down to maximizing $\langle s\rangle.$ Because our interest will be in the limit of a large QMC projection time $\beta$ our focus will be in optimizing the component of $\langle s\rangle=A\overline{s}^{\beta}$ which causes $\langle s\rangle$ to decay exponentially with $\beta$ - i.e. $\overline{s}$.

Naively, though, this appears to be difficult to optimize as computing $\langle s\rangle$ itself has a sign problem. It turns out, though, that locally optimizing $\langle s(R)\rangle$ over a class of unitary rotations $R$ is significantly easier then computing it. In particular, $\overline{s(R)}=\exp(-\Delta E(R))$ Kolodrubetz et al. (2013); Spencer et al. (2012); Ceperley (1981) where $\Delta E(R)\equiv E(R)-\overline{E(R)}$ is the gap between the ground state $E(R)$ of the fermionic Hamiltonian $H(R)$ and the ground state $\overline{E(R)}$ of an effective bosonic version of the Hamiltonian $|H(R)|$ Loh et al. (1990); Hatano and Suzuki (1992); Troyer and Wiese (2005) defined as

[TABLE]

with $\overline{E(R)}=\langle\overline{\Psi}|\overline{H(R)}|\overline{\Psi}\rangle$

where $|\overline{\Psi}\rangle$ is the ground state of $\overline{H(R)}$. Note this definition ensures that the propagator used in QMC, $|G|\equiv 1-\tau|H|$ (where $\tau$ is a small constant) contains strictly positive matrix elements. While $E(R)$ is hard to compute (there is a sign problem in evaluating it), its value is independent of the choice of the unitary rotation - i.e. $E(R)=E$. On the other hand $\overline{E(R)}$ is sign-free and hence computable by Monte Carlo on large systems. Maximizing $\overline{s(R)}$ then comes down to maximizing $\overline{E(R)}$.

Note this metric is closely related to, but not identical to having a ground state with only positive amplitudes or a Hamiltonian with no positive off-diagonal matrix elements. For example, the Heisenberg model on a bipartite lattice has no sign problem (i.e. no exponentially large relative variance) in QMC but does have positive off-diagonal Hamiltonian matrix elements and a non-positive ground state (although both these can be removed by applying the Marshall sign rule Marshall and Peierls (1955)). $\Delta E(I)$ where $I$ is the identity rotation is correctly zero in this case.

II Unitary Rotation

Here we describe the parameterization of our unitaries. While a unitary rotation exists which removes the sign problem (i.e. the eigenstate basis), it is often harder to find then the ground state itself. Instead, we focus on unitary rotations which have simple representations. Two such classes are unitary rotations from shallow-depth quantum circuits and basis rotations on the set of single-particle orbitals used to represent the Hamiltonian: i.e. given a Hamiltonian written as $H=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}+\sum_{ijkl}V_{ijkl}c_{i}^{\dagger}c_{j}^{\dagger}c_{k}c_{l}$ we can write $b_{j}\equiv\sum_{k}U_{jk}c_{k}$ where $U$ is a $N\times N$ unitary matrix for a system of $N$ sites. $U$ can be parameterized in various ways including as the eigenvalues of an orthogonal matrix $H$ or as $U=e^{A}$ where $A$ is a skew-symmetric matrix, i.e. $A=-A^{T}$. In our numerical examples, we focus on single particle basis rotations although the approach we describe works also with quantum circuits.

III Optimization

Our goal now is to maximize $\overline{E(R)}$. In this section we demonstrate a series of techniques to optimize this quantity which we benchmark on single particle basis rotations of the Hubbard model under periodic boundary conditions,

[TABLE]

with $\langle i,j\rangle$ denoting nearest neighbors.

We optimize $\overline{E(R)}$ via 3 separate techniques. First, we optimize $\overline{E(R)}$ via exact diagonalization using finite differences. When diagonalization is too costly, we turn to the second technique, projector quantum Monte Carlo Trivedi and Ceperley (1990); Ceperley and Alder (1980) (PQMC). PQMC is a form of QMC which propagates walkers for a projection time $\beta$ giving the energy

[TABLE]

where walkers are sampled from an initial wave-function $|\Psi_{\textrm{init}}\rangle$. In PQMC, by taking large $\beta$, we can compute $\overline{E(R)}\equiv\lim_{\beta\rightarrow\infty}\overline{E(R;\beta)}$ in polynomial time; derivatives can be computed using finite differences giving an algorithm which scales linearly in the number of parameters and otherwise is similar in cost to a sign-free QMC simulation of the same rotated Hamiltonian (i.e. polynomial in system size).

Parameters are then updated using the optimization scheme described in Luo and Clark (2019) inspired by Lou and Sandvik (2007). For each of the unique entries $v_{i}$ of our unitary parameterization, the next parameter is determined by

[TABLE]

where $\alpha$ is a random number chosen between 0 and 1 and $\gamma$ controls the size of the random step.

Although the cost of the direct PQMC approach is polynomial it can still be expensive especially since this rotation can result in a quadratic number of terms per (matrix) row (which is a significant increase from the linear number of terms in the Hubbard model but standard for QMC on materials Booth et al. (2009)). To overcome this expense, we develop and benchmark a further approximation to our algorithm which significantly accelerates the optimization. Instead of directly optimizing $\overline{E(R)}$, we optimize

[TABLE]

where $\Psi_{\textrm{u(R)}}\equiv 1/2^{N}$ is the uniform wave-function in the basis $R$. This is a strict upper bound in energy to $\overline{H(R)}$ and is equivalent to $\overline{E(R;\beta=0)}$ in the PQMC when the initial wave-function is chosen as $|\Psi_{\textrm{u(R)}}\rangle$. This approximation turns out to be reasonable as the variational energy $\overline{E_{V}(R)}$ tracks the ground state energy $\overline{E(R)}$ and so pushes the parameters in the correct direction (see fig. 2). This approximation can be made even faster by choosing a relatively fixed few random states ($O(10)$) and using these states to optimize $\overline{E_{V}(R)}$. We further utilize the Jax library Bradbury et al. (2018) and implement the derivative $d\overline{E_{V}(R)}/dR$ using automatic differentiation. Because these derivatives (and the energy) are computed from only $O(10)$ Monte Carlo samples, the statistical error is large, but the correlations introduced by using the same configurations better controls the path generated by the gradients.

III.1 Results

III.1.1 1D

We perform initial tests optimizing one-dimensional Hubbard chains for $N=\{3,4,8\}$. We evaluate $\Delta E(R)$ via exact diagonalization and perform derivative free optimization using Nelder–Mead optimization Johnson ; Nelder and Mead (1965). In these runs, $R$ is parameterized by the eigenvectors of a real Hermitian matrix. Our optimization completely removes the sign-problem at $N=3$ at 2/3 doping and $N=4$ Hubbard model at half filling for all values of $U$ (see supplementary material for the specific matrices); a posteri, we discovered these rotations are known in the CTQMC literature Shinaoka et al. (2015). Interestingly at $N=4$ (see supplementary), our optimization procedure finds a basis which does not make all the off-diagonal signs of $|G|$ positive despite still removing the sign problem (although by hand $G$ can be made positive). We also include an example for $N=8$ where we improved but were unable to completely remove the sign problem (see fig. 1).

III.1.2 2D

Next we consider a variety of width-4 Hubbard models. First we implement projector optimization for the $n=0.5$, $U/t=1$, $4\times 4$ Hubbard model. We parameterize $R$ as the matrix exponential of a real skew-symmetric matrix and start from a random unitary which has a bad sign problem (worse than not rotating at all).

Fig. 2a shows the improvement in $\overline{E(R)}$ during optimization. By rounding the final unitary and removing noise, we obtain a sparse unitary rotation which rotates each site operator into linear combinations of 4-operator plaquettes (see inset of fig. 2a). The rotated operators don’t map to their counterpart in the original basis, i.e. operator $b_{(i,j)}$ doesn’t contain $c_{(i,j)}$. We find that $\Delta E(R)=0.137$ and achieve an increase in the average sign which goes as $\exp[6.27\beta]$ at large $\beta$.

We test our optimized rotation for a range of doping $n$ and interaction strengths $U$ in fig. 3. While only optimized for a single $U$ and $n$, we find this rotation mitigates the sign problem over the real-space basis for the $4\times 4$ Hubbard model for all $0.125<n<1.0$ when $U/t\leq 4$. The ability for a single optimized rotation to be applicable to a larger regime of phase space shows the versatility of our approach.

We then turn to width-4 cylinders with periodic boundary conditions. We perform optimization on the $16\times 4$ Hubbard model with $n=0.5$ and $U/t=2$, using automatic differentiation of $\overline{E_{V}(R)}$ and varying step sizes (see supplementary), shown in fig. 2b. By rounding and removing noise of the final unitary we obtain a sparse rotation which improves the sign by $\exp[12.86(5)\beta]$. We observe this rotation also works for a variety of doping and $U$ values, e.g. $n=0.5$ $U/t=6$ and $n=0.75$ $U/t=4$ have improvements of the average sign by $\exp[3.78(6)\beta]$ and $\exp[8.64(6)\beta]$ respectively.

Upon analyzing the unitary, we find a clear structure that can be used for any width-4 cylinder. The single particle orbitals in the bulk of the lattice become a linear combination over each column of width four, leading to 4 creation/annihilation operators (with implicit spin indices)

[TABLE]

At the left and right edges of the system, despite periodic boundaries, two consecutive columns are rotated as

[TABLE]

We illustrate this rotation in the inset of fig. 2b. Similar to the structure of the $4\times 4$ rotation, the new operators are linear combinations of nearby sites but lack the original site with respect to the old basis. However, the bulk is made of 2-site pairs of sites rather than 4 site plaquettes, and the edge plaquettes are not spaced apart in the $\hat{x}$ direction.

Testing this for a larger system, we observe an average sign reduction over not rotating on a $32\times 4$ cylinder at $n=0.5$ doping of $\exp[22.5(4)\beta]$ and $\exp[16.0(3)\beta]$ for $U/t=4$ and $U/t=6$ respectively. These energy differences $\overline{E(I)}-\overline{E(R)}$ per site at $U=4$ are increasing with respect to system size (from $L=8$ to $L=32$), suggesting in the thermodynamic limit that our rotation provides a better lower bound to the true fermionic energy than identity rotation (details can be found in the supplementary).

Finally, we can take the edge rotation found in eq. 8 and apply it both to the $4\times 4$ and $8\times 8$ lattices. For the $4\times 4$ case, we find that the rotation improves the sign more not than rotating, but not as well as our optimized $4\times 4$ unitary, i.e. $\exp[5.40\beta]$ compared to $\exp[6.27\beta]$ for $n=1$ $U/t=1$ respectively. When using this rotation to tile an $8\times 8$ unitary, we find that we improve the sign only for $U/t\lesssim 1$, e.g. at $U/t=1$ we find an improvement of $\exp[-3.0(1)\beta]$ and $\exp[-3.9(1)\beta]$ for $n=0.5$ and $n=0.875$ respectively.

In addition to the above results, we include in the supplementary an optimization of a 6x6 lattice to show that the original PQMC optimization of $\overline{E(R)}$ scales beyond what is accessible by exact diagonalization. We also include a study on $4\times 4$ systems of how the rotation affects the average sign in PQMC with respect to $\beta$ and the benefit of introducing annihilation Booth et al. (2009); Kalos and Pederiva (2000) on top of the rotation.

III.2 Conclusion

We have introduced a new approach for finding a basis which reduces the sign problem in bulk systems of arbitrary dimensions. For fixed accuracy and projection time, this improves the calculation speed of a ground-state QMC by many orders of magnitude; a similar approach will also improve finite-temperature calculations. Our method uses QMC to maximize the energy of the “bosonic” version of the rotated Hamiltonian $|H(R)|$; this maximization corresponds to maximizing the average sign of a QMC calculation in that rotated basis. In many ways, our sign-problem mitigation approach is similar in flavor to variational Monte CarloUmrigar et al. (2007); Neuscamman et al. (2012); Kochkov and Clark (2018); Otis and Neuscamman (2019); Zaletel and Pollmann (2019); Haghshenas et al. (2019); Carleo and Troyer (2017); Kochkov and Clark (2018) and many of the tricks from that community will be applicable here. Although we used PQMC in this Letter, any QMC which works in an orthogonal basis (world-line Monte Carlo, path integral ground state, etc.) have the same sign problem and are equally improved by the rotation. In addition, we developed a variational optimization approach which instead optimizes the variational energy of the uniform state $|\Psi_{u(R)}\rangle$ with respect to $|H(R)|$ over a small number of configurations, for a significant speedup and at polynomial scaling. This heuristic works because of the empirical observation that this variational energy and the projector energy track each other closely up to some (slowly varying) fixed offset.

To illustrate the method, we optimized the single-particle basis of $4\times 4$ and $16\times 4$ Hubbard model at $n=0.5$ doping and $U/t=1$ and $U/t=2$ respectively. This resulted in an exponential increase in average sign, by approximately $\exp\left[6.278\beta\right]$ and $\exp\left[12.86(5)\beta\right]$ leading to significant speedups for a typical $\beta\in[10,100]$. The rotations we found are sparse, generalize to a wider range of $U$ and $n$, and the latter rotation generalizes to all width-4 cylinders.

A uniting theme among these optimized rotations is their structure of rotating into linear combinations of a small number of nearby sites. The new basis in both cases seem relatively insensitive to both doping and $U/t$, with each $b_{(i,j)}$ not including the matching $c_{(i,j)}$. In addition, a 4-operator plaquette is apparent in both geometries, although in slightly different forms. Other techniques beyond QMC simulations may also benefit from this insight, although its unclear if it extends beyond the Hubbard model itself. The automatic optimization of our method should be useful for other sign-plagued models to potentially turn up other sparse bases.

Acknowledgements.

Acknowledgments:— We acknowledge Luiz Santos for helping to initiate this work trying to find a sign-problem mitigating unitary, early collaboration on this project, and useful discussions and insights related to this work. BKC thanks David Ceperley with whom he collaborated on a related unpublished project as a graduate student that attempted to mitigate the sign problem by minimizing the sum of the badly-signed off-diagonal matrix elements. BKC also thanks Shiwei Zhang, Miles Stoudenmire, and Giuseppe Carleo for conversations. RL thanks Shivesh Pathak for useful conversations and Yubo Yang for visualization assistance. We acknowledge support from the Department of Energy (DOE) Award No. DE-SC0020165. This project is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This work also made use of the Illinois Campus Cluster, a computing resource that is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with the National Center for Supercomputing Applications (NCSA) and which is supported by funds from the University of Illinois at Urbana-Champaign. Note: During the final stages of this work, a number of related works appeared on the arXiv Torlai et al. (2020), Hangleiter et al. (2020), and Klassen et al. (2020).

Appendix A Supplementary Material

A.1 Hubbard Chain Rotations

A.1.1 Sign Free Basis Rotations

Our algorithm developed two basis rotations that completely removed the sign problem for all values of $U$ but did not diagonalize the Hamiltonian directly. For the $N=3$, $n=2/3$ unitary matrix is

[TABLE]

and for $N=4$, $n=1$ the unitary matrix is

[TABLE]

We illustrate the $N=4$ rotation and optimization in fig. S1, showing the Hamiltonian under the rotation of the above unitary.

A.1.2 8 Site Hubbard Model Optimization

The final unitary for the 8 site Hubbard chain is

[TABLE]

A.2 Exploring the $4\times 4$ Hubbard Model

A.2.1 Bulk Measurement of the Sign problem

Using PQMC, we can evaluate the sign problem directly for the 4x4 Hubbard model. We measure $\Delta E=E-\overline{E(I)}$ in fig. S2 using energies from Dagotto et al. (1992) with fixed $S_{z}=0$. As expected, for intermediate to high $U$ for any number of electrons results in a sign problem. For $n=0.875$, there is a maximum in the sign problem for $U/t=4$, a doping that has been noted for its difficulty in numeric simulations LeBlanc et al. (2015). Even at half-filling, where other techniques do not have a sign problem Wang et al. (2015), PQMC simulations would suffer strongly.

A.2.2 Optimized Unitary Rotation

During the optimization, shown in fig. 2a, we can track the ratio of $\Delta E$ to compare with the true ground state fermionic energy known from exact diagonalization. This is tracked in fig. S3. At the optimal point, step 24, we reach a $\Delta E(R)/\Delta E(I)=0.176(1)$, while the rounded sparse unitary reaches $\Delta E(R)/\Delta E(I)=0.02$, a factor of 8.8 improvement.

The rounded sparse unitary can described as

[TABLE]

where

[TABLE]

and $s_{1}=s_{2}=-1$; $s_{3}=s_{4}=1$ and the indices are taken $\mod 4$. To create this unitary, we start with the original optimized unitary $R_{i}j$ at step 24 of our optimization, we then set to zero any values $\left|R_{ij}\right|<0.1$ and equally weight the remaining values, keeping the original signs intact. This new unitary has 64 non-zero entries, out of a possible 256 entries.

In fig S4 we compare the average sign of the identity rotation against our optimized basis finding significant improvement and allowing projection out to much larger $\beta$ albeit at the cost of a denser basis. We also compare against another explicit approach to mitigate sign-problems, the use of annihilation Booth et al. (2009); Kalos and Pederiva (2000) of opposite signed walkers as is used in FCIQMC (while the annihilation technique of FCIQMC is used here, we do not apply the initiator approximation). While the sign-problem is insensitive to the number of walkers used in PQMC, the number of walkers does matter when using annihilation. Using a number of walkers which is approximately $1\%$ of the Hilbert space size, we find two interesting things. First, in this example, basis rotation improves the sign problem significantly more then improvements seen using annihilation. Second, annihilation and basis rotation can be combined and, in this situation, we find annihilation actually helps more in the case where the basis rotation has already partially mitigated the sign problem.

A.3 Scaling the optimization - $6\times 6$ Hubbard Model

Because the scaling or our algorithm with system size is polynomial, one can perform projector optimization of $\overline{E(R)}$ over larger system sizes; we optimize the 630-parameter space of a $6\times 6$ Hubbard model at $U/t=1$ and $n=0.5$ over a random starting basis using finite differences. We successfully use this method over 25 steps to improve the average sign from a random rotation by $\exp\left[11.09(4)\beta\right]$ (see the inset of fig. S5).

Similarly, we can perform optimization on $\overline{E_{V}(R)}$ using finite differences, shown in fig. S5. We update optimization parameters (hyperparameters) every 100 steps or so, corresponding color changes. Steps 0-549 use 11 states to calculate derivatives, steps 550-749 use 21 states, and steps 750-850 use 41 states. The step size $\gamma$ also varies: steps 0-249 use $\gamma=0.01$, steps 250-649 use $\gamma=0.005$, steps 650-749 use $\gamma=0.001$, and steps 750-849 use $\gamma=0.0001$. The final rotation found improves the average sign by $\approx\exp\left[2.95(3)\beta\right]$ over not rotating ($R=I$).

A.4 Hyperparameters and Optimization of the $16\times 4$ Optimization

During optimization of the $16\times 4$ Hubbard model at $n=0.5$ and $U/t=2$, we use 20 walkers fixed throughout the 18800 steps, 19 randomly chosen alongside the completely doubly occupied state. For the first 100 steps we use $\gamma=0.01$ then $\gamma=0.001$ from steps 100-500, followed by $\gamma=0.0001$ for 500 steps. Afterward, we apply 8 sequences of 100 steps of $\gamma=0.001$ followed by $\gamma=0.0001$ for 1000 steps, until step 8800 which uses $\gamma=10^{-5}$

To round the final noisy unitary, we set to zero any values $\left|R_{ij}\right|<0.1$ and equally weight the remaining values, keeping the original signs intact. This new unitary has 160 non-zero entries, out of a possible 4096 entries.

A.5 Energy Scaling of the Width-4 Hubbard Models

First, we measure for various doping $\overline{E_{V}(R)}$ for the $8\times 4$, $16\times 4$ and $32\times 4$ lattices using using 500, 500, and 250 sampled states respectively, shown in fig. S6. Using the unitary construction defined by eqs. 7,8 (shown in the first row of fig. S6), as system sizes increase, the larger $U/t$ values begin to have less of a sign problem. For smaller $U/t$, the energy gap per site remains relatively constant, suggesting that the rotation is mitigated compared to rotating for $L\gg 1$.

To test this, we look at $\overline{E(I)}-\overline{E(R)}$ at $U/t=4$ and $n=0.5$ for $L=4,8,16,32$ $L\times 4$ Hubbard models, shown in fig. S7. Fitting the energies to a line (excluding $L=4$) we find a new a variational lower bound to the ground state energy of $-1.23544(5)$/site, greater than the lower bound of -1.429490(4)/site from not rotating.