What is an intuitive explanation for why group velocity is [math]\frac{\mathrm{d} \omega}{\mathrm{d}k}[/math] and phase velocity is [math]\frac{\omega}{k}[/math]?

That's pretty easy to answer. If you assume the wave is monochromatic, then if you want to travel with the wave so that the phase stays constant, you must travel at a speed where the argument of the cos function is constant. [math]cos(\omega t - kx) = const\rightarrow \omega t - kx = constvelocity = \frac{dx}{dt}[/math] [math]kx = \omega t - constx = \frac{(\omega t - const)}{k}\frac{dx}{dt} = \frac{w}{k}[/math] For group velocity, it is a little tougher.  But if you assume that instead of being monochromatic, the wave contains multiple frequencies that travel at different speeds, then you can compose the wave as a sum of a bunch of plane waves: [math]A(x,t) = \[ \int_{-\infty}^{\infty} A(k) e^{ i (kx - \omega(k)t }\,dk\][/math] where [math]\omega[/math] is a function of [math]k[/math]. If we assume that the wave is narrowband (that is, the frequencies aren't spread out that much), then we can approximate the function w(k) with its linearization: [math]w(k) ≈ \omega_0 + (k-k_0)\frac{d \omega}{dk}[/math] Moving this out, we get [math]A(x,t) = \[ \int_{-\infty}^{\infty} A(k) e^{ i (kx - \omega_0 t + k_0 \frac{d \omega}{dk} t - k \frac{d \omega}{dk} t }\,dk\][/math] extracting the terms that are not functions of k, you get [math]A(x,t) = \exp{it(\frac{d \omega}{dk} k_0 - w_0)} \[ \int_{-\infty}^{\infty} A(k) e^{ i (k(x - \frac{d \omega}{dk} t }\,dk\][/math] From there, you can see that the outside multiple gives it a wavy component, and it has a magnitude of 1. From there, you can take the magnitude of the whole function, which will give you just the amplitude -- or the height of the wave, which height is what you want to follow everywhere. Inside the integral, to maintain a constant height, you must keep the argument of [math]\exp{i k(x - \frac{d \omega}{dk} t)}[/math] constant.  To do that, keep [math]x - \frac{d \omega}{dk} t = constx = const + \frac{d \omega}{dk} t[/math] [math]\frac{dx}{dt} = \frac{d \omega}{dk}[/math] And there is your derivation right there. (If anyone wants to edit this and put it in LaTex to make it look better, feel free to do so, I just don't consider it worth the effort at the moment).

This is going to be a somewhat nonstandard answer that is the result of me trying to find a neat and intuitively satisfying explanation for what phase and group velocities are in the most general context possible. So, apologies if this turns out to fly in the face of what people working on optics typically deal with. If there are blatant factual inaccuracies in what follows, please correct me. First of all, the notion of velocity as defined for particles is naturally associated with real space. This is not surprising at all as particles are after all objects most naturally defined on real space. By contrast, waves are objects that are more naturally defined on the Fourier space, so it makes sense to say that the notion corresponding to velocity in this case should be naturally associated with the Fourier space as opposed to real space. Now, velocity itself is not a well-defined notion here; what we will instead look at is the inverse velocity, which I shall denote by [math]\boldsymbol\upsilon[/math] (I don't know if it's discernible, but that's supposed to be an upsilon as opposed to a [math]\boldsymbol v[/math] which LaTeX renders as somewhat thinner). This shouldn't be too alien; after all lattice spacings in the reciprocal lattice of a crystal (which is nothing but its Fourier space) are measured in units of inverse length. Now, let's talk about a wave which may be defined in any number of dimensions (say [math]N[/math]). We can express this in general as a complex valued function [math]\psi(t,\boldsymbol x)[/math]. This can be written in terms of real-valued functions [math]A(t,\boldsymbol x)[/math] and [math]\phi(t,\boldsymbol x)[/math] as [math]\psi = A e^{i\phi}[/math]. Of course, [math]A[/math] is the amplitude while [math]\phi[/math] is the phase. We think of the wave as being formed by two families of surfaces, one composed of surfaces of constant [math]A[/math] and the other of surfaces of constant [math]\phi[/math], billowing and propagating through space with time. We may thus assign two inverse velocity fields [math]\boldsymbol \upsilon_g(t,\boldsymbol x)[/math] and [math]\boldsymbol\upsilon_p(t,\boldsymbol x)[/math] which basically measure how slowly these surfaces are moving through space. More precisely, we have [math]\boldsymbol\upsilon_g \partial_t A + \nabla_{\boldsymbol x} A=0,[/math] [math]\boldsymbol\upsilon_p \partial_t \phi + \nabla_{\boldsymbol x} \phi=0.[/math] As the subscripts may indicate, I'm eventually going to be identifying these inverse velocities as the 'reciprocals' of the group and phase velocities respectively. Of course, to be able to distill two specific vectors out of this, we need to first stipulate that the inverse velocity fields are uniform vector fields. So, setting [math]\boldsymbol\upsilon_g[/math] and [math]\boldsymbol\upsilon_p[/math] as constant in the two differential equations above yields the following general solution: [math]A = A(t-\boldsymbol\upsilon_g\cdot \boldsymbol x),[/math] [math]\phi = \phi(t-\boldsymbol\upsilon_p\cdot \boldsymbol x).[/math] From now onwards we shall thus treat [math]A[/math] and [math]\phi[/math] as single variable functions, so that [math]A'[/math] and [math]\phi'[/math] denote [math]A'(t-\boldsymbol\upsilon_g\cdot \boldsymbol x)[/math] and [math]\phi'(t-\boldsymbol\upsilon_p\cdot \boldsymbol x)[/math] respectively. Also, we'll use the shorthands [math]\phi_0=\phi(0)[/math] and [math]\phi'_0=\phi'(0)[/math]. We know from quantum mechanics that given that the functions are square integrable or alternatively, that we integrate over some finite region and subtract off boundary terms, the expectation value for [math]\omega[/math] and [math]\boldsymbol k[/math] are given by [math]\langle \omega\rangle= -i\int\mathrm d^N\boldsymbol x\,\psi^\star \partial_t\psi,[/math] [math]\langle \boldsymbol k\rangle= i\int\mathrm d^N\boldsymbol x\,\psi^\star \nabla_{\boldsymbol x}\psi.[/math] You may now readily verify that the inverse phase velocity is given by [math]\boldsymbol \upsilon_p = \langle \boldsymbol k\rangle/\langle \omega\rangle[/math]. To obtain an expression for [math]\boldsymbol \upsilon_g[/math] in terms of [math]\omega[/math] and [math]\boldsymbol k[/math], we'll have to make the assumption that [math]\phi=\phi_0+\phi'_0(t-\boldsymbol\upsilon_p\cdot \boldsymbol x).[/math] This essentially has to do with the fact that we require the wave to be of narrow bandwidth. We also note that [math]A[/math] may be Fourier decomposed as [math]A = \int\mathrm d\omega\, \hat A(\omega) \exp[i\omega(t-\boldsymbol\upsilon_g\cdot \boldsymbol x)].[/math] Plugging these expressions for [math]A[/math] and [math]\phi[/math] back into [math]\psi = A e^{i\phi}[/math] gives us [math]\psi = e^{i\phi_0}\int\mathrm d\omega\, \hat A(\omega) \exp[i(\omega + \phi'_0)t[/math] [math]- \,i(\omega\boldsymbol \upsilon_g + \phi'_0\boldsymbol\upsilon_p)\cdot \boldsymbol x].[/math] Furthermore, carrying out a change of variables [math]\omega \mapsto \omega-\phi'_0[/math] yields [math]\psi = e^{i\phi_0}\int  \mathrm d\omega\, \hat A(\omega-\phi'_0) \exp[i\omega t-i\boldsymbol k\cdot \boldsymbol x],[/math] where [math]\boldsymbol k = \omega\boldsymbol  \upsilon_g + \phi'_0(\boldsymbol\upsilon_p-\boldsymbol\upsilon_g)[/math]. We thus obtain a sort of inverse dispersion relation that immediately gives us the result [math]\boldsymbol \upsilon_g = d\boldsymbol k/d\omega[/math]. It may be readily seen that when [math]N=1[/math], the inverse phase and group velocities are actually reciprocals of the phase and group velocity respectively, just as expected. A well-defined notion of taking the reciprocal of a vector doesn't carry over into higher dimensions (unless you specify an entire basis of vectors), but for what it's worth, [math]\boldsymbol \upsilon_p\cdot \boldsymbol v_p=\boldsymbol \upsilon_g\cdot \boldsymbol v_g=1[/math].

I don't really think you can do a lot better than the derivation: https://en.wikipedia.org/wiki/Group_velocity#Derivation