# Exponential Derivative Trick

The matrix exponential can appear in a variety of computational problems. Unfortunately the standard finite difference routine to obtain its derivative is often inaccurate. Using complex numbers we can do better!

## Complex Step

From Ref. [1]
For normal finite difference schemes we can approximate a derivative with

$F^\prime(x) = \frac{F(x+h)-F(x)}{h}+O(h^2)$

with $$h\ll 1$$ . When $$F(x+h)\approx F(x)$$ we run into subtraction errors. A useful trick to mitigate this, when dealing with analytic functions, is to use the complex step which is

$\boxed{ F^\prime (x) = \frac{\text{Im}\left[F(x+ih)\right]}{h}} +O(h^2)$

### Short proof:

Expand $$F(x+ih)$$ about $$x$$ using a taylor expansion

$F(x+ih) = F(x)+ihF^\prime(x) +O(h^2)$

Asking what the imaginary part is to first order yields the derivative. Note that $$\text{Re}\left[F(x+ih)\right]=F(x)$$, so we can obtain the function at a point and its derivative in one calculation.

## Exponential Complex Step

From Ref. [2]
The same trick can be applied for matrices, particularly the exponential matrix function $$\exp(A)$$. Consider the Fréchet derivative (really a directional derivative) of a function $$f(A)=\exp(A)$$. The derivative of $$f$$ can be approximated as

$\frac{df}{dA_{ij}}=\frac{\text{Im}\left[ \exp(A+ihE_{ij}) \right]}{h}$

where $$(E_{ij})_{kl}=\delta_{ik}\delta_{jl}$$ is the single-entry matrix. Note: Other matrices can be used for $$E$$.

## Python Example

Here is an example for getting the derivative at the identity. The derivative is analytic here, and is simply $$\delta_{ij}$$

import scipy.linalg as spsl

A = np.zeros([4,4],dtype=np.complex)
#derivative at i,j
i = 1
j = 2
h = 1e-6

Ah = A.copy()
Ah[i,j]+=1j*h
mA = spsl.expm(Ah)
print("exp(A)=\n",np.real(mA),"\n d/dA_ij exp(A)=\n",np.imag(mA)/h)
#normal finite difference
Ah =A.copy()
Ah[i,j]+=h
finiteDiff = np.real(spsl.expm(Ah)-spsl.expm(A))/h
print("Finite differences:\n",finiteDiff)


Output:

exp(A)=
[[ 1.  0.  0.  0.]
[ 0.  1.  0.  0.]
[ 0.  0.  1.  0.]
[ 0.  0.  0.  1.]]
d/dA_ij exp(A)=
[[ 0.  0.  0.  0.]
[ 0.  0.  1.  0.]
[ 0.  0.  0.  0.]
[ 0.  0.  0.  0.]]
Finite differences:
[[ 0.  0.  0.  0.]
[ 0.  0.  1.  0.]
[ 0.  0.  0.  0.]
[ 0.  0.  0.  0.]]


While the finite difference scheme here worked, it required two matrix exponential calculations rather than one; in terms of operations the complex step method is better for larger matrices. See Ref [2] for cases in which the complex step has a lower relative error than finite difference schemes.

# References

[1] Squire, William, and Trapp, George, Using complex variables to estimate derivatives of real functions, SIAM Review 40, 1998, pp. 110-112. DOI:10.1137/S003614459631241X

[2] Al-Mohy, Awad H., and Higham, Nicholas J., The complex step approximation to the Fréchet derivative of a matrix function, Numer Algor (2010) 53:133–148 DOI:10.1007/s11075-009-9323-y

Tags:

Categories:

Updated: