I'd do it like so:

1. For arbitrary H and epsilon sufficiently close to 0,

2. Taking a derivative of both sides with respect to epsilon,

.

3. We can evaluate this equation at epsilon = 0 to get

where the LHS is conveniently the directional derivative of f:

.

4. Since this last equation holds for arbitrary H, we have the desired

Of course, there are some technical details left to fill in: for instance in step 1 we assumed that the set of invertible matrices is open, and in step 3, that matrix inversion is smooth (follows from Cramer's Rule).