How to understand magnetization work? What is the best way to understand that the magnetization work has the form
$$\delta W=HdM?$$

On a side note: Is this even the correct formula, or should one rather use
$$\delta W=BdM?$$
I'm guessing the difference comes from assuming only weak magnetization, so that $$B=H+M\propto H$$ but I'm not sure.
 A: The work done $\delta W$ by the system in a time $\delta t$, by an electric field $\mathbf{E}$ acting on a current density is $\mathbf{J}$:
$$ \delta W = \delta t \int_V\mathbf{J}\cdot\mathbf{E}\,\mathrm{d}V,$$
where we assume a quasi-static and reversible process so as to ignore hysteresis.
You can get the charge density (due to free, not bound, charges) $\mathbf{J}$ from Maxwell's equations:
$$ \nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}, $$
so that:
$$ \delta W = \delta t \left \{ \int_V (\nabla\times \mathbf{H})\cdot\mathbf{E}\,\mathrm{d}V - \int_V \frac{\partial \mathbf{D}}{\partial t}\cdot\mathbf{E}\,\mathrm{d}V  \right \} = \\  = \delta t \left \{ \left [ \int_V \nabla\times (\mathbf{H}\times\mathbf{E})\,\mathrm{d}V  + \int_V \mathbf{H}\cdot(\nabla\times \mathbf{E})\,\mathrm{d}V \right ] - \int_V \frac{\partial \mathbf{D}}{\partial t}\cdot\mathbf{E}\,\mathrm{d}V  \right \} . $$
The $\nabla \cdot (\mathbf{H}\times\mathbf{E})$ can be cast as a surface integral at the surface $S$ around $V$, via Gauss' theorem, and set to zero via "physical" boundary conditions, that is $\mathbf{E}, \mathbf{H} \rightarrow 0$ away from the (localised) charges and currents.
The $\nabla \times \mathbf{E}$ can be re-cast via another of Maxwell's equations:
$$ \nabla\times\mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}. $$
SO, the work done by the system is:
$$  \delta W =  \int_V\,\mathrm{d}V \,\mathbf{H}\cdot \mathrm{d}\mathbf{B} + \int_V\,\mathrm{d}V \,\mathbf{E}\cdot \mathrm{d}\mathbf{D}. $$
Just focussing on magnetic fields:
$$  \delta W =  \int_V\,\mathrm{d}V \,\mathbf{H}\cdot \mathrm{d}\mathbf{B}. $$
Using $\mathbf{B} = \mu_0 (\mathbf{H} + \mathbf{M})$, one can also re-write the above as:
$$  \delta W =  \mu_0\int_V\,\mathrm{d}V \,\mathbf{H}\cdot \mathrm{d}\mathbf{H} + \mu_0\int_V\,\mathrm{d}V \,\mathbf{H}\cdot \mathrm{d}\mathbf{M}, $$
whose physical interpretation is work done when changing external field $\mathbf{H}$ (first term) and when inducing the internal system's response ($\mathbf{M}$).
A: The answer depends on what you mean by "$H$", because there are several equally important H-fields one can define.
 To proceed let $\mathbf {\mathfrak{H}}$ denote the external biasing magnetic field, say, one that is generated by a permanent magnet or by a coil energized by a constant current source, and place some magnetizable material in the field where, of course, $\nabla \times \mathbf {\mathfrak{H}} = 0$ and $\mathbf {\mathfrak{B}}=\mu_0\mathbf {\mathfrak{H}}$ is the associated induction field in vacuum, and let $\mathbf {M}$ denote the magnetic moment density of the polarizable material we place in this field.
 In its most $practical$ form, the volume density of infinitesimal magnetic energy/work required to increase the polarization from $\mathbf {M}$ by $\delta \mathbf {M}$ is $$\delta W_0 = \mu_0 \mathbf {\mathfrak{H}}\cdot\delta \mathbf {M} \tag{0}\label{0}$$ I say this is the most practical because this work is accessible to direct measurement by being dependent on the $external$ field that is under the control of the experimenter. If the measurement of $\delta W_0$ is done isothermally then $$\delta \mathfrak{F} = \delta W_0-SdT\\=\mu_0\mathbf {\mathfrak{H}}\cdot\delta \mathbf {M}-SdT$$ represents the change in the $magnetic free energy$ of the polarizable material.
