# Matrix derivative

I am looking to compute the derivative of the following expression: $$\frac{\partial}{\partial X}\mathrm{tr}\left[A\exp(X)\right]$$ where $$A$$ is both a symmetrical and positive-definite matrix and $$X$$ is a symmetrical matrix. As a result, $$\exp(X)$$ is symmetrical and definite-positive.

Just as a reminder, the matrix exponential $$\exp(X)$$ is defined as: $$\exp(X)=I+X+\frac12X^2+\dots$$ so if we diagonalize it ($$X=UDU^T$$), then we have that $$\exp(X)=U\exp(D)U^T$$.

I have been looking at this for hours and I am unable to calculate a simple solution as soon as the matrix is not a trivial $$1\times1$$ matrix.

** Edit **

There are a few answers and hints that relate to $$\mathrm{tr}(A^TB)=\sum A\odot B$$, but this doesn't work in general for non-linear function, as far as I understand. For example: $$\frac{\partial}{\partial X}\mathrm{tr}\left[A^TX^2\right]=AX^T+X^TA$$ for generic square matrices $$A$$ and $$X$$. As these two matrices do not commute in general, we can't simplify the expression.

** Edit 2 **

I actually found that the expression is fairly more complicated than I expected. Part of the answer is described here: https://mathoverflow.net/a/403404/468380

** Final edit **

I believe that the answer from Lutz Lehmann is the closest to what it is possible to do given my original question. Basically, this is a very challenging problem, with no closed form solution. Fortunately, I was able to rewrite my original problem using matrix logarithms and in this case, I do get a commuting derivative which simplifies greatly the work.

• The trace of the matrix product $tr(A^TB) = A ⊙ B$, so don't worry about $A$. The trace is additive in the sense that $tr(I+X + \frac{1}{2} X^2 ...) = tr(I) + tr(X) + \frac{1}{2}tr(X^2)...$. I suggest you start by deriving $\frac{∂}{∂X} tr(X^k) = k X^{k-1}$. Commented Nov 29, 2022 at 20:43
• @CharlieS, is $\odot$ the Frobenius inner product? Commented Nov 29, 2022 at 20:53
• In any case, both the trace and multiplication with $A$ are linear operations, so applying the chain rule shouldn't be too difficult. Commented Nov 29, 2022 at 21:30
• OP: please consider writing an answer to your own question to detail the solution that you have found. Commented Nov 30, 2022 at 7:35
• SIde note: the derivative of the matrix exponential is described in more detail in Higham's book Functions of matrices. Commented Nov 30, 2022 at 7:35

Per the very helpful wikipedia page "derivative of the exponential map", based around the "Baker-Campbell-Hausdorff formula", the relevant formula for a directional derivative is $$\newcommand{\ad}{\operatorname{ad}}$$ $$\frac{d}{dt}e^{X+tH}\Big|_{t=0}=e^X\phi_1(-\ad_X)[H]$$ where $$\phi_1(z)=\frac{e^z-1}z=1+\frac{z}2+\frac{z^2}{3!}+…+\frac{z^n}{(n+1)!}+…$$ and $$\ad_X[Y]=[X,Y]$$, so that $$(\ad_X)^2[Y]=[X,[X,Y]]$$ etc.

Now including the trace we get $$\newcommand{\Tr}{\operatorname{Tr}} \Tr(A·[X,Y])=\Tr(AXY-AYX)=-\Tr([X,A]·Y) \\ \Tr(A·[X,[X,Y]])=-\Tr([X,A]·[X,Y])=\Tr([X,[X,A]]·Y) \\ \Tr(Ae^X·\phi_1(-\ad_X)[H])=\Tr(\phi_1(+\ad_X)[Ae^X]·Y)$$

So that finally, how helpful it may ever be, $$\frac{\partial \Tr(Ae^X)}{\partial X}=\phi_1(\ad_X)[Ae^X]$$

I like the derivation where some very large $$N$$ is picked and $$\exp(X)$$ approximated as $$(I+\frac{X}{N})^N+O(N^{-2})$$. Then \begin{align} \exp(X+tH)&=\left(I+\frac{X}{N}+t\frac{H}{N}\right)^N+O(N^{-2}) \\ &=\left(I+\frac{X}{N}\right)^N +t\sum_{k=0}^{N-1}\left(I+\frac{X}{N}\right)^{N-k}\frac{H}{N}\left(I+\frac{X}{N}\right)^{k-1} +O(N^{-2},t^2) \end{align}

