This page was generated from docs/tutorials/MatrixRepresentationsOfGeometricFunctions.ipynb. Interactive online version: Binder badge.

Matrix Representations of Geometric Functions

More info can be found in Chapter 5 of “New foundations for classical mechanics”, by David Hestenes.

This notebook shows how some matrix groups can be represented in geometric algebra. Not as spinors in CGA, or PGA, just as functions in plain old Geometric Algebra.
This is done by:
  1. Creating a geometric function

  2. Apply it to an orthonormal frame

  3. Convert the resultant frame into a matrix

The matrix is defined as the inner product of each basis element of original and transformed frame.

\[M_{ij} = a_i\cdot b_j = a_i\cdot f(a)_j\]

(or vice-versa with the i,j, you get the point). Since we are going to do this repeatedly, define a func2Mat()

[1]:
import numpy as np

def func2Mat(f,I):
    '''
    Convert a function acting on a vector into a matrix, given
    the space defined by psuedoscalar I
    '''
    A = I.basis()
    B = [f(a) for a in A]
    M = [float(b | a) for a in A for b in B]
    return np.array(M).reshape(len(B), len(B))

Start with initializing a euclidean N-dimensional algebra and assign our pseudoscalar to \(I\), pretty standard.

[2]:
from clifford import Cl
from math import *

l,b = Cl(3)        # returns (layout,blades). you can change dimesion here
I = l.pseudoScalar

Anti-symmetric

This is so easy,

\[x \rightarrow x\cdot B\]
[3]:
B = l.randomIntMV()(2) # we use randIntMV because its easier to read
f = lambda x:x|B
func2Mat(f,I=I)
[3]:
array([[ 0., -3.,  9.],
       [ 3.,  0., -6.],
       [-9.,  6.,  0.]])

Whats the B? you can read its values straight off the matrix.

[4]:
B
[4]:
(3^e12) - (9^e13) + (6^e23)

Diagonal ( Directional Scaling)

A bit awkward this one, but its made by projection onto each basis vector, then scaling the component by some amount.

\[x \rightarrow \sum{\lambda_i (x\cdot e_i) e_i}\]
[5]:

ls = range(1,len(I.basis())+1) # some dilation values (eigenvalues)
A = I.basis()


d = lambda x: sum([(x|a)/a*l for a,l in zip(A,ls)])
func2Mat(d,I=I)
[5]:
array([[1., 0., 0.],
       [0., 2., 0.],
       [0., 0., 3.]])

Orthgonal, Rotation

\[x\rightarrow Rx\tilde{R}\]

where

\[R=e^{B/2}\]
[6]:
B = l.randomMV()(2)
R = e**(B/2)
r = lambda x: R*x*~R
func2Mat(r,I=I)
[6]:
array([[-0.02877844,  0.83278728,  0.5528446 ],
       [-0.2119728 , -0.54557927,  0.81080873],
       [ 0.97685175, -0.09385421,  0.19222917]])

The inverse of this is ,

\[x\rightarrow \tilde{R}xR\]
[7]:
rinv = lambda x: ~R*x*R # the inverse rotation
func2Mat(rinv,I=I)
[7]:
array([[-0.02877844, -0.2119728 ,  0.97685175],
       [ 0.83278728, -0.54557927, -0.09385421],
       [ 0.5528446 ,  0.81080873,  0.19222917]])

Orthogonal, Reflection

\[x \rightarrow -axa^{-1}\]
[8]:
a = l.randomIntMV()(1)
n = lambda x: -a*x/a
func2Mat(n,I=I)
[8]:
array([[ 0.97530864, -0.09876543, -0.19753086],
       [-0.09876543,  0.60493827, -0.79012346],
       [-0.19753086, -0.79012346, -0.58024691]])
[9]:
a
[9]:
-(1^e1) - (4^e2) - (8^e3)

Notice the determinant for reflection is -1, and for rotation is +1.

[10]:
from numpy.linalg import det
det(func2Mat(n,I=I)), det(func2Mat(r,I=I))
[10]:
(-0.9999999999999999, 1.0000000000411782)

Symmetric

This can be built up from the functions we just defined ie Rotation*Dilation/Rotation

