The derivative of matrix inversion

Idea. Reduce to inverting a matrix close to $I$, where the Neumann series gives an explicit expansion.

Step 1: factorisation. For $t$ small, $A+tH=A(I+tA^{-1}H)$. We have $\|tA^{-1}H\|\le|t|\|A^{-1}\|\|H\|<1$ for $|t|<\frac1{\|A^{-1}\|\|H\|}$.

Step 2: Neumann series. For $\|X\|<1$, $(I+X)^{-1}=\sum_{k\ge0}(-X)^k=I-X+X^2(I+X)^{-1}$, hence
$$\|(I+X)^{-1}-I+X\|\le\frac{\|X\|^2}{1-\|X\|}.$$
(Indeed $(I+X)^{-1}-I+X=X^2(I+X)^{-1}$, as one checks by multiplying by $I+X$.) With $X=tA^{-1}H$ this is $O(t^2)$.

Step 3: expansion.
$$(A+tH)^{-1}=(I+tA^{-1}H)^{-1}A^{-1}=\big(I-tA^{-1}H+O(t^2)\big)A^{-1}=A^{-1}-tA^{-1}HA^{-1}+O(t^2).$$

Step 4: conclusion. So $\frac{d}{dt}(A+tH)^{-1}|_{t=0}=-A^{-1}HA^{-1}$. Since $H\mapsto-A^{-1}HA^{-1}$ is linear and the remainder is $O(\|H\|^2)$ uniformly, this is the Fréchet derivative.

Consistency check. Differentiating $AA^{-1}=I$ gives $HA^{-1}+A\,D=0$, i.e. $D=-A^{-1}HA^{-1}$. For $n=1$: $\frac{d}{da}\frac1a=-\frac1{a^2}$.

Remark. Matrices do not commute, so one must keep the order $A^{-1}\,H\,A^{-1}$; for $n\ge2$ it is *not* $-HA^{-2}$ in general. The map $\iota$ is $C^\infty$ (even analytic) on $GL_n$: iterating, the $k$-th derivative involves products with $k+1$ copies of $A^{-1}$.