Lie Groups and Lie Algebra — Visual Experiments¶

Goal: Build geometric intuition for SO(3), SE(3), exponential maps, and tangent-space reasoning through small rigid-motion experiments.

What we'll explore¶

#	Experiment	Question answered
1	Rotation order matters	Why is SO(3) non-commutative in practice?
2	Local linearization vs exact exponential	When does the Lie algebra approximation stop being trustworthy?
3	Geodesic interpolation vs Euler-angle lerp	Why is interpolation on the manifold different from coordinate-wise interpolation?
4	Twists generate rigid motion in SE(3)	How does a single algebra element produce a coupled rotation-translation motion?

---

Linked theory: ../research/05-lie-groups-and-lie-algebra.md

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm
from scipy.spatial.transform import Rotation

np.random.seed(42)
plt.rcParams.update({
    'figure.dpi': 110,
    'font.size': 11,
    'axes.titlesize': 12,
    'axes.labelsize': 11,
    'figure.facecolor': 'white',
    'axes.facecolor': '#f8f9fa',
    'axes.grid': True,
    'grid.alpha': 0.35,
})

import numpy as np import matplotlib.pyplot as plt from scipy.linalg import expm from scipy.spatial.transform import Rotation np.random.seed(42) plt.rcParams.update({ 'figure.dpi': 110, 'font.size': 11, 'axes.titlesize': 12, 'axes.labelsize': 11, 'figure.facecolor': 'white', 'axes.facecolor': '#f8f9fa', 'axes.grid': True, 'grid.alpha': 0.35, })

In [2]:

def hat(omega):
    wx, wy, wz = omega
    return np.array([
        [0.0, -wz, wy],
        [wz, 0.0, -wx],
        [-wy, wx, 0.0],
    ])


def so3_exp(rotvec):
    theta = np.linalg.norm(rotvec)
    if theta < 1e-10:
        return np.eye(3) + hat(rotvec)
    K = hat(rotvec / theta)
    return np.eye(3) + np.sin(theta) * K + (1.0 - np.cos(theta)) * (K @ K)


def so3_log(R):
    return Rotation.from_matrix(R).as_rotvec()


def se3_hat(xi):
    omega = xi[:3]
    v = xi[3:]
    mat = np.zeros((4, 4))
    mat[:3, :3] = hat(omega)
    mat[:3, 3] = v
    return mat


def se3_exp(xi, t=1.0):
    return expm(se3_hat(xi) * t)


def apply_rotation(R, points):
    return points @ R.T


def apply_transform(T, points):
    homog = np.hstack([points, np.ones((points.shape[0], 1))])
    moved = homog @ T.T
    return moved[:, :3]


def draw_frame(ax, R, title, origin=None):
    if origin is None:
        origin = np.zeros(3)
    colors = ['#d62728', '#2ca02c', '#1f77b4']
    labels = ['x-axis', 'y-axis', 'z-axis']
    for idx in range(3):
        axis = R[:, idx]
        ax.quiver(
            origin[0], origin[1], origin[2],
            axis[0], axis[1], axis[2],
            color=colors[idx], linewidth=2.5, arrow_length_ratio=0.12,
        )
        ax.text(*(origin + 1.12 * axis), labels[idx], color=colors[idx])
    ax.set_title(title)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-1.2, 1.2)
    ax.set_ylim(-1.2, 1.2)
    ax.set_zlim(-1.2, 1.2)


def set_equal_3d(ax, limit=1.2):
    ax.set_xlim(-limit, limit)
    ax.set_ylim(-limit, limit)
    ax.set_zlim(-limit, limit)

