Notebook

Topic 21: The Matrix Exponential¶

Author: Thane Downing thanedowning@gmail.com

Introduction¶

The matrix exponential is a useful tool in solving systems of linear differential equations, helping to describe the behavior of linear systems, as well as having several other intersting uses. The matrix exponential works much like the scalar exponential and has many of the same properties. While perhaps initially unfamiliar, the matrix exponential developes readily from very familiar concepts and contributes to a straightforward understanding of how matrix math fits into a larger framework.

Explanation of the theory¶

The eponential of any square matrix A can be writen as $e^A$ , and can be understood in a similar way that we understand the scalar equivalent exponential $e^x$ . From calculus, we know that

$\frac{d}{dx}e^x = xe^x$ The same is true of the matrix exponential

$\frac{d}{dx}e^A = Ae^A$ The scalar exponential has the property

$(e^x)^{-1} = e^{-x}$ which follows directly from algebra. The matrix exponential shares this property such that

$(e^A)^{-1} = e^{-A}$ . The Taylor series can be use to define

$e^x$ .

$e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \ldots + \frac{n^2}{n!} + \ldots$ The matrix exponential can also be defined this way.

$e^A = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \ldots + \frac{A^2}{n!} + \ldots$ In fact, this is the most general way to define the matrix exponential, and it works for any square matrix. Therefore, in general we will define the matrix exponential as

$e^A = \sum_{n=0}^{\infty} \frac{1}{n!}A^n$

One classic property of scalar exponentials that doesn't translate generally to matrix exponentials is the property

$e^ae^b = e^{a+b}$ For the matrix exponential case, in general

$e^Ae^B \neq e^{A+B}$ This is because

$A$ and

$B$ must be able to commute in order for these to be equal, i.e.

$AB = BA$ . This can be seen using the definition given above;

$e^Ae^B$ can only be rearranged into

$e^{A+B}$ (by multiplying to instances of the Taylor series definition above) if, in the following product

$(A+B)(A+B) = A^2 + AB + BA + B^2$ , the sum of

$AB + BA$ can be writen

$2AB$ or

$2BA$ . In general this is not the case, therefore in general

$e^Ae^B \neq e^{A+B}$ .

An expanded form to represent the matrix exponential, that only works for diagonalizable matrices, is given as $e^A = X e^{\Lambda} X^{-1}$ , where $\Lambda$ is the diagonal matrix of the eigenvalues of $A$ and $X$ is an invertable matrix made up of the eigenvectors of $A$ . Taking this one step futher we can find that

$e^\Lambda = \begin{pmatrix} e^{\lambda_{1}} & & & \\ & e^{\lambda_{2}} & & \\ & & \ddots & \\ & & & e^{\lambda_{n}} \end{pmatrix}$ which makes computation quite straightforward.

The matrix exponential is handy in solving systems of linear differential equations, such as those encountered in continuous time state space equations. Given some linear ordinary differential equation $\frac{dx}{dt} = ax$ , we can determine that the solution is $x(t) = e^{at}$ . Now suppose we have a system of linear differential equations such that $\frac{d\textbf{x}}{dt} = A\textbf{x}$ . The solution to this system is $\textbf{x}(t) = e^{At}\textbf{x}(0)$ . Using the expanded form given above, this solution can be rewriten $\textbf{x}(t) = e^{At}\textbf{x}(0) = Xe^{\Lambda t}X^{-1}\textbf{x}(0)$ . We can use the eigenvalues of A to determine the behavior of the solution similar to the way we can determine system dynamics in other settings. The solution will decay if Re $(\lambda_i) < 0$ , it will just oscillate if Re $(\lambda_i) = 0$ ( $\lambda_i$ is purley imaginary), and it blow up if Re $(\lambda_i) > 0$ , where $\lambda_i$ is the eigen value associated with a given eigenvector in the system.

Occasionally, this solution using the matrix exponential will be represented in terms of the Laplace transform as in

