Compression of correlation matrices and an efficient method for forming matrix product states of fermionic Gaussian states

Here we present an efficient and numerically stable procedure for compressing a correlation matrix into a set of local unitary single-particle gates, which leads to a very efficient way of forming the matrix product state (MPS) approximation of a pure fermionic Gaussian state, such as the ground state of a quadratic Hamiltonian. The procedure involves successively diagonalizing subblocks of the correlation matrix to isolate local states which are purely occupied or unoccupied. A small number of nearest-neighbor unitary gates isolate each local state. The MPS of this state is formed by applying the many-body version of these gates to a product state. We treat the simple case of compressing the correlation matrix of spinless free fermions with definite particle number in detail, although the procedure is easily extended to fermions with spin and more general BCS states (utilizing the formalism of Majorana modes). We also present a density matrix renormalization group–like algorithm to obtain the compressed correlation matrix directly from a hopping Hamiltonian. In addition, we discuss a slight variation of the procedure which leads to a simple construction of the multiscale entanglement renormalization ansatz of a fermionic Gaussian state, and present a simple picture of orthogonal wavelet transforms in terms of the gate structure we present in this paper. As a simple demonstration, we analyze the Su-Schrieffer-Heeger model (free fermions on a one-dimensional lattice with staggered hopping amplitudes).

Figure 1

(a) Shows the occupations $n_b$ and corresponding entanglement $S_1\left(n_b\right)$ from diagonalizing a block of $B=16$ sites in the middle of a system of free gapless fermions on $N=1000$ sites at half-filling. The minimum and maximum eigenvalues $n_1$ and $n_{16}$ differ from 0 and 1 by $\approx 1.74\times {10}^{-11}$. The eigenvalues closest to $\frac{1}{2},\frac{1}{2}-n_8=n_9-\frac{1}{2}\approx 0.21$, have entropies $S_1\left(n_8\right)=S_1\left(n_9\right)\approx 0.60$, which are close to the maximum of $S_1\left(\frac{1}{2}\right)=ln\left(2\right)\approx 0.69$. (b) Shows examples of eigenvectors from the same diagonalization. The eigenvectors with eigenvalues near 0 and 1, which contribute very little to the entanglement, are localized in the middle of the block, while the eigenvectors with eigenvalues closer to $\frac{1}{2}$, which contribute most to the entanglement, have large support on the edges of the block.

Figure 2

Examples of approximate occupied and unoccupied eigenvectors of $\mathrm{\Lambda }$ obtained from diagonalizing ${\mathrm{\Lambda }}_{\mathcal{B}}$ where subblock $\mathcal{B}$ are sites $1,...,B$. $\mathrm{\Lambda }$ is formed from the ground state of $\widehat{H}=-t{\sum }_{i=1}^{N-1}({\widehat{a}}_i^{†}{\widehat{a}}_{i+1}+\text{H.c.})$ for $N=1024$ at half-filling $\left(N_F=N/2\right)$. A block size of $B=12$ is used. Eigenvectors with highest $\left(n_{\text{occ}}\right)$ and lowest $\left(n_{\text{unocc}}\right)$ eigenvalues found from diagonalizing subblock $\mathcal{B}$ are shown. We find $1-n_{\text{occ}}=2.4\times {10}^{-15}$ and $n_{\text{unocc}}=7.3\times {10}^{-16}$, so the occupations are accurate to nearly machine double precision. $1-n_{\text{occ}}$ and $n_{\text{unocc}}$ should be equal at half-filling (because of particle-hole symmetry), but are different in this case as a result of roundoff errors.

Figure 3

Definition of a gate used throughout the paper. Example for $N=6$ sites for a gate at site $i=3$. Unless specified otherwise, circuits are in a direct sum space. We take the convention that multiplying a matrix from the top by a vector corresponds to multiplying the matrix on the right by a column vector.

Figure 4

In (a) we show schematically the procedure to obtain, given an approximate eigenvector $\vec{v}$ of the correlation matrix $\mathrm{\Lambda }$, the set of local rotation gates that make up our compressed correlation matrix. The example shown is for a block size $B=4$ and system size $N=8$. (b) Shows that, by conjugating the correlation matrix by the gates obtained, the correlation matrix is approximately partially diagonalized.

Figure 5

(a) Shows the overall gate structure obtained by the diagonalization procedure. These gates form the total $N\times N$ unitary $V$ which approximately diagonalizes our correlation matrix $\mathrm{\Lambda }$. By conjugating a diagonal matrix with the appropriate occupations of 0 or 1 found in the diagonalization procedure by this set of gates, we get an approximation for the correlation matrix. (b) Shows an example of the correlations allowed by representing the correlation matrix $\mathrm{\Lambda }$ with a diagonal matrix conjugated by a finite depth circuit of depth $<N/2$. The gray area (the “light cone”) represents sites where there can be nonzero correlations with the first site. The circles in the middle represent a diagonal matrix with 1's and 0's on the diagonal, which is conjugated by a unitary change of basis approximated here by a finite depth circuit. For the circuit depth shown, there can not be correlations with the last two sites. A circuit of depth $\ge N/2$ is required to allow for the possibility of nontrivial correlations across the entire system.

Figure 6

An example of an alternative diagonalization scheme resulting in a MERA-like gate structure. Here, we show a section of the first two renormalization steps, with 12 sites shown in the first layer and 6 renormalized sites shown in the second. A block size of $B=4$ is used. For this block size, there are two layers of disentanglers and one layer of isometries per level of the MERA. Open legs at the top of each layer correspond to diagonal modes of the correlation matrix (with eigenvalues 0 or 1) and are ignored at the next layer.

