# Geometric Algebra¶

See Wikipedia for an idea of what this is.

Also see Example usage.

## Spaces¶

class pymbolic.geometric_algebra.Space(basis=None, metric_matrix=None)[source]
basis_names

A sequence of names of basis vectors.

metric_matrix

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

dimensions
is_orthogonal
is_euclidean
pymbolic.geometric_algebra.get_euclidean_space(n)[source]

Return the canonical n-dimensional Euclidean Space.

## Multivectors¶

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

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.

data

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

space

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 galgebra.mv.Mv object in galgebra, 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, also read as ‘$$A$$ removed from $$B$$’.

A>>B

Right contraction $$A\llcorner B$$ (|_)of multivectors, also read as ‘$$A$$ without $$B$$’.

Warning

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

scalar_product(other)[source]

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

Often written $$A*B$$.

x(other)[source]

Return the commutator product.

See (1.1.55) in [HS].

Often written $$A\times B$$.

__pow__(other)[source]

Return self to the integer power other.

Unary operators

inv()[source]

Return the multiplicative inverse of the blade self.

Often written $$A^{-1}$$.

rev()[source]

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

Often written $$A^\dagger$$.

invol()[source]

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$$.

dual()[source]

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

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

__inv__()[source]

Return the dual of self, see dual().

norm_squared()[source]
__abs__()[source]
I

Return the pseudoscalar associated with this object’s Space.

Comparisons

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

zap_near_zeros(tol=None)[source]

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

close_to(other, tol=None)[source]

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

project(r)[source]

Often written $$\langle A\rangle_r$$.

xproject(r, dtype=None)[source]

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.

odd()[source]

even()[source]

as_scalar()[source]
as_vector(dtype=None)[source]

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

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

Helper functions

map(f)[source]

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

## 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(
e0 * -1.8499999999999999
+ e1 * 2.9879999999999995
+ e2 * -5.201600000000001)


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)
MV(1)
>>> print(one.I)
MV(e0 * 1)
>>> print(one.I ** 2)
MV(-1.0)

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


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

def test_geometric_algebra(dims):
pytest.importorskip("numpy")

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))