Seeking guidance in algorithmically implementing GA operations

Hi @tza,
Thanks for your question. Since you indicated that you didn’t understand my treatment of the Poincare duality in the course notes, I thought I would wait and see how the discussion developed before piping up preferring to see how others might explain things better than I had.
If you’ve watched the SIGGRAPH course video then you’ll know that I see my task as helping you move from the stance of “PGA as magic” to “PGA as understanding” – if you want to. By magic I mean that you are satisfied when you have code that runs and does what you want it to, but you might now be able to explain why it works to your colleague. I consider that a valid position and am not interested in forcing anyone to go further that they want to in the direction of understanding.
If you do want however to understand how euclidean PGA works then I think there are some remaining issues that I’d like to address in this and my next post.

In this post I’d like to deal with the “right complement” and in a separate post I’ll address the issue of why euclidean PGA uses the dual construction (planes as 1-vectors).

The right complement is defined above as follows:

The right complement is basically what thing X^* do we have to right multiply X by to get the pseudoscalar?

This definition is well-defined as long as X^{-1} exists, then X^*=X^{-1}I where I is the pseudoscalar.

But when the metric is degenerate, as with euclidean PGA, and X^2=0, then it’s not well-defined. For example, let X = e_{012}, the ideal point in the z-direction. Then the z-plane e_3 =: X^* certainly satisfies XX^*=I but it’s not unique. In general let X_\alpha^*=e_3 + \alpha e_0, then X^* also satisfies XX^*=I. since e_0^2=0. Geometrically, X_\alpha^* represents the plane parallel to the z-plane at a signed distance of \alpha from it. So if one plane is a right complement then all planes parallel to it are also a right complement according to the above definition. Hence the “right complement” isn’t well-defined.

So if we want to have a well-defined map we need to specify it otherwise. There’s a standard way to do this called the canonical dual basis. [More on this in this post here.] It’s basically the algorithm that the previous posts show: given a simple blade X as the wedge of k basis 1-vectors, define X^* to be the blade formed by wedging together the (n-k) basis 1-vectors that do NOT occur in X (with a sign adjustment if necessary). That certainly will satisfy the previous “right complement” condition XX^{*}=I but due to the restrictive definition, it is now unique. The final step is to use linearity to extend the definition to a general k-vector.