$e^{At} = \mathcal{L}^{-1}(sI - A)^{-1}$ However, using the Neumann expansion

$(sI - A)^{-1} = \frac{1}{s}(I + A/s + A^2/s^2 + \ldots)$ , then taking the inverse Laplace transform produces the same Taylor series form of the matrix exponential as that given above.

In the context of continuous time LTI system defined by the state space equations

$\dot{\textbf{x}} = A\textbf{x}(t) + B\textbf{f}(t)\\ \textbf{y}(t) = C\textbf{x}(t)+D\textbf{f}(t)$ the matrix exponential is used in the solution

$\textbf{x}(t) = e^{At}\textbf{x}(0) + \int_{0}^{t} e^{A(t-\lambda)}B\textbf{f}(t)d\lambda$ If instead of starting at

$t = 0$ , we start at some arbitrary time

$\tau$ , then the solution to the state at

$t$ will be of the form

$\textbf{x}(t) = e^{A(t-\tau)}\textbf{x}(\tau) + \int_{\tau}^{t} e^{A(t-\tau)}B\textbf{f}(t)d\lambda$ The matrix

$e^A(t-\tau)$ can be seen as a method of moving from some arbitrary state

$\textbf{x}(\tau)$ to another

$\textbf{x}(t)$ , and is therefore sometimes refered to as the propagation matrix or the state-transition matrix denoted

$\Phi(t,\tau)$ . This matrix has a couple interesting propeties:

$\begin{matrix} \Phi(t,t) = I \\ \frac{\partial \Phi(t,\tau)}{\partial t} = A \Phi(t,\tau) \\ \Phi(t,\tau) = (\Phi(\tau,t))^{-1} \\ \end{matrix}$

Although there are contrived examples where the solution to the matrix exponential is nice and tidy, in general the solution will be rather complex. For instance, the general solution to a $e^A$ , where $A$ is a 2x2 matrix

$A = \begin{pmatrix} a & b\\ c & d \end{pmatrix}$ is given by

$\begin{matrix} e^A = \frac{1}{\Delta} \begin{pmatrix} e^{\frac{a+d}{2}}(\Delta cosh(\frac{1}{2}\Delta)+(a-d)sinh(\frac{1}{2}\Delta)) & 2b e^{\frac{a+d}{2}}sinh(\frac{1}{2}\Delta)\\ 2c e^{\frac{a+d}{2}}sinh(\frac{1}{2}\Delta) & e^{\frac{a+d}{2}}(\Delta cosh(\frac{1}{2}\Delta)+(d-a)sinh(\frac{1}{2}\Delta)) \end{pmatrix} \\ \Delta = \sqrt{(a-d)^2+4bc} \end{matrix}$

Simple Numerical Examples¶

There are several examples we can use to verify the results of the preceeding section. First, we generate a random square matrix and use the matrix exonential function from the scipy.linalg library to calculate the matrix exponential.

In [1]:

import numpy as np
from scipy.linalg import expm
from math import factorial, sinh, cosh, exp, sqrt, cos, ceil
import random
import matplotlib.pyplot as plt

In [2]:

A = np.array([[random.random(),random.random(),random.random(),random.random()],
              [random.random(),random.random(),random.random(),random.random()],
              [random.random(),random.random(),random.random(),random.random()],
              [random.random(),random.random(),random.random(),random.random()]])
print "A =\n", A
print "e^A =\n", expm(A)

A =
[[ 0.41777503  0.53103038  0.04252436  0.41216531]
 [ 0.69812543  0.79625546  0.76295725  0.41021156]
 [ 0.53157639  0.82114634  0.18914389  0.5079495 ]
 [ 0.16940129  0.18442986  0.58999039  0.67698445]]
e^A =
[[ 2.11211077  1.38598808  0.72110255  1.16504324]
 [ 2.09831048  3.57898164  1.9539366   1.83430856]
 [ 1.59212573  2.1362369   2.19697908  1.61725546]
 [ 0.88873006  1.12036052  1.30709399  2.58273717]]

