The auxiliary magnetic field is not zero in your example. The reason is that $\mathbf H$ must also obey continuity conditions at the boundary and these act as source terms. These source terms are the surface equivalents of the 'magnetization charge' $\nabla\cdot\mathbf M$, and they're present at surfaces where the magnetization has a discontinuity $\Delta\mathbf M$ with a nonzero component across the surface normal $\hat{\mathbf n}$.
As always, the continuity conditions at the boundary the limiting form of Maxwell's equations at the boundary. In this case they come from the magnetic Gauss law and the magnetostatic Ampère law, which for the $\mathbf B$ field read
$$
\oint_S\mathbf B\cdot\mathrm d\mathbf S=0
\quad\text{and}\quad
\oint_C \mathbf B\cdot \mathrm d\mathbf l=\mu_0I_\mathrm{enc},
\tag 1
$$
and in their boundary form
$$
\hat{\mathbf n}\cdot\Delta \mathbf B=0
\quad\text{and}\quad
\hat{\mathbf n}\times\Delta \mathbf B=\mu_0\mathbf K.
\tag 2
$$
This means that the normal component of $\mathbf B$ should be continuous across any surface. It should then be clear that the same cannot be true for the auxiliary $\mathbf H$ field, because $\mathbf H=\tfrac{1}{\mu_0}\mathbf B-\mathbf M$ is the sum of a continuous $\mathbf B$ field and a discontinuous magnetization $\mathbf M$.
Moreover, this discontinuity can be directly quantified from the boundary-form Maxwell equation by simply plugging in the definition of $\mathbf H$, to get
$$
\hat{\mathbf n}\cdot\Delta \mathbf H=-\frac{1}{\mu_0}\hat{\mathbf n}\cdot\Delta \mathbf M.
\tag 3
$$
This acts as a source term for $\mathbf H$ and it's what makes it nonzero throughout.
On the other hand, the specific case you mention, in which the entire space is filled with the magnetized material, is trickier to deal with. In general, infinite bulk sources (as opposed to infinite planes or lines) are iffy propositions in electromagnetism, and they don't always give physical results. (For an extreme example, a constant charge density throughout all space, there's simply no acceptable resolution.)
In such cases what one needs to do is to take a finite volume and then take the limit as it becomes bigger and bigger. One then needs to check that the answer is (i) independent of the shape of the volume, and (ii) independent of the relative position of the probe point with respect to the centre of expansion; if the answer is yes on both counts then that is your answer.
For the specific case when all of space is filled with a magnetized material, this seems to be doable. The easiest way to approach this is to use slabs of infinite area and finite width along the magnetization direction, and then take the width to infinity. Here in every finite case $\mathbf H=\mathbf M$ is easily verified from the boundary conditions, so $\mathbf B=2\mu_0\mathbf M$ is independent of the centre of expansion and therefore the correct answer.
As an independent check, it is relatively easy to verify that a constantly magnetized sphere of material produces a uniform magnetic field inside it, so you could take the limit as the radius goes to infinity and you will (presumably) get the same answer.
To round things up, I'll refer you to Jackson's Classical Electrodynamics §5.8 (pp. 191-194 on the 3rd edition, Wiley 1999), which provides a more thorough treatment. In particular, to get a complete formulation in terms of $\mathbf H$ you also need to provide a condition for its tangential component. To do this, you simply rephrase Ampère's law as
$$
\oint_C \mathbf H\cdot \mathrm d\mathbf l=I_\mathrm{free},
\tag 4
$$
which then boils down to
$$
\hat{\mathbf n}\times\Delta \mathbf H=\mathbf K_\mathrm{free}.
\tag 5
$$
In your case there are no free currents, so the tangential component of $\mathbf H$ is continuous everywhere (i.e. it behaves exactly like a static electric field!), but there is no escaping the discontinuity in the normal component (which is again analogous to that of an electric field).
Finally, note that Jackson includes some additional simplifications which can be brought to bear in linear media, in which case $\mathbf B=\mu\mathbf H$ and you can rephrase $(3)$ as $\hat{\mathbf n}\cdot \mathbf H_1=\frac{\mu_2}{\mu_1}\hat{\mathbf n}\cdot \mathbf H_1$ (assuming isotropic media, though that can be relaxed easily). However, this does not apply in your case since you have a fixed magnetization, so your medium is not linear.
