Frechet derivative of a composition of functions over matrices In control theory, the discrete Lyapunov equation is defined as
\begin{align*}
A^T X A + Q  = X,
\end{align*} 
where $A \in \mathcal{M}(n \times n; \mathbb R)$ and $Q \in \mathbb {S}_{++}$ ( positive definite matrices). There is a theorem stating if the spectral radius of $A$ satisfies $\rho(A) < 1$ and for fixed $Q > 0$, there exists a unique $X \in \mathbb {S}_{++}$ which solves above equation.
Let $D = \{A \in \mathcal{M}(n \times n; \mathbb R): \rho(A) < 1\}$ and fix $Q$. Suppose we define some scalar valued function $f$ over $X$ which are solutions of Lyapunov equation over $D$. To make it more concrete, let us define this scalar valued function to be $f(X) = \text{tr}(X)$. This function can be also viewed as a function $g$ over $D$, i.e., it is a composition
\begin{align*}
g \colon A \xrightarrow{h} X \xrightarrow{f} \text{tr}(X).
\end{align*}
 Now I would like to differentiate $g$ with respect to $A$. Is it possible to find an explict formula for this Frechet derivative? The difficulty is the first function $h$ is not explicitly defined. Another question to ask is whether this $h : A \mapsto X$ is continuous.
 A: It would be really messy but you can utilize the "naive"(in the numerical sense) solution of the Lyapunov equation which is $\textrm{vec}(X) = (I-A^T\otimes A)^{-1}\textrm{vec}(Q)$ and trace condition is a row vector of $1$s and $0$s whose $1$s hitting every diagonal element on $\textrm{vec}(X)$.
Hence the explicit (again theoretical) expression for $g:D\to \mathbb{R}$ is 
$$
g(A) = \begin{bmatrix}1&0&\cdots&0&\color{red}{0}&\color{red}{1}&\color{red}{0}&\cdots&\color{red}{0}&\color{blue}{0}&\color{blue}{0}&\color{blue}{1}&\cdots\end{bmatrix}(I-A^T\otimes A)^{-1}\textrm{vec}(Q)
$$
with colors, trying to encode the entries multiplying each row group, resembling the log det problems.
A: Assume a small variation $\Delta A$ of the elements of $A$. Then, for the the new solution $X+\Delta X$ of the Lyapunov equation we have
$$(A+\Delta A)^T(X+\Delta X)(A+\Delta A)+Q=X+\Delta X \qquad \qquad(1)$$
Taking into account the unperturbed equation $A^TXA+Q=X$ for small variations $\Delta A$ (we ignore second order terms) we obtain
$$(\Delta A)^T X A+A^T X(\Delta A)=\Delta X-A^T(\Delta X)A\qquad \qquad (2)$$
Consider a variation $\Delta a_{ij}$ of the $(i,j)$-element in $A$. Then, this variation will induce a variation $\Delta_{i,j} X$ (this is a slight abuse of notation to differentiate on the effects of the different element variations) on $X$  that should satisfy
$$\Delta a_{ij}(e_j e_i^T X A+A^T Xe_ie_j^T)=\Delta_{i,j} X-A^T(\Delta_{i,j} X)A\qquad \qquad(3)$$
where $e_i$ is the $i$-th column of the identity matrix. Since $\Delta [tr(X)]=tr(\Delta X)$ the desired matrix $$S=\frac{\partial [tr(X)]}{\partial A}$$ 
will have elements given by

$$S_{ij}=\lim_{\Delta a_{ij}\rightarrow 0}\frac{tr(\Delta_{i,j}X)}{\Delta a_{ij}}$$

Applying the vec operator in (3) we obtain
$$vec(\Delta_{i,j}X)=(\mathbb{I}-A^T\otimes A^T)^{-1}vec(A^TXe_ie_j^T+e_je_i^TXA)\Delta a_{ij}$$
For the trace we have
$$tr(\Delta_{i,j}X)=vec^T(\mathbb{I})vec(\Delta_{i,j}X)=vec^T(\mathbb{I})(\mathbb{I}-A^T\otimes A^T)^{-1}vec(A^TXe_ie_j^T+e_je_i^TXA)\Delta a_{ij}$$
and therefore 

