intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Geometric Algebra\n",
    "Below I will try to motivate the case for Geometric Algebra (GA) using rotation of vectors in various spaces and compare to the other most common approaches.\n",
    "\n",
    "![](media/Counterclockwise_rotation.png)\n",
    "\n",
    "(Source: https://commons.wikimedia.org/wiki/File:Counterclockwise_rotation.png)\n",
    "\n",
    "Prerequisites\n",
    "- Rotation matrices: good to know and to understand the motivation but not absolutely necessary\n",
    "- Complex numbers: used to introduce Geometric Algebra rotors\n",
    "\n",
    "TODO List\n",
    "- Add graphics\n",
    "- Add practical examples for PGA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rotating 2D vectors\n",
    "Given a 2D vector $v = \\begin{bmatrix} x \\\\ y \\end{bmatrix} = x e_x + y e_y$ how do we rotate the vector given a rotation angle $\\phi$? Let's look at some different solutions to this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using 2x2 Rotation Matrices\n",
    "A common solution is to use rotation matrices. In the 2D case it will be a 2x2 matrix of the following form\n",
    "\n",
    "\\begin{equation}\n",
    "R(\\phi) = \n",
    "\\begin{bmatrix}\n",
    "cos(\\phi) & -sin(\\phi)\\\\\n",
    "sin(\\phi) & cos(\\phi)\n",
    "\\end{bmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Multiplying the vector and the matrix together we arrive at\n",
    "\n",
    "\\begin{equation}\n",
    "v' = R(\\phi) v = (x cos(\\phi) - y sin(\\phi)) e_x + (x sin(\\phi) + y cos(\\phi)) e_y\n",
    "\\end{equation}\n",
    "\n",
    "Rotations can also be composed by multiplying two or more matrices together to obtain a single matrix\n",
    "\n",
    "\\begin{equation}\n",
    "v' = R_a R_b v = R v$\n",
    "\\end{equation}\n",
    "\n",
    "The downside of this method is that we need to calculate the `2x2` matrix which contains four numbers which seems more than actually necessary. Also while rotation matrices\n",
    "do exist for other dimensions it is usually not easy to come up with them. They are often either memorized or require a lookup."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Complex Numbers\n",
    "Another common tool for 2D rotation are the complex numbers where we introduce an imaginary unit $i$ that squares to $-1$. First note that if we multiply a complex number with the imaginary unit $i$ the complex number will get rotated by 90°. To verify this we multiply $1 + i$ by $i$ the result will be $(1 + i) i = -1 + i$.\n",
    "\n",
    "TODO: Add image of 90° rotation\n",
    "\n",
    "By multiplying by integer powers of $i$ we can rotate by integer multiples of 90°, for example multiplying by $i^3$ will rotate by 270°.\n",
    "\n",
    "TODO: Add image of 90° 180° 270° 360° rotations\n",
    "\n",
    "However in order to create arbitrary rotations phi, we need to form the exponential $e^{i \\phi}$ and multiply by it.\n",
    "The resulting complex number of the exponential in terms of real and imaginary parts is given by [Euler's formula](https://en.wikipedia.org/wiki/Euler%27s_formula)\n",
    "\n",
    "\\begin{equation}\n",
    "e^{i \\phi} = cos(\\phi) + i sin(\\phi)\n",
    "\\end{equation}\n",
    "\n",
    "Now in order to achieve what we wanted, we need to represent our vector $v$ as a complex number $c = x + i y$ (ie. the x coordinate becomes the real part and the y coordinate becomes the imaginary part).\n",
    "Then we can apply the formula from above to rotate by an arbitrary angle $\\phi$:\n",
    "\n",
    "\\begin{equation}\n",
    "c' = c e^{i \\phi} = (x + i y) e^{i \\phi} = (x + i y) (cos(\\phi) + i sin(\\phi)) = (x cos(\\phi) - y sin(\\phi)) + i (x sin(\\phi) + y cos(\\phi))\n",
    "\\end{equation}\n",
    "\n",
    "TODO: Show image of complex plane general rotation\n",
    "\n",
    "If we compare this to what we arrived at using rotation matrices the result is identical. Similarly we can also compose rotations by multiplying two or more complex numbers together\n",
    "to arrive at a single complex number that does both rotations. We now also only need to store two numbers for a rotation. Nevertheless the approach still leaves some questions. The conversion between vectors and complex numbers seems arbitrary. Furthermore this approach only works in 2 dimensions and it's not clear how to generalize it yet."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using 2D Geometric Algebra\n",
    "Usually introductions to GA begin by defining various rules and going over derivations before doing anything useful with them. I will also define some rules but try to get to the interesting stuff more quickly.\n",
    "\n",
    "Like for the standard 2D vector algebra in geometric algebra we have two basis vectors $e_x, e_y$ using which arbitrary vectors $v = x e_x + y e_y$ can be formed. Multiplying two same basis vectors together will result in $+1$ if they are the same (ie. $e_x e_x = 1, e_y e_y = 1$) similarly to how the dot product in standard vector algebra works.\n",
    "\n",
    "What is new is that we can also multiply two different basis vectors and the result will not be zero, but can't be simplified further $e_x e_y = e_{xy}$. $e_xy$ here is not something new, it just stands for the two basis vectors multiplied together as a shorthand way of writing. These elements made up of two basis vectors are called bivectors.\n",
    "\n",
    "Importantly the order of the product matters. A rule is that when you swap the factors of a product of basis vectors you pick up a minus sign, for example. $e_{xy} = e_x e_y = -e_y e_x = -e_{yx}$.\n",
    "\n",
    "This kind of product is called the Geometric Product and it is the central idea of Geometric Algebra. Let's now use these three rules and see what some results are when we use them\n",
    "\n",
    "\\begin{aligned}\n",
    "e_x e_y e_x & = & \\text{(rewrite as shorthand)} \\\\\n",
    "e_{xyx} & = & \\text{(Swap neighbours on the right, pick up minus sign)} \\\\\n",
    "-e_{xxy} & = & \\text{(Multiplying same basis vectors results in 1, e_xx = e_x e_x = 1)} \\\\\n",
    "-e_y &\n",
    "\\end{aligned}\n",
    "\n",
    "Hopefully this gives some idea of how the rules work and how to use them. Now for something more interesting, let's see what happens if we multiply of these together, that is, squaring it:\n",
    "\n",
    "\\begin{equation}\n",
    "e_{xy}^2 = e_{xy} e_{xy} = e_{xyxy} = -e_{xyyx} = -e_{xx} = -1\n",
    "\\end{equation}\n",
    "\n",
    "We can see the square of the bivector $e_{xy}$ is $-1$, the same as the square of the imaginary unit of the complex numbers.\n",
    "\n",
    "\n",
    "Back to the rotation problem, to resolve the disconnect of vectors and complex numbers we make use of the things we just learnt. Vectors are represented as before as $v = x e_x + y e_y$. As mentioned before $e_{xy}$ can be identified as the imaginary unit $i$ hence we present complex numbers as $a + b e_{xy}$ and a rotation by an arbitrary angle $\\phi$ can be performed in the same way as before using Euler's formula\n",
    "\n",
    "\\begin{equation}\n",
    "R(\\phi) = e^{\\phi e_{xy}} = cos(\\phi) + e_{xy} sin(\\phi)\n",
    "\\end{equation}\n",
    "\n",
    "The entity $R$ you get after exponentiating is called a rotor (because it rotates when you multiply by it, d'uh). With this we don't need to perform the arbitrary conversion between vectors and complex numbers, we can simply multiply the vector and the rotor together as is\n",
    "\n",
    "\\begin{equation}\n",
    "v R(\\phi) = (x e_x + y e_y) (cos(\\phi) + e_{xy} sin(\\phi)) = e_x (x cos(\\phi) - y sin(\\phi)) + e_y (x sin(\\phi) + y cos(\\phi))\n",
    "\\end{equation}\n",
    "\n",
    "This result is the same as the rotation formula we got for the rotation matrix and the procedure is pretty much identical to what we did with complex numbers. We can also see that the imaginary unit is in fact a rotation in the XY plane as the bivector $e_{xy}$ is formed by multiplying the two basis vectors together, so in some sense it represents the XY plane.\n",
    "\n",
    "So far it seems that at best we gained some clarity, which is good, but might not justify using it over the well known complex numbers and vectors yet. The big advantage of this approach as we will see in the next\n",
    "sections is that this generalizes to any dimension and does not require any memorization besides the basic geometric algebra rules we learnt above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rotating 3D vectors\n",
    "Now we go up one dimension to three dimensional space. We want to rotate a 3D vector $v = x e_x + y e_y + e_z$ by the angles $\\alpha, \\beta, \\gamma$\n",
    "around the three axes (or planes) of rotation. As we will see things will get a bit more complicated in the usual approaches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using 3x3 Rotation Matrices\n",
    "As in the 2D case we can look up the 3D rotation matrix formula on Wikipedia or some other place\n",
    "\n",
    "\\begin{equation}\n",
    "R(\\alpha, \\beta, \\gamma) =\n",
    "\\begin{bmatrix}\n",
    "cos(\\alpha) cos(\\beta) & cos(\\alpha) sin(\\beta) sin(\\gamma) - sin(\\alpha) cos(\\gamma) & cos(\\alpha) sin(\\beta) cos(\\gamma) + sin(\\alpha) sin(\\gamma) \\\\\n",
    "sin(\\alpha) cos(\\beta) & sin(\\alpha) sin(\\beta) sin(\\gamma) + cos(\\alpha) cos(\\gamma) & sin(\\alpha) sin(\\beta) cos(\\gamma) - cos(\\alpha) sin(\\gamma) \\\\\n",
    "-sin(\\beta) & cos(\\beta) sin(\\gamma) & cos(\\beta) cos(\\gamma)\n",
    "\\end{bmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "Multiplying a 3D vector with it will give us the desired result $v' = R(\\alpha, \\beta, \\gamma) v$. Some people might be able to memorize this formula but I think it is obvious to see that the rotation matrices for higher dimensions\n",
    "will be even harder to come up with. Additionally the `3x3` matrix now contains nine numbers which is way more than the original three we started with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using 3D Geometric Algebra\n",
    "We now have three basis vector $e_x, e_y, e_z$ and as a result three independent bivectors $e_{xy}, e_{xz}, e_{yz}$.\n",
    "Let's see what the bivectors square to using the rules from before:\n",
    "\n",
    "\\begin{aligned}\n",
    "e_{xy}^2 = e_{xyxy} = -e_{xyyx} = -1 \\\\\n",
    "e_{xz}^2 = e_{xzxz} = -e_{xzzx} = -1 \\\\\n",
    "e_{yz}^2 = e_{yzyz} = -e_{yzzy} = -1\n",
    "\\end{aligned}\n",
    "\n",
    "All three of our bivectors square to $-1$. Some readers might be familiar with quaternions which have the same property for their three imaginary units $i, j, k$ and are like the 3D versions of the 2D complex numbers. Previously we demonstrated that $e_{xy}$ can be used to perform a rotation in the XY plane, so one might already suspect what $e_{xz}$ and $e_{yz}$ can be used for.\n",
    "\n",
    "One might expect that the rotor $R_{xz} = e^{\\pi e_{xz}}$ would produce a rotation by 90° ($\\pi$ in degrees) in the XZ plane when multiplying a 3D vector by it.\n",
    "This is almost correct. It turns out that the 2D case was special and the general formula for applying rotors requires a two sided product, also called the sandwich product\n",
    "\n",
    "\\begin{equation}\n",
    "v' = R v \\widetilde{R}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\widetilde{R}$ means reversing the order of all basis vectors appearing in it (eg. $\\widetilde{\\phi e_{xz}} = \\phi e_{zx} = -\\phi e_{xz})$.\n",
    "\n",
    "Another result of this is that the angle in the exponent needs to be half of the desired angle as we are now multiplying by it twice. As stated before the 2D case was special and one could then verify that this sandwich product yields the same results as the one sided product in 2D.\n",
    "\n",
    "So, to put this all together, a 3D rotation in the XZ plane can be generated using\n",
    "\n",
    "\\begin{equation}\n",
    "R_{xz} = e^{\\frac{\\phi}{2} e_{xz}} = cos(\\frac{\\phi}{2}) + e_{xz} sin(\\frac{\\phi}{2})\n",
    "\\end{equation}\n",
    "\n",
    "by sandwiching a 3D vector with it\n",
    "\n",
    "\\begin{equation}\n",
    "v' = R_{xz} v \\widetilde{R}_{xz} = (cos(\\frac{\\phi}{2}) + e_{xz} sin(\\frac{\\phi}{2})) v \\widetilde{(cos(\\frac{\\phi}{2}) + e_{xz} sin(\\frac{\\phi}{2}))} = (cos(\\frac{\\phi}{2}) + e_{xz} sin(\\frac{\\phi}{2})) x (cos(\\frac{\\phi}{2}) - e_{xz} sin(\\frac{\\phi}{2}))\n",
    "\\end{equation}\n",
    "\n",
    "Now to perform arbitrary 3D rotations in multiple planes we can just apply multiple of these rotors, for example a rotation in XY plane then the XZ plane then YZ plane\n",
    "by angles $\\alpha, \\beta, \\gamma$ respectively can be performed using three sandwich products in the following way:\n",
    "\n",
    "\\begin{equation}\n",
    "v' = R_{yz}(\\gamma) R_{xz}(\\beta) R_{xy}(\\alpha) v \\widetilde{R}_{xy}(\\alpha) \\widetilde{R}_{xz}(\\beta) \\widetilde{R}_{yz}(\\gamma)\n",
    "\\end{equation}\n",
    "\n",
    "The three rotors $R_{xy}(\\alpha), R_{xz}(\\beta), R_{yz}(\\gamma)$ can also be combined into a single rotor $R(\\alpha, \\beta, \\gamma)$ by multiplying them together. The same rotation around three axes can then be achieved\n",
    "using a single sandwich product\n",
    "\n",
    "\\begin{equation}\n",
    "v' = R(\\alpha, \\beta, \\gamma) v \\widetilde{R}(\\alpha, \\beta, \\gamma)\n",
    "\\end{equation}\n",
    "\n",
    "As we have three rotors and each rotor contains a unique bivector and a scalar part it turns out the combined rotor will contain four numbers, much less than the nine of the `3x3` rotation matrix.\n",
    "\n",
    "Hopefully it is now possible to see that the Geometric Algebra approach will generalize to any dimension. For example in 4D we would have six planes of rotation and thus six bivectors, and one can generate rotations in them using exactly the same approach with the exponential and sandwich product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Translating and rotating 2D vectors\n",
    "\n",
    "Besides rotation another common geometric operation is translation. In the following sections we will see where Geometric Algebra becomes really useful and provides more than just being insightful, minimal and generalizing well.\n",
    "We're going back to two dimensions to keep things less verbose.We're going back to two dimensions to keep things less verbose. The goal now is to translate a vector $v = x e_x + y e_y$ by an offset $T = t_x e_x + t_y e_y$ and also rotate is as before by an angle $\\phi$ in the three axes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Standard Vectors and Addition\n",
    "\n",
    "The most straightforward way of achieving translation is simply adding the translation vector to the original vector\n",
    "\n",
    "\\begin{equation}\n",
    "v' = v + T = (x + t_x) e_x + (y + t_y) e_y\n",
    "\\end{equation}\n",
    "\n",
    "Now we also want to rotate the vector to arrive at the result. We can use the 2x2 rotation matrix from before to do this\n",
    "\n",
    "\\begin{equation}\n",
    "v'' = R(\\phi) v'\n",
    "\\end{equation}\n",
    "\n",
    "This works, but has several issues. For one, we are using two different operations to represent translations and rotations. Furthermore it also makes combining these transformations (combined rotation and translation) harder. We can't just compose multiple transforms into a single one as we could do for rotation (eg. by multiplying two rotation matrices together we could get a single rotation matrix)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Standard Vectors with homogenous coordinates and 3x3 matrices\n",
    "\n",
    "To fix the issue of translation and rotation being represented in different ways using different operations (vector addition for translation, matrix multiplication for rotation) we want a way to unify them. For this the traditional solution is to use homogenous coordinates. Here we simply pad the 2D vector using a constant $1$ to get a vector with three elements\n",
    "\n",
    "\\begin{equation}\n",
    "v_h = \\begin{bmatrix} v \\\\ 1 \\end{bmatrix} = \\begin{bmatrix} x \\\\ y \\\\ 1 \\end{bmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "Translation can then be done using `3x3` matrices as the factors in the third column will add a term proportional to the constant $1$\n",
    "\n",
    "\\begin{aligned}\n",
    "x' = R_{11} x + R_{12} y + R_{13} 1 \\\\\n",
    "y' = R_{21} x + R_{22} y + R_{23} 1\n",
    "\\end{aligned}\n",
    "\n",
    "The last row of the `3x3` matrix is usually chosen so that the resulting vector will also have a $1$ in the fourth position (ie. $R_{31} = R_{32} = 0, R_{33} = 1$).\n",
    "\n",
    "Unlike the vector addition solution, we can now compose the transformations by multiplying `3x3` matrices together to get a single `3x3` matrix that does the combined rotation and translation. However now we have even more numbers.\n",
    "A `3x3` matrix will contain nine numbers, far more than the minimal amount of three (two for translation, one for rotation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using 2D Projective Geometric Algebra\n",
    "In addition to the previous two basis vectors $e_x, e_y$ we now add a third basis vector $e_0$. The interesting thing about this basis vector is that instead of squaring to $+1$ it will square to $0$. This might seem strange at first but this allows performing translations in the same way we performed rotations using the exponential function and rotors. Remember, when we started with rotations we multiplied by $i$ (or the bivector $e_{xy}$ in GA) which generated a rotation by 90° and later we introduced the exponential with Euler's formula to produce a rotation of an arbitrary angle. Here we will do the same for translation.\n",
    "\n",
    "Another change now is that previously we used $e_x, e_y$ as the two basis vectors to represent our points. Instead, we now use the bivectors $e_{y0}, e_{0x}$ instead and we will see that this seemingly unnatural choice will help us achieve the translation we wanted. We also add a third bivector $e_xy$ to the point representation with coefficient $1$. A point with coordinates $x, y$ is now represented as\n",
    "\n",
    "\\begin{equation}\n",
    "P = x e_{y0} + y e_{0x} + e_{xy}\n",
    "\\end{equation}\n",
    "\n",
    "Now to do the equivalent of the 90° rotation for translation, let's notice that right-multiplying the point with $(1 + e_{0x})$ generates a translation of $1$ in the x direction (ie. increments the coefficient of $e_{0y}$ by $1$)\n",
    "\n",
    "\\begin{aligned}\n",
    "& (x e_{y0} + y e_{0x} + e_{xy}) (1 + e_{0x}) & = \\\\\n",
    "& (x e_{y0} + y e_{0x} + e_{xy}) + e_{xy0x} & = \\\\\n",
    "& (x e_{y0} + y e_{0x} + e_{xy}) + e_{y0} & = \\\\\n",
    "& (x + 1) e_{y0} + y e_{0x} + e_{xy}\n",
    "\\end{aligned}\n",
    "\n",
    "So now we can translate by $1$ in both the x and the y direction (a translation in the y direction would be done by right-multiplying with $(1 + e_{y0})$). As with rotation, to generate a translation by an arbitrary amount we now\n",
    "form the exponential to get a rotor $R = e^{\\frac{t_x}{2} e_{0x} + \\frac{t_y}{2} e_{y0}}$. It turns out here we can't use Euler's formula to get the result of the exponential as it only applies for values that square to $-1$ such as our bivector (or the imaginary unit $i$). Instead the solution is a much simpler formula\n",
    "\n",
    "\\begin{equation}\n",
    "e^{\\frac{t_x}{2} e_{0x} + \\frac{t_y}{2} e_{y0}} = 1 + (\\frac{t_x}{2} e_{0x} + \\frac{t_y}{2} e_{y0})\n",
    "\\end{equation}\n",
    "\n",
    "so we just add a scalar $1$. This is also why there was a one plus in our x-translation by $1$ earlier as $e^{e_{0x}} = 1 + e_{0x}$. Now by multiplying a rotor that does rotation and a rotor that does translation, we can have\n",
    "more general rotors that do both rotation and translation. These are often also called motors. In general they are applied as before using the sandwich product. A 2D motor contains four numbers (a scalar, two translation bivectors $e_{y0}, e_{0x}$ and a rotation bivector $e_{xy}$), very close to the minimum of three required for 2D translation and rotation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications of (Projective) Geometric Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolating motors\n",
    "A cool thing we can now do that we couldn't easily do before is interpolating between two motors.\n",
    "\n",
    "![](media/interpolate_motors.png)\n",
    "\n",
    "Given two $m_1$ and $m_2$ we can interpolate them using the logarithm. Taking the logarithm of a motor allows us to get the term that was in the exponent before we took the exponential. Then we can linearly interpolate the two exponents and afterwards take the exponential again which should give us a nice and smooth interpolation that blends both translation and rotation. Putting those thoughts into a formula we arrive at\n",
    "\n",
    "\\begin{equation}\n",
    "interpolate(m_1, m_2, \\alpha) = e^{\\alpha log(m_1) + (1 - \\alpha) log(m_2)}\n",
    "\\end{equation}\n",
    "\n",
    "with interpolation factor $\\alpha \\in [0, 1]$. This requires us to compute the logarithm of a motor. Fortunately for the 2D case there is a simple formula\n",
    "\n",
    "\\begin{equation}\n",
    "log(m) = \\frac{m}{||m||} = \\frac{m}{\\sqrt{m \\widetilde{m}}}\n",
    "\\end{equation}\n",
    "\n",
    "While in other dimensions we have [a way of computing the exponential of (most interesting) GA multivectors](https://www.euclideanspace.com/maths/algebra/clifford/algebra/functions/exponent/index.htm) the logarithm unfortunately doesn't generalize as easily. Formulas for the most important cases are known though including the 3D case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other important things that were not explicitly mentioned\n",
    "- Basis vectors can square to $[-1, 0, +1]$\n",
    "- Geometric product can be decomposed into an inner an exterior product $a b = a \\cdot b + a \\wedge b$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further study\n",
    "- [bivector.net](https://bivector.net/): Contains GA tools and links to various resources\n",
    "- [CoffeeShop](https://enkimute.github.io/ganja.js/examples/coffeeshop.html): Interactive PGA examples\n",
    "- [Geometric Algebra talk at SIGGRAPH2019](https://youtu.be/tX4H_ctggYo): Great introduction to PGA many motivating examples"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": 3
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}