A bstract 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 bonanza shown by these variational 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 ϕ 4 scalar field theory in 2D. Finally, we analyze the rôle played by higher order correlators in our results and show that strong subadditivity is satisfied.