\[x \rightarrow r(d(r^{-1}(x))\]

which if you write it out, looks kind of dumb

\[x \rightarrow R[\sum{\lambda_i (\tilde{R}x R\cdot e_i) e_i]R}\]

So, the antisymmetric matrix is interpreted as a set dilations about a some orthogonal frame rotated from the basis (what basis,eh? exactly what basis!).

More generally we could include reflection in the \(R\) too.

[11]:

g = lambda x: r(d(rinv(x)))
func2Mat(g,I=I)
[11]:
array([[2.30480895, 0.44215098, 0.13438513],
       [0.44215098, 2.61247833, 0.36292709],
       [0.13438513, 0.36292709, 1.08271272]])

Eigen stuffs

By definition the eigen-stuff is the invariants of the transformation, sometimes this is a vector, and other times it is a plane.

Rotation

The eigen blades of a rotation are really the axis and plane of rotation.

[12]:
from numpy.linalg import eig

vals, vecs = eig(func2Mat(r,I=I))
np.round(vecs,3)
[12]:
array([[ 0.626+0.j   , -0.136-0.535j, -0.136+0.535j],
       [ 0.293+0.j   ,  0.676+0.j   ,  0.676-0.j   ],
       [ 0.723+0.j   , -0.157+0.463j, -0.157-0.463j]])

If you checkout the real column, and compare this to the bivector which generated this rotation (aka the generator), after its been normalized

[13]:
B/(abs(B))
[13]:
(0.72272^e12) - (0.29331^e13) + (0.62581^e23)
[14]:
B
[14]:
(1.68666^e12) - (0.68452^e13) + (1.46049^e23)
[15]:
vals
[15]:
array([ 1.        +0.j        , -0.69106427+0.72279332j,
       -0.69106427-0.72279332j])
[16]:

cos(abs(B)), sin(abs(B))
[16]:
(-0.6910642689158123, 0.7227933150132574)

Symmetric

For the symmetric matrix, the invariant thing the orthonormal frame along which the dilations take place

[17]:
vals, vecs = eig(func2Mat(g,I=I))
np.round(vecs,5).T
[17]:
array([[ 0.55284,  0.81081,  0.19223],
       [ 0.83279, -0.54558, -0.09385],
       [-0.02878, -0.21197,  0.97685]])

This is easily found by using the rotation part of the symmetric operator,

[18]:
[R*a*~R  for a in I.basis()]
[18]:
[-(0.02878^e1) - (0.21197^e2) + (0.97685^e3),
 (0.83279^e1) - (0.54558^e2) - (0.09385^e3),
 (0.55284^e1) + (0.81081^e2) + (0.19223^e3)]

Primitive Visualization in 2D

[19]:
from pylab import linspace, plot,axis,legend

def plot_ps(ps,**kw):
    x = [p[e1] for p in ps]
    y = [p[e2] for p in ps]
    plot(x,y, marker='o',ls='',**kw)


l,b = Cl(2)
locals().update(b)
I = l.pseudoScalar

## define function of interest
B = l.randomMV()(2)
R = e**(B/2)
f = lambda x:  R*x*~R

## loop though cartesian grid and apply f,
ps,qs=[],[]
for x in linspace(-1,1,11):
    for y in linspace(-1,1,11):
        p = x*e1 + y*e2
        q = f(p)
        ps.append(p)
        qs.append(q)

plot_ps(ps,label='before')
plot_ps(qs,label='after')
axis('equal')
legend()
/home/docs/checkouts/readthedocs.org/user_builds/clifford-alex/envs/latest/lib/python3.8/site-packages/traitlets/traitlets.py:3030: FutureWarning: --rc={'figure.dpi': 96, 'savefig.transparent': True} for dict-traits is deprecated in traitlets 5.0. You can pass --rc <key=value> ... multiple times to add items to a dict.
  warn(
[19]:
<matplotlib.legend.Legend at 0x7f6a014c8a30>
../_images/tutorials_MatrixRepresentationsOfGeometricFunctions_35_2.svg
[20]:
func2Mat(f,I =I )
[20]:
array([[ 0.82037719,  0.57182276],
       [-0.57182276,  0.82037719]])