  A variational method for studying the ground state of strongly interacting
quantum many-body bosonic systems is presented. Our approach constructs a class
of extensive variational non-Gaussian wavefunctions which extend Gaussian
states by means of nonlinear canonical transformations (NLCT) on the fields of
the theory under consideration. We illustrate this method with the one
dimensional Bose-Hubbard model for which the proposal presented here, provides
a family of approximate ground states at arbitrarily large values of the
interaction strength. We find that, for different values of the interaction,
the non-Gaussian NLCT-trial states sensibly improve the ground state energy
estimation when the system is in the Mott phase.
