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.