gasparse 0.0.6a0

A python library written entirely in C for Geometric Algebras to deal with sparse multivector arrays

A python library to do computations on compiled C code using sparse representation of multivector arrays.

Installation

Building and installing from source

cd gasparse pip install .

• For usage of this package go to examples

• To understand about compatibility between geometric algebras see compatibility

• An overview of the overloaded operators for geometric algebras is presented in Table of Operators

Why this package?

gasparse is a Geometric Algebra (GA) library which combines low level code evaluation with the python C API to deploy very efficient operations on multivector arrays, operations such that are easily expressed using the python syntax.

This package implements Geometric Algebras. Implements operations between multivector arrays. The fundamental principle behind this package is the use of overloaded operators that call optimized C functions on multivector arrays. The entire algebras are generated by computing the signs and bitmaps maps (the cayley table) for the binary operations. On calling gasparse.GA(p,q,r) the specified algebra \(\mathcal{G}_{p,q,r}\) is generated.

There are three alternatives for dealing with Geometric Algebra kingdon, Clifford and numba they are implemented fully in python. Unlike this package they can implement symbolic computations. Because they are lazy or something, they use numba or JIT compilers and symbolic computations to do the heavy lifting for them.

List of some of the disadvantages of the packages kingdon, Clifford and numba:

• Overhead. Because generating code is expensive. Because JIT compiling is also expensive.

• Code cannot be optimized further (limited by the JIT compilers, JIT compiling only gets you so far…)

• Compare my approach with theirs… (show when my approach is better/worse, benchmarking)

Our approach however is able to compute the entire cayley table for big algebras very efficiently, reducing the single overhead that usually condemns code generation. Our approach leverages the relationship between subalgebras to very efficiently compute signs for the cayley table. Having this option makes the use of bigger geometric algebras more appealing. For this purpose we wrote specialized software that computes operations for this type of multivectors. See Large Algebras to take real advantage of the different available computation modes.

The key features of this package are

• Do not need to optimize code symbolically since operations are implemented in a low level language.

• Leverage sparseness of the multivectors.

• Do not need to use numba or other JIT compilers to speed-up numerical computations.

• Numerical computations on multivector arrays are computed in compiled C code.

• Multiple multivector types for which operations can be dispatched the most efficiently way possible.

• There is compatibility between geometric algebras (can do computations between multivectors of different algebras)

• Very fast computation of cayley tables.

• Can efficiently deal with big geometric algebras without big overhead.

• Can easily scale for other types of data structures (multivector type).

• In futures iterations improve performance by leveraging knowlegde of the C language. (vector calls, assembly)

Another cool feature of this package is that for someone who is already familiar with programming with kingdon Clifford and numba is going to easily adapt to this package’s “syntax”, since most of the python operators are similar.

Examples

A simple example

Operations between multivector arrays are dispatched to compiled C code

>>> import gasparse
>>> vga = gasparse.GA(3) # Intialize a 3D vanilla geometric algebra
>>> x = vga.mvarray([1,2,3],grades=1)
>>> y = vga.mvarray([4,5,6],grades=1)
>>> print(x*y) # The product is executed in low level C code
32 + -3*e12 + -6*e13 + -3*e23

Creating geometric algebras:

import gasparse
vga = gasparse.GA(3)
cga = gasparse.GA(4,1)
cga = gasparse.GA(metric=[-1,1,1,1,1])
ga = gasparse.GA(q=3,r=4)

Creating multivector arrays:

import gasparse
vga = gasparse.GA(3)
locals().update(vga.basis()) # Update the global variables e, e1, e2, e3, e12, e13, e23, e123.
values = [[0.1,1,2,3],[0.4,4,5,6]]
x = vga.mvarray(values,grades=[0,1])
x = vga.mvarray(values,basis=['e1','e3','e123','e12'])
x = vga.mvarray([1,2,3,4,5,6,7,8]) # Consider all basis elements
x = vga.mvarray(values,basis=[1, e2, e123, e23]) # Use the variables to create the multivector

Note that for the last line the basis can be any list of gasparse.mvarray with the restriction that the multivector array must be 0-dimensional.

