Projective Geometric Algebra

In the conformal model, null vectors are instrumental in throwing off the shackles of coordinates. We cleverly built null vectors starting from a wacky inner product that induced negative squared norms.

But why take the scenic route? Why not simply define a basis vector \(\newcommand{\vv}[1]{\mathbf{#1}} \newcommand{\a}{\vv{a}} \newcommand{\b}{\vv{b}} \newcommand{\P}{\vv{P}} \newcommand{\Q}{\vv{Q}} \newcommand{\I}{\vv{I}} \newcommand{\x}{\vv{x}} \newcommand{\y}{\vv{y}} \newcommand{\z}{\vv{z}} \newcommand{\E}{\vv{E}} \newcommand{\e}{\vv{e}} \vv{v}\) to be a null vector? What stops us demanding \(\vv{v}^2 = 0\) by fiat?

Nothing. However, it results in a degenerate inner product, leading to a trickier algebra that historically deterred researchers. For example, \(\I\) is no longer invertible, which makes duals less accessible.

Happily, other intrepid researchers overcame the obstacles, leading to projective geometric algebra (PGA). Simpler and more efficient than the conformal model, projective geometric algebra still handles rotations and translations beautifully.

Look familiar? Let \(f(X) = e^{X/2}\) and:

\[ V = f(12 \z\infty) f(\y\x \gamma) f(8\x\infty) f(\z\y \beta) f(2\y\infty) f(\x\z \alpha) \]

Again, the map \(p \mapsto V p V^{-1}\) as \(\alpha, \beta, \gamma\) range over all angles produces the above animation. But this time we started from a 4-dimensional vector space, leading to a 16-dimensional geometric algebra. Thus the number of entries in a multivector is at most the number of entries in a traditional workhorse of computer graphics, a 4x4 matrix.

There are drawbacks. We can no longer wedge 4 points to describe a circle. We lose succinct inversion and dilation. Perhaps for proving advanced geometry theorems we may be better off with the conformal model. But for many applications projective geometric algebra may be ideal.

PGA Tour

We copy over definitions from older articles:

{-# LANGUAGE CPP #-}
{-# LANGUAGE OverloadedLists #-}
import Control.Monad
import Data.Bool
import Data.List
import Data.Map (Map)
import qualified Data.Map as M
import Data.Maybe
#ifdef __HASTE__
import Data.IORef
import Haste.DOM
import Haste.Events
import Haste.Graphics.Canvas
import Haste.Graphics.AnimationFrame
#endif

type Multivector = Map String Double

m ! s = M.findWithDefault 0 s m

a ~== b = abs (a - b) < 10**(-6)

wedge :: Multivector -> Multivector -> Multivector
wedge m n = M.filter (not . (~== 0)) $
  M.fromListWith (+) $ catMaybes $ f <$> M.assocs m <*> M.assocs n where
  f a@(s, _) b@(t, _)
    | length u == length s + length t = Just c
    | otherwise                       = Nothing
    where c@(u, _) = basisMul a b

dot :: Multivector -> Multivector -> Multivector
dot m n = M.filter (not . (~== 0)) $
  M.fromListWith (+) $ catMaybes $ f <$> M.assocs m <*> M.assocs n where
  f a@(s, _) b@(t, _)
    | length u == length t - length s = Just c
    | otherwise                       = Nothing
    where c@(u, _) = basisMul a b

mul :: Multivector -> Multivector -> Multivector
mul x y = M.filter (not . (~== 0)) $
  M.fromListWith (+) $ basisMul <$> M.assocs x <*> M.assocs y

The bluntest approach might be to define one null vector for each dimension in our space. Points will certainly be represented by null vectors, but this takes degeneracy too far. In particular, if we want dot products to give useful information, then we need \(n\) non-null basis vectors.

Therefore we start with the standard inner product of \(\mathbb{R}^n\), then add a null vector \(\e_0\) giving us an \((n+1)\)-dimensional space known as \(\mathbb{R}_{n,0,1}\). In our code, we write 0 for \(\e_0\) and handle this case specially, defaulting to the standard behaviour for all other basis vectors.

basisMul (s, m) (t, n) = gnome [] (s ++ t) (m * n) where
  gnome pre (a:b:rest) c
    | a < b  = gnome (pre ++ [a]) (b:rest) c
    | a == b = case a of
      '0' ->   ([], 0)
      _   ->   back pre rest c
    | 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)

Next, we move to the dual space \(\mathbb{R}^*_{n,0,1}\) to spread \(\e_0\) around. Recall the dual of a \(k\)-blade is an \((n-k)\)-blade representing its orthogonal complement.

For \(\mathbb{R}^*_{2,0,1}\), if \(\e_1\) and \(\e_2\) are non-null basis vectors (which show up as 1 and 2 in our code), then we represent 1-vectors with the duals of \(\e_0, \e_1, \e_2\).

Already we are hampered by degeneracy. We usually multiply by \(\I\) to compute a dual, but this fails for deriving \(\E_0\) from \(\e_0\), since \(\e_0\) is a null vector. Furthermore, without this portal gun to teleport us between a space and its dual, we are forced to think "backwards" all the time: wedge products describe intersections and meet products describe spans. (If an orthogonal complement grows, then the original subspace shrinks, and vice versa.)

It so happens we can temporarily pretend to have a standard inner product and multiply \(\e_0, \e_1, \e_2\) by \(\I\) to get their duals:

\[ \E_0 = \e_1 \e_2, \E_1 = \e_2 \e_0, \E_2 = \e_0 \e_1 \]

The \(\E_1\) and \(\E_2\) represent the \(x\) and \(y\) directions, and \(\E_0\) is the origin. All but \(\E_0\) are null.

We represent euclidean points and ideal points as follows:

one :: Multivector
one = [("", 1)]

canon = mul one

point2 x y = canon [("20", x), ("01", y), ("12", 1)]
ideal2 x y = canon [("20", x), ("01", y)]

With these representations, for a finite point \(\P\) we have \(\P^2 = -1\), and for an ideal point \(V\) we have \(V^2 = 0\):

square x = mul x x
prop_finitePointSquaredIsMinusOne = square (point2 123 456) == [("",-1)]
prop_idealPointSquaredIsZero = square (ideal2 123 456) == []

As we’re in a projective space, \(\lambda \P\) also represents \(\P\) for any nonzero scalar \(\lambda\); a finite point is normalized when \(\P^2 = -1\). Ideal points always square to zero.

Like before, we avoid geometrically meaningless squares thanks to null vectors. The square of a (normalized) point is a constant, so is independent of the origin.

We test our code by reproducing the multiplication table of the geometric algebra basis vectors from one of Gunn’s papers:

basis2 :: [Multivector]
basis2 = one : map (flip M.singleton 1) (words "0 1 2 12 20 01 012")

pseudo2 = [("012", 1)] :: Multivector

table1 :: [[Multivector]]
table1 = [[mul r c | c <- basis2] | r <- basis2]
crudeTable1 = putStr $ unlines $ map (unwords . map (show . M.toList)) table1
[("",1.0)] [("0",1.0)] [("1",1.0)] [("2",1.0)] [("12",1.0)] [("02",-1.0)] [("01",1.0)] [("012",1.0)]
[("0",1.0)] [] [("01",1.0)] [("02",1.0)] [("012",1.0)] [] [] []
[("1",1.0)] [("01",-1.0)] [("",1.0)] [("12",1.0)] [("2",1.0)] [("012",1.0)] [("0",-1.0)] [("02",-1.0)]
[("2",1.0)] [("02",-1.0)] [("12",-1.0)] [("",1.0)] [("1",-1.0)] [("0",1.0)] [("012",1.0)] [("01",1.0)]
[("12",1.0)] [("012",1.0)] [("2",-1.0)] [("1",1.0)] [("",-1.0)] [("01",-1.0)] [("02",-1.0)] [("0",-1.0)]
[("02",-1.0)] [] [("012",1.0)] [("0",-1.0)] [("01",1.0)] [] [] []
[("01",1.0)] [] [("0",1.0)] [("012",1.0)] [("02",1.0)] [] [] []
[("012",1.0)] [] [("02",-1.0)] [("01",1.0)] [("0",-1.0)] [] [] []

The line \(a x + b y + c = 0\) is represented by:

line2 a b c = [("1", a), ("2", b), ("0", c)] :: Multivector

A euclidean line \(\a\) is normalized when \(a^2 = 1\).

We provide similar routines for 3D geometry:

point3 x y z = canon [("032", x), ("013", y), ("021", z), ("123", 1)]
ideal3 x y z = canon [("032", x), ("013", y), ("021", z)]
plane3 a b c d = [("1", a), ("2", b), ("3", c), ("0", d)] :: Multivector
pseudo3 = [("0123", 1)] :: Multivector

Our demo needs to retrieve the coordinates of a euclidean point for the 3D case:

-- Assumes input is a euclidean point.
coords3 :: Multivector -> [Double]
coords3 m = [-get "023", get "013", -get "012"] where
  d     = m!"123"
  get s = m!s / d

Squared away

Define the squared norm of a euclidean \(k\)-vector \(\vv{X}\) by:

\[ || \vv{X} ||^2 = \vv{X}^2 \]

(This formula does not apply to all \(k\)-vectors, but rather only those representing points, lines, planes, and so on.)

The square of an ideal element is always zero, so we also define the square of the ideal norm \(||\vv{X}||^2_\infty\) to be the sum of the squares of the coefficients, and say an ideal \(k\)-vector \(\vv{X}\) is normalized when \(||\vv{X}||^2_\infty = 1\).

This can be formulated without coordinates:

\[ ||\vv{X}||^2_\infty = ||\vv{X}\vee\P||^2 \]

where \(\P\) is any normalized finite point. However, our implementation simply reads off the coefficients:

scalarAdd = M.insertWith (+) ""
scalarMul a v = (*a) <$> v

normSquared :: Multivector -> Double
normSquared v = square v ! ""

idealNormSquared :: Multivector -> Double
idealNormSquared = sum . map (^2) . M.elems

normalize :: Multivector -> Multivector
normalize v
  | d ~== 0 = scalarMul (sqrt (1/idealNormSquared v)) v
  | otherwise = scalarMul (sqrt (-1/d)) v
  where d = normSquared v

The meet product

Recall in the conformal model we computed wedge products to find the span of subspaces, such as the line defined by two points, or the plane defined by three points. Because PGA uses the dual space, we need meet products instead.

Usually we take advantage of \(A \vee B = (A^* \wedge B^*)^*\). But without an easy dual operation, we resort to a shuffle product to compute the meet, a sort of geometric version of the inclusion-exclusion principle.

We need a helper function that returns all partitions of a sequence into two subsequences, along with the sign of the permutation obtained by concatening the subsequences.

signedPartitions sz xs = f sz xs where
  f n xs
    | n > length xs = []
    | n == 0 = [(1, [], xs)]
    | (x:xt) <- xs = [(s, x:as, bs) | (s, as, bs) <- f (n - 1) xt]
                  ++ [(parity n s,as, x:bs) | (s, as, bs) <- f n xt]
  parity n | odd n = negate
           | otherwise = id

prop_signedPartitionsExample = signedPartitions 2 [1..5] ==
  [(1,[1,2],[3,4,5])
  ,(-1,[1,3],[2,4,5])
  ,(1,[1,4],[2,3,5])
  ,(-1,[1,5],[2,3,4])
  ,(1,[2,3],[1,4,5])
  ,(-1,[2,4],[1,3,5])
  ,(1,[2,5],[1,3,4])
  ,(1,[3,4],[1,2,5])
  ,(-1,[3,5],[1,2,4])
  ,(1,[4,5],[1,2,3])
  ]

To compute meets with a shuffle, we need the dimension of the space, which in a cleaner library might be held somewhere by each Multivector value. (Such a library could also prevent mixing up multivectors from different spaces.) Instead, we require the caller to supply the dimension.

vee :: Int -> Multivector -> Multivector -> Multivector
vee dim m n = M.filter (not . (~== 0)) $
  M.fromListWith (+) $ concat $ f <$> M.assocs m <*> M.assocs n where
  f a@(s, m) b@(t, n)
    | dim > length s + length t = []
    | otherwise = [(bs, fromIntegral sgn*m*n) |
      (sgn, as, bs) <- signedPartitions (dim - length t) s,
      intersect as t == []]

vee2 = vee 3
vee3 = vee 4

Products in 2D

Let’s start with the geometric product in 2D on normalized inputs.

Multiplication by \(\I\) is commutative, and computes the orthogonal complement in the Euclidean 2D space. (Not the geometric algebra itself, otherwise we would have easy access to the dual!)

For a euclidean line \(\a\), the product \(\a\I\) is the ideal point in the direction perpendicular to the line.

For a euclidean point \(\P\), the product \(\P\I\) is \(\e_0\), the ideal line.

prop_linePerp = line2 1 2 3 `mul` pseudo2 == ideal2 1 2
prop_pointPerp = point2 42 9000 `mul` pseudo2 == [("0",-1)]

The product of two euclidean lines is:

\[ \a\b = \a \cdot \b + \a \wedge \b \]

The dot product \(\a \cdot \b\) is the cosine of the angle \(\theta\) between them.

Let \(\P\) denote their normalized point of intersection. If \(\P\) is finite, that is, if the lines intersect, then \(\a \wedge \b = (\sin \theta) \P\). If \(\P\) is ideal, that is, the lines are parallel, then \(\a \wedge \b = d_{ab} \P\) where \(d_{ab}\) denotes the distance between the lines.

The product of two euclidean points is:

\[ \P \Q = -1 + d_{PQ} \vv{V} \]

where \(d_{PQ}\) is the distance between them and \(\vv{V}\) is an ideal line representing the direction perpendicular to the line joining them. We define

\[ \P \times \Q = d_{PQ} \vv{V} \]

hence \(d_{PQ} = || \P \times \Q ||_\infty\).

The product of a euclidean line and a euclidean point is:

\[ \a \P = \a \cdot \P + \a \wedge \P = \a^\perp_\P + d_{\a\P}\I \]

where \(\a^\perp_\P = \a \cdot \P\) is the line through \(\P\) perpendicular to \(\a\) and \(d_{\a\P}\) is the distance between the point and the line.

It turns out we also have meaningful products when one or both of the inputs are ideal.

Transformations

There are more geometric interpretations to explore and we have yet to touch the 3-dimensional case. But let’s fast-forward to the part we’ve all been waiting for: rotors!

The reflection of \(\vv{X}\) in the \((n-1)\)-hyperplane \(\a\) is:

\[ \a\vv{X}\a \]

This formula reflects points, lines, planes, and so on.

We could finish here, as all Euclidean isometries can be generated from reflections. Composing reflections through \(\a\) then \(\b\) is simply:

\[ \b\a\vv{X}\a\b \]

If \(\vv{R} = \b \a\) is normalized then we call it a rotor, and transform a subspace \(\vv{X}\) by computing \(\vv{R} \vv{X} \tilde{\vv{R}}\), where \(\tilde{\vv{R}}\) is the reversal of \(\vv{R}\), meaning we multiply the factors in reverse order. Calculation shows \(\vv{R}\tilde{\vv{R}} = 1\).

Rotations have an exponential form:

\[ e^{\frac{\theta}{2} \P} \]

where \(\P\) is the normalized point (in 2D) or axis (in 3D) at the center of a rotation by the angle \(\theta\). For euclidean points, \(\P^2 = -1\) so \(\P\) takes the role played by \(\vv{i}\) in other models.

For example, a rotation around the \(y\)-axis by 120 degrees is given by:

approxEq v w = all (~== 0) $ M.elems $ M.unionWith (-) v w

prop_rotateExample =
  axis == [("13",-1)] && rotor `approxEq` [("",0.5),("13",-sqrt(3)/2)]
  where
  rotor = scalarAdd (cos (a/2)) $ scalarMul (sin (a/2)) axis
  axis = normalize $ point3 0 0 0 `vee3` point3 0 1 0
  a = 2*pi/3

Observe we computed the meet product of two points to find the multivector describing the \(y\)-axis. It turns out to have a simple form, though in general lines in 3D are less straightforward in PGA:

prop_lineExample =
  line == [("01",-23),("02",-14),("03",17),("12",-7),("13",-3),("23",-7)]
  where
  line = point3 1 2 3 `vee3` point3 6 5 4

In 2D, translation by \(d\) units can be viewed as a rotation of size \(d\) with the center at an ideal point \(\P\) in the perpendicular direction. By considering the Taylor series expansion, and recalling \(\P^2 = 0\), we find:

\[ e^{\frac{d}{2} \P} = 1 + \frac{d}{2} \P \]

For a 3D translation, the axis is an ideal line, which Gunn states can be found by computing \((\E_0 \vee \vv{V})\I\) where \(\vv{V}\) is the direction of travel. It seems we find the line joining the origin and \(\vv{V}\) then take its orthogonal complement to yield the axis.

\[ e^{\frac{d}{2} (\E_0 \vee \vv{V})\I} \]

For example, the following rotor translates 12 units along the \(z\)-axis:

prop_translateExample =
  zdir == [("03",-1)] && rotor == [("",1),("03",-6)]
  where
  zdir = mul (M.singleton "123" 1 `vee3` ideal3 0 0 1) pseudo3
  rotor = scalarAdd 1 $ scalarMul (12/2) zdir

Rotations and translations involving the coordinate axes have a simple form so are great for demos.

moon a b c p = foldl1 mul $ rs ++ [p] ++ map inv (reverse rs) where
  inv = M.mapWithKey f where
    f [] v = v
    f _  v = -v
  rs =
    [ tra "03" 12
    , rot "21" c
    , tra "01" 8
    , rot "32" b
    , tra "02" 2
    , rot "13" a
    ]
  rot axis ang = [("", cos (ang/2)), (axis, sin (ang/2))] :: Multivector
  tra dir d = [("", 1), (dir, d/2)] :: Multivector

Above, our definition for \(V\) uses notation from our adventure in the conformal model: the symbols \(\infty, \x, \y, \z\) correspond to \(\e_0,\e_1,\e_2,\e_3\) which are 0 1 2 3 in our code.

Into orbit

Lastly, we slap together some code to animate a cube on this webpage.

#ifdef __HASTE__
-- Cube faces and their colours.
cubeFC = zip
  [[0, 1, 3, 2], [4, 6, 7, 5],
   [0, 4, 5, 1], [2, 3, 7, 6],
   [1, 5, 7, 3], [0, 2, 6, 4]]
-- http://the-rubiks-cube.deviantart.com/journal/Using-Official-Rubik-s-Cube-Colors-268760351
  [RGB 0 0x51 0xba, RGB 0 0x9e 0x60,
   RGB 0xff 0xd5 0, RGB 0xff 0xff 0xff,
   RGB 0xff 0x58 0, RGB 0xc4 0x1e 0x3a]

paint canvas tRaw = sequence_ $ renderOnTop canvas . drawFace <$> filter vis cubeFC
  where
  ms = round tRaw :: Int
  ang s = fromIntegral (mod ms (s*1000)) * 2 * pi / fromIntegral (s*1000)
  cubeVertices = (\[a,b,c] -> point3 a b c) <$> replicateM 3 [-1, 1]
  ws = coords3 . moon (ang 2) (ang 5) (ang 23) <$> cubeVertices
  drawFace (face, rgb) = color rgb . fill . path $ project . (ws!!) <$> face
  project :: [Double] -> (Double, Double)
  project [x, y, z] = (x/z*160 + 160, y/z*160 + 160)
  -- Cull faces pointing away from us.
  vis (face, _) = d!"" <= 0 where
    (u:v:_) = zipWith (zipWith (-)) t ts
    d = foldl1' mul [toGA u, toGA v, M.singleton "zyx" 1, toGA h]
    ts@(h:t) = (ws!!) <$> face
    toGA = M.fromList . zip (pure <$> "xyz")

main = withElems ["canvas"] $ \[cElem] -> do
  Just canvas <- fromElem cElem
  paused <- newIORef False
  let
    anim t = do
      render canvas $ fill $ rect (0, 0) (320, 320)
      paint canvas t
      b <- readIORef paused
      unless b $ void $ requestAnimationFrame anim
    pause = do
      b <- readIORef paused
      writeIORef paused $ not b
      when b $ void $ requestAnimationFrame anim
  _ <- cElem `onEvent` MouseDown $ const pause
  requestAnimationFrame anim
#endif

Ben Lynn blynn@cs.stanford.edu 💡