def hat(omega): wx, wy, wz = omega return np.array([ [0.0, -wz, wy], [wz, 0.0, -wx], [-wy, wx, 0.0], ]) def so3_exp(rotvec): theta = np.linalg.norm(rotvec) if theta < 1e-10: return np.eye(3) + hat(rotvec) K = hat(rotvec / theta) return np.eye(3) + np.sin(theta) * K + (1.0 - np.cos(theta)) * (K @ K) def so3_log(R): return Rotation.from_matrix(R).as_rotvec() def se3_hat(xi): omega = xi[:3] v = xi[3:] mat = np.zeros((4, 4)) mat[:3, :3] = hat(omega) mat[:3, 3] = v return mat def se3_exp(xi, t=1.0): return expm(se3_hat(xi) * t) def apply_rotation(R, points): return points @ R.T def apply_transform(T, points): homog = np.hstack([points, np.ones((points.shape[0], 1))]) moved = homog @ T.T return moved[:, :3] def draw_frame(ax, R, title, origin=None): if origin is None: origin = np.zeros(3) colors = ['#d62728', '#2ca02c', '#1f77b4'] labels = ['x-axis', 'y-axis', 'z-axis'] for idx in range(3): axis = R[:, idx] ax.quiver( origin[0], origin[1], origin[2], axis[0], axis[1], axis[2], color=colors[idx], linewidth=2.5, arrow_length_ratio=0.12, ) ax.text(*(origin + 1.12 * axis), labels[idx], color=colors[idx]) ax.set_title(title) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.set_xlim(-1.2, 1.2) ax.set_ylim(-1.2, 1.2) ax.set_zlim(-1.2, 1.2) def set_equal_3d(ax, limit=1.2): ax.set_xlim(-limit, limit) ax.set_ylim(-limit, limit) ax.set_zlim(-limit, limit)

---

Experiment 1: Rotation order matters¶

Question: Why can two 90-degree rotations give different answers depending on their order?

This is the fastest way to feel that SO(3) is a group but not a commutative one. Composition is the natural operation, and the order of composition matters.

In [3]:

np.random.seed(1)
Rx = so3_exp(np.array([np.pi / 2, 0.0, 0.0]))
Ry = so3_exp(np.array([0.0, np.pi / 2, 0.0]))
R_xy = Rx @ Ry
R_yx = Ry @ Rx

np.random.seed(1) Rx = so3_exp(np.array([np.pi / 2, 0.0, 0.0])) Ry = so3_exp(np.array([0.0, np.pi / 2, 0.0])) R_xy = Rx @ Ry R_yx = Ry @ Rx

In [4]:

fig = plt.figure(figsize=(12, 5))
fig.suptitle('Experiment 1: $R_x R_y$ and $R_y R_x$ are different frames', fontsize=13, fontweight='bold')

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
draw_frame(ax1, R_xy, 'Rotate about y, then x')
set_equal_3d(ax1)

ax2 = fig.add_subplot(1, 2, 2, projection='3d')
draw_frame(ax2, R_yx, 'Rotate about x, then y')
set_equal_3d(ax2)

plt.tight_layout()
plt.show()

fig = plt.figure(figsize=(12, 5)) fig.suptitle('Experiment 1: $R_x R_y$ and $R_y R_x$ are different frames', fontsize=13, fontweight='bold') ax1 = fig.add_subplot(1, 2, 1, projection='3d') draw_frame(ax1, R_xy, 'Rotate about y, then x') set_equal_3d(ax1) ax2 = fig.add_subplot(1, 2, 2, projection='3d') draw_frame(ax2, R_yx, 'Rotate about x, then y') set_equal_3d(ax2) plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The final coordinate frames are not the same, even though both use the same two 90-degree rotations.
2. Why it looks this way: Rotations in 3D do not commute. The second rotation always acts in the coordinate system produced by the first one.
3. Key takeaway: Treating rotations like ordinary vectors is dangerous because the natural operation is group composition, not vector addition.

---

Experiment 2: Local linearization vs exact exponential¶

Question: How far can we trust the first-order approximation R ≈ I + hat(omega)?

Lie algebra updates are linear and convenient, but they only approximate the real manifold locally. This experiment shows the approximation error as the rotation angle grows.

In [5]:

np.random.seed(2)
axis = np.array([1.0, 1.0, 0.5])
axis = axis / np.linalg.norm(axis)
angles = np.linspace(0.0, np.pi, 80)
test_vector = np.array([1.0, 0.2, -0.4])

orthogonality_errors = []
vector_errors = []
for theta in angles:
    rotvec = axis * theta
    R_exact = so3_exp(rotvec)
    R_linear = np.eye(3) + hat(rotvec)
    orthogonality_errors.append(np.linalg.norm(R_linear.T @ R_linear - np.eye(3), ord='fro'))
    vector_errors.append(np.linalg.norm(R_exact @ test_vector - R_linear @ test_vector))