Next, we can test our Taylor series expansion method (making sure to use enough iterations to get a good approximation) to verify that this approach will give us the same thing as the library function does.

In [3]:

taylor_estimate = np.eye(4) + A
add = A
for i in range(2,200):
    add = np.matmul(add,A)/i
    taylor_estimate = taylor_estimate + add
print taylor_estimate

[[ 2.11211077  1.38598808  0.72110255  1.16504324]
 [ 2.09831048  3.57898164  1.9539366   1.83430856]
 [ 1.59212573  2.1362369   2.19697908  1.61725546]
 [ 0.88873006  1.12036052  1.30709399  2.58273717]]

We can also test the expanded form of the matrix exponential by using the scipy.linalg library to find the eigenvalues and eigenvectors of $A$ to verify that this method also produces what we expect to see.

In [4]:

lam, X = np.linalg.eig(A)
eigen_decomp = np.matmul(np.matmul(X,expm(np.diag(lam))),np.linalg.inv(X))
print eigen_decomp

[[ 2.11211077  1.38598808  0.72110255  1.16504324]
 [ 2.09831048  3.57898164  1.9539366   1.83430856]
 [ 1.59212573  2.1362369   2.19697908  1.61725546]
 [ 0.88873006  1.12036052  1.30709399  2.58273717]]

We can check the basic property $(e^A)^{-1} = e^{-A}$

In [5]:

print "e^-A =\n", expm(-A)
print "(e^A)^-1 =\n", np.linalg.inv(expm(A))

e^-A =
[[ 0.76307929 -0.33076492  0.17330506 -0.21782085]
 [-0.30979663  0.71708154 -0.50396241 -0.05396847]
 [-0.25085094 -0.497817    1.11972926 -0.23443581]
 [-0.0012397   0.05469565 -0.40770436  0.60419532]]
(e^A)^-1 =
[[ 0.76307929 -0.33076492  0.17330506 -0.21782085]
 [-0.30979663  0.71708154 -0.50396241 -0.05396847]
 [-0.25085094 -0.497817    1.11972926 -0.23443581]
 [-0.0012397   0.05469565 -0.40770436  0.60419532]]

A quick check will show that indeed, in general $e^Ae^B \neq e^{A+B}$

In [6]:

B = np.array([[random.random(),random.random(),random.random(),random.random()],
              [random.random(),random.random(),random.random(),random.random()],
              [random.random(),random.random(),random.random(),random.random()],
              [random.random(),random.random(),random.random(),random.random()]])
C = expm(A)*expm(B)
D = expm(A+B)
print "(e^A)*(e^B) =\n",C 
print "e^(A+B) =\n",D

(e^A)*(e^B) =
[[  2.89313406   0.6620096    0.70024447   0.35380883]
 [  4.2540861   12.28011182   5.30642275   3.5365327 ]
 [  2.72511048   4.51514568   6.84983879   2.1878806 ]
 [  1.759763     2.07105209   2.90330633   5.14254724]]
e^(A+B) =
[[  7.49914959   7.82533417   7.16019721   5.87058551]
 [ 17.40289883  21.92329833  19.12141699  15.75151267]
 [ 14.0086614   17.02401286  15.89451394  12.54361878]
 [ 12.57003962  14.49902883  13.63391563  12.01335549]]

We can even test the formula for computing a 2x2 matrix $\begin{pmatrix} a & b\\ c & d \end{pmatrix}$ using the terms a, b, c, and d.

In [7]:

a = random.random()
b = random.random()
c = random.random()
d = random.random()

