My rough understanding about lattice simulations of bosonic quantum field theories is that the partition function can be approximated by explicitly summing over a large number of field configurations, chosen by some Monte Carlo algorithm for instance. But is there a similar way to simulate fermions in order to get at nonperturbative effects?
If an interaction is quadratic in the fermionic fields, like the Yukawa interaction $ \phi \bar{\psi} \psi$, we can evaluate the fermionic path integral analytically for each particular field configuration $\phi$ in our Monte Carlo algorithm. But what if there is something not quadratic like a four-fermion interaction $(\bar{\psi}\psi)^2$?
I know that if the number of lattice sites is finite the perturbation series must end at some finite order (since there are a finite number of Grassmann variables), so in principle this can be solved `nonperturbatively' by simply summing over all terms. But of course that would be extremely inefficient, so how is this treated in practice?