Using numpy to generate random multivector arrays

We can convert between numpy arrays to multivector arrays and vice versa by using lists as intermidiate data structures. To show an example where we convert a numpy array to a multivector array we generate a random numpy array and then convert it back to a multivector array. The user has to make sure that the innermost dimension has size compatible with the specified grade in ga.mvarray. To get the sizes of the grades the user can use ga.size(grades), as is exemplified in the script bellow. The following script generates 5 random multivectors of grade zero and two of the three dimensional vanilla geometric algebra.

>>> import gasparse
>>> import numpy as np
>>> ga = gasparse.GA(3)
>>> arr = np.random.rand(5,ga.size(1,2)) # innermost dimension must be the the size of grades 1 and 2
>>> print(arr)
[[0.90962674 0.84695676 0.62962863 0.69754318 0.32404308 0.66473111]
 [0.66384851 0.74067395 0.62313971 0.40263883 0.85645313 0.06053186]
 [0.62515404 0.33892925 0.92988035 0.26066636 0.51058016 0.52560483]
 [0.71055042 0.68262854 0.40054357 0.62849844 0.56987662 0.60513613]
 [0.5360391  0.88132078 0.55923661 0.45492674 0.67648109 0.52545563]]
>>> x = ga.mvarray(arr.tolist(),grades=[1,2]) # only accepts lists as input
>>> print(x)
[[0.90962674*e1 + 0.84695676*e2 + 0.62962863*e12 + 0.69754318*e3 + 0.32404308*e13 + 0.66473111*e23],
 [0.66384851*e1 + 0.74067395*e2 + 0.62313971*e12 + 0.40263883*e3 + 0.85645313*e13 + 0.060531857*e23],
 [0.62515404*e1 + 0.33892925*e2 + 0.92988035*e12 + 0.26066636*e3 + 0.51058016*e13 + 0.52560483*e23],
 [0.71055042*e1 + 0.68262854*e2 + 0.40054357*e12 + 0.62849844*e3 + 0.56987662*e13 + 0.60513613*e23],
 [0.5360391*e1 + 0.88132078*e2 + 0.55923661*e12 + 0.45492674*e3 + 0.67648109*e13 + 0.52545563*e23]]

Note that in the above example the basis elements of the multivectors are ordered by bitmaps. In the context of generating random multivectors this is irrelevant. But in other situations it may not be helpfull to have this mapping between lists/numpy arrays and multivector arrays as such we advise to either separate the lists into values of grade one and values of grade two.

import gasparse
import numpy as np
ga = gasparse.GA(3)
arr1 = np.random.rand(5,ga.size(1))
arr2 = np.random.rand(5,ga.size(2))
x = ga.mvarray(arr1.tolist(),grades=1) + ga.mvarray(arr1.tolist(),grades=2)

or using ga.basis()

import gasparse
import numpy as np
ga = gasparse.GA(3)
arr = np.random.rand(5,ga.size(1,2))
basis1 = list(ga.basis(grades=1).values())
basis2 = list(ga.basis(grades=2).values())
x = ga.mvarray(arr.tolist(),basis=basis1+basis2)

Converting gasparse.mvarray to lists

To get a list with the values of the multivectors use the function x.tolist(grades) where grades can be an integer or a list of integers <=p+q+r. If no arguments are given then all grades are considered. Attention: If multivectors have values in grades that are ommited in the arguments then information will be lost. Example of getting lists

>>> import gasparse
>>> ga = gasparse.GA(3)
>>> x = ga.mvarray([[1,1,2,3],[1,4,5,6]],grades=[0,2])
>>> print(x)
[[1 + 1*e12 + 2*e13 + 3*e23],
 [1 + 4*e12 + 5*e13 + 6*e23]]
>>> values,basis = x.tolist(0,2) # returns only grades zero and two
>>> values,basis = x.tolist([0,2]) # returns only grades zero and two
>>> print(values,basis,sep='\n')
[[1.0, 1.0, 2.0, 3.0], [1.0, 4.0, 5.0, 6.0]]
[1, 1*e12, 1*e13, 1*e23]
>>> values,basis = x.tolist() # returns a list for all grades
>>> print(values,basis,sep='\n')
[[1.0, 1.0, 2.0, 3.0], [1.0, 4.0, 5.0, 6.0]]
[1, 1*e12, 1*e13, 1*e23]

