There is a good way to make the harmonic product spectrum more robust to the effects of noise. You generate a synthetic spectrum starting from a histogram of zero crossing intervals that has been "smeared" back into a continuous signal with Gaussian peaks by applying a kernel density estimate. The result can then be passed through the HPS to find the fundamental. The decorrelated zero crossings caused by noise are much less problematic than working with the HPS of the original signal.
Sounds interesting, do you an example of the performance you might expect from such a system? I must say that out out of the estimators I tried, HPS is the most difficult to deal with. You need information about how the harmonics are related to the fundamental and how many harmonics you expect before applying the algorithm. Maybe there are some adaptive methods that exists for this job but HPS is intrinsically ad-hoc, which is not very nice.