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?
Set up Lagrange multiplier:
Differentiation and setting equal to 0 gives an eigen-point problem:
How do we solve this? If anyone knows a nice GA way to do it let me know afterwards :)
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
With our new notation we can construct the reflection as a linear matrix operation:
where P here is a 32x1 vector of coefficients.
We therefore need to find the eigenvectors of the matrix:
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
To enforce P is a valid conformal point we will optimise in the euclidean domain and map into conformal
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?
Defining a cost function:
Leads to another eigenproblem, we want the minimal non-negative eigenvalue of:
Bringing in the linear algebra again, we define the matrix that performs projection to grade g as:
Which leads to a final matrix as:
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