np.random.seed(2) axis = np.array([1.0, 1.0, 0.5]) axis = axis / np.linalg.norm(axis) angles = np.linspace(0.0, np.pi, 80) test_vector = np.array([1.0, 0.2, -0.4]) orthogonality_errors = [] vector_errors = [] for theta in angles: rotvec = axis * theta R_exact = so3_exp(rotvec) R_linear = np.eye(3) + hat(rotvec) orthogonality_errors.append(np.linalg.norm(R_linear.T @ R_linear - np.eye(3), ord='fro')) vector_errors.append(np.linalg.norm(R_exact @ test_vector - R_linear @ test_vector))

In [6]:

fig, axes = plt.subplots(1, 2, figsize=(13.5, 4.8))
fig.suptitle('Experiment 2: Tangent-space linearization is only local', fontsize=13, fontweight='bold')

axes[0].plot(angles, orthogonality_errors, color='#d62728', linewidth=2)
axes[0].set_title('Orthogonality error of $I + \hat{\omega}$')
axes[0].set_xlabel('Rotation angle (radians)')
axes[0].set_ylabel('Frobenius norm error')

axes[1].plot(angles, vector_errors, color='#1f77b4', linewidth=2)
axes[1].set_title('Action error on a test vector')
axes[1].set_xlabel('Rotation angle (radians)')
axes[1].set_ylabel('Vector error norm')

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(13.5, 4.8)) fig.suptitle('Experiment 2: Tangent-space linearization is only local', fontsize=13, fontweight='bold') axes[0].plot(angles, orthogonality_errors, color='#d62728', linewidth=2) axes[0].set_title('Orthogonality error of $I + \hat{\omega}$') axes[0].set_xlabel('Rotation angle (radians)') axes[0].set_ylabel('Frobenius norm error') axes[1].plot(angles, vector_errors, color='#1f77b4', linewidth=2) axes[1].set_title('Action error on a test vector') axes[1].set_xlabel('Rotation angle (radians)') axes[1].set_ylabel('Vector error norm') plt.tight_layout() plt.show()

<>:5: SyntaxWarning: invalid escape sequence '\h'
<>:5: SyntaxWarning: invalid escape sequence '\h'
/tmp/ipykernel_2366/1322818632.py:5: SyntaxWarning: invalid escape sequence '\h'
  axes[0].set_title('Orthogonality error of $I + \hat{\omega}$')

What you see¶

1. What the plot shows: For tiny angles, the linear approximation behaves well. As the angle grows, it rapidly stops being orthogonal and its action on vectors becomes inaccurate.
2. Why it looks this way: The Lie algebra is a tangent space, not the manifold itself. First-order expansions are good near the identity but do not preserve the group constraints globally.
3. Key takeaway: Use algebra updates locally, then map back to the group with the exponential map instead of staying in the linear approximation.

---

Experiment 3: Geodesic interpolation vs Euler-angle lerp¶

Question: Why is interpolation on SO(3) different from interpolating angle coordinates directly?

A geodesic interpolation uses the logarithm and exponential maps. Euler-angle interpolation works in coordinates, not on the manifold itself.

In [7]:

np.random.seed(3)
R0 = np.eye(3)
R1 = so3_exp(np.array([1.2, -0.8, 0.9]))
final_euler = Rotation.from_matrix(R1).as_euler('xyz', degrees=False)
base_vector = np.array([1.0, 0.0, 0.0])
steps = np.linspace(0.0, 1.0, 50)

geodesic_points = []
euler_points = []
rotvec_total = so3_log(R1)
for t in steps:
    R_geo = so3_exp(t * rotvec_total)
    R_euler = Rotation.from_euler('xyz', t * final_euler).as_matrix()
    geodesic_points.append(R_geo @ base_vector)
    euler_points.append(R_euler @ base_vector)

geodesic_points = np.array(geodesic_points)
euler_points = np.array(euler_points)

