Restriction of sparse operator representations to fixed particle manifold

[babbush]

Very often we would like to analyze properties of operators on a restricted particle number manifold. For instance, many (in fact nearly all) simulation algorithms depend in some way on the norm of ||H||. But if you know that your simulation will occur on the fixed electron number manifold (almost always the case), then you are really only concerned with the highest eigenvalue of H that has an eigenvector confined to that manifold. Another example is that you'd like to simulate a particular state, perhaps the ground state of a jellium or Hubbard model Hamiltonian. But the ground states of these models are vacuum - what you really want is the ground state in a particular particle sector. You want to use Lanczos algorithm for efficiency but that's not easy to unless you first project the operator into the desired particle number manifold.

Usually, the task is to obtain a state on a fixed particle number manifold and then take expectation values of that state with other operators. A very concrete example is that we'd like to know the effect that Trotter error has on the eigenvalues estimated by phase estimation. So what we can do is to take the expectation value of the Trotter error operator with the ground state of our Hamiltonian. But again, perhaps we are looking at a Hubbard model and want the ground state on that particular manifold.

Since we usually use the sparse representation of states and operators for such purposes, that is where this code will be most relevant. But what is the best strategy for performing this task? For instance, one could delete all the rows and columns of the matrix on the wrong sector to make the operator more compact. At least in the occupation-number basis (the one associated with Jordan-Wigner), I think this is not too complicated. I suppose you need to look at the Hamming weight of the binary representation of the row/column number. Right?

@jarrodmcc any thoughts on this? What strategy would you suggest?

[idk3]

At least in the occupation-number basis (the one associated with Jordan-Wigner), I think this is not too complicated. I suppose you need to look at the Hamming weight of the binary representation of the row/column number. Right?

Can you give an example? Both $Z_2$ and $Z_1 Z_2$ have their nonzero entries in the same positions, so I'm not completely sure what this means.

[jarrodmcc]

I think removing select rows and columns from the Hamiltonian might be a bit more difficult than expected for this. Two (human-time) simpler approaches I would attempt are

1. Lanzcos and similar sparse diagonalization methods won't deviate from the symmetry subspace the initial state starts in. Change the current sparse routines to take an initial vector, and give it one with the correct particle number. This is also a good idea from the point of view of number of iterations and efficiency.

2. Construct a sparse representation of the symmetry operator you want to maintain. For example the number operator N. Add a penalty term to the Hamiltonian that is of the form C (N - n * I)^2 where C is some positive constant, n is the number of electrons you would like, and I is the identity. For modest values of C, this will shift the eigenspectrum of undesired states upwards without contaminating the values you are interested in.

A third, perhaps more time-consuming option, is to explore the symmetry reductions done in Bravyi's latest paper, but this is an alternative to JW and may introduce other complications depending on what you are looking for.

[jarrodmcc]

Now that I think about it a little more, I suppose constructing an operator only in the space of the particle number you are interested in wouldn't be too hard either. This is what's usually done in CI codes for quantum chemistry.

I guess my thinking it would be a little tricky came from the way we currently construct operators in the computational basis that might not generalize to representations that are not JW. That is, we currently build the full operator in the computational basis, which could then be sliced down using something like numpy slice operations. This works as long as the computational basis corresponds to eigenstates of the operator we want to slice with respect to, which is the case for JW but probably not other representations.

A more general solution would enumerate states using all possible (a_{i0}^\dagger ... a{i_n}^\dagger), where n is the number of particles, and build the Hamiltonian in that basis. Our representation would be quite slow for doing this, but is possible. In the future, if both our QubitOperator and FermionOperator had matrix-free actions on these spaces, it would be much more efficient.

[babbush]

Thanks Jarrod!