# Least Squares Ellipsoid Fit

## September 20, 2020

In one of my previous posts, I demonstrated how to fit a sphere using the least squares method. In this post I’ll show how you can also fit an ellipsoid using a least squares fit. It ends up being a bit simpler than the sphere. I’ll also include a simple Python example to perform least square ellipsoid fits.

# Mathematical background

The equation of a three dimensional ellipsoid can be described as

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$$

where $$x$$, $$y$$, and $$z$$ are the cartesian coordinates of some dataset. Now let’s rearrange the equation to become

$$\beta_a x^2 + \beta_b y^2 + \beta_c z^2 = 1$$

where $$\beta_a = 1 / a^2$$, $$\beta_b = 1 / b^2$$, and $$\beta_c = 1/c^2$$. Now we can express the problem in matrix form as the following linear system of equations as

$$\mathbf{A} \mathbf{B} = \mathbf{O}$$

where $$\mathbf{A}$$ is a matrix of our squared data components, $$\mathbf{B}$$ is a vector of $$\beta$$ parameters of the ellipsoid that we intend to solve for, and $$\mathbf{O}$$ is a vector of ones. To be explicit we have

$$\begin{bmatrix} x_{1}^2 & y_{1}^2 & z_{1}^2 \\ x_{2}^2 & y_{2}^2 & z_{2}^2 \\ \vdots & \vdots & \vdots \\ x_{n}^2 & y_{n}^2 & z_{n}^2 \\ \end{bmatrix} \begin{bmatrix} \beta_a \\ \beta_b \\ \beta_c \\ \end{bmatrix} = \begin{bmatrix} 1.0 \\ 1.0 \\ \vdots \\ 1.0 \end{bmatrix}$$

for $$n$$ number of data points.

Now a least squares problem is a simple optimization which finds the model parameters which minimize the L2 norm between the model and the data. This can be expressed as

$$\mathrm{arg\,min}_{\mathbf{B}} (\mathbf{A} \mathbf{B} - \mathbf{O})^2$$

and it happens to have a solution of

$$\mathbf{B} = (\mathbf{A}^\mathrm{T}\mathbf{A})^{-1} \mathbf{A}^\mathrm{T} \mathbf{O}$$

which is found by setting the first derivative of the linear equation to zero.

Once $$\mathbf{B}$$ is solved for, it is only a simple rearrangement to solve for $$a$$, $$b$$, and $$c$$.

$$a = \sqrt{\frac{1}{\beta_a}} \\ b = \sqrt{\frac{1}{\beta_b}} \\ c = \sqrt{\frac{1}{\beta_c}} \\$$

# Example ellipsoid data

Let’s generate data for an ellipsoid where $$a = 1.2$$, $$b = 0.2$$, and $$c = 0.9$$. This is easiest using parametric equation as

$$x = a \cos u \sin v \\ y = b \sin u \sin v \\ z = c \cos v$$

where $$u \in [0, 2\pi]$$ and $$v \in [0, \pi]$$.

We can quickly generate the cartesian points in Python for this ellipsoid using the following code.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

a = 1.2
b = 0.2
c = 0.9

u = np.linspace(0., np.pi*2., 20)
v = np.linspace(0., np.pi, 20)
u, v = np.meshgrid(u,v)

x = a*np.cos(u)*np.sin(v)
y = b*np.sin(u)*np.sin(v)
z = c*np.cos(v)

# turn this data into 1d arrays
x = x.flatten()
y = y.flatten()
z = z.flatten()

#   3D plot of ellipsoid
fig = plt.figure()
ax.scatter(x, y, z, zdir='z', s=20, c='b',rasterized=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()


Which gives us the following data and plot in cartesian coordinates. # Python code to fit ellipsoid

Now in order to fit an ellipsoid to the above data, we need to create the matrix $$\mathbf{A}$$ and vector $$\mathbf{O}$$. Then we’ll use numpy’s least squares solver to find the parameters of the ellipsoid. Lastly we’ll adjust the parameters such that our ellipse is in standard form.

# build our regression matrix
A = np.array([x**2, y**2, z**2]).T

# vector of ones
O = np.ones(len(x))

# least squares solver
B, resids, rank, s = np.linalg.lstsq(A, O)

# solving for a, b, c
a_ls = np.sqrt(1.0/B)
b_ls = np.sqrt(1.0/B)
c_ls = np.sqrt(1.0/B)

print(a_ls, b_ls, c_ls)


Running all of the code together will give you a_ls=1.1999999999999995, b_ls=0.19999999999999993 and c_ls=0.9000000000000006. This is double precision error on our original ellipsoid of $$a = 1.2$$, $$b = 0.2$$, and $$c = 0.9$$. Not too bad! Hopefully you can use this to formulate other least squares problems, or ellipsoids in higher dimensions.