Exponential

Overview

This module allows you to compute the exponential of a matrix.

f(A) = e^A.

We will take advantage of the following relationship to scale the matrix:

e^A = (e^{\frac{A}{\sigma}})^{\sigma}.

Method

The matrix exponential is computed using Chebyshev polynomials. First, we need to scale the matrix so that the spectrum of the matrix is contained in [-1,1]. Then we compute the exponential, and rescale it by successive squaring.

We also take advantage of polynomial factoring to speed up the computation of the Chebyshev polynomials. This can all be hardcoded (as in the below example) because we only need 32 coefficients to accurately compute the exponential.

matrix_dimension = InputMatrix.shape[0]
min_val, max_val = Helper.EigenCircle(InputMatrix)
spectral_radius = max_val
sigma_val = 1
sigma_counter = 1
print spectral_radius
while spectral_radius/sigma_val > 1.0:
  sigma_val = sigma_val*2
  sigma_counter = sigma_counter + 1
print sigma_val, sigma_counter
ScaledMatrix = InputMatrix/sigma_val

x = numpy.linspace(-1,1,1000)
y = [math.exp(i) for i in x]
coefficient_list = numpy.polynomial.chebyshev.chebfit(x,y,32)
numpy.set_printoptions(precision=15)
print coefficient_list

# Basic
T0 = scipy.sparse.identity(matrix_dimension)
T1 = ScaledMatrix.copy()
T2 = 2.0*ScaledMatrix.dot(T1) - T0
T3 = 2.0*ScaledMatrix.dot(T2) - T1

# First Half
T4 = 2.0*ScaledMatrix.dot(T3) - T2
T5 = 2.0*ScaledMatrix.dot(T4) - T3
T6 = 2.0*ScaledMatrix.dot(T5) - T4
T7 = 2.0*ScaledMatrix.dot(T6) - T5
T8 = 2.0*ScaledMatrix.dot(T7) - T6

# Compute the second half
OutputMatrix = coefficient_list[16]*T8
OutputMatrix = OutputMatrix + coefficient_list[15]*T7
OutputMatrix = OutputMatrix + coefficient_list[14]*T6
OutputMatrix = OutputMatrix + coefficient_list[13]*T5
OutputMatrix = OutputMatrix + coefficient_list[12]*T4
OutputMatrix = OutputMatrix + coefficient_list[11]*T3
OutputMatrix = OutputMatrix + coefficient_list[10]*T2
OutputMatrix = OutputMatrix + coefficient_list[9]*T1
OutputMatrix = 0.5*T8.dot(OutputMatrix)

# Add in the first half
OutputMatrix2 = (coefficient_list[8])*T8
OutputMatrix2 = OutputMatrix2 + (coefficient_list[7] + coefficient_list[9]/2.0)*T7
OutputMatrix2 = OutputMatrix2 + (coefficient_list[6] + coefficient_list[10]/2.0)*T6
OutputMatrix2 = OutputMatrix2 + (coefficient_list[5] + coefficient_list[11]/2.0)*T5
OutputMatrix2 = OutputMatrix2 + (coefficient_list[4] + coefficient_list[12]/2.0)*T4
OutputMatrix2 = OutputMatrix2 + (coefficient_list[3] + coefficient_list[13]/2.0)*T3
OutputMatrix2 = OutputMatrix2 + (coefficient_list[2] + coefficient_list[14]/2.0)*T2
OutputMatrix2 = OutputMatrix2 + (coefficient_list[1] + coefficient_list[15]/2.0)*T1
OutputMatrix2 = OutputMatrix2 + (coefficient_list[0] + coefficient_list[16]/2.0)*T0
OutputMatrix = OutputMatrix + OutputMatrix2

OutputMatrix2 = Chebyshev.Compute(ScaledMatrix,coefficient_list)
print scipy.sparse.linalg.norm(OutputMatrix2 - OutputMatrix)

identity_mat = scipy.sparse.identity(matrix_dimension)
for i in range(0,sigma_counter-1):
  OutputMatrix = OutputMatrix.dot(OutputMatrix)
return OutputMatrix

How to Cite

From what I could tell, people usually cite the "19 dubious ways" paper when they use this method.

@article{moler1978nineteen,
  title={Nineteen dubious ways to compute the exponential of a matrix},
  author={Moler, Cleve and Van Loan, Charles},
  journal={SIAM review},
  volume={20},
  number={4},
  pages={801--836},
  year={1978},
  publisher={SIAM}
}