# Conformal Geometric Algebra with Python

Part 2

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

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

# 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