np.random.seed(3) R0 = np.eye(3) R1 = so3_exp(np.array([1.2, -0.8, 0.9])) final_euler = Rotation.from_matrix(R1).as_euler('xyz', degrees=False) base_vector = np.array([1.0, 0.0, 0.0]) steps = np.linspace(0.0, 1.0, 50) geodesic_points = [] euler_points = [] rotvec_total = so3_log(R1) for t in steps: R_geo = so3_exp(t * rotvec_total) R_euler = Rotation.from_euler('xyz', t * final_euler).as_matrix() geodesic_points.append(R_geo @ base_vector) euler_points.append(R_euler @ base_vector) geodesic_points = np.array(geodesic_points) euler_points = np.array(euler_points)

In [8]:

u = np.linspace(0, 2 * np.pi, 40)
v = np.linspace(0, np.pi, 20)
xs = np.outer(np.cos(u), np.sin(v))
ys = np.outer(np.sin(u), np.sin(v))
zs = np.outer(np.ones_like(u), np.cos(v))

fig = plt.figure(figsize=(12.5, 5.4))
fig.suptitle('Experiment 3: The manifold path differs from coordinate-wise interpolation', fontsize=13, fontweight='bold')
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.plot_surface(xs, ys, zs, color='#dfe7fd', alpha=0.12, linewidth=0)
ax.plot(geodesic_points[:, 0], geodesic_points[:, 1], geodesic_points[:, 2], color='#2ca02c', linewidth=2.5, label='Geodesic via exp/log')
ax.plot(euler_points[:, 0], euler_points[:, 1], euler_points[:, 2], color='#d62728', linewidth=2.0, linestyle='--', label='Euler-angle lerp')
ax.scatter(*base_vector, color='black', s=40, label='Start vector')
ax.scatter(*geodesic_points[-1], color='#1f77b4', s=45, label='End vector')
ax.set_title('Path traced by a rotated unit vector on the sphere')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.legend(loc='upper left')
set_equal_3d(ax)
plt.tight_layout()
plt.show()

u = np.linspace(0, 2 * np.pi, 40) v = np.linspace(0, np.pi, 20) xs = np.outer(np.cos(u), np.sin(v)) ys = np.outer(np.sin(u), np.sin(v)) zs = np.outer(np.ones_like(u), np.cos(v)) fig = plt.figure(figsize=(12.5, 5.4)) fig.suptitle('Experiment 3: The manifold path differs from coordinate-wise interpolation', fontsize=13, fontweight='bold') ax = fig.add_subplot(1, 1, 1, projection='3d') ax.plot_surface(xs, ys, zs, color='#dfe7fd', alpha=0.12, linewidth=0) ax.plot(geodesic_points[:, 0], geodesic_points[:, 1], geodesic_points[:, 2], color='#2ca02c', linewidth=2.5, label='Geodesic via exp/log') ax.plot(euler_points[:, 0], euler_points[:, 1], euler_points[:, 2], color='#d62728', linewidth=2.0, linestyle='--', label='Euler-angle lerp') ax.scatter(*base_vector, color='black', s=40, label='Start vector') ax.scatter(*geodesic_points[-1], color='#1f77b4', s=45, label='End vector') ax.set_title('Path traced by a rotated unit vector on the sphere') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.legend(loc='upper left') set_equal_3d(ax) plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The green geodesic path and the red Euler-lerp path generally do not coincide on the unit sphere.
2. Why it looks this way: Geodesic interpolation uses a shortest-path construction in the group manifold, while Euler interpolation moves in a chosen coordinate chart and inherits that chart's distortions.
3. Key takeaway: Interpolation on rotations should respect the manifold geometry, especially when smoothness and shortest-path behavior matter.

---

Experiment 4: Twists generate rigid motion in SE(3)¶

Question: What does the exponential of a twist look like as an actual rigid-body motion?

A twist is one element of se(3). Exponentiating it produces a one-parameter motion in SE(3), coupling rotation and translation into a single group action.

In [9]:

np.random.seed(4)
body_points = np.array([
    [0.0, 0.0, 0.0],
    [0.6, 0.0, 0.0],
    [1.2, 0.0, 0.0],
    [1.2, 0.5, 0.0],
    [1.2, 1.0, 0.0],
    [0.0, 0.0, 0.4],
    [1.2, 0.0, 0.4],
    [1.2, 1.0, 0.4],
], dtype=float)

