Fast update algorithm for the quantum Monte Carlo simulation of the Hubbard model

This paper presents an efficient algorithm for computing the transition probability in auxiliary field quantum Monte Carlo simulations of strongly correlated electron systems using a Hubbard model. This algorithm is based on a low rank updating of the underlying linear algebra problem, and results in significant computational savings. The computational complexity of computing the transition probability and Green’s function update reduces to $\mathcal{O}\left(k^2\right)$ during the $k\text{th}$ step, where $k$ is the number of accepted spin flips, and results in an algorithm that is faster than the competing delayed update algorithm. Moreover, this algorithm is orders of magnitude faster than traditional algorithms that use naive updating of the Green’s function matrix.

Figure 1

(Color online) (a) CPU times taken by delayed and submatrix update algorithms for one complete lattice sweep for a variety of reset sizes $q$. System size is $N=4800$ and the total number of lattice site spin flips is $m=4800$. The Green’s function matrix is updated every $q$ steps. It is clear that the computational cost of submatrix algorithm is almost constant with the reset size $q$, whereas the cost of delayed algorithm increases linearly with reset size $q$. (b) A detailed representation of CPU times taken by delayed and submatrix algorithms. Figure shows the total CPU time taken for spin flip operations and the Green’s function updates for a given reset size $q$. The time taken for spin-flip operations increases linearly for delayed algorithm while it remains constant and is extremely cheap for submatrix algorithm. On the contrary, the computational cost of the Green’s function update is comparable for both algorithms, albeit delayed algorithm is faster due to less number of operations.

Figure 2

(Color online) Further breakdown of CPU times taken by the submatrix algorithm in updating the Green’s function update [Eq. (30)]. In the legend, LHS solve and RHS solve refer to ${\mathbf{L}}_k^{-1}{\mathbf{G}}_0\left(\mathbf{p},:\right)$ and ${\mathbf{G}}_0\left(:,\mathbf{p}\right){\mathbf{U}}_k^{-1}$ operations, respectively, and diagonal scaling refers to multiplication by ${\mathbf{D}}_{1+\gamma }^{-1}$. These are the additional costs associated with submatrix algorithm. The computational cost associated with DGEMM would be same as that associated with delayed algorithm.