Matrix Exponential: Padé vs. Closed Form

time.py

import time
import numpy as np

from scipy.linalg import expm

def angvel2skewsym(w):
  """Converts an angular velocity to a skew-symmetric representation.
  """
  return np.array([
    [0, -w[2], w[1]],
    [w[2], 0, -w[0]],
    [-w[1], w[0], 0],
  ])

theta = 2.0
twist = np.asarray([0.1, -0.02, 0.0, 0.0, 0.09, 0.1])
v, w = twist[:3], twist[3:]

# ensure w is unit norm
w /= np.linalg.norm(w)
w_skew = angvel2skewsym(w)

times1 = []
for i in range(10000):
  tic = time.time()
  ori1 = np.eye(3) + (np.sin(theta)*w_skew) + (1 - np.cos(theta))*(w_skew @ w_skew)
  toc = time.time()
  times1.append(toc-tic)

times2 = []
for i in range(10000):
  tic = time.time()
  ori2 = expm(theta * w_skew)
  toc = time.time()
  times2.append(toc-tic)

print("closed form: {}s".format(np.mean(times1)))
print("pade approximation: {}s".format(np.mean(times2)))
print("speedup: {}".format(np.mean(times2) / np.mean(times1)))