The method is based on sampling matchings, not uniformly, but with stationary distribution.
Recall that for any state function f(x), its time average value is also given by å_x p_x f (x), so that we would have the following.
A function of the type f (M) := 1 if |M| = 0 and f (M) := 0 if |M| >0 (just a counting of the times the zero matching is sampled) would provide the following, and therefore a direct computation of Z(l). However this would result in a slow converging method.