Derivative of the matrix exponential with respect to its matrix argument I was trying to find the Frechet derivative of $f = \exp(X)$, where $X \in \mathbb{R}^{n\times n}$ is positive definite. I thought it ought to be $\exp(X)$.

I see results where the derivative is with respect to a scalar argument, but this question has not been asked before.

I tried to see if I could find $Df_X$ starting with 
$Df_X[H] = \exp(X+H) - \exp(X)$.  If I can show that the right hand side evaluates to $I + XH + X^2H/2 + \ldots = \exp(X)H$, I am done.

After I use the series definition, however,  I am lost  because I see no reason to  assume that $A$ and $H$ commute.

Please help.

EDIT
Following the suggestion in the comment, I try to compute the Gateaux derivative as $\exp(X + tH)$ by writing down the first few terms.

$\exp(X+tH) = I + (X + tH) + (X^2 + tXH + tHX + t^2H^2)/2 + \cdots$

$\dfrac{d}{dt}\exp(X+tH)\Big|_{t=0} = H + (XH+HX)/2 + \cdots$

And now am stuck again. It seems the expression on the right cannot be rearranged to give what I want.
I think it  is the derivative of the trace of the exponential, not the exponential itself that yields $\exp(X)$
Let $f : \mathbb R^{n \times n} \to \mathbb R^{n \times n}$ be defined by $f (X) = \exp(X)$. Hence,

$\begin{array}{rl} f (X + h V) &amp;= \exp(X + h V)\\ &amp;= I_n + (X + h V) + \frac{1}{2!} (X + h V)^2 + \frac{1}{3!} (X + h V)^3 + \cdots\\ &amp;= I_n + X + h V + \frac{1}{2!} X^2 + \frac{h}{2!} (X V + V X) + \frac{h^2}{2!} V^2 + \frac{1}{3!} X^3 + \\ &amp;\,\,\,\,\,+ \frac{h}{3!} (X^2 V + X V X + V X^2) + \frac{h^2}{3!} (X V^2 + V X V + V^2 X) + \frac{h^3}{3!} V^3 + \cdots\\ &amp;= f (X) + h \left( V + \frac{1}{2!} (X V + V X) + \frac{1}{3!} (X^2 V + X V X + V X^2) + \cdots\right) + \cdots\end{array}$

Thus, the directional derivative of $f$ in the direction of $V$ at $X$ is given by

$$\begin{array}{rl} D_V f (X) &amp;= \displaystyle\lim_{h \to 0} \frac{1}{h} \left( f (X + h V) - f (X) \right)\\ &amp;= V + \frac{1}{2!} (X V + V X) + \frac{1}{3!} (X^2 V + X V X + V X^2) + \cdots\end{array}$$

We then write

$$D_V f (X) = M_0 + \frac{1}{2!} M_1 + \frac{1}{3!} M_2 + \frac{1}{4!} M_3 + \cdots$$

where

$$\begin{array}{rl} M_0 &amp;= V\\ M_1 &amp;= X V + V X =: \{X,V\}\\ M_2 &amp;= X^2 V + X V X + V X^2\\ &amp;\vdots\\ M_k &amp;= \displaystyle\sum_{i=0}^k X^{k-i} V X^i\end{array}$$
The formula given by @ user251257 is closed and not so ugly, at least when $X$ is diagonalizable and we know the (eventually numerical) decomposition  $X=PDP^{-1}$ where $D=diag((\lambda_i))$. Let $f(X)=e^X$; then $Df_X(H)=\int_0^1 e^{\alpha X}He^{(1-\alpha)X} d\alpha=\int_0^1 B(\alpha) d\alpha$ is a $n\times n$ matrix; since $B(\alpha)= Pe^{\alpha D}P^{-1}HPe^De^{-\alpha D}P^{-1}$, its $(k,l)$ entry $b_{k,l}$ is in the form $\sum_{ij}u_{ij}e^{\alpha(\lambda_i-\lambda_j)}$ where $u_{ij}$ is a linear function of the $(h_{ij})$. Of course, for every $(k,l)$,  $\int_0^1 b_{kl} d \alpha$ is easy to calculate (from a numerical point of view).

Example. Let $X=\begin{pmatrix} 5.6&amp; 0.4&amp; 0.2\\ 0.4&amp; 8.3&amp; -8.5\\ 0.2&amp; -8.5&amp; 13.3\end{pmatrix}$ (a random symmetric $&gt;0$ matrix) and $H=E_{1,1}$. Then 

$b_{1,1}(t)\approx (.99302845389679037258*(.97923083022132240582*exp(5.6517287919679228980*t)+0.013765826093233907222*exp(1.8877987424462300666*t)+0.000031797582234059493609*exp(19.660472465585847035*t)))*exp(5.6517287919679228980-5.6517287919679228980*t)+(.11773898646801403107*(.11610306332719986462*exp(5.6517287919679228980*t)+0.0016321530422941489699*exp(1.8877987424462300666*t)+0.0000037700985200174380332*exp(19.660472465585847035*t)))*exp(1.8877987424462300666-1.8877987424462300666*t)-(0.005658693915992128895*(-0.0055800692513703492171*exp(5.6517287919679228980*t)-0.000078443468620371001349*exp(1.8877987424462300666*t)-1.8119600140866963913*10^{-7}*exp(19.660472465585847035*t)))*exp(19.660472465585847035-19.660472465585847035*t)$

and $\int_0^1b_{1,1}(t)dt\approx 1854.0328800470185607$.

Let  $\Delta=e^{A+0.001E_{1,1}}-e^A$; then $\Delta_{1,1}\approx  1.854$.