We present a class of correlated wavefunctions generated by exponentials of two-body on-site Hermitian operators that can be evaluated with polynomial computational cost via a Hamiltonian similarity transformation. Wavefunctions of this form have been studied with variational Monte Carlo methods, but we present a formalism to perform non-stochastic calculations. The Hausdorff series generated by these Jastrow factors can be summed exactly without truncation resulting in a set of equations with polynomial computational cost. The correlators include the density-density, collinear spin-spin, spin-density cross terms, and on-site double occupancy operators. The resulting non-Hermitian many-body Hamiltonian can be solved in a biorthogonal mean-field approach with only a small set of correlation terms required for accurate calculations in systems with local interactions. Although the energy of the model is unbound, projective equations in the spirit of coupled cluster theory lead to well-defined solutions. The theory is tested on the one and two-dimensional repulsive Hubbard model where it yields accurate results for large systems with low computational cost. Symmetry projection methods are included to further improve the reference wavefunction and results under strong correlation without sacrificing good quantum numbers resulting in very accurate energies for small systems and producing a better ground state for the calculation of other properties.