Solution (source code)

= Solution

Differentiability at $a$ means that there is a <linear map> $L:V\to V_1$ with
$$
\frac{\|f(a+h)-f(a)-Lh\|_1}{\|h\|}\longrightarrow0\quad(h\to0).
$$
This map is the <total derivative> $Df_a$. In the <Frobenius norm>, $|AB|\le|A||B|$, by the <Cauchy-Schwarz inequality> applied to <matrix> entries. Expanding with the order of all factors preserved gives
$$
(A+H)^3=A^3+A^2H+AHA+HA^2+AH^2+HAH+H^2A+H^3.
$$
The remainder after the terms linear in $H$ has norm at most $3|A||H|^2+|H|^3=o(|H|)$. Hence
$$
\boxed{Df_A(H)=A^2H+AHA+HA^2.}
$$
For $|C|<1$, the <Neumann series> $I-C+C^2-\cdots$ converges absolutely in the finite-dimensional space of <matrices>. Multiplying its partial sums by $I+C$ and passing to the limit shows that it is $(I+C)^{-1}$. Its remainder after $I-C$ has norm at most $|C|^2/(1-|C|)=o(|C|)$, proving the required expansion.

If $A$ is an <invertible matrix> and $|H|<1/|A^{-1}|$, then $A+H=A(I+A^{-1}H)$ is an <invertible matrix>. Thus the <general linear group> is an open subset of the <matrix> space. The preceding expansion also gives
$$
(A+H)^{-1}=A^{-1}-A^{-1}HA^{-1}+o(|H|).
$$
Apply the <chain rule> to $g(X)=(X^{-1})^3$, using the cubic <total derivative> at $A^{-1}$ and this <total derivative> of inversion:
$$
\boxed{Dg_A(H)=-A^{-3}HA^{-1}-A^{-2}HA^{-2}-A^{-1}HA^{-3}.}
$$
The three terms generally cannot be combined by commuting $A$ and $H$.