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).