The term linear in $$t$$ now is approximately equal to $$\left(I+\frac{X}{N}\right)^N·\frac1N\sum_{k=0}^{N-1} \left(I-\frac{X}{N}\right)^kH\left(I+\frac{X}{N}\right)^{k} \\ =\left(I+\frac{X}{N}\right)^N·\frac1N\sum_{k=0}^{N-1} \left(I-\frac{\ad_X}{N}\right)^kH+O(N^{-2})$$ with $$\left(I-\frac{X}{N}\right)H\left(I+\frac{X}{N}\right)=H-\frac{[X,H]}N+O(N^{-2})$$. This expression can be further approximately recognized as Riemann sum for the integral $$e^X·\int_0^1\exp(-s\ad_X)[H]\,ds=e^X·\phi_1(-\ad_X)[H]$$ as claimed, as the integration now proceeds completely in the commutative operator algebra generated by $$\ad_X$$.

• This is indeed what I think is the correct approach in general. That being said, it's problematic as it involves an infinite sum that's not straightforward to compute numerically. I will check but I think that I can rewrite the problem as a perturbative one, so I can truncate the sum early (and it becomes tractable).
– PC1
Commented Nov 30, 2022 at 17:58
• I'm really not sure how helpful that eventually is. $\phi_1(M)$ can be evaluated with the same efficiency as $\exp(M)$, but notice that if $X$ is $n\times n$, then $ad_X$ acts on its argument as vector, that is, in a suitable basis $ad_X$ is a $n^2\times n^2$ matrix. Commented Nov 30, 2022 at 20:21
• This gets very complicated to compute in practice, I went back to my original problem and reformulated it as a matrix logarithm. In this form, I get something like $[\log(X),I]=0$ so it's much more tractable.
– PC1
Commented Dec 1, 2022 at 4:51

Observing that both the trace and the multiplication with a matrix $$A$$ are linear operations, it is easy to apply the chain rule. For this, first see that $$\frac{\partial}{\partial X_{ij}} [\text{tr} X ] = \frac{\partial}{\partial X_{ij}} \left[\sum_{l} X_{ll}\right] = \delta_{ij}.$$ From there, you can get the rest by applying the chain rule as long as you know what the derivative of $$e^X$$ is.

By you can also write it out via index notation: $$\frac{\partial}{\partial X_{ij}} [\text{tr} Ae^X ] = \frac{\partial}{\partial X_{ij}} \left[\sum_{l} (Ae^X)_{ll}\right] = \frac{\partial}{\partial X_{ij}} \left[\sum_{lk} A_{lk}(e^X)_{kl}\right] = \sum_{lk} A_{lk} \frac{\partial}{\partial X_{ij}} \left[(e^X)_{kl}\right].$$ This comes down to the same: As long as you know what the derivative of $$e^X$$ is, you can compute the derivative of the original expression.

$$\def\o{{\tt1}} \def\LR#1{\left(#1\right)} \def\op#1{\operatorname{#1}} \def\vc#1{\op{vec}\LR{#1}} \def\rs#1{\op{Unvec}\LR{#1}} \def\Diag#1{\op{Diag}\LR{#1}} \def\diag#1{\op{diag}\LR{#1}} \def\trace#1{\op{Tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\c#1{\color{red}{#1}} \def\CLR#1{\c{\LR{#1}}} \def\fracLR#1#2{\LR{\frac{#1}{#2}}} \def\Sk#1{\LR{\sum_{k=0}^\infty #1}} \def\Sj{\sum_{j=\o}^k}$$Use a colon as a convenient product notation for the trace, i.e. \eqalign{ A:B &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij} \;=\; \trace{A^TB} \\ A:A &= \|A\|^2_F \\ } The properties of the underlying trace function allow the terms in such a product to be rearranged in numerous ways, e.g. \eqalign{ A:B &= B:A \\ A:B &= A^T:B^T \\ \LR{AB}:C &= A:\LR{CB^T} \;=\; B:\LR{A^TC} \\ } Now we can use a power series expansion of the exponential to calculate the differential and gradient of your function \eqalign{ \phi &= A^T:e^{X} \\ &= A^T:\Sk{\frac{X^k}{k!}} \\ d\phi &= A^T:\Sk{\frac{\c{dX^k}}{k!}} \\ &= A^T:\Sk{\;\frac{\c{\Sj X^{k-j}\:dX\:X^{j-\o}}}{k!}} \\ &= \Sk{\Sj \frac{X^{j-\o}\:A\:X^{k-j}}{k!}}^T:dX \\ \grad{\phi}{X} &= \Sk{\Sj \frac{X^{k-j}\:A\:X^{j-\o}}{k!}} \\ } where the last line takes advantage of the fact that both of your matrices are symmetric.

## Daleckii-Krein

Since $$X$$ is symmetric it can be diagonalized which permits the use of the Daleckii-Krein Theorem \eqalign{ X &= QBQ^T,\quad B = \Diag{b},\quad Q^TQ=I \\ F &= f(X) \\ dF &= Q\LR{R\odot\LR{Q^TdX\,Q}}Q^T \\ d\phi &= A^T:dF \;=\; Q\LR{R\odot\LR{Q^TA^TQ}}Q^T:dX \\ \grad{\phi}{X} &= Q\LR{R\odot\LR{Q^TA^TQ}}Q^T \\ } where $$\odot$$ denotes the Hadamard product.

