We study the variational solution of generic interacting fermionic lattice systems using fermionic Gaussian states and show that the process of gaussification, leading to a nonlinear closed equation of motion for the covariance matrix, is locally optimal in time by relating it to the time-dependent variational principle. By linearizing our nonlinear equation of motion around the ground-state fixed point we describe a method to study low-lying excited states leading to a variational method to determine the dispersion relations of generic interacting fermionic lattice systems. This procedure is applied to study the attractive and repulsive Hubbard model on a two-dimensional lattice, as well as the stability of the Hofstadter butterfly structure in the presence of interactions.