Geometric Algebra

See Wikipedia for an idea of what this is.

New in version 2013.2.

Also see Example usage.


class pymbolic.geometric_algebra.Space(basis=None, metric_matrix=None)

A sequence of names of basis vectors.


A (dims,dims)-shaped matrix, whose (i,j)-th entry represents the inner product of basis vector i and basis vector j.

  • basis – A sequence of names of basis vectors, or an integer (the number of dimensions) to use the default names e0 through eN.
  • metric_matrix – See metric_matrix. If None, the Euclidean metric is assumed.

Return the canonical n-dimensional Euclidean Space.


class pymbolic.geometric_algebra.MultiVector(data, space=None)

An immutable multivector type. Its implementation follows [DFM]. It is pickleable, and not picky about what data is used as coefficients. It supports pymbolic.primitives.Expression objects of course, but it can take just about any other scalar-ish coefficients.


A mapping from a basis vector bitmap (indicating blades) to coefficients. (see [DFM], Chapter 19 for idea and rationale)


See the following literature:

[DFM] L. Dorst, D. Fontijne, and S. Mann, Geometric Algebra for Computer Science: An Object-Oriented Approach to Geometry. Morgan Kaufmann, 2010.

[HS] D. Hestenes and G. Sobczyk, Clifford Algebra to Geometric Calculus: A Unified Language for Mathematics and Physics. Springer, 1987.

The object behaves much like the corresponding sympy.galgebra.GA.MV object in sympy, especially with respect to the supported operators:

Operation Result
A+B Sum of multivectors
A-B Difference of multivectors
A*B Geometric product \(AB\)
A^B Outer product \(A\wedge B\) of multivectors
A|B Inner product \(A\cdot B\) of multivectors
A<<B Left contraction \(A\lrcorner B\) of multivectors
A>>B Right contraction \(A\llcorner B\) of multivectors


Many of the multiplicative operators bind more weakly than even addition. Python’s operator precedence further does not match geometric algebra, which customarily evaluates outer, inner, and then geometric.

In other words: Use parentheses everywhere.

mapper_method = 'map_multivector'

More products


Return the scalar product, as a scalar, not a MultiVector.

Often written \(A*B\).


Return the commutator product.

See (1.1.55) in [HS].

Often written \(A\times B\).


Return self to the integer power other.

Unary operators


Return the multiplicative inverse of the blade self.

Often written \(A^{-1}\).


Return the reverse of self, i.e. the multivector obtained by reversing the order of all component blades.

Often written \(A^\dagger\).


Return the grade involution (see Section 2.9.5 of [DFM]), i.e. all odd-grade blades have their signs flipped.

Often written \(\widehat A\).


Return the dual of self, see (1.2.26) in [HS].

Written \(\widetilde A\) by [HS] and \(A^\ast\) by [DFW].


Return the dual of self, see dual().


Return the pseudoscalar associated with this object’s Space.


MultiVector objects have a truth value corresponding to whether they have any blades with non-zero coefficients. They support testing for (exact) equality.


Remove blades whose coefficient is close to zero relative to the norm of self.

close_to(other, tol=None)

Grade manipulation


Generate all blades in self, optionally only those of a specific grade.


Return a new multivector containing only the blades of grade r.

Often written \(\langle A\rangle_r\).

xproject(r, dtype=None)

If r == 0, return self.project(0).as_scalar(). If r == 1, return self.project(1).as_vector(dtype). Otherwise, return self.project(r).


Return a set of grades occurring in self.


If self only has components of a single grade, return that as an integer. Otherwise, return None.


Extract the odd-grade blades.


Extract the even-grade blades.


New in version 2014.2.


New in version 2014.2.


Return a numpy vector corresponding to the grade-1 MultiVector self.

If self is not grade-1, ValueError is raised.

Helper functions


Return a new MultiVector with coefficients mapped by function f, which takes a single coefficient as input and returns the new coefficient.

  • data

    This may be one of the following:

    • a numpy.ndarray, which will be turned into a grade-1 multivector,
    • a mapping from tuples of basis indices (together indicating a blade, order matters and will be mapped to ‘normalized’ blades) to coefficients,
    • an array as described in data,
    • a scalar–where everything that doesn’t fall into the above cases is viewed as a scalar.
  • space – A Space instance. If None or an integer, get_euclidean_space() is called to obtain a default space with the right number of dimensions for data. Note: dimension guessing only works when a numpy.ndarray is being passed for data.

Example usage

This first example demonstrates how to compute a cross product using MultiVector:

>>> import numpy as np
>>> import pymbolic.geometric_algebra as ga
>>> MV = ga.MultiVector

>>> a = np.array([3.344, 1.2, -0.5])
>>> b = np.array([7.4, 1.1, -2.0])
>>> np.cross(a, b)
array([-1.85  ,  2.988 , -5.2016])