$$S_{ij}=vec^T(\mathbb{I})(\mathbb{I}-A^T\otimes A^T)^{-1}vec(A^TXe_ie_j^T+e_je_i^TXA)$$

A: Define the variables
$$\eqalign{
 M &= (I\otimes I-A\otimes A) \in {\mathbb R}^{n^2\times n^2} \cr
 x &= {\rm vec}(X),\,\,q={\rm vec}(Q),\,\,\,y = {\rm vec}(I)\,\in {\mathbb R}^{n^2} \cr
}$$
Then we can rearrange and vectorized the Lyapunov equation 
$$\eqalign{
 Q &= X - A^TXA \cr
 q &= M^Tx \cr
}$$
Taking the differential yields the relationship between $dx$ and $dA$
$$\eqalign{
M^Tdx &= -dM^Tx \cr
dx &= M^{-T}(dA\otimes A+A\otimes dA)^Tx \cr
}$$
The function we are actually interested in is $$\phi={\rm tr}(X)=I:X$$ where the colon denotes the trace/Frobenius product, i.e. $\,\,A:B\equiv{\rm tr}(A^TB)$.
Take the differential of this function
$$\eqalign{
d\phi &= I:dX = y:dx = y^T:dx^T \cr
 &= y^T:x^T(dA\otimes A+A\otimes dA)M^{-1} \cr
 &= xy^TM^{-T}:(dA\otimes A+A\otimes dA) \cr
}$$
Now we need to decompose the LHS of the product into a sum of Kronecker factors
$$\eqalign{
xy^TM^{-1} &= \sum_{k=1}^r B_k\otimes C_k \cr
B_k,C_k &\in {\mathbb R}^{n\times n}
}$$
We also need to know the rule for a Kronecker-Frobenius mixed product
$$(A\otimes B\otimes C):(X\otimes Y\otimes Z)=(A:X)\,(B:Y)\,(C:Z)$$
Substitute the Kronecker factorization into the differential to obtain our final result
$$\eqalign{
d\phi
 &= \sum_{k=1}^r B_k\otimes C_k:(dA\otimes A+A\otimes dA) \cr
 &= \bigg(\sum_{k=1}^r (A:B_k)C_k + (A:C_k)B_k\bigg):dA \cr\cr
S &= \frac{\partial\,{\rm tr}(X)}{\partial A} \cr
 &= \sum_{k=1}^r (A:B_k)C_k + (A:C_k)B_k \cr
 &= \sum_{k=1}^r {\rm tr}(A^TB_k)C_k + {\rm tr}(A^TC_k)B_k \cr\cr
}$$
For more information about the Kronecker product factorization, look for papers by Pitsianis and vanLoan. It turns out to be yet another (albeit clever) application of the SVD.
A: The Kronecker commutation matrix, $K,\,$ provides the most direct route to a solution. 
Using it we can write the differential of the Lyapunov equation, vectorize it, and solve for $dx$
$$\eqalign{
dX - A^TdXA &= dA^TXA+A^TXdA \cr
\big(I\otimes I-A\otimes A\big)^T\,dx &= \Big((A^TX^T\otimes I)K+(I\otimes A^TX)\Big)\,da \cr
 C^Tdx &= B\,da \cr
 dx &= C^{-T}B\,da \cr
}$$
The trace function can also be put into vector form and solved directly
$$\eqalign{
 {\rm tr}(X) &= {\rm vec}(I)^T{\rm vec}(X) = i^Tx \cr
 d\,{\rm tr}(X) &= i^Tdx = i^TC^{-T}B\,da = (B^TC^{-1}i)^Tda \cr
 \frac{\partial\,{\rm tr}(X)}{\partial a} &= B^TC^{-1}i \cr
 \frac{\partial\,{\rm tr}(X)}{\partial A} &= {\rm Mat}(B^TC^{-1}i) \cr
}$$
