problem6
Four points are chosen independently and at random on the surface of a sphere (using the uniform distribution). What is the probability that the center of the sphere lies inside the resulting tetrahedron?
This is question A6 on the 1992 Putnam collegiate math competition and is explored brilliantly in Grant Sanderson's awesome 3blue1brown video. I find this question to be a lot of fun for two reasons:
- Despite being one of the more difficult questions on one of the hardest math tests in the world, it's easy to understand what the question is asking.
- At first glance, I don't have the faintest idea of what the answer might be or how I'd figure it out.
This post is completely inspired by Grant's video which contains fascinating insight as to how one might have reached the solution. Before reading any further, I highly recommend watching the video because, as Grant points out, the problem solving process is far more interesting than the result and Grant does an incredible job of illuminating this.
In the problem solving spirit, I wanted to take a stab at the question myself with a slightly different approach, by running a sufficiently large number of random simulations and counting the results.
Note: this has been done before in this excellent Medium post, I'll do things slightly differently in a few places so both posts are well worth the read!
How to check if a point is enclosed within a tetrahedron?¶
With a bit of linear algebra, this turns out to be relatively straightforward. The three dimensional case, in which a point in 3D space is enclosed within a tetrahedron, is described here, but the same method extends just as well to any dimension.
The slightly simpler two dimensional case looks like the following:
Given a triangle's vertices and a point $p$
$v1 = (x_1, y_1)$
$v2 = (x_2, y_2)$
$v3 = (x_3, y_3)$
$p = (x_p, y_p)$
The point $p$ is contained within the triangle if the sign of the determinant of all of the following matrices are the same:
$ D_0 = \begin{pmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{pmatrix} $
$ D_1 = \begin{pmatrix} x_p & y_p & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{pmatrix} $
$ D_2 = \begin{pmatrix} x_1 & y_1 & 1 \\ x_p & y_p & 1 \\ x_3 & y_3 & 1 \end{pmatrix} $
$ D_3 = \begin{pmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_p & y_p & 1 \end{pmatrix} $
With a little help from numpy, this is easy enough to implement in a function that takes a point and a list of vertices and returns whether the sign of the determinant of all of the above matrices are the same.
from copy import copy, deepcopy
from typing import List
import numpy as np
def contains(point: np.ndarray, vertices: List[np.ndarray]) -> bool:
"""Returns whether point is contained within the given vertices
Uses a generalized version of this method
http://steve.hollasch.net/cgindex/geometry/ptintet.html
"""
if len(point) + 1 != len(vertices):
raise ValueError("Required that num vertices = dimension + 1")
for vertex in vertices:
if len(vertex) != len(point):
raise ValueError("Dimension of point and vertices do not match")
# create the matrix D0 above with each vertex stacked vertically
# and the right-most column filled with 1s
d0 = deepcopy(vertices)
d0 = [np.append(vertex, 1) for vertex in d0]
dimension = len(point)
# generate a list of the determinants of the matrices described above
determinants = [np.linalg.det(d0)]
# create each matrix d1...dn and store the determinant
for i in range(0, dimension+1):
d_matrix = copy(d0)
d_matrix[i] = np.append(point, 1)
determinants.append(np.linalg.det(d_matrix))
# return whether the sign of all determinants are the same
return all(np.sign(det) == np.sign(determinants[0]) for det in determinants)
A simple test¶
Let's test this implementation by creating a right triangle with the vertices $(0, 0)$, $(1, 0)$, and $(0, 1)$ and verify that the point $(0.2, 0.2)$ is contained within the triangle while the point $(0.6, 0.6)$ is not.
p1 = np.array([0, 0])
p2 = np.array([1, 0])
p3 = np.array([0, 1])
vertices = [p1, p2, p3]
point = np.array([0.2, 0.2])
assert contains(point, vertices)
point = np.array([0.6, 0.6])
assert not contains(point, vertices)
Generating a random point on the unit sphere¶
Looks good! Next comes the "four points are chosen independently and at random on the surface of a sphere" part of the problem. Random points on the unit circle or unit sphere can be represented as vectors with unit length1. Generating a random vector with unit length simply involves generating a vector of random numbers, then dividing each of the vector's elements by its length so that the length becomes one.
1. More precisely, when talking about length we mean the vector's Euclidean norm.
def random_unit_point(dimension: int) -> np.ndarray:
"""Generate a random point with length 1 for the given dimension
random_unit_point(dimension=2) returns a random point on the unit circle
random_unit_point(dimension=3) returns a random point on the unit sphere
"""
vec = np.array([np.random.randn() for _ in range(dimension)])
return vec / np.linalg.norm(vec)
Helper classes¶
Two classes will help in writing cleaner, more pythonic code for running the simulation. First, a Point class that subclasses numpy.ndarray, but adds the attributes x, y, and z.
Additionally, a Triangle class that takes the three vertices that form the triangle and implements the magic __contains__ method, providing the ability to write statements like:
if point in triangle:
num_points_enclosed += 1
class Point(np.ndarray):
def __new__(cls, input_array):
obj = np.asarray(input_array).view(cls)
return obj
@property
def x(self):
return self[0]
@property
def y(self):
return self[1]
@property
def z(self):
return self[2]
class Triangle:
def __init__(self, p1: np.ndarray, p2: np.ndarray, p3: np.ndarray):
self.p1 = p1
self.p2 = p2
self.p3 = p3
def __contains__(self, point: np.ndarray) -> bool:
return contains(point, [self.p1, self.p2, self.p3])
p1 = Point([0, 0])
p2 = Point([1, 0])
p3 = Point([0, 1])
assert Point([0.2, 0.2]) in Triangle(p1, p2 , p3)
assert Point([0.6, 0.6]) not in Triangle(p1, p2 , p3)
p2 = Point([-1, 0])
p3 = Point([0, -1])
assert Point([-0.2, -0.2]) in Triangle(p1, p2 , p3)
assert Point([-0.2, 0.2]) not in Triangle(p1, p2 , p3)
Visualizing the 2D implementation¶
Here, it'll help to stop and visualize what's been implemented so far to verify that all looks good. It's useful to visualize the problem in 2D by generating three random points on the unit circle, drawing the triangle that it forms, and coloring the triangle blue if the origin is enclosed in the triangle or red in the case that it's not. Visualizing this should provide some confidence that the Point class, Triangle class, and the contains function all look correct.
%matplotlib inline
from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
origin = Point([0, 0])
for _ in range(5):
p1 = random_unit_point(dimension=2)
p2 = random_unit_point(dimension=2)
p3 = random_unit_point(dimension=2)
t = np.linspace(0, np.pi*2,100)
plt.plot(np.cos(t), np.sin(t), linewidth=1)
plt.plot(origin.x, origin.y, 'ro')
color = 'blue' if origin in Triangle(p1, p2, p3) else 'red'
t1 = plt.Polygon([p1, p2, p3], color=color)
plt.gca().add_patch(t1)
plt.show()
Running the simulation¶
So far so good! We're ready to run the simulation by generating 10,000 triangles and counting the number of times that the origin is enclosed within our random triangle.
origin = Point([0, 0])
num_runs = 10**5
times_origin_in_triangle = 0
for _ in range(num_runs):
p1 = random_unit_point(dimension=2)
p2 = random_unit_point(dimension=2)
p3 = random_unit_point(dimension=2)
t = Triangle(p1, p2, p3)
times_origin_in_triangle += origin in t
print(
f"Origin lies within a random triangle on the unit circle "
f"{100 * times_origin_in_triangle / num_runs}% of the time."
)
Simulating the 3D case¶
The 3D case follows closely from the 2D case with a Tetrahedron class that mirrors the Triangle class but with four vertices in three dimensions. With this, we're ready to simulate the creation of 10,000 tetrahedrons and check the number of times that the origin is enclosed.
class Tetrahedron:
def __init__(self, p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray):
self.p1 = p1
self.p2 = p2
self.p3 = p3
self.p4 = p4
def __contains__(self, point: np.ndarray) -> bool:
"""Returns whether point is contained within tetrahedron"""
return contains(point, [self.p1, self.p2, self.p3, self.p4])
p1 = Point([0, 0, 0])
p2 = Point([1, 0, 0])
p3 = Point([0, 1, 0])
p4 = Point([0, 0, 1])
t = Tetrahedron(p1, p2, p3, p4)
assert Point([0.01, 0.01, 0.01]) in t
assert Point([0.5, 0.5, 0.5]) not in t
origin = Point([0, 0, 0])
num_runs = 10**5
num_times_origin_contained = 0
for _ in range(num_runs):
p1 = random_unit_point(dimension=3)
p2 = random_unit_point(dimension=3)
p3 = random_unit_point(dimension=3)
p4 = random_unit_point(dimension=3)
t = Tetrahedron(p1, p2, p3, p4)
num_times_origin_contained += origin in t
print(
f"Origin lies within random tetrahedron on the unit sphere "
f"{100 * num_times_origin_contained / num_runs}% of the time."
)
Additional resources¶
Tada! If this post was interesting at all, I can hardly say enough good things about the entire 3blue1brown channel which does an amazing job of providing intuitive visual explanations for challenging concepts. Of particular relevance to this post is a great explanation on the determinant.
I've also recently discovered Jeremy Kun's awesome blog which has nice primers on different areas of mathematics. Jeremy's book, A Programmer's Introduction to Mathematics, is on my reading list.
Finally, I struggled with linear algebra throughout school until a friend of mine in grad school recommended Gilbert Strang's books and lectures which helped me a ton with gaining an inuitive understanding of linear algebra concepts.