Grade projections to the scalar grade

When multivectors are grade projected to the scalar grade (grade zero) the resulting multivector is going to be of type 'scalar'. This enables us to dispatch operations that are way more efficient e.g.

>>> import gasparse
>>> from gasparse import mvarray as mv
>>> ga = gasparse.GA(3)
>>> x = ga.mvarray([[1,1,2,3],[1,4,5,6]],grades=[0,1])
>>> y = x/mv.sqrt(abs((x*~x)(0)))  # normalize the mvarray
>>> y = ~x/(x*~x)(0) # Take the inverse of the mvarray
>>> norm_sq = (x*~x)(0) # Compute the norm square of the mvarray
>>> print(norm_sq.type())
GA(3).mvarray.scalar
>>> print(norm_sq.tolist(0)[0]) # print the values as a list
[[15.0], [78.0]]

Generating Large Geometric Algebras

For large geometric algebras we recomend the user to chose the computation mode 'large'. This computation mode disable the computation of bitmaps (this is done online) and only generate the cayley table for the geometric product, the other products use bitmap comparison to discard certain products. Another reason to use the 'large' computation mode is that computing cayley tables for big algebras while using the default computation mode ('generic') will result in the process to be killed.

>>> import gasparse
>>> import timeit
>>> timeit.timeit(lambda: gasparse.GA(10,compute_mode='large'), number=5)/5
0.001923231399996439
>>> timeit.timeit(lambda: gasparse.GA(10), number=5)/5
0.03866826760004187
>>> timeit.timeit(lambda: gasparse.GA(12,compute_mode='large'), number=5)/5
0.03538742340006138
>>> timeit.timeit(lambda: gasparse.GA(12), number=5)/5
0.47024899120006014
>>> timeit.timeit(lambda: gasparse.GA(15,compute_mode='large'), number=5)/5
2.3766122478000398
>>> timeit.timeit(lambda: gasparse.GA(15), number=5)/5
Killed

NOTE: For algebras with n>=10 the subscripts that correspond to the basis vectors of index 10 and above are represented by symbols rather than numbers. This happens because bitmaps get converted to characters via (char)value + '1'. Thus for indices bigger than 9 the corresponding symbols are the ones followed by 9 in the ASCII table. Concretely 10, 11, 12, 13, 14 and 15 are represented by the symbols ':', ';', '<', '=', '>' and '?' respectively. In a subsquent revision we might consider printing multivectors differently. Also note that the representation in this form makes it impossible to define elements via their basis since e=> or e2: is not valid sintax for variables. However we can use ga.mvarray([1],basis='e=>') and ga.mvarray([1],basis='e2:') to create valid variables for the basis bivectors \(e_{12}\wedge e_{13}\) and \(e_2\wedge e_{10}\) respectively.

Computing with big geometric algebras using compute_mode='large' gives us huge performance benefits with respect to the 'generic' mode concretely we show the difference in performance

>>> import gasparse
>>> import timeit
>>> import numpy as np
>>> gal = gasparse.GA(12,compute_mode='large')
>>> ga = gasparse.GA(12,compute_mode='generic')
>>> arr1 = np.random.rand(10,ga.size(1)).tolist()
>>> arr2 = np.random.rand(10,ga.size(1)).tolist()
>>> xl1 = gal.mvarray(arr1,grades=1,dtype='sparse')
>>> x1 = ga.mvarray(arr1,grades=1,dtype='sparse')
>>> xl2 = gal.mvarray(arr2,grades=1,dtype='sparse')
>>> x2 = ga.mvarray(arr2,grades=1,dtype='sparse')
>>> time_generic = timeit.timeit(lambda: x1*x2, number=5)/5
>>> time_large = timeit.timeit(lambda: xl1*xl2, number=5)/5
>>> print("generic is ", time_generic/time_large, " times slower then large",sep='')
generic is 5.6724215181259625 times slower then large
>>> time_generic = timeit.timeit(lambda: x1+x2, number=5)/5
>>> time_large = timeit.timeit(lambda: xl1+xl2, number=5)/5
>>> print("generic is ", time_generic/time_large, " times slower then large",sep='')
generic is 23.555375285590276 times slower then large
>>> time_generic = timeit.timeit(lambda: x1.prod(), number=5)/5
>>> time_large = timeit.timeit(lambda: xl1.prod(), number=5)/5
>>> print("generic is ", time_generic/time_large, " times slower then large",sep='')
generic is 2.2492217006908293 times slower then large
>>> time_generic = timeit.timeit(lambda: x1.sum(), number=5)/5
>>> time_large = timeit.timeit(lambda: xl1.sum(), number=5)/5
>>> print("generic is ", time_generic/time_large, " times slower then large",sep='')
generic is 12.2753188093821 times slower then large

