  We propose a numerical method which embeds the variational non-Gaussian
wavefunction approach within exact diagonalization, allowing for efficient
treatment of correlated systems with both electron-electron and electron-phonon
interactions. Using a generalized polaron transformation, we construct a
variational wavefunction that absorbs entanglement between electrons and
phonons into a variational non-Gaussian transformation; exact diagonalization
is then used to treat the electronic part of the wavefunction exactly, thus
taking into account high-order correlation effects beyond the Gaussian level.
Keeping the full electronic Hilbert space, the complexity is increased only by
a polynomial scaling factor relative to the exact diagonalization calculation
for pure electrons. As an example, we use this method to study ground-state
properties of the two-dimensional Hubbard-Holstein model, providing evidence
for the existence of intervening phases between the spin and charge-ordered
states. In particular, we find one of the intervening phases has strong charge
susceptibility and binding energy, but is distinct from a charge-density-wave
ordered state, while the other intervening phase displays superconductivity at
weak couplings. This method, as a general framework, can be extended to treat
excited states and dynamics, as well as a wide range of systems with both
electron-electron and electron-boson interactions.