>>> mv_a = MV(a)
>>> mv_b = MV(b)
>>> print -mv_a.I*(mv_a^mv_b)
MV(-1.85*e0 + 2.988*e1 + -5.2016*e2)

This simple example demonstrates how a complex number is simply a special case of a MultiVector:

>>> import numpy as np
>>> import pymbolic.geometric_algebra as ga
>>> MV = ga.MultiVector
>>> sp = ga.Space(metric_matrix=-np.eye(1))
>>> sp
Space(['e0'], array([[-1.]]))

>>> one = MV(1, sp)
>>> one
MultiVector({0: 1}, Space(['e0'], array([[-1.]])))
>>> print one
>>> print one.I
>>> print one.I ** 2

>>> print (3+5j)*(2+3j)/(3j)
>>> print (3+5*one.I)*(2+3*one.I)/(3*one.I)
MV(6.33333333333 + 3.0*e0)

The following test demonstrates the use of the object and shows many useful properties:

def test_geometric_algebra(dims):

    import numpy as np
    from pymbolic.geometric_algebra import MultiVector as MV  # noqa

    vec1 = MV(np.random.randn(dims))
    vec2 = MV(np.random.randn(dims))
    vec3 = MV(np.random.randn(dims))
    vec4 = MV(np.random.randn(dims))
    vec5 = MV(np.random.randn(dims))

    # Fundamental identity
    assert ((vec1 ^ vec2) + (vec1 | vec2)).close_to(vec1*vec2)

    # Antisymmetry
    assert (vec1 ^ vec2 ^ vec3).close_to(- vec2 ^ vec1 ^ vec3)

    vecs = [vec1, vec2, vec3, vec4, vec5]

    if len(vecs) > dims:
        from operator import xor as outer
        assert reduce(outer, vecs).close_to(0)

    assert (vec1.inv()*vec1).close_to(1)
    assert (vec1*vec1.inv()).close_to(1)
    assert ((1/vec1)*vec1).close_to(1)
    assert (vec1/vec1).close_to(1)

    for a, b, c in [
            (vec1, vec2, vec3),
            (vec1*vec2, vec3, vec4),
            (vec1, vec2*vec3, vec4),
            (vec1, vec2, vec3*vec4),
            (vec1, vec2, vec3*vec4*vec5),
            (vec1, vec2*vec1, vec3*vec4*vec5),

        # Associativity
        assert ((a*b)*c).close_to(a*(b*c))
        assert ((a ^ b) ^ c).close_to(a ^ (b ^ c))
        # The inner product is not associative.

        # scalar product
        assert ((c*b).project(0)) .close_to(b.scalar_product(c))
        assert ((c.rev()*b).project(0)) .close_to(b.rev().scalar_product(c))
        assert ((b.rev()*b).project(0)) .close_to(b.norm_squared())

        assert b.norm_squared() >= 0
        assert c.norm_squared() >= 0

        # Cauchy's inequality
        assert b.scalar_product(c) <= abs(b)*abs(c) + 1e-13

        # contractions

        # (3.18) in [DFM]
        assert abs(b.scalar_product(a ^ c) - (b >> a).scalar_product(c)) < 1e-13

        # duality, (3.20) in [DFM]
        assert ((a ^ b) << c) .close_to(a << (b << c))

        # two definitions of the dual agree: (1.2.26) in [HS]
        # and (sec 3.5.3) in [DFW]
        assert (c << c.I.rev()).close_to(c | c.I.rev())

        # inverse
        for div in list(b.gen_blades()) + [vec1, vec1.I]:
            assert (div.inv()*div).close_to(1)
            assert (div*div.inv()).close_to(1)
            assert ((1/div)*div).close_to(1)
            assert (div/div).close_to(1)
            assert ((c/div)*div).close_to(c)
            assert ((c*div)/div).close_to(c)

        # reverse properties (Sec 2.9.5 [DFM])
        assert c.rev().rev() == c
        assert (b ^ c).rev() .close_to((c.rev() ^ b.rev()))

        # dual properties
        # (1.2.26) in [HS]
        assert c.dual() .close_to(c | c.I.rev())
        assert c.dual() .close_to(c*c.I.rev())

        # involution properties (Sec 2.9.5 DFW)
        assert c.invol().invol() == c
        assert (b ^ c).invol() .close_to((b.invol() ^ c.invol()))

        # commutator properties

        # Jacobi identity (1.1.56c) in [HS] or (8.2) in [DFW]
        assert (a.x(b.x(c)) + b.x(c.x(a)) + c.x(a.x(b))).close_to(0)

        # (1.57) in [HS]
        assert a.x(b*c) .close_to(a.x(b)*c + b*a.x(c))