jsEval "curl_module('../compiler/Map.ob')"
All Hail Geometric Algebra!
Rotate a cube by 90 degrees around the \(x\)-axis, then rotate it by 90 degrees around the \(z\)-axis. This is equivalent to a single rotation. What is the axis and angle of this single rotation?
If we hold a physical cube in our hands, we can turn it a couple of times and intuitively find the answer. Is there a similarly straightforward mathematical method that rigorously solves this problem?
What if we generalize? That is, given any two rotations that fix the origin, what is the single rotation that is equivalent to their composition? Can we solve this just as easily?
Yes! With geometric algebra, we can solve the first problem comfortably with pen and paper, and apply the same method to its generalization.
We’ll take a whirlwind tour of A Survey of Geometric Algebra and Geometric Calculus by Alan Macdonald. See also:
-
The history of geometric algebra according to me: how principles familiar to functional programmers led to a single elegant powerful framework that unifies well-known tools for studying geometry: Cartesian coordinates, Plücker coordinates, complex numbers, vector algebra, quaternions.
-
A paper comparing geometric algebra and vector algebra for ray tracing by Daniel Fontijne and Leo Dorst.
Geometric algebra for (almost) free
Let \(V\) be a vector space over \(\mathbb{R}\), and consider the free algebra on the set \(V\). In other words, we compute with multinomials where every vector is an indeterminate, except the order of the indeterminates matters. For example:
\[ \newcommand{\n}{\mathbf{n}} \newcommand{\u}{\mathbf{u}} \newcommand{\v}{\mathbf{v}} \newcommand{\w}{\mathbf{w}} \newcommand{\i}{\mathbf{i}} \newcommand{\e}{\mathbf{e}} \newcommand{\f}{\mathbf{f}} \newcommand{\vv}[1]{\mathbf{#1}} (2 \vv{vuvw} + 1)(3\v - 2) = 6\vv{vuvwv} - 4\vv{vuvw} + 3\v - 2 \]
We call these things multivectors, and they are easy to implement with maps.
import Map
type Multivector = Map String Double
freeMul :: Multivector -> Multivector -> Multivector
freeMul x y = mfilter (/= 0) $
fromListWith (+) $ f <$> assocs x <*> assocs y where
f (s, m) (t, n) = (s ++ t, m * n)
freeMul (fromList [("vuvw", 2), ("", 1)]) (fromList [("v", 3), ("", -2)])
What a pleasant algebra! There is an identity element. The product is associative, and distributes over addition. What a shame it’s also useless: how do we geometrically interpret \(\vv{uvwuuv}\)?
There’s no way. But for finite-dimensional \(V\), by constraining our free algebra with one simple relation, we obtain the vaunted geometric algebra, by requiring that, for any vector \(\v\):
\[ \v\v = \v\cdot\v \]
where the right-hand side denotes the standard dot product. We could generalize by replacing the dot product with any quadratic form, and by considering any vector space over any field, but we focus on the standard dot product on \(\mathbb{R}^n\), where the resulting geometric algebra is denoted \(\mathbb{G}^n\).
We have \(\v^2 = |\v|^2 \), and by the polarization identity:
\[ \u \cdot \v = \frac{1}{2} ((\u + \v) \cdot (\u + \v) - \u \cdot \u - \v \cdot \v) \]
From our constraint, this is:
\[ \u \cdot \v = \frac{1}{2} ((\u + \v)^2 - \u^2 - \v^2) \]
By elementary algebra, this implies:
\[ \u \cdot \v = \frac{1}{2} (\u\v + \v\u) \]
Hence when \(\u\) and \(\v\) are orthogonal, we have \(\u\v = -\v\u\).
As its name suggests, geometric algebra is useful for geometry yet possesses nice algebraic properties. Even division is possible sometimes: the inverse of a nonzero vector \(\v\) is \(\v / |\v|^2\).
[By the way, this section’s title is unoriginal.]
The canonical basis
To recap:
-
We can replace \(\v\v\) with \(|\v|^2\).
-
We can replace \(\u\v\) with \(-\v\u\) when \(\v\) and \(\u\) are orthogonal.
Thus if \(\e_1, …, \e_n\) is an orthonormal basis for \(\mathbb{R}^n\), then every multivector of \(\mathbb{G}^n\) can be written as a linear combination of products of the form:
\[ \e_{i_1}…\e_{i_k} \]
where \(i_1,…,i_k\) is a strictly increasing subsequence of \([1..n]\). These products are a canonical basis for \(\mathbb{G}^n\), which implies the basis for an \(n\)-dimensional vector space yields a canonical basis of size \(2^n\) for its geometric algebra.
Here’s a function showing the details. It computes the geometric product in canonical form, assuming each character that appears in the input represents an orthonormal basis vector of the space. For clarity, we employ the inefficient gnome sort:
-- Computes geometric product of multivectors written in terms of some
-- orthonormal basis. Returns result in canonical form.
mul :: Multivector -> Multivector -> Multivector
mul x y = mfilter (/= 0) $
fromListWith (+) $ basisMul <$> assocs x <*> assocs y
basisMul (s, m) (t, n) = gnome [] (s ++ t) (m * n) where
gnome pre (a:b:rest) c
-- Vectors are in order: move on.
| a < b = gnome (pre ++ [a]) (b:rest) c
-- Same vector: remove both copies since magnitude = 1.
| a == b = back pre rest c
-- Wrong order: flip and back up one step. Since the vectors are
-- orthogonal, we flip the sign of the geometric product.
| a > b = back pre (b:a:rest) (-c)
where
back [] rest c = gnome [] rest c
back pre rest c = gnome (init pre) (last pre:rest) c
gnome pre rest c = (pre ++ rest, c)
one = singleton "" 1
canon = mul one . fromList
prop_exampleOrthonormal = canon [("uvwuuv", 1)] == fromList [("uw", -1)]
-- Returns canonical basis given orthonormal basis.
canonicalBasis :: String -> [String]
canonicalBasis = subsequences
If the \(k\) is the same for every term when an multivector is written canonically, then we say the multivector is a \(k\)-vector. A 0-vector is also called a scalar, a 1-vector a vector, a 2-vector a bi-vector, and a 3-vector a tri-vector:
-- Assumes input is in canonical form.
name :: Multivector -> String
name x
| size x == 0 = "zero multivector"
| all ((== k) . length) ss = show k ++ "-vector" ++ t
| otherwise = "mixed multivector"
where
s:ss = keys x
k = length s
t | null (aka k) = ""
| otherwise = " or " ++ aka k
aka 0 = "scalar"
aka 1 = "vector"
aka 2 = "bivector"
aka 3 = "trivector"
aka _ = ""
checkName l s = name (canon l) == s
checkName [] "zero multivector"
checkName [("", 42)] "0-vector or scalar"
checkName [("x", -5), ("y", 98)] "1-vector or vector"
checkName [("ab", 2), ("bc", 3)] "2-vector or bivector"
checkName [("abc", 4)] "3-vector or trivector"
checkName [("abcde", 0.01)] "5-vector"
checkName [("", 2), ("bc", 3)] "mixed multivector"
Lastly, \(n\)-vectors are also called pseudoscalars. We write \(\vv{I}\) for the pseudoscalar \(\e_1 … \e_n\), which turns out to be unique up to sign for any choice of basis. Then \(\vv{I}^{-1} = \e_n … \e_1 = (-1)^{n-1}\vv{I}\).
The fundamental identity
Let \(\e, \f\) be an orthonormal basis of a plane.
sqr x = mul x x
sqr (canon [("e", 3), ("f", 4)]) == canon [("", 25)]
mul (canon [("e", 3), ("f", 4)]) (canon [("e", 5), ("f", 6)])
== canon [("", 39), ("ef", -2)]
mul (canon [("e", 3), ("f", 4)]) (canon [("e", -4), ("f", 3)])
== canon [("ef", 25)]
The geometric product of two vectors is a scalar and a bivector. Furthermore, appearances of \(0\) and \(25 = 5^2 = 3^2 + 4^2\) hints at a deeper truth. Expanding \((a \e + b \f)(c \e + d \f)\) proves the fundamental identity: for any vectors \(\u\) and \(\v\) in our plane:
\[ \u\v = \u\cdot\v + \u\wedge\v \]
where \(\u\wedge\v\) is the outer product or wedge product, the bivector component of the geometric product of two vectors, and satisfies:
\[ \u\wedge\v = |\u| |\v| \sin \theta \i \]
where \(\theta\) is the angle between \(\u\) and \(\v\), and \(\i\) is the bivector \(\e \f\).
We interpret \(\i\) to mean the plane generated by \(\e\) and \(\f\). With vector algebra, we represent planes with normal vectors, which only makes sense in exactly three dimensions. With geometric algebra, we represent a plane by multiplying two orthonormal vectors that generate them; this works in any dimension.
We find \(\i^2 = -1\), which explains the choice of notation. In fact, for any given plane we can construct a complex number system in this manner, that is, multivectors of the form \(a + b\i\). Such numbers are called geometric algebra complex numbers.
We can also write:
\[ \u\v = |\u| |\v| e^{\i\theta} \]
which resembles a standard complex number in polar form.
Suppose \(|\u| = |\v|\). Then left multiplying by \(u\) gives:
\[ \v = \u e^{\i\theta} \]
Thus the geometric algebra interpretation of \(e^{\i\theta}\) is a rotation by \(\theta\) in the plane \(\i\), rather than a particular point on a unit circle measured from some predefined zero bearing. This contrast is a rotational analogue of the difference between vectors and Cartesian points.
Rotations, in any dimension
Complex numbers excel at describing rotations in two dimensions, and quaternions in three. Geometric algebra complex numbers excel at describing rotations in any dimension.
Let \(\i\) be the product of two orthonormal vectors representing the plane in which we wish to rotate. Let \(\theta\) be the angle about the origin we wish to rotate. Let \(\u\) be a vector.
Decompose \(\u\) with respect to \(\i\):
\[ \u = \u_{\perp} + \u_{\parallel} \]
That is, we find vectors \(\u_{\perp} \cdot \i = 0\) and \(\u_{\parallel} \wedge \i = 0\) satisfying the above summation. These are unique and readily computed, but for this section we only need their existence.
Rotation only affects the second term, so \(\u\) rotates to:
\[ \begin{array}{ccl} \u_\perp + \u_\parallel e^{\i\theta} &=& \u_\perp e^{-\i\theta/2} e^{\i\theta/2} + \u_\parallel e^{\i\theta/2} e^{\i\theta/2} \\ &=& e^{-\i\theta/2} \u_\perp e^{\i\theta/2} + e^{-\i\theta/2} \u_\parallel e^{\i\theta/2} \\ &=& e^{-\i\theta/2} \u e^{\i\theta/2} \end{array} \]
where we have used \(\u_\perp \i = \i \u_\perp\) and \(\u_\parallel \i = - \i \u_\parallel\), identities that can be verified by setting \(\i = \e\f\) for orthonormal vectors \(\e,\f\) and a touch of algebra.
In 3D, define \(\n\) by \(\e\f\n =\vv{I}\), which means \(\n\) is a unit normal to \(\i\). Then we can phrase a rotation of a vector \(\u\) about the axis \(\n\) by the angle \(\theta\) as:
\[ e^{-\n\vv{I}\theta/2} \u e^{\n\vv{I}\theta/2} \]
This formula implies Euler’s rotation theorem, as composing two rotations results in an expression of the same form.
The multivector \(e^{-\n\vv{I}\theta/2}\) is called a rotor.
Problem solved
We now solve our opening problem with a minimum of paperwork. Let \(\newcommand{\x}{\vv{x}} \newcommand{\y}{\vv{y}} \newcommand{\z}{\vv{z}} \x, \y, \z\) be an orthonormal basis of \(\mathbb{R}^3\), so \(\vv{I} = \x\y\z\). The following code takes coordinates to a multivector and back:
toGA :: [Double] -> Multivector toGA = fromList . zip (pure <$> "xyz") fromGA :: Multivector -> [Double] fromGA m = (m!) . pure <$> "xyz"
We focus on the right half of the rotation formula. For a rotation of 90 degrees around the \(x\)-axis followed by a rotation of 90-degrees around the \(z\)-axis, it reads:
\[ \begin{array}{ccc} e^{\x\vv{I}45^\circ} e^{\z\vv{I}45^\circ} &=& (\cos 45^\circ + \y\z \sin 45^\circ)(\cos 45^\circ + \x\y \sin 45^\circ) \\ &=& \frac{1}{2} + \frac{\y\z + \x\y - \x\z}{2} \\ &=& \cos 60^\circ + c(\x + \z + \y)\vv{I} \sin 60^\circ \end{array} \]
for some constant \(c > 0\). Therefore the combined rotation is 120 degrees around the axis \(\x + \y + \z\), or in Cartesian coordinates, the axis through the origin and (1, 1, 1).
We can easily write code to compute the axis and angle of the composition of any two 3D rotations with axes through the origin, and use it to check our work:
pi = 3.141592654
atanTaylor 1 = pi * 0.25
atanTaylor x = foldr (+) 0 $ reverse $ take 25 $
zipWith (/) (iterate (*(-x*x)) x) (iterate (+2) 1)
acos x = 2*atanTaylor(sqrt (1-x*x) / (1+x))
cossinSmall theta = cordic lim tab theta 1 0 where
lim = 25
cordic n ((p2, a):rest) beta x y
| n == 0 = (kn * x, kn * y)
| otherwise = cordic (n - 1) rest (beta - sig a) x' y'
where
sig = if beta < 0 then negate else id
x' = x - sig p2 * y
y' = sig p2 * x + y
kn = kvalues!!lim
kvalues = scanl1 (*) $ (\k -> 1/sqrt(1+0.5^(2*k))) <$> [0..]
tab = zip pows $ atanTaylor <$> pows where pows = iterate (*0.5) 1
cossin theta
| theta < 0 = second negate $ cossin (-theta)
| theta <= pi/4 = cossinSmall theta
| theta <= pi/2 = (\(a,b) -> (b,a)) $ cossin (pi/2 - theta)
| theta <= pi = first negate $ cossin (pi - theta)
| theta <= 2*pi = (\(a,b) -> (-a, -b)) $ cossin $ theta - pi
| otherwise = cossin $ theta - 2*pi
cos = fst . cossin
sin = snd . cossin
rotr :: [Double] -> Double -> Multivector
rotr axis theta = canon [("", c), ("yz", nx*s), ("zx", ny*s), ("xy", nz*s)]
where
c = cos (theta / 2)
s = sin (theta / 2)
[nx, ny, nz] = map (/ sqrt(sum $ map (^2) axis)) axis
m # s = maybe 0 id $ mlookup s m
axisAngle :: Multivector -> ([Double], Double)
axisAngle m = ([m#"yz", -(m#"xz"), m#"xy"], acos (m#"") * 2)
prop_rotationExample =
(theta * 180 / pi) ~== 120 && all (~== 1) (map (/ head axis) axis)
where
a ~== b = abs (a - b) < 0.1^5
x90 = rotr [1, 0, 0] $ pi/2
z90 = rotr [0, 0, 1] $ pi/2
(axis, theta) = axisAngle $ mul x90 z90
prop_rotationExample
Let’s go for a spin!
We animate the rotations in question on this very webpage using our geometric algebra routines:
cubeVertices = replicateM 3 [-1, 1] :: [[Double]]
faces = [[0, 1, 3, 2], [4, 6, 7, 5],
[0, 4, 5, 1], [2, 3, 7, 6],
[1, 5, 7, 3], [0, 2, 6, 4]] :: [[Int]]
conjugate r u = rotInv r `mul` u `mul` r
rotInv = mapWithKey $ \k v -> if null k then v else -v
-- Darken, so the white side is visible.
rgbs :: [[Int]]
rgbs = map (\c -> c * 220 `div` 255) <$>
-- http://the-rubiks-cube.deviantart.com/journal/Using-Official-Rubik-s-Cube-Colors-268760351
[[0, 0x51, 0xba], [0, 0x9e, 0x60],
[0xff, 0xd5, 0], [0xff, 0xff, 0xff],
[0xff, 0x58, 0], [0xc4, 0x1e, 0x3a]]
rx t = rotr [1, 0, 0] $ t * pi / 2
rz t = rotr [0, 0, 1] $ t * pi / 2
even n = n .&. 1 == 0
odd n = n .&. 1 == 1
vshadeQ tRaw = conjugate r where
r | even sec = rbase
| odd sec = rbase `mul` ([rx, rz]!!(div sec 2 `mod` 2)) t
(sec, msec) = (tRaw `mod` 12000) `divMod` 1000
rs = scanl mul one $ cycle [rx 1, rz 1]
rbase = rs!!div sec 2
t = fromIntegral msec / 1000
r111 t = rotr [1, 1, 1] $ t * 2 * pi / 3
vshadeA tRaw = conjugate r where
r | even sec2 = rbase
| otherwise = rbase `mul` r111 t
(sec2, msec2) = (tRaw `mod` 12000) `divMod` 2000
rs = scanl mul one $ repeat $ r111 1
rbase = rs!!div sec2 2
t = fromIntegral msec2 / 2000
paint vshade ctx t = do
jsEval $ "cls(" ++ ctx ++ ");"
mapM_ render $ filter vis (zip faces rgbs)
where
render (face, c) = jsEval_ $ ctx ++ ".fillStyle = \"" ++ encodeColor c ++ "\";"
++ ctx ++ ".fill(" ++ encodePath(project . (ws!!) <$> face) ++ ");"
-- f (face, c) = color (toColor c) . fill . path $ project . (ws!!) <$> face
ws = nudge . fromGA . vshade (t) . toGA <$> cubeVertices
-- Arrange signs so that x-axis goes right, y-axis goes up,
-- and z-axis goes out of the screen.
project :: [Double] -> (Double, Double)
project [x, y, z] = ((-x/z + 0.5) * 320 + 160, (y/z - 0.5) * 320 + 160)
-- Move the up, right, and into the screen.
nudge = zipWith (+) [2, 2, 5]
-- Cull faces pointing away from us.
vis (face, _) = d!"" <= 0 where
-- Use GA to compute a cross product and dot product.
dual = (`mul` singleton "zyx" 1)
(u:v:_) = zipWith (zipWith (-)) (tail ts) ts
d = mul (dual (mul (toGA u) (toGA v))) $ toGA (head ts)
ts = (ws!!) <$> face
encodeColor rgb = "rgb(" ++ unwords (show <$> rgb) ++ ")"
encodePath (h:t) = "new Path2D(\"M " ++ showPoint h
++ concatMap ((" L "++) . showPoint) t ++ " Z\")"
where
showPoint (x, y) = show x ++ " " ++ show y
go t = do
paint vshadeQ "ctxQ" t
paint vshadeA "ctxA" t
jsEval_ "maybeAnim();"
jsEval_ [r|maybeAnim = () => {
if (!paused) requestAnimationFrame(t => repl.run("chat", ["Main"], "go " + Math.round(t)));
};
maybeAnim();
|]
We’ve barely scratched the surface. With the geometric algebra homogeneous model of \(\mathbb{R}^n\), we get:
-
4x4 matrices for manipulating points in 3D, just like with vector algebra: rotations, translations, perspective projection, and so on.
-
Simple formulas for subspaces (lines, planes, and beyond). for example, a line is defined by two points \(p, q\), so its equation is \(p \wedge q\); a plane is defined by three points \(p, q, r\) , so its equation is \(p \wedge q \wedge r\).
-
Simple formulas for transforming subspaces, aka outermorphisms: \(f(p \wedge q) = f(p) \wedge f(q)\).
-
Simple formulas for finding the intersection of subspaces: \(a^* \cdot b\).
Did you notice gaps in our reasoning throughout? If not, good! Otherwise, see my linear algebra notes.