Note that however mixed algebras operations are still using old technology similar to what is done with 'generic', so don’t expect any performance benefits for mixed algebras operations.

Here are how computation time increases when we construct bigger geometric algebras:

>>> import gasparse
>>> import timeit
>>> arr = [0]*15
>>> for i in range(1,16):
>>>     arr[i-1] = timeit.timeit(lambda: gasparse.GA(i,compute_mode='large'), number=5)/5
>>> print(arr)
[2.722999852267094e-06,
 1.558800067869015e-06,
 1.9256000086897985e-06,
 2.8394002583809198e-06,
 7.089399878168478e-06,
 1.9103799786535092e-05,
 6.513800035463646e-05,
 0.00023918759980006142,
 0.0009205148002365604,
 0.0036732804001076147,
 0.01586245879989292,
 0.06812388460020884,
 0.22016883979995328,
 0.8322477006000554,
 3.5297154857998976]

Compatibility between Geometric Algebras

The user has to be carefull when computing operations between multivectors of different algebras. Two algebras of \(n\) and \(m\) dimension \(n<m\) are compatible if the first \(n\) elements of the metric array/tensor are equal. Similarly we can say that two geometric algebras are compatible if the metric tensors of both geometric algebras fully overlap with one another. To illustrate a context where two algebras are imcompatible consider generating a 3D geometric algebra and an algebra of 4 dimensions where the first basis vector is negative and the other positives. The following scripts shows the error obtained after attempting an operation between multivectors of incompatible algebras

>>> import gasparse
>>> ga1 = gasparse.GA(metric=[1,1,1])
>>> ga2 = gasparse.GA(metric=[-1,1,1,1])
>>> x = ga1.mvarray([1,3,4],grades=1)
>>> y = ga2.mvarray([2,4,7,6],grades=1)
>>> x+y
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: Probably Incompatible Algebras!

Overloaded Operators

Table of Operators

Operation	Expression	Python
Geometric product	\(ab\)	a*b
Inner product	\(a \cdot b\)	a|b
Outer product	\(a \wedge b\)	a^b
Regressive product	\(a \vee b\)	a&b
Divide a by b	\(a/b\)	a/b
Sum a with b	\(a+b\)	a+b
Subtract b from a	\(a-b\)	a-b
Reverse of a	\(a^\dagger\)	~a
Grade projections	\(\langle a\rangle_{1,3}\)	a(1,3)
Dual of a	\(a^* = aI\)	a.dual()
Undual of a	\(a^{-*} = aI^{-1}\)	a.undual()

Division is only available when the second argument is either a ‘scalar’ type multivector array, ‘float’ or ‘int’. The scalar product can be computed using the inner or gemetric product and projection to scalars \(a*b=\langle ab\rangle=\langle a\cdot b \rangle\rightarrow\) (a|b)(0) or (a*b)(0). We can also use lists to project to specified grades a([1,3]). Note that dualization when the pseudoscalar is null, that is \(I^2=0\), is defined via the relationship between the basis vectors as \(e_J^\dagger e_J^* = I\) where \(e_J\) are basis multivectors that span all the geometric algebra. The undual operation is defined as the operation that gives back the initial multivector \((a^*)^{-*} = a\).