How to find the derivative of the adjugate operator? How can I find the general derivative of the function $\mathbf{A} \to \mathrm{adj}(\mathbf{A})$ ?
Where $\mathbf{A}$ is invertible I got the following:
$\mathrm{D}f_{\mathrm{adj}(\mathbf{A})(H)} = \mathrm{tr}(\mathrm{adj}(\mathbf{A}) \cdot H)\cdot \mathrm{adj}(\mathbf{A})^{-1}   -   (\mathrm{adj}(\mathbf{A})) \cdot  H \cdot\mathrm{adj}(\mathbf{A})^{-1}$
Is it right? (And if it is, it's still not enough for every $\mathbf{A}$..)
 A: If $A$ is invertible, then $f(A)=adj(A)=\det(A)A^{-1}$. Then $Df_A:H\rightarrow tr(adj(A)H)A^{-1}-\det(A)A^{-1}HA^{-1}=tr(adj(A)H)A^{-1}-adj(A)HA^{-1}$.
EDIT. There was a mistake in this part.
If $A$ is not invertible, it is more complicated and I give only an idea. $H$ is given and $K=Df_A(H)$ is unknown. One has $A.adj(A)=adj(A).A=\det(A).I=0$ but, when $\det(A)=0$, these relations do not characterize $adj(A)$; moreover , the derivative of RHS is not $0$.
We obtain $(*)$ $H.adj(A)+AK=KA+adj(A).H=tr(adj(A)H).I$ but these relations do not characterize $K$. 
Note that, when $rank(A)\leq n-2$, $adj(A)=0$ and the previous relations reduce to $AK=KA=0$.
When $rank(A)=n-1$, $rank(adj(A))=1$ and the solutions $K$ of $(*)$ are in an affine space of dimension $1$, that is, $K=u+tv$ where $u,v$ are known vectors (dependent on $H$) and $t$ is an unknown scalar.
A: I will try two elementary approaches.
In the first approach, define $\delta:A\mapsto\det(A)$ and $f:A\mapsto\operatorname{adj}(A)$. Using Jacobi's formula $D\delta_A(H)=\operatorname{tr}(\operatorname{adj}(A)H)$, we immediately obtain
$$
\bbox[10px,border:2px solid red]{
Df_A(H)=\big((-1)^{i+j}\operatorname{tr}(\operatorname{adj}(A_{ji})H_{ji})\big)_{i,j\in\{1,2,\ldots,n\}}
}\tag{0}
$$
where $A_{ji}$ and $H_{ji}$ are respectively the submatrices of $A$ and $H$ obtained by deleting the $j$-th rows and the $i$-th columns.
In the second approach, we continue from the answer by user91684. If $A\in M_n(\mathbb C)$ is invertible and $B=\operatorname{adj}(A)$, since $f(A)=\delta(A)A^{-1}$, by the product rule we obtain
\begin{align}
Df_A(H)
&=\operatorname{tr}\left(BH\right)A^{-1}-BHA^{-1}\\
&=\det(A)\left(\operatorname{tr}\left(A^{-1}H\right)A^{-1}-A^{-1}HA^{-1}\right).\tag{1}
\end{align}
If $A$ is singular, consider the special case $A=S$ first, where
$$
S=\Sigma\oplus0_{k\times k}
$$
for some nonsingular diagonal matrix $\Sigma$. Here $k\ge1$ is the nullity of $A$. Let
\begin{aligned}
S'&=0\oplus I_k,\\
S_t&=\Sigma\oplus tI_k=S+tS'.\\
\end{aligned}
When $t\ne0$ is small, $S_t$ is invertible. Let $s=\det(\Sigma)$. Hence we can apply $(1)$ to obtain
\begin{align}
Df_{S_t}(K)
&=\det(S_t)\left(\operatorname{tr}(S_t^{-1}K)S_t^{-1}-S_t^{-1}KS_t^{-1}\right)\\
&=t^ks\left[\operatorname{tr}\left((S^++t^{-1}S')K\right)(S^++t^{-1}S')
-(S^++t^{-1}S')K(S^++t^{-1}S')\right].\tag{2}
\end{align}
In view of $(0)$, the value of $Df_A(H)$ is continuous in the elements of $A$. Therefore $Df_S(K)=\lim_{t\to0}Df_{S_t}(K)$. We consider three possibilities according to the nullity $k$ of $A$:

*

*When $k=1$, the expression in $(2)$ is in the form of $ts[M_0+t^{-1}M_1+t^{-2}M_2]$, where $M_0,M_1,M_2$ are some constant matrices. Since $\lim_{t\to0}Df_{S_t}(H)$ must exist, we must have $M_0=0$ and $\lim_{t\to0}Df_{S_t}(H)=sM_1$. Thus
$$
Df_S(K)=sM_1
=s\left(\operatorname{tr}(S'K)S^+ + \operatorname{tr}(S^+K)S' - S'KS^+ - S^+KS'\right).\tag{3}
$$

*When $k=2$, we obtain from $(2)$
$$
Df_S(K)=\lim_{t\to0}Df_{S_t}(K)
=s(\operatorname{tr}(S'K)S'-S'KS')
=s\pmatrix{0\\ &\operatorname{adj}(K')},\tag{4}
$$
where $K'$ denotes the trailing principal $2\times2$ submatrix of $K$.

*When $k\ge3$, the expression in $(2)$ is $t^{k-2}$ times a polynomial in $t$. Therefore
$$
Df_S(K)=\lim_{t\to0}Df_{S_t}(K)=0.\tag{5}
$$
This can be justified by considering the directional derivatives of $f_S$. Since the nullity of $S$ is at least $3$, the nullity of a rank-one update in the form of $S+hE_{ij}$ is at least $2$. Therefore $\operatorname{adj}(S+hE_{ij})$ is identically zero for all indices $i,j$ and scalar $h$. In other words, the directional derivative of $f$ at $S$ along every direction $E_{ij}$ vanishes. Since $f$ is also continuously differentiable, $Df_S$ must be zero.

Now, return to the general case. Since $A$ is singular, $A=USV^\ast$ for some diagonal matrix $S=\Sigma\oplus0$ and invertible matrices $U,V$. One can pick a singular value decomposition if one wants to, but other choices of $U,S,V$, such as those obtained from a rank decomposition, are also legitimate. At any rate, $A^g=(V^\ast)^{-1}S^+U^{-1}$ is a reflexive generalised inverse of $A$ in the sense that it satisfies the first two Penrose conditions $AA^gA=A$ and $A^gAA^g=A^g$.
For any matrix $H$, let $K=U^{-1}H(V^\ast)^{-1}$. Since
$$
\operatorname{adj}(A+H)
=\operatorname{adj}\left(U(S+K)V^\ast\right)
=\operatorname{adj}(V^\ast)\operatorname{adj}(S+K)\operatorname{adj}(U),
$$
we have
$$
Df_A(H)=\operatorname{adj}(V^\ast)Df_S(K)\operatorname{adj}(U).
$$
It follows from $(1)$ and $(3)-(5)$ that
$$
\bbox[10px,border:2px solid red]{
Df_A(H)=\begin{cases}
\operatorname{tr}\left(BH\right)A^{-1}-BHA^{-1}
&\text{ when } k=0,\\
\operatorname{tr}(BH)A^g + \operatorname{tr}(A^gH)B - BHA^g - A^gHB
&\text{ when } k=1,\\
\det(\Sigma)\det(UV^\ast)\,Y^\ast\operatorname{adj}\left(XHY^\ast\right)X
&\text{ when } k=2,\\
0&\text{ when } k\ge3,\\
\end{cases}}\tag{6}
$$
where $B=\operatorname{adj}(A)$ and $X,Y\in M_{2,n}(\mathbb C)$ are the submatrices taken from the last two rows of $U^{-1}$ and $V^{-1}$ respectively.
Computationally, if we choose a SVD $A=USV^\ast$, then $A^g$ is just $A^+$ and $X^\ast, Y^\ast$ are identical to the submatrices taken from the last two columns of $U$ and $V$ respectively.
One mystery remains. It isn't clear why the values of the expressions in $(6)$ for the cases $k=1$ or $k=2$ are independent of the choice of the decomposition $A=USV^\ast$, but such independence must be true because the value of $Df_A(H)$ is unique.