The problem with this approach is that any constitutive $functional$ $relationship$ $\mathbf{M}=\mathbf{M}(\mathfrak{H})$  between the induced polarization $\mathbf{M}$ and the external bias field $\mathfrak{H}$ depends also on the shape of the magnetizable sample and is dependent not just on the material itself. The reason for the shape dependence is the magnetic charges (poles) that are induced within and of the surface on the sample and the accompanying so-called $demagnetizing$ $field$ here denoted by $\mathbf{\tilde{H}}$. The volume pole and surface pole density is $p=-\nabla \cdot \mathbf{M}$ and $\mathbf{M}\cdot \mathbf{\hat{n}}$, resp. where $\mathbf{\hat{n}}$ is the surface normal. 
The integral
$$U_p=\frac{\mu_0}{2}\int \int \frac{dp_1 dp_2}{r_{12}}$$ represents the total magnetostatic interaction energy of an arbitrary configuration of infinitesimal poles $dp$. (The factor $\frac{1}{2}$ is there because we run through the poles $twice$ during the double integration.) 
With some vector analytic gymnastics you can show that $U_p=-\frac{\mu_0}{2} \int \mathbf {\tilde {H}}\cdot\mathbf{M}d\mathcal{v}$ where $\mathbf {\tilde {H}}=\int {dp}\frac{\mathbf{r}^0}{r^2}$; and with a little more gymnastics you can also show that $\int \mathbf {\tilde {H}}\cdot \delta \mathbf{M} d\mathcal{v} =\int \delta \mathbf{M} \cdot \mathbf {\tilde {H}}d\mathcal{v}$ and therefore $$\delta U_p= -\mu_0\int \mathbf {\tilde {H}}\cdot \delta \mathbf{M} d\mathcal{v} $$
 The total field inside the sample is the $sum$ of the external bias field and the field of the induced poles 
$$\mathbf{H}=\mathbf{\mathfrak{H}}+\mathbf {\tilde {H}} \tag{1}\label{1}$$
A relationship $\mathbf{M}=\mathbf{M}(\mathbf{H})$ is now presumably $independent$ of the sample shape and depends only on the material, so it would be constitutive. Then instead of $\delta W_0$ we can calculate the work involved creating both the polarization and the field within:
$$\delta W_1 = \mu_0\mathbf{H}\cdot \delta\mathbf{M} \\
= \mu_0\mathbf{\mathfrak{H}}\cdot \delta\mathbf{M}+\mu_0\mathbf {\tilde {H}}\cdot \delta\mathbf{M} \\$$ and integrating over space we get
$$\delta w_1 = \int \delta W_1 d\mathcal{v} 
\\= \delta w_0 - \delta u_p 
=  \int \delta W_0 d\mathcal{v} - \int \delta U_p d\mathcal{v} \tag{2}\label{2}$$
Let the sample be an ellipsoid of volume $\mathcal{V}$ and place it in a homogeneous external field $\mathfrak{H}$, then it is known that the internal field $\mathbf{H}$ will also be homogeneous and its components are related to the demagnetizing tensor components so that $M_k = D_k H_k$ (k=1,2,3) and $u_p= \mathcal{V} \frac{\mu_0}{2} \sum_k D_k H_k^2$.
Since $\mathfrak{H}$ is the $actual$ control variable in an experiment it is more convenient and practical to represent the work/energy relations with a kind of magnetic "free enthalpy" $\mathfrak{G}$ rather than using $\delta \mathfrak{F}$ as in $\eqref{0}$:
$$\delta W_2 = \delta (W_0 - \mathbf {\mathfrak{H}}\cdot \mathbf {M} )\\
= - \mathbf {M}\cdot\delta \mathbf {\mathfrak{H}} \\
\delta \mathfrak{G} = - \mathbf {M}\cdot\delta \mathbf {\mathfrak{H}} - SdT
\tag{3}\label{3} $$.
The usual method of finding $\mathbf{H}$ and $\mathbf{M}=\mathbf{M}(\mathbf{H})$ from knowing $\mathfrak{H}$ is either by employing a sample that is of the form of a long and thin needle where the surface poles are $negligible$ and thus $\mathbf {\tilde {H}} \approx 0$ and $\mathbf{H} \approx \mathfrak{H}$, or a toroid where there are no surface poles and thus $\mathbf {\tilde {H}} = 0$ and $\mathbf{H}= \mathfrak{H}$.
