  In this work we provide a method to study the entanglement entropy for
non-Gaussian states that minimize the energy functional of interacting quantum
field theories at arbitrary coupling. To this end, we build a class of
non-Gaussian variational trial wavefunctionals with the help of exact nonlinear
canonical transformations. The calculability \emph{bonanza} shown by these
variational \emph{ansatze} allows us to compute the entanglement entropy using
the prescription for the ground state of free theories. In free theories, the
entanglement entropy is determined by the two-point correlation functions. For
the interacting case, we show that these two-point correlators can be replaced
by their nonperturbatively corrected counterparts. Upon giving some general
formulae for general interacting models we calculate the entanglement entropy
of half space and compact regions for the $\phi^4$ scalar field theory in 2D.
Finally, we analyze the r\^ole played by higher order correlators in our
results and show that strong subadditivity is satisfied.