The right complement defined using the canonical dual basis in this way is not however coordinate-free! To see this, in the 2D Grassmann algebra \bigwedge(\mathbb{R}^2) consider replacing e_1 with e_0+e_1 so the new basis \{e_0', e_1'\} = \{e_0, e_0+e_1\}. Then I' = e'_0 \wedge e'_1 = I. Construct the canonical dual bases: for the first basis \{e_0^*, e_1^*\} = \{e_1, -e_0\}, while WRT to the second basis \{{e'}_0^*, {e'}_1^*\} = \{{e'}_1, -{e'}_0\} = \{e_0+e_1, -e_0\} where the second equality follows from substituting. So expressing the dual of e_0 = e_0' WRT the first basis is e_1 while WRT the second basis it is e_1' = e_0+e_1. So this is not coordinate-free.

Exercise: Let \{e_0, e_1, e_2\} be a basis of of the 3D real Grassmann algebra \bigwedge(\mathbb{R}^3). Define a second basis \{e_0', e_1', e_2'\} = \{e_0, e_1, e_0+e_2\}. Show that the canonical dual of e_0 WRT the first basis is e_1\wedge e_2 = e_{12} while WRT the second basis is e_1\wedge (e_0 + e_2) = -e_{01}+e_{12}.

Takeaway: One of the big selling points of geometric algebra is that it is coordinate-free. A danger of coordinate-dependent code in such a system (even when it works for a given coordinate system) is that once that coordinate system changes, the code will probably not work any longer and the idea of looking at the coordinates won’t occur to people since they have been told “GA is coordinate-free.”

Everyone has to decide for himself how comfortable they are with non-coordinate-free code. Read on if you want to learn how one further step yields a coordinate-free “dual” that works in all metrics.

It’s called “Poincare duality” or the “dual coordinate map” and it’s described in various forms in my articles for the past 10 years.

Here’s the short form: take the canonical dual map X^* obtained above. We first write it somewhat more precisely than we did above. Let e_I be some basis k-blade where I is the index set consisting of k of the first n+1 integers starting with 0. Let I^\perp be the index set obtained by taking the complement of I within \{0,1,...n\}. Then e_I \wedge e_{I^\perp} = \pm I.

To continue with our focus on n-D euclidean PGA: It can be considered as built on the dual Grassmann algebra G* in which the 1-vectors are points and the (n-1)-vectors are planes, and the wedge operation is meet. Consider the standard Grassmann algebra G in which the 1-vectors are points and the (n-1)-vectors are planes, and the wedge operation is join. Let the elements of G be written with superscripts instead of subscripts. So the 1-vectors are written as e^i, etc. Then the dual coordinate map J: G^* \rightarrow G is defined on basis blades by J(e_I) = e^{(I^\perp)}. Extend it by linearity to the whole algebra.

Basic result: J maps a geometric element of G^* to the same geometric element of G, and is coordinate-free. (Hint: to understand how this works, note that the map D: G^* \rightarrow G given by D(e_I) = e^I is essentially the canonical dual basis map that was previously defined above as a self-mapping G^*\rightarrow G^*. So J is essentially obtained by applying the canonical dual basis map twice: e_I \rightarrow e_{I^\perp} \rightarrow e^{I^\perp}, which yields the identity, since the dual of the dual is the identity map. This also explains why J is coordinate-free since the second step undoes any coordinate-dependence introduced by the first.)

Example: if P \in \bigwedge^3 is a point in G^*, then J(P) \in \bigwedge^1 represents the same point in G, etc.

Takeaway: if you’re already calculating I^\perp for the canonical dual basis approach (above), all you need to do to make it coordinate-free is to interpret it as an index in the “other” algebra, i. e., instead of e_{I^\perp}, think e^{(I^\perp)}, and, to be safe, keep track that you’ve moved over to G. Since the main application of this dual map (in my experience) is to calculate the regressive product, this “visit” to the other algebra should be short-lived: a\vee b := J^{-1}(J(a) \wedge J(b)).

Implementation: A single bit in the data structure for your GA element is sufficient to keep track of whether the element lives in G or G^*, and protects against invalid operations (like trying to take the wedge of an element in G and an element in G^*).

Wrap-up: Essentially I have described three versions of the “dual” map. I would also (continue to) recommend that more care be taken with the notation/terminology used here in order to know exactly which of these three maps are meant in a given context. Note that in the below I use G^* as the source algebra since this post is in the context of euclidean PGA for which that is the case (why that is, I’ll address in my next post.). There’s nothing magic about it, if your GA is based on a standard Grassmann algebra you should flip the roles of G and G^* below.

  1. “dual of X”: The grade-reversing dual map G^*\rightarrow G^* that I described above based on the “canonical dual basis”, that is coordinate-dependent. It’s also referred to in previous post as “right complement” (same goes obviously in principle for the “left complement”).
  2. “dual coordinates of X”: Poincare duality, also grade reversing G\rightarrow G^*, coordinate-free, mapping an element in one Grassmann element to the same geometric element in the other Grassmann algebra, and
  3. “orthogonal complement of X” or “polar of X”: P: G^* \rightarrow G^*, P(X) = X I. I don’t use the pseudoscalar inverse here since in a degenerate metric it doesn’t exist. This operation is essential for doing metric geometry in both degenerate and non-degenerate metrics. It’s sometimes called the pole-polar pairing or metric polarity (a term from projective geometry). I’ll have more to say about its usage in euclidean PGA in my next post.

The regressive product: Note that 1) and 2) can be used to calculate the regressive product for any metric but only 2) is coordinate-free. 3) can be used when the metric is non-degenerate. In all cases the regressive product looks like a \vee b = M^{-1}(M(a) \wedge M(b)) where M is one of these 3 maps and \wedge is the native wedge product in G or G^* as the case may be.

Observant readers of my other PGA writings will notice that the “canonical dual basis” in the only new theme here, throwing new light on Poincare duality (see the above discussion). Hopefully now that there’s more people actually programming and using euclidean PGA (hurrah!), these comments can be read with fresh eyes and help clarify some of the new and unfamiliar features of PGA. And perhaps it can contribute more conscious use of the proper terminology:-)

Disclaimer: in the interests of readability I have been quite cavalier about the various \pm sign adjustments that have to be made to make the whole thing work as advertised. I trust that the reader can take it as an exercise to figure out where to insert them.

@tza: if you still have questions on what I’ve written please don’t hesitate to post or send a message – it will naturally help if you can specify as exactly as possible where you lost the thread.

I welcome general feedback – I’ll be updating the SIGGRAPH course notes with this information so I want to make it as clear as I can without sacrificing accuracy.

3 Likes