Figure 7

Here, we show an example of a discrete wavelet transform written in the gate notation introduced in this paper. We show the D4 wavelet, which corresponds to a fermionic Gaussian MERA with one layer of disentanglers and one layer of isometries per layer. ${\mathrm{w}}_1$ and ${\mathrm{s}}_1 ({\mathrm{w}}_2$ and ${\mathrm{s}}_2$) label wavelet and scaling functions for the first (second) layer. Taking ${\theta }_1=\pi /6$ and ${\theta }_2=5\pi /12$, we reproduce the conventional scaling coefficients for the D4 WT, ${\vec{a}}^T=(a_1,a_2,a_3,a_4)=(1+\sqrt{3},3+\sqrt{3},3-\sqrt{3},1-\sqrt{3})/\left(4\sqrt{2}\right)$.

Figure 8

Here, we show explicitly how to obtain the scaling and wavelet coefficients of the D4 WT from the circuit construction. Taking ${\theta }_1=\pi /6$ and ${\theta }_2=5\pi /12$, in (a) and (b) we reproduce the conventional scaling coefficients for the D4 WT, ${\vec{a}}^T=(a_1,a_2,a_3,a_4)=(1+\sqrt{3},3+\sqrt{3},3-\sqrt{3},1-\sqrt{3})/\left(4\sqrt{2}\right)$, up to a trivial reversal in the order. In (c) with the same choice of angles we reproduce the conventional wavelet coefficients $(a_4,-a_3,a_2,-a_1)$, again up to a trivial reversal and sign.

Figure 9

Tensor diagram showing the structure of gates $\left\{{\widehat{V}}_{{\mathcal{B}}_i}\right\}$ for $i=1,...,N-1$ obtained in our procedure and how they contract to form an MPS. The white (black) triangles represent projectors onto the occupied (unoccupied) state. The ordering of the occupied and unoccupied states is determined by the ordering of the occupations found in the diagonalization procedure, one particular example at half-filling is shown here. Here, we show a system with $N=8$ sites and a block size $B=4$. The diagram on the right shows that once the sites are rotated into a basis where one of the modes is occupied or unoccupied (generally with some alternating pattern), the fully occupied or unoccupied modes can be projected out. The transformations $\left\{{\widehat{V}}_{{\mathcal{B}}_i}\right\}$, including the projections, can be directly interpreted as the tensors composing the MPS representation of our many-body state if we do an exact contraction, or we can apply them as a set of gates as explained in the text.

Figure 10

Block size required to obtain a relative error in the total energy of less than ${10}^{-6}$ as a function of the calculated energy gap (in units of $t$) for $N=128$ sites.

Figure 11

Relative error in the total energy as a function of the block size $B$ for $N=128$ sites and $\delta =0$.

Figure 12

Examples of occupied and unoccupied modes found in the diagonalization process. (a) Shows occupied/unoccupied modes for $\delta =0.4$ (energy gap $\approx 0.806135t$). (b) Shows occupied/unoccupied modes for $\delta =0$ (energy gap $\approx 0.146088t$).

Figure 13

Examples of deviations in occupations at the end of the diagonalization procedure for $N=128$ sites. (a) Shows errors in the occupations for $\delta =0.4$ (energy gap $\approx 0.806135t$). (b) Shows errors in the occupations for $\delta =0.0$ (energy gap $\approx 0.146088t$).

Figure 14

Block size $B$ needed for a relative error in the energy of $<{10}^{-6}$ as a function of number of sites $N$ for spinless, gapless fermions with open boundary conditions at half-filling. As expected from arguments about the entanglement of a critical system, we find $B\sim ln\left(N\right)$, tested up to $N=2^{16}=65536$ sites (note the $log$ scale on the $x$ axis). To study systems of this size and avoid the $O\left(N^3\right)$ diagonalization of the hopping Hamiltonian, we obtain the correlation matrix using the GDMRG algorithm as explained in Appendix pp2.

Figure 15

Relative error in the energy for the proposed GMERA construction for increasing number of sites for a block size $B=10$. The system analyzed is the ground state of free fermions hopping on a lattice with open boundary conditions. All errors are below ${10}^{-6}$. As expected for a MERA, the error is seen to saturate for large $N$, indicating a fixed block size is sufficient to obtain an accuracy $<{10}^{-6}$ up to very large system sizes.

Figure 16

A plot of the time to form the MPS approximation of gapped and gapless free-fermion ground states at half-filling as a function of sites $N$ using gates from a GMPS. The bond dimensions are chosen large enough such that the relative errors in the energy of the MPS are below ${10}^{-6}$. The block size of the GMPS used to form the MPS are the minimum required to obtain a GMPS with a relative error in the energy of ${10}^{-6}$. A cutoff in the singular values of the SVD of ${10}^{-11}$ was used when applying the gates to form the MPS using the method described in Sec. 3d. For the gapped case, the SSH model with $\delta =0.1$ is used, corresponding to an energy gap of $\mathrm{\Delta }\approx 0.2t$ (exact as $N\to \infty$).

Figure 17

(a) Shows an example of an effective Hamiltonian centered at site 5 for the GDMRG algorithm. The effective Hamiltonian is the diagonal square block of the transformed Hamiltonian matrix covering the indicated sites. The example is for $N=10$ sites and a block size of $B=3$. Here, the center is only one site, but could be more to improve convergence just like in the standard DMRG algorithm. (b) Shows, for a sweep to the right, how to obtain the new left block from the previous effective Hamiltonian by rotating with the appropriate gates found and taking the submatrix of the indicated sites.