Derivative of the inverse of a symmetric matrix w.r.t itself I'm trying take the derivative of a symmetrix matrix $\mathbf{C}$ with respect to itself.
$$
\begin{equation}
\frac{\partial \mathbf{C}^{-1}}{\partial \mathbf{C}}
\end{equation}
$$
Using the indicial notation, above equation can be rewritten as follows
$$
\begin{equation}
\frac{\partial C_{ij}^{-1}}{\partial C_{kl}}
\end{equation}
$$
At first I've used the following formula,
$$
\begin{equation}
\frac{\partial C_{ij}^{-1}}{\partial C_{kl}} = -C^{-1}_{ik}C^{-1}_{lj}
\end{equation}
$$
But I quickly realized that we've lost the symmetry of the problem now.
I read The Matrix Cookbook and the other posts about the same problem but unfortunately, I couldn't understand the things they've done.
For example in this article, at Eq.(100) authors have used the property below when taking the derivative of Eq.(99)
$$
\begin{equation}
\frac{\partial \mathbf{C}^{-1}}{\partial \mathbf{C}} = -\mathbf{C}^{-1} \boxtimes \mathbf{C}^{-T} \mathbf{I}_s
\end{equation}
$$
Where $\boxtimes$ is the square product, $\mathbf{I}_s$ is the symmetric fourth-order identity tensor and they are defined as follows
$$
\begin{align}
(\mathbf{A} \boxtimes \mathbf{B})_{ijkl} &= \mathbf{A}_{ik}\mathbf{B}_{jl} \\
(\mathbf{I}_s)_{ijkl} &= \frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})
\end{align}
$$
I couldn't understand how did they achieve this result and how can I derive it myself.

In the single-variable case, we have that
$$\dfrac{d}{dt}C(t)^{-1}=-C(t)^{-1}\dfrac{dC(t)}{dt}C(t)^{-1}.$$
This can obtained by differentiating the expression $C(t)C(t)^{-1}=I$ on both sides with some simple algebra. This directly generalizes to the multivariable case by expressing
$$C=\sum_{i,j}c_{ij}e_ie_j^T.$$
Then, we have that
$$\dfrac{d}{dc_{kl}}C^{-1}=-C^{-1}\dfrac{dC}{dc_{kl}}C^{-1}=-C^{-1}e_ke_l^TC^{-1},$$
from which we get
$$\dfrac{d}{dc_{kl}}C_{ij}^{-1}=-e_i^TC^{-1}e_ke_l^TC^{-1}e_j=-C_{ik}^{-1}C_{lj}^{-1}.$$
When $C$ is symmetric, then it can be written as
$$C=\sum_{i}c_{ii}e_ie_i^T+\sum_{i>j}c_{ij}(e_ie_j^T+e_je_i^T).$$
$$\begin{array}{rcl}
\dfrac{d}{dc_{kl}}C^{-1}&=&-C^{-1}\dfrac{dC}{dc_{kl}}C^{-1}=-C^{-1}(e_ke_l^T+e_le_k^T)C^{-1},\ \mathrm{for}\ k\ne l\\
\dfrac{d}{dc_{kk}}C^{-1}&=&-C^{-1}\dfrac{dC}{dc_{kk}}C^{-1}=-C^{-1}e_ke_k^TC^{-1}
\end{array}$$
then we have that
$$\begin{array}{rcl}
\dfrac{d}{dc_{kl}}C_{ij}^{-1}&=&-e_i^TC^{-1}(e_ke_l^T+e_le_k^T)C^{-1}e_j=-C_{ik}^{-1}C_{lj}^{-1}-C_{il}^{-1}C_{kj}^{-1},\ \mathrm{for}\ k\ne l\\
\dfrac{d}{dc_{kk}}C_{ij}^{-1}&=&-e_i^TC^{-1}e_ke_k^TC^{-1}e_j=-C_{ik}^{-1}C_{kj}^{-1}.\end{array}$$
