# 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.

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)$$

• Could you please explain more on how you get to equation $(3)$? More specifically, I don't understand how a variation of $i,j$ entry in $A$ will only lead a variation of $i,j$ entry in $X$? Also I am not quite familiar with the vectorization operation. Could you point some references? Thanks. May 4, 2018 at 16:53
• @jing007 For the derivation of (3), I have replaced in (2) $\Delta A=(\Delta a_{ij})e_ie_j^T$ since only a variation of the $(i,j)$ element of $A$ is considered. The coresponding variation in $X$ is not only in the $(i,j)$ element of $X$. $\Delta_{i,j} X$ denotes a matrix (the variation in $X$) that is associated with the variation of the $(i,j)$ element of $A$ not a specific matrix entry. The necessary properties for the vec operation and kronecker products can be found in ime.unicamp.br/~cnaber/Kronecker.pdf
– RTJ
May 4, 2018 at 17:18
• Thanks. Could you also explain why $I - A^T \otimes A^T$ is invertible? May 4, 2018 at 17:35
• @jing007 The following general property is true: If $\lambda_i$ are the eigenvalues of $A$ and $\mu_j$ are the eigenvalues of $B$ then the eigenvalues of $A\otimes B$ are $\lambda_i\mu_j$ (see en.wikipedia.org/wiki/Kronecker_product ). Thus, $\lambda(I-A^T\otimes A^T)=1-\lambda (A^T\otimes A^T)=1-\lambda_i\lambda_j\neq 0$ if all eigenvalues of $A$ are lying on the interior of the unit circle. Thus, matrix $I-A^T\otimes A^T$ has no zero eigenvalue and is therefore invertible.
– RTJ
May 4, 2018 at 17:48
• Thanks again for your clarification. Do you think the mapping $h : A \mapsto X$ is continuous? May 4, 2018 at 18:09

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.

• I assume you mean $\mathrm{vec}(X) = (I - A^\top \otimes A)^{-1} \mathrm{vec}(Q)$? May 4, 2018 at 6:20
• @KwinvanderVeen Oops. Right, thanks. May 4, 2018 at 6:22
• You could use a vec expression for $0$-$1$ vector $$g(A)={\rm vec}(I_n)^T(I_n\otimes I_n-A\otimes A)^{-T}{\rm vec}(Q)$$
– greg
May 4, 2018 at 13:21
• @greg Yes I was about to fix it but RTJ's answer already figured that part so I left it as is. May 4, 2018 at 13:28

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.

• Thanks for your answer. Your derivation is very nice. But I could not make sense of your taking differential step. I know intuitively or loosely you can do that. What is the exact math meaning? May 5, 2018 at 20:28

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 }

• How do we make sense of the taking 'differential' step? Is it differential form? May 5, 2018 at 19:42
• @jing007 If you're not familiar with differentials, you can treat it similarly to taking the derivative with respect to a scalar, e.g. time. The product rule for differentials (or time derivatives) is straightforward $$d(ABC)=dA\,BC+A\,dB\,C+AB\,dC$$ since $(A,B,C)$ are matrices, and therefore do not commute, their relative order cannot be changed. The differential of a sum is also very simple $$d(A+B)=dA+dB$$ Finally, the differential (time derivative) of a constant is obviously zero.
– greg
May 5, 2018 at 21:23
• Thanks for your clarification. But unfortunately this still does not make much sense to me. I know by 'hand-waving' this is perfectly fine. But I would like some rigor on taking this 'differential'. The only way to make sense of it on my head is to put this equation into some manifold space. Yet, I have not identified them. May 5, 2018 at 21:34
• More specifically, I am a little uncomfortable of freely using the symbols "$dX, dA$" before assigning some mathematical meaning to them. May 5, 2018 at 21:50