What is the permeability of a permanent magnet? The following is a diagram that shows the B- and H-fields of a permanent magnet.
Inside the magnet, the H-field is in the opposite direction to the B-field, because of the magnetisation M-field and the relationship
$$ \vec{B} = \mu_0 (\vec{H} + \vec{M})$$
(in SI units).
It is also commonly stated that the relationship between $\vec{B}$ and $\vec{H}$ can be written as $\vec{B} = \mu \vec{H}$, where in general, $\mu$ can be a tensor that depends on the H-field, position, time, temperature etc. (i.e. not necessarily restricted to linear media).
Yet when I look up the relative permeability for magnets (for example here for Neodymium magnets), it is given as a positive scalar $>1$. But surely if the B-field and H-field are in opposite directions then $\mu$ is negative?
What is the source of my confusion?

 A: The book that the wikipedia article refers to is "Design of Rotating Electrical Machines" by Pyrhonen et al, one can read some pages via Google. It is an engineering book, and it uses engineering-type concepts like "reluctance" that are useful to calculate "magnetic circuits" and things like that. 
In those kind of calculations, a  strong permanent magnet acts like an air gap, because its magnetization hardly changes in an applied field (as long as its value does not get too close to the coercive field in the wrong direction). So it can be modeled with $\mu_r \approx 1$. 
A: The relation between $B$ and $H$ (or for the electric case between $E$ and $D$) are completely material dependent and have no closed form expression. They are called constitutive relations. There are only models for these relations. The permeability model is just the simplest one.
Even when one would say $μ$ is tensorial and allowed to depend on space, time, $H$ and $B$, this (extended permeability) model is still not fully general. This is because there may be dynamic effects ($M$ may depend on past values of $H$) and/or dependencies on spatially neighboring points of $H$. Dynamic effects are very common (a.k.a. frequency dependence)
If one really wants to write down an equation capturing all this, one has probably something like this (this is from the Wikipedia link above):
$M = \frac{1}{\mu_0} \iiint d^3r' dt' \chi_m(r',t',r,t,H) H(r',t')$
However, I think this equation is not very useful. In my opinion, it is equally helpful to say $M$ is formed by applying a causal operator to $H$ (potentially nonlinear and time variant). But even that is not most general:
Magnetization can also be caused by other physical quantities such as temperature (Seebeck effect), electric field (Hall Effect), stress (Piezomagnetic effect), radio frequency waves (Nuclear magnetic resonance -> e.g. Bloch equaitions) with all these effects having different equations to describe magnetization.
Concluding, it is probably best to see magnetization (and its electric equivalent polarization) as a quantity which is given by many physical effects and has to be determined by the equations specific for the effects to be taken into account. A permeability based model is only valid for isotropic material, low fields and slow temporal variations.
For permanent magnet the first breaks down (fields are constantly there and are not "low") and the second very quickly with increasing frequency.
So for a magnetized magnet a permeability based model makes no sense.
One can specify a permeability for the (unmagnetized) material the magnet is made of, though. This is probably what you saw.
A good starting point is probably an advanced book on electromagnetism, such as Jacksons. More specific topics are likely handled in special book on magnetics, superconduction, spectroscopy (NMR or else methods), nonlinear optics or optoelectronics.
A: To your question "So is the concept of a relative permeability invalid for any material with a permanent magnetisation?" let me quote [1] 

(page 26) The formal theory of magnetostatics, as presented in elementary textbooks on electricity, is usually based on the linear approximation $\mathbf{M} = \chi \mathbf{H}$ or $\mathbf{B} = (1+4\pi\chi) \mathbf{H}$, where the susceptibility $\chi$ and the permeability $\mu = 1+4\pi\chi$ are assumed to be constants for a given material at a given temperature. For ferromagnetic materials, such a relation is never more than a very crude approximation usable over a limited range; it has some usefulness in practical applications, for example in the conventional "magnetic circuit" formulas, but it has almost no value in basic physical theory. In ferromagnetism, permeability - so far as it warrants introduction at all- is not a quantity to be used in the theoretical formulas
  but a quantity to be calculated by means of them. Furthermore, on the basis of the definitions that we have adopted for the flux density B and magnetizing force H at points of a magnetic specimen, there is no reason for even expecting that the magnetization will
  be a function of either of these field vectors. And in fact the physical theory of ferromagnetism justifies this scepticism: for sufficiently small bodies, the magnetizatiou depends, among other things, on the size and shape of the body, so that no M vs H relation
  exists.
(page 86) From these considerations it is clear that for ferromagnetics a "permeability" $\mu$, must never be defined simply as the value of $B/H$ but must be defined as a value of $B/H$, of $dB/dH$, or of $\Delta B/\Delta H$ under specified conditions. Such special definitions lead to such concepts as initial permeability, reversible permeability,
  incremental permeability, and ideal permeability, whose uses are limited and mostly of engineering nature.
  For calculation of the field in a gap, or at a more general point outside a magnetized body, it is often a good approximation to assume that B ~ H everywhere inside the material and therefore to replace it with a linear material for which $\mu = \infty$. The electrostatic analogue is a perfect conductor. The internal H is zero, and the flux lines emerge from the specimen surface along its unit normal; the specimen surface is an equipotential for the total scalar potential $\phi$, of which, in any region not occupied by conduction currents, H is the negative gradient.
  The specification $B =\mu H$, with $\mu$ assumed constant, fixes the magnetization distribution and therefore the demagnetizing factor of a non-ellipsoidal specimen. Even in this approximation, however, the calculation of the demagnetizing factor (now a function of $\mu$) is an unsolved problem of potential theory except in the following
  cases: (1) the ellipsoid and its degenerate forms; (2) the limit $\mu \to 1$, where the magnetization becomes uniform, the theorems of Ch. 3, § 3 hold, and the calculation reduces to a calculation of the volume-average H due to uniform M; (3) certain two-dimensional problems in the limit $\mu \to \infty$. Examples of (2) are the circular cylinder in a longitudinal or transverse field and the infinitely long rectangular bar in a transverse field. An example of (3) is the infinitely long rectangular bar in a transverse field. Since numerical values for these cases are not given in the standard
  literature, they are given in the Appendix, with some corresponding ellipsoid values for comparison.

Brown: Magnetostatic Principles in Ferromagnetism, North Holland, 1962
