I understand the matrix multiplication behind Grover's algorithm, but I'd like to get an intuitive grasp on why sequence of gates Hadamard-Phase-Hadamard does inversion about the mean. Can anyone help?

I'm not sure whether it's completely intuitive if it still has some formulas in it, but here's a try. (I prefer term "reflection", since it makes the geometrical interpretation a bit simpler, but I think they are used interchangeably.)

1. First, let's convince ourselves that Controlled Z does inversion about the $$|1...1\rangle$$ state.
Controlled Z flips the phase of $$|1...1\rangle$$ basis state and keeps the phases of other states unchanged, so it can be written as $$\mathcal{I} - 2|1...1\rangle \langle 1...1|$$, which is the reflection about the $$|1...1\rangle$$ (with an extra global phase of -1, which can be ignored).

2. Second, let's see how to represent reflection about any state $$|\psi\rangle$$ in terms of reflection about the $$|1...1\rangle$$ state, given that we know a unitary $$U$$ that prepares state $$|\psi\rangle$$ from the $$|1...1\rangle$$, i.e., $$|\psi\rangle = U|1...1\rangle$$.
Reflection about $$|\psi\rangle$$ is

$$2|\psi\rangle \langle \psi| - \mathcal{I} = 2 U|1...1\rangle \langle 1...1|U^\dagger - \mathcal{I} = U(2|1...1\rangle \langle 1...1| - \mathcal{I})U^\dagger$$

So you can represent that reflection by applying the following sequence of steps:

• $$U^\dagger$$
• reflection about the $$|1...1\rangle$$
• $$U$$
3. Finally, let's see how to prepare the mean state $$\sum_j |j\rangle$$, starting with the $$|1...1\rangle$$ state.
The easiest way is to apply X gate to each qubit to get $$|0...0\rangle$$, and then to apply H gate to each qubit to get an equal superposition of all basis states.

Putting this all together, we get the following procedure:

• apply H to each qubit
• apply X to each qubit
• apply Controlled Z with most qubits as control and last qubit as target
• apply X to each qubit
• apply H to each qubit