xi = np.array([0.0, 0.0, 1.1, 0.8, 0.0, 0.35])
times = [0.0, 0.4, 0.8, 1.2]
se3_snapshots = [apply_transform(se3_exp(xi, t=t), body_points) for t in times]

np.random.seed(4) body_points = np.array([ [0.0, 0.0, 0.0], [0.6, 0.0, 0.0], [1.2, 0.0, 0.0], [1.2, 0.5, 0.0], [1.2, 1.0, 0.0], [0.0, 0.0, 0.4], [1.2, 0.0, 0.4], [1.2, 1.0, 0.4], ], dtype=float) xi = np.array([0.0, 0.0, 1.1, 0.8, 0.0, 0.35]) times = [0.0, 0.4, 0.8, 1.2] se3_snapshots = [apply_transform(se3_exp(xi, t=t), body_points) for t in times]

In [10]:

fig = plt.figure(figsize=(14, 5.4))
fig.suptitle('Experiment 4: A twist traces a rigid motion curve in $SE(3)$', fontsize=13, fontweight='bold')

for idx, (t, points) in enumerate(zip(times, se3_snapshots), start=1):
    ax = fig.add_subplot(1, 4, idx, projection='3d')
    ax.scatter(points[:, 0], points[:, 1], points[:, 2], c=np.linspace(0.2, 0.9, len(points)), cmap='viridis', s=40)
    for segment in [(0, 2), (2, 4), (0, 5), (2, 6), (4, 7), (5, 6), (6, 7)]:
        p, q = segment
        ax.plot(
            [points[p, 0], points[q, 0]],
            [points[p, 1], points[q, 1]],
            [points[p, 2], points[q, 2]],
            color='#6c757d',
            linewidth=1.8,
        )
    ax.set_title(f't = {t}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-0.5, 1.8)
    ax.set_ylim(-0.6, 1.6)
    ax.set_zlim(-0.1, 1.0)

plt.tight_layout()
plt.show()

fig = plt.figure(figsize=(14, 5.4)) fig.suptitle('Experiment 4: A twist traces a rigid motion curve in $SE(3)$', fontsize=13, fontweight='bold') for idx, (t, points) in enumerate(zip(times, se3_snapshots), start=1): ax = fig.add_subplot(1, 4, idx, projection='3d') ax.scatter(points[:, 0], points[:, 1], points[:, 2], c=np.linspace(0.2, 0.9, len(points)), cmap='viridis', s=40) for segment in [(0, 2), (2, 4), (0, 5), (2, 6), (4, 7), (5, 6), (6, 7)]: p, q = segment ax.plot( [points[p, 0], points[q, 0]], [points[p, 1], points[q, 1]], [points[p, 2], points[q, 2]], color='#6c757d', linewidth=1.8, ) ax.set_title(f't = {t}') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.set_xlim(-0.5, 1.8) ax.set_ylim(-0.6, 1.6) ax.set_zlim(-0.1, 1.0) plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The same rigid body moves through a sequence of poses that combine turning and translating, while preserving the body's internal distances.
2. Why it looks this way: exp(t * hat(xi)) stays inside SE(3) for every t, so each snapshot is a valid rigid motion rather than an arbitrary affine transformation.
3. Key takeaway: A twist is a local coordinate in the algebra, but its exponential is a globally valid pose update on the motion group.

---

Summary¶

Experiment	Key insight
1. Rotation order matters	SO(3) composition is non-commutative.
2. Local linearization vs exact exponential	Tangent-space approximations are only reliable near the identity.
3. Geodesic interpolation vs Euler lerp	Coordinate interpolation and manifold interpolation are different operations.
4. Twists generate rigid motion	se(3) provides local coordinates whose exponential stays on SE(3).

Connection to the broader theory¶

These experiments capture the practical reason Lie groups matter: rotations and rigid motions are not ordinary vectors, yet we still need local coordinates for optimization and inference. The Lie algebra gives that local linear space, and the exponential map returns us to the correct geometric object.