Hugo Hadfield
Cambridge University PhD student, Signal Processing and Communications Laboratory
Conformal Geometric Algebra with Python
Part 2
Hugo Hadfield 2018
Under the hood
- Clifford is built as flat, dense arrays of blade coefficients
- The sign flips for products are calculated in advance
- The various product functions are then exposed as JITed closures with the Numba JIT (Just in Time) compiler
- It is pretty fast!
Mixing GA and LA
- Dense, flat representation is not the most efficient
- Does make interfacing with linear algebra systems very easy
- Currently many GA algorithms still terminate in a linear algebra solution step
- Need to be able to easily merge between GA and LA
Toy Problem
- Given a set of lines in the world, what point lies closest to all of them in a least squares sense?
C= \min_P \left( -\sum^N_i P \cdot(L_iPL_i) \right)
P^2 = 0
C= \min_P \left( -\sum^N_i P \cdot(L_iPL_i) \right)
P^2 = 0
\therefore L = -\sum^N_i P \cdot(L_iPL_i) + \lambda P^2
Set up Lagrange multiplier:
Differentiation and setting equal to 0 gives an eigen-point problem:
\sum^N_iL_iPL_i = 2\lambda P
How do we solve this? If anyone knows a nice GA way to do it let me know afterwards :)
\sum^N_iL_iPL_i = 2\lambda P
Bring in the linear algebra
Define: as the matrix that when multiplied with a 32x1 coefficient vector is equivalent to the left geometric product with Li
Define: is equivalent to the right geometric product with Li
\underrightarrow{L_i}
\underleftarrow{L_i}
With our new notation we can construct the reflection as a linear matrix operation:
\left( \sum^N_i \underleftarrow{L_i} \hspace{0.1cm} \underrightarrow{L_i} \right) P = 2\lambda P
where P here is a 32x1 vector of coefficients.
We therefore need to find the eigenvectors of the matrix:
M = \left( \sum^N_i \underleftarrow{L_i} \hspace{0.1cm} \underrightarrow{L_i} \right)
The clifford Python package provides us with the tools to solve this easily
from clifford.g3c import *
from clifford.tools.g3c import *
from pyganja import *
# Make a load of lines
line_list = generate_random_object_cluster(10, random_line)
# Allocate some space for our finished matrix
accumulator_matrix = np.zeros((32,32))
# Loop over our lines and construct the matrix
for i in range(len(line_list)):
# Get the line as a left gmt matrix
L_i_l = get_left_gmt_matrix(line_list[i].value)
# Get the line as a right gmt matrix
L_i_r = get_right_gmt_matrix(line_list[i].value)
# Multiply and add
accumulator_matrix += L_i_r@L_i_l
# Raise the matrix to a very high power
power_mat = np.linalg.matrix_power(accumulator_matrix/10,256)
# Get a point that lies on the first line as an approximation to the e-vector
p_start = normalise_n_minus_1(((line_list[0]|eo)*einf*(line_list[0]|eo))(1))
# Apply the matrix
p_end = layout.MultiVector(value=power_mat@p_start)(1)
# Remove any junk that has come along with it
final_point = normalise_n_minus_1((p_end*einf*p_end)(1))
# Visualise
draw(line_list+[final_point], scale=0.1)Non linear optimisation
- Python has a well developed scientific ecosystem
- SciPy (Scientific Python) contains a set of non-linear optimisers
- Clifford is designed to work well with these optimisers
Lets take our toy problem and try non-linear optimisation to solve it this time
C= \min_P \left( -\sum^N_i P \cdot(L_iPL_i) \right)
To enforce P is a valid conformal point we will optimise in the euclidean domain and map into conformal
P = \frac{1}{2}p^2 e_\infty + p + e_0
= \text{up}(p)
from clifford.g3c import *
from clifford.tools.g3 import *
from clifford.tools.g3c import *
from pyganja import *
# Import scipy minimize
from scipy.optimize import minimize
# Define our cost function
def cost_function(p_numpy, line_list):
p = np_to_euc_mv(p_numpy)
P = up(p)
return -sum((P|(l*P*l))[0] for l in line_list)
# Make a load of lines
line_list = generate_random_object_cluster(10,
random_line)
# We will start at the origin
x0 = np.zeros(3)
# Run a non-linear optimiser to find the result
opt_result = minimize(cost_function, x0,
args=tuple([line_list]),
method='L-BFGS-B')
# Get our final conformal point
final_point = up(np_to_euc_mv(opt_result.x))
# Visualise
draw(line_list+[final_point], scale=0.1)Optimal Sphere Fitting
- Another fun toy problem, this is a shameless knock off of Leo Dorst's AGACSE 2018 tutorial
- Given a set of points in 3D space what is the best fitting sphere?
\sum_i^N(P_i\cdot\sigma)^2/\sigma^2
Defining a cost function:
Leads to another eigenproblem, we want the minimal non-negative eigenvalue of:
\frac{1}{N}\sum_i^NP_i\left(P_i\cdot [\hspace{0.2cm}]\right)
Bringing in the linear algebra again, we define the matrix that performs projection to grade g as:
M_g
Which leads to a final matrix as:
\left( \frac{1}{N}\sum_i^N \underrightarrow{P_i} M_0 \underrightarrow{P_i } \right)
from clifford.g3c import *
from clifford.tools.g3 import *
from clifford.tools.g3c import *
from pyganja import *
# Sample a load of points on a sphere with a bit of noise
radius = 2.5
center = 3*e1 + 2.3*e2 + 5*e3
sigma = 0.1
sphere_points = [up((random_euc_mv().normal()*radius + random_euc_mv()*sigma) + center) for i in range(50)]
# Loop over our points and construct the matrix
accumulator_matrix = np.zeros((32,32))
for i in range(len(sphere_points)):
# Get the point as a left gmt matrix
P_i_l = get_left_gmt_matrix(sphere_points[i].value)
# Multiply and add
accumulator_matrix += P_i_l@mask0@P_i_l
# Find the eigenvalues of the matrix
e_vals, e_vecs = np.linalg.eig(accumulator_matrix)
# Find the smallest non negative eigenvalues
min_eval = np.inf
min_eval_index = -1
for i in range(len(e_vals)):
if e_vals[i] < min_eval and e_vals[i] > 0:
min_eval = e_vals[i]
min_eval_index = i
# Get the associated eigenvector and visualise it
best_sphere = layout.MultiVector(value=e_vecs[:,min_eval_index]).normal()
draw([best_sphere]+sphere_points,static=False,scale=0.3)Best fitting circle
- The best fitting sphere minimises the spherical radial error in the data set
- The second best fitting eigenvector is also a sphere that minimises the remaining error
- Best sphere meets seconds best sphere is the best circle?
from clifford.tools.g3 import *
from clifford.tools.g3c import *
from pyganja import *
# Sample a load of points on a circle with a bit of noise
center = 3*e1 + 2.3*e2 + 5*e3
sphere_points = [up(((np.random.randn()*e1 + np.random.randn()*e2).normal()*2.5 + random_euc_mv()*0.01) + center) for i in range(50)]
# Loop over our points and construct the matrix
accumulator_matrix = np.zeros((32,32))
for i in range(len(sphere_points)):
# Get the point as a left gmt matrix
P_i_l = get_left_gmt_matrix(sphere_points[i].value)
# Multiply and add
accumulator_matrix += P_i_l@mask0@P_i_l
# Find the eigenvalues of the matrix
e_vals, e_vecs = np.linalg.eig(accumulator_matrix)
# Find the smallest and second smallest non negative eigenvalues
min_eval = np.inf
min_eval2 = np.inf
min_eval_index = -1
min_eval_index2 = -1
for i in range(len(e_vals)):
this_e_val = e_vals[i]
if this_e_val > 0:
if this_e_val < min_eval:
min_eval2 = min_eval
min_eval = this_e_val
min_eval_index2 = min_eval_index
min_eval_index = i
elif this_e_val < min_eval2:
min_eval2 = this_e_val
min_eval_index2 = i
# Get the associated eigenvectors and visualise
best_sphere = layout.MultiVector(value=mask1@np.real(e_vecs[:, min_eval_index])).normal()
second_best_sphere = layout.MultiVector(value=mask1@np.real(e_vecs[:, min_eval_index2])).normal()
best_circle = meet(best_sphere*I5,second_best_sphere*I5)(3).normal()
draw([best_circle] + sphere_points, browser_window=True, scale=0.1)
Summary
- Clifford allows CGA to mix freely with traditional linear algebra and non-linear optimisation
- CGA can be used to solve interesting and difficult problems in an intuitive way
CGAPython2
By Hugo Hadfield
CGAPython2
A tutorial presentation on Conformal Geometric Algebra and the python clifford package for geometric algebra. Originally presented to faculty and students from Brno University of Technology and Masaryk University. The clifford python package can be found here: https://github.com/pygae/clifford