All we need is the symmetric $$R$$ matrix, which lies at the heart of the theorem \eqalign{ R_{ij} &= \begin{cases} {\Large\frac{f(b_i)\,-\,f(b_j)}{b_i\,-\,b_j}} \quad{\rm if}\; b_i\ne b_j \\ \\ \qquad f'(b_j) \qquad {\rm otherwise} \\ \end{cases} }

• The vectorization approach is quite interesting as I have other similar equations. Thank you for pointing that out.
– PC1
Commented Dec 8, 2022 at 17:47
• Unfortunately, the vectorization approach won't work because the $M$ matrix is singular. I think it can be salvaged by including a vector from the nullspace of $M$ and enforcing a symmetric constraint on $X$. While attempting that, I recalled the Daleckii-Krein theorem, which seems like an even better approach.
– greg
Commented Dec 9, 2022 at 6:28
• @PC1 The vectorization approach was removed from the post since it could not be salvaged. Daleckii-Krein seems like the simplest solution.
– greg
Commented Dec 13, 2022 at 19:51
• Thank you for the help, this is indeed the most appropriate approach I believe too.
– PC1
Commented Dec 13, 2022 at 21:23

This is an unpleasant derivative, except in special cases such as when $$A$$ is a power of $$X$$, or where you only care about the derivative of your function in direction of $$X$$.

You can brute-force the derivative using the Zassenhaus formula; another approach, assuming that the eigenvalues of $$X$$ are distinct, is to take the spectral decomposition $$X = VDV^T$$ and observe that $$\operatorname{tr}(A\exp(X)) = \operatorname{tr}(AV\exp(D)V^T) = \exp(D) : V^TAV = \sum_{i} e^{\lambda_i} v_i^T A v_i,$$ where $$(\lambda_i, v_i)$$ are $$X$$'s eigenpairs.

Now you can apply standard formulas for the derivative of eigenvectors and eigenvalues with respect to variations $$\delta X$$ in $$X$$:

$$d\left[\operatorname{tr}(A\exp(X))\right]\delta X = \sum_i\left[ e^{\lambda_i} (v_i^T\delta X v_i) (v_i^T A v_i) + 2e^{\lambda_i} \sum_{j\neq i} \left( \frac{1}{\lambda_i-\lambda_j}[v_i^T \delta X v_j] [v_i^T A v_j] \right)\right].$$

• That's what I also see, if $A$ and $X$ commute then the solution is straightforward to express, unfortunately this is not true in general. I can certainly assume that the eigenvalues of $X$ are all distinct, this is certainly a good way to approach the problem. I will check if I can use any structure on the eigenvalue distribution, maybe that expression can be computed approximately without having to consider all the terms.
– PC1
Commented Nov 30, 2022 at 17:43

What about Matrix calculus identities which states that
\begin{align} \text{Differentials}: \qquad & \text{d}(\operatorname {tr} (\mathbf {X} ))= \displaystyle \operatorname {tr} (\text{d}(\mathbf {X}) ), \\[1em] \text{Trace of \mathbf{X}}: \qquad & \frac {\partial \operatorname {tr} (\mathbf {X} )}{\partial \mathbf {X} }= \mathbf {I}, \\[1em] \text{Matrix multiply with \mathbf{X}}: \qquad & \frac {\partial \operatorname {tr} (\mathbf {AX} )}{\partial \mathbf {X} } = \frac {\partial \operatorname {tr} (\mathbf {XA} )}{\partial \mathbf {X}} \equiv \mathbf {A}, \\[1em] \text{Trace of exponential}: \qquad & \frac {\partial \operatorname {tr} \left(e^{\mathbf {X} }\right)}{\partial \mathbf{X}} = e^{\mathbf {X}}. \end{align}
\begin{align} \text{Expression in question}: \qquad & \frac {\partial \operatorname {tr} \left(\mathbf {A} e^{\mathbf {X} }\right)}{\partial \mathbf{X}} = ~... \text{work in progress} ... \end{align}
• We must have that $\mathrm{tr}(AB)=\mathrm{tr}(BA)$ for any (compatible) matrices $A$ and $B$. The last expression in your answer is incorrect. The main problem is that $A$ and $X$ do not commute in general.
• In other words, $\exp(X)A\ne\partial_X\mathrm{tr}(A\exp(X))\ne A\exp(X)$ in general.
• After working on this for some time, I don't believe that simple matrix identities can be used to solve this problem. It gets very intricate as soon as $[A,X]\ne0$ or the dimension of the square matrices is larger than $1\times1$.