Determining whether a large series of 3D points all line on a plane

Question

TL;DR: For a large series of precise 3D coordinates that describe a real-world orbit, how can we determine if they all lie exactly on a plane?

The Problem

I've used NASA's highly precise SPICE toolkit to calculate the hourly position of the moon relative to the Earth in Cartesian coordinates. The frame of reference for these points is the J2000 frame: Earth's equator (x and y), with the z-axis running orthogonally through the poles.

We know the moon's orbit does not lie on the equatorial plane. The Earth-Moon orbit is typically described as a plane that is 5.14 degrees inclined from the equatorial plane.

(This is a 28-day period from Jan. 2018, hence the small overlap since the sidereal month is 27.32 days.)

What I want to determine is whether this orbit is perfectly "flat" -- whether the Moon orbits the Earth on a perfect plane rotated from the terrestrial equator -- or whether there is even a slight variance on the z axis of this approximate Earth-Moon orbital plane.

It's certainly very close to a plane, as we'd expect, as you can see from a year's worth of z coordinates from the Earth's frame of reference. But orbit's are hugely complicated given all the possible aberrations, which the SPICE toolkit usually accounts for.

One idea I had was to choose three points in the orbit, spaced about 9 days apart, calculate the plane for those three points, and then calculate the distance from that plane for each point. This seems inelegant and possibly biased by which three points I choose, though in theory all I need to know is whether the distance of each point to the 3-point plane is zero, regardless of which three points define the plane.

It seems to me there is a more elegant solution than this brute-force approach?

Answer 1

Let's say that you have position vectors $\{p_i:i\in I\}\subset \mathbb R^3$. If they were coplanar, you would have a vector $v\in\mathbb R^3$ such that $v\cdot p_i=0$ for every $i\in I$. Since this does not necessarily exist, what you can do instead is minimize the quadratic function $$ f(v) = \sum_{i\in I} (v\cdot p_i)^2 $$ over $v\in S^2$.

This gives a good candidate for the "average orbit normal".

Edit 1: you have to subtract the barycenter of the points $p_i$ first, if you want the problem to be translation invariant.

Edit 2: the gradient of $f$ can be computed easily and is $$ \nabla f(v) = 2\sum_{i\in I} (v\cdot p_i)p_i = 2\left(\sum_{i\in I}p_i\otimes p_i\right)v. $$ With Lagrange multipliers, you want to find $v\in S^2$ and $\lambda\in\mathbb R$ such that $\nabla f(v)=2\lambda v$, because the constraint is $v\cdot v=1$. This is equivalent to finding the eigenvalues of the $3\times3$ matrix $\sum_{i\in I}p_i\otimes p_i$.

Edit 3: alternatively, call $P=\sum_{i\in I}p_i\otimes p_i$ and notice that $f(v)=v^TPv$. Therefore the minimum on $S^2$ is achieved at the eigenvectors corresponding to the lowest eigenvalue of $P$.

Edit 4: here is a picture of some points, the plane, the unit normal, and the segments indicating the distances of the points to the plane:

Edit 5: here is some Python code performing the computation

import numpy as np

# 120 random points in 3D
p = np.random.randn(120, 3)

# compute the svd
u, s, vh = np.linalg.svd(p, full_matrices=False)

# the best normal is the last row of vh
n = vh[-1]

# the sum of the squared distances is the square of the last singular value
sum_of_squared_distances = s[-1]**2

# verify the last statement
np.allclose((p@n) @ (p@n), s[-1]**2)