H = expm(np.array([[a,b],[c,d]]))
G = (1/sqrt((a-d)**2+4*b*c))*np.array([
    [exp((a+d)/2)*(sqrt((a-d)**2+4*b*c)*cosh(.5*sqrt((a-d)**2+4*b*c))+(a-d)*sinh(.5*sqrt((a-d)**2+4*b*c))),
     2*b*exp(((a+d)/2))*sinh(.5*sqrt((a-d)**2+4*b*c))],
    [2*c*exp(((a+d)/2))*sinh(.5*sqrt((a-d)**2+4*b*c)),
     exp(((a+d)/2))*(sqrt((a-d)**2+4*b*c)*cosh(.5*sqrt((a-d)**2+4*b*c))+(d-a)*sinh(.5*sqrt((a-d)**2+4*b*c)))]])
print "H =\n", H
print "G =\n", G

H =
[[ 1.76853665  1.72718204]
 [ 1.81716025  3.22362205]]
G =
[[ 1.76853665  1.72718204]
 [ 1.81716025  3.22362205]]

Engineering Application¶

One of the most applicable uses of the matrix exponential is in propagating the states of a system described by state space equations. We can choose from thousands of examples of physical systems that have well defined state space representations to use as the basis for this exercise. The following code is flexible enough to implement any LTI constant coefficient system. Only the A, B, C, and D matrices need to be changed (along with the corresponding dimensions of u, x, and y).

The following example models an RLC circuit with two capacitors in parallel, where the input is an AC voltage source with an amplitude of five volts, the intermidiary states are the voltages over the two capcitors, and the output is the current through the first resistor, as shown in the picture below. RLC_Circuit

In [8]:

c1 = 10*10**-6 # In Farads
c2 = 20*10**-6 # In Farads
r1 = 200 # In Ohms
r2 = 100 # In Ohms

A = np.array([[-1/(c1*r1)-1/(c1*r2),1/(c1*r2)],
              [1/(c2*r2),-1/(c2*r2)]])
B = np.array([[1/(c1*r1)],[0]])
C = np.array([[-1/r1,0]])
D = np.array([[1/r1]])

x_0 = np.array([[0],[4]])

t_step = 0.001
t_start = 0
t_end = 15

bins = int(ceil((t_end-t_start)/t_step))

u_vals = np.zeros((bins,1))
x_vals = np.zeros((bins,2))
y_vals = np.zeros((bins,1))


P = expm(A*t_step) # This is the propagation matrix

x_val = x_0

for i in range(bins):
    u_val = 5*cos(i*t_step)
    y_val = np.matmul(C,x_val)+D*u_val
    u_vals[i] = u_val
    x_vals[i,0] = x_val[0]
    x_vals[i,1] = x_val[1]
    y_vals[i] = y_val
    x_val = np.matmul(P,x_val)+np.matmul(P,B)*u_val*t_step 

In the following plots, the red line shows the voltage of the input signal u(t) (AC voltage source), the blue and black dashed lines show the voltage on capacitors 1 and 2 respectively, and the green line shows the current through resistor 1 in amps. The x axis is the number of time samples.

In [9]:

plt.figure(1)
plt.title("System Dynamics of an RLC Circuit")

plt.subplot(211)
plt.plot(u_vals, 'r-', x_vals[:,0], 'b--', x_vals[:,1])
plt.title("System Dynamics of an RLC Circuit")
plt.ylabel('Voltage (in Volts)')
plt.axis([0, 15000, -5.5, 5.5])

plt.subplot(212)
plt.plot( y_vals, 'g-')
plt.ylabel('Current (in Amps)')
plt.axis([0, 15000, -5.5, 5.5])

plt.show

Out[9]:

<function matplotlib.pyplot.show>

Sources¶

"Mathematical Methods and Algorithms for Signal Processing" by Todd K. Moon and Wynn C. Stirling

http://mathworld.wolfram.com/MatrixExponential.html

http://web.mit.edu/18.06/www/Spring17/Matrix-Exponentials.pdf

http://web.mit.edu/16.unified/www/FALL/signalssystems/Lecture11_12.pdf