  The numerical computation of many hadronic correlation functions is
exceedingly difficult due to the exponentially decreasing signal-to-noise ratio
with the distance between source and sink. Multilevel integration methods,
using independent updates of separate regions in space-time, are known to be
able to solve such problems but have so far been available only for pure gauge
theory. We present first steps into the direction of making such integration
schemes amenable to theories with fermions, by factorizing a given observable
via an approximated domain decomposition of the quark propagator. This allows
for multilevel integration of the (large) factorized contribution to the
observable, while its (small) correction can be computed in the standard way.
