Derivative of Matrix Power with resepect to entries Let's consider a matrix $A = \mathbb R^{d\times d}$. I'm interested in the entry wise derivative of $A^n$ that is if
$$B = A^n$$
I'd like to find
$$ c_{ij} := \frac{\partial}{\partial a_{ij}} b_{ij}.$$
If we now consider $C = (c_{ij})_{ij}$ as a matrix, is there a simple way to express $c_{ij}$ or could we even express it in terms of the matrix $A$?
First we obviously have $c_{ij} = 0$ for $n=0$ and $c_{ij} = 1$ for $n=1$, but after that I found it hard to come up with a general formula.
 A: $
\def\bbR#1{{\mathbb R}^{#1}}
\def\LR#1{\left(#1\right)}
\def\p{\partial}
\def\grad#1#2{\frac{\p #1}{\p #2}}
\def\trace#1{\operatorname{Tr}\LR{#1}}
\def\gradLR#1#2{\LR{\grad{#1}{#2}}}
$Let $\{e_k\}$ denote the standard vector basis for $\bbR{d}$ and use them to construct the standard matrix basis for $\bbR{d\times d}$
$$\eqalign{
E_{jk} &= e_j\,e_k^T = E_{kj}^T \\
}$$
We'll also need the Frobenius product, which is a concise notation for the trace
$$\eqalign{
A:B &= \sum_{i=1}^d\sum_{j=1}^d A_{ij}B_{ij} \;=\; \trace{A^TB} \\
A:A &= \|A\|^2_F \\
}$$
The properties of the underlying trace function allow the terms in a
Frobenius product to be rearranged in many different ways, e.g.
$$\eqalign{
A:B &= B:A \\
A:B &= A^T:B^T \\
C:\LR{AB} &= \LR{CB^T}:A \;=\; \LR{A^TC}:B \\
}$$
The standard matrix basis is orthonormal with respect to the Frobenius product
$$E_{ij}:E_{k\ell} = \delta_{ik}\,\delta_{j\ell}$$
and this provides a convenient way to extract the $(i,j)$ component of a matrix
$$\eqalign{
A_{ij} &= E_{ij}:A \\
}$$

Okay, now we are ready to try the $n=2$ case
$$\eqalign{
B &= A^2 \\
dB &= dA\,A + A\,dA  \\
\grad{B}{A_{ij}} &= E_{ij}\,A + A\,E_{ij}  \\
\grad{B_{ij}}{A_{ij}} &= E_{ij}:\gradLR{B}{A_{ij}} \\
 &= E_{ij}:E_{ij}\,A \;+\; E_{ij}:A\,E_{ij} \\
 &= \LR{E_{ji}\,E_{ij} + E_{ij}E_{ji}}:A \\
 &= \LR{E_{jj} + E_{ii}}:A \\
}$$
Let's try $n=3$
$$\eqalign{
B &= A^3 \\
dB &= dA\,A^2 + A\,dA\,A + A^2\,dA \\
\grad{B}{A_{ij}} &= E_{ij}\,A^2 + A\,E_{ij}\,A + A^2\,E_{ij} \\
\grad{B_{ij}}{A_{ij}} &= E_{ij}:\gradLR{B}{A_{ij}} \\
 &= E_{ij}:E_{ij}\,A^2 \;+\; E_{ij}:A\,E_{ij}\,A \;+\; E_{ij}:A^2\,E_{ij} \\
 &= E_{ji}E_{ij}:A^2 \;+\; E_{ij}:AE_{ij}A \;+\; E_{ij}E_{ji}:A^2 \\
 &= \LR{E_{ii}+E_{jj}}:A^2 \;+\; E_{ij}:AE_{ij}A \\\\
}$$
So the pattern that I see emerging is
$$\eqalign{
B &= A^n \\
\grad{B_{ij}}{A_{ij}}
 &=  E_{ij}:\LR{\sum_{k=0}^{n-1} A^{k}E_{ij}A^{n-k-1}} \\
 &= \LR{E_{ii}+E_{jj}}:A^{(n-1)}
  \;+\; E_{ij}:\LR{\sum_{k=1}^{n-2} A^{k}E_{ij}A^{n-k-1}} \\
}$$
where one could interpret the term
$\LR{A^kE_{ij}A^{n-k-1} = A^ke_i\;e_j^TA^{n-k-1}}$
as the matrix product of the $i^{th}$ column of $A^k$ and the $j^{th}$ row of $A^{n-k-1}$
