Random Point in Ball.tex

\documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage{import}
\usepackage{verbatim}
\import{.}{preamble.tex}

\title{How to Find a Random Point in a High-Dimensional Ball}
\author{James Zhang}
\date{July 23, 2022}

\begin{document}

\maketitle

\section{Introduction}

Hey everyone, I'm James. Last year, I watched some Summer of Math Exposition videos, and I really enjoyed them, so this year, I'm going to make one myself.

One of the videos from last year that caught my attention was \href{https://www.youtube.com/watch?v=4y_nmpv-9lI&list=PLnQX-jgAF5pTkwtUuVpqS5tuWmJ-6ZM-Z&index=6&t=3s}{``The BEST Way to Find a Random Point in a Circle"} by Justin. Coincidentally, two of the courses that I took at university this year mentioned higher-dimensional volumes and geometry. So in this video, I'm going to extend Justin's results and talk about how to find a uniform random point in a high-dimensional ball. As you will see, the results in higher dimensions are quite different from those in 2 dimensions.

First of all, let me clarify what I mean by a high-dimensional ball. A ball is a generalization of a circle. A circle (and its interior) with radius $r$ is the set of all points $(x, y)$ such that $x^2 + y^2 \leq r^2$. An $n$-dimensional ball with radius $r$ is the set of all points in $\R^n$ such that $x_1^2 + \ldots + x_n^2 \leq r^2$.

\section{Rejection Sampling}

The first method that Justin mentioned in his video was rejection sampling. The idea of rejection sampling is that we select a random point from a square, where each coordinate is selected uniformly at random between $-1$ and 1. If the point selected happened to be in the unit circle, we're done. If not, we reject that point and select another one. This generalizes easily to $n$ dimensions, so I'm going to implement it. Here's my code:

\begin{verbatim}
def rejection_sampling(n):
    trials = 0
    while True:
        trials += 1
        vec = np.empty(n)
        for i in range(n):
            vec[i] = (random.random() - 0.5) * 2
        if (np.dot(vec, vec) <= 1):
            return (vec, trials)
\end{verbatim}

The code first enters an infinite loop so that we can repeat the selection process until we get a valid point. For each trial, we generate an array with $n$ entries and initialize each entry with a uniform random number between $-1$ and 1. Since the dot product of a vector with itself is the square of its length, we can use the dot function of numpy to check whether the point is in the ball, and if it is, we return the point. Notice that we also track the number of trials it takes to get a valid point and we return it along with the vector.

\subsection{Python Simulation}

I selected 3141 points for dimensions 1 to 14 and recorded the average number of trials. Here are the results:

\begin{center}
    \includegraphics[scale=0.3]{manim/project/images/Trials vs. Dimension Linear.png}
    \includegraphics[scale=0.3]{manim/project/images/Trials vs. Dimension Log.png}
\end{center}

As you can see, as the dimension increases, the number of trials increases very rapidly. In just 14 dimensions, each point on average already requires more than 10,000 trials. We can imagine that in higher dimensions, it takes even more trials to get a volid point. Even when we switch to log scale for the number of trials, the graph still appears to be concave up. This indicates that the number of trials may grow faster than exponentially with dimension. How did that happen?

\subsection{Volume of the Unit Ball}

Let $X$ be the number of trials that we need to select a valid point in the ball. Since each trial is independent, and we repeat until we succeed, $X$ is a geometric random variable.

Let $\beta_n$ be the $n$-dimensional volume of the $n$-dimensional unit ball. Let $\kappa_n$ be the volume of the $n$-dimensional box that we sample from, where $\kappa$ stands for ``cube''. Since each component is randomly selected between $-1$ and 1, the box has side length 2, so its $n$-dimensional volume is $2^n$.

For any particular trial, the probability that the point we selected is in the ball is equal to the ratio between the volume of the ball and the volume of the box. Therefore, to find this probability, we need to find $\beta_n$, the volume of the $n$-dimensional unit ball.

% Reference: Vector Calculus, Linear Algebra, and Differential Forms: A Unified Approach by Hubbard & Hubbard

We will set up an integral to figure out $\beta_n$. First, we divide the ball into small slices, where each slice is approximately a small cylinder. We vary the $x_1$ coordinate from $-1$ to 1 and add up the volumes of these small cylinders. The height of each cylinder is $dx_1$. The base is an $(n-1)$-dimensional ball. In the special case of three dimensions, the base of each slice is a circle.

Since the radius of the $n$-dimensional ball is 1, and this side has length $x_1$, the radius of the base is $r = \sqrt{1 - x_1^2}$ by the Pythagorean Theorem. Since volume scales appropriately according to dimension, the $(n-1)$-dimensional volume of the base is $r^{n-1} \beta_{n-1}$, where $\beta_{n-1}$ is the volume of the $(n-1)$-dimensional unit ball. Thus, the $n$-dimensional volume of this small slice is base times height, or $dV = r^{n-1} \beta_{n-1} \, dx_1$.

Setting up the integral, we have
\[
    \beta_n = \int dV = \int_{-1}^1 r^{n-1} \beta_{n-1} \, dx_1 = \beta_{n-1} \int_{-1}^1 (1 - x_1^2)^{\frac{n-1}{2}} \, dx_1.
\]

We pulled out $\beta_{n-1}$ from the integral because it is just a constant. It remains to figure out the integral
\[
    c_n = \int_{-1}^1 (1 - t^2)^{\frac{n-1}{2}} \, dt
\]
because then we have the recurrence relation $\beta_n = c_n \beta_{n - 1}$.

Trig substitution with $t = \sin(\theta)$ gives
\begin{align*}
    c_n     & = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos^{n}(\theta) \, d\theta   \\
    c_{n-2} & = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos^{n-2}(\theta) \, d\theta
\end{align*}

Integration by parts for $c_n$, with $u = \cos^{n-1}(\theta)$ and $dv = \cos(\theta) \, d\theta$ gives

\[
    c_n = \cos^{n-1}(\theta) \sin(\theta) \bigg|_{-\frac{\pi}{2}}^{\frac{\pi}{2}} + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \sin^2(\theta) (n-1) \cos^{n-2}(\theta) \, d\theta
\]

The first term on the right hand side is 0, so
\[
    c_n = (n-1) \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \sin^2(\theta) \cos^{n-2}(\theta) \, d\theta
\]

Recall that $\sin^2(\theta) = 1 - \cos^2(\theta)$, so
\[
    c_n = (n-1) \bracks{\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos^{n-2}(\theta) \, d\theta - \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos^n(\theta) \, d\theta} = (n-1) c_{n-2} - (n-1) c_n
\]

Rearranging, we get
\[
    c_n = \frac{n-1}{n}c_{n-2}.
\]

The base cases are $c_1 = \int_{-1}^1 (1 - t^2)^0 \, dt = 2$ (by simple integration) and $c_2 = \int_{-1}^1 (1 - t^2)^{\frac{1}{2}} \, dt = \frac{\pi}{2}$ (area of a semicircle). A quick check from geometry: we have $\beta_0 = 1$ (a single point has 0-dimensional volume 1), $\beta_1 = 2$, and $\beta_2 = \pi$.

\ignore{For video: Recall that $\beta_n = c_n \beta_{n - 1}$. If we go from $n-1$ dimensions to $n$ dimensions, $\kappa$ increases by a factor of 2, whereas $\beta$ is multiplied by a factor of $c_n$. Therefore, the ratio between $\beta$ and $\kappa$ gets multiplied by a factor of $\frac{c_n}{2}$ each time we raise a dimension.}

Recall that $\beta_n = c_n \beta_{n - 1}$. In other words, the ratio between $\beta_n$ and $\beta_{n - 1}$ is $c_n$. If $c_n$ were a constant that is less than 2, say 1.5, then every time we raise a dimension, we double the volume of the box ($\kappa_{n} = 2 \kappa_{n - 1}$), but the volume of the ball only increases by a factor of 1.5. The ratio between the volume of the ball and the volume of the box ($\frac{\beta_n}{\kappa_n}$) will thus decrease by a factor of $\frac{3}{4}$. Since the probability that a particular trial succeeds is equal to that ratio, the probability will decay exponentially by a factor of $\frac{3}{4}$ each time.

However, in reality, $c_n = \frac{n-1}{n}c_{n-2}$, and since $\frac{n-1}{n} < 1$, this means that the multiplier $c_n$ is not a constant: it decreases every time! In fact, there will be an $n$ such that $c_n < 1$, and from that point on, the volume of the ball decreases as we increase dimension. (Here, we're talking about the numerical value of the volume, ignoring the units.)

Therefore, the ratio and thus the probability decays faster than exponentially. Since the expected value of a $\mathrm{Geom}(p)$ random variable is $\frac{1}{p}$, this implies that the expected number of trials needed to sample a valid point increases faster than exponentially.

\section{Sphere and Distance Method}

Another method to select a random point from a circle that Justin mentioned was to turn to polar coordinates. We first choose $\theta$ uniformly, then choose $r$ so that the points $(r, \theta)$ are distributed uniformly in the circle.

We can try to apply this idea to higher dimensions. In three dimensions, we can try to choose $r, \theta, \vphi$ in some way such that the points represented by $(r, \theta, \vphi)$ in spherical coordinates are distributed uniformly in a three-dimensional ball. However, we can't choose $\vphi$ uniformly because when $\vphi$ is close to $\frac{\pi}{2}$, that is, we're choosing points near the equator, we expect to see way more points than when $\vphi$ is close to 0 or $\pi$ (or near the poles). There are higher-dimensional analogs of the spherical coordinates in higher dimensions, but we can imagine that as the dimensions go up, it would be increasingly difficult to keep track of the distributions of each variable.

Instead, we take a similar approach: First choose a direction uniformly at random, and then choose a random number between 0 and 1 to be the distance from the origin.

Choosing a random direction in $n$-dimensional space is the same as choosing a random point on the $n$-dimensional unit sphere. A sphere is the boundary of a ball, and points on the unit sphere satisfy $x_1^2 + \ldots + x_n^2 = 1$.

So now, our task is divided into two parts: Choosing a random point on the unit sphere and choosing a distance from the origin. We'll focus on the first part for now and worry about choosing a distance later.

\subsection{Random Point on Sphere}

\subsubsection{Failed Attempt}

\ignore{For video: We choose a random vector from a box with side length 2 where each component is chosen uniformly from -1 to 1.}

To select a random point on a sphere, we might try the following method: We choose a random vector from the box $[-1, 1]^n$ where each component is chosen uniformly in $[-1, 1]$. Then we normalize this vector by dividing away its length. The resulting vector is a unit vector, so it lies on the unit sphere.

This seems like a good approach, but after some thought, we can see that the random points generated in this way aren't distributed uniformly on the sphere. Intuitively, It is more likely to choose a point near the diagonals of the box than near the center of the faces. Let's look at a two-dimensional example:

% Draw picture of square

The distance from the origin to the upper-right corner is $\sqrt{2}$, whereas the distance from the origin to the upper edge is 1. Consider these two cones, which are just triangles in 2D. If a random vector falls into the green cone, it will be scaled to a point on this green arc. If a random vector falls into the blue cone, it will be scaled to a point on this blue arc. The cones are drawn such that the blue and the green arcs have the same length. Then these two cones are similar, and the ratio between their areas is $\parens{\frac{\sqrt{2}}{1}}^2$.

For $n$ dimensions, the diagonal distance will be $\sqrt{n}$ times larger, and the ratio between the volumes of the cones will be $\parens{\frac{\sqrt{n}}{1}}^n$. This again grows faster than exponentially, which means that the probability that a random point lands on one region of the sphere is more than exponentially bigger than the probability that a random point lands on some other region of the sphere with the same area.

\subsubsection{Rotational Invariance and the Multivariate Normal Distribution}

The problem with the box is that it is not rotationally invariant: If we rotate the box by a certain amount, the rotated box doesn't look the same as the original. Meanwhile, the ball and the sphere are all rotationally invariant: However we rotate them, they look the same.

Every rotation in $\R^n$ can be characterized by an orthogonal matrix. An orthogonal matrix is a matrix whose columns (and rows) form an orthonormal basis. (Note: All rotations are isometries, but not all isometries are rotations.) For instance, in $\R^2$, we can express a rotation like this:
\[
    Q = \begin{bmatrix}
        \cos(\theta) & -\sin(\theta) \\
        \sin(\theta) & \cos(\theta)
    \end{bmatrix}
\]

The columns $\xymatrix{\cos(\theta)}{\sin(\theta)}$ and $\xymatrix{-\sin(\theta)}{\cos(\theta)}$ are orthonormal because the length of both vectors are
\[
    \sqrt{\cos^2(\theta) + \sin^2(\theta)} = 1
\]
and the dot product between the two vectors is 0. Similarly, the rows of $Q$ are also orthonormal.

% Manim rotation animation

To figure out a valid way to choose a point uniformly on the sphere, we need to think outside the box. Literally. Introducing the multivariate normal distribution. As you might know, the probability density function of a normal random variable is shaped like a bell curve. Here, we're using the standard multivariate normal distribution, denoted $\mathcal{N}(\mathbf{0}, I)$, which is a vector of independent random variables, each of which follows the standard normal distribution $\mathcal{N}(0, 1)$. (Here I'm using ``normal'' and ``Gaussian'' interchangeably.)

You might also know that any linear combination of independent normal random variables is still a normal random variable. We can use this fact to show that the standard multivariate normal distribution is rotationally invariant:

Suppose $Q$ is the orthogonal matrix representing an arbitrary rotation in $\R^n$, and
$X = \begin{bmatrix}
        X_1 \\ \vdots \\ X_n
    \end{bmatrix}$ is a standard multivariate normal. Then each $X_i$ is a Gaussian with mean 0 and variance 1.

Let $q_{1,1}, \ldots, q_{1,n}$ be the first row of $Q$. Then the first element of $Y = QX$ is $Y_1 = q_{1,1} X_1 + \ldots + q_{1,n} X_n$. Since $Q$ is orthogonal, $q_{1,1}^2 + \ldots + q_{1,n}^2 = 1$.

Since the $X_i$'s are independent Gaussian random variables, the linear combination is still a Gaussian random variable. Moreover, by linearity of expectation, the mean of $Y_1$ is $E[Y_1] = q_{1,1} E[X_1] + \ldots + q_{1,n} E[X_n] = 0$, and by independence, the variance of $Y_1$ is $\Var(Y_1) = q_{1,1}^2 \Var(X_1) + \ldots + q_{1,n}^2 \Var(X_n) = q_{1,1}^2 + \ldots + q_{1,n}^2 = 1$.

We have shown that $Y_1$ is a standard normal random variable. This is the same for each component of $Y = QX$. A quick calculation of covariance shows that the components of $Y$ are independent:
\begin{align*}
    \mathrm{Cov}(Y_1, Y_2)
     & = E[Y_1 Y_2] - E[Y_1] E[Y_2]                                                                                            \\
     & = E[(q_{1,1} X_1 + \ldots + q_{1,n} X_n)(q_{2,1} X_1 + \ldots + q_{2,n} X_n)]                                           \\
     & = q_{1,1} q_{2,1} E[X_1^2] + \ldots + q_{1,n} q_{2,n} E[X_n^2] + \sum_{i,j \leq n, i \neq j} q_{1,i} q_{2,j} E[X_i X_j] \\
     & = E[X_1^2] (q_{1,1} q_{2,1} + \ldots + q_{1,n} q_{2,n}) + \sum_{i,j \leq n, i \neq j} q_{1,i} q_{2,j} E[X_i] E[X_j]     \\
     & = 0
\end{align*}

So $Y$ is still a standard multivariate normal. Thus $X$ is rotationally invariant.

Here's a graph of the joint density function of a 2-dimensional standard multivariate normal random variable. We can see from the density function that the distribution is indeed symmetric around the origin.

\begin{center}
    \includegraphics[scale=0.3]{manim/project/images/Multivariate Normal PDF.png}
\end{center}

Now, we can first generate a random vector that follows the multivariate normal distribution and normalize it, and it will give a point on the unit sphere selected uniformly at random. To see why it is uniform, choose two different small patches on the sphere, $A$ and $B$, that are congruent to each other. Suppose that $p$ is the probability that a random unit vector falls in $A$. Then we can rotate the sphere so that $A$ coincides with $B$. Thus the probability that a random vector falls in $B$ after the rotation is $p$. Because the multivariate normal distribution is rotationally invariant, the distribution of points after the rotation is the same as the distribution before the rotation, so the probability that a random unit vector falls in $A$ after the rotation is also $p$. Since every region on the sphere can be considered as being made up by these small patches, the probability that a chosen point falls in that region is proportional to the area of the region.

There is still the problem of how to generate a Gaussian random variable from a uniform random variable, since the CDF of a Gaussian random variable does not have a closed elementary form. This is a good topic for a separate video, and I'm just going to rely on the python library to do this, but for those who are interested, there is an algorithm called the Box-Muller Transform which can do that. The key feature of that algorithm is that it runs in polynomial time.

\subsection{Choosing the Distance}

After finding a uniform random point on the sphere, it's now time to choose the distance from the origin. Here, we use the same method as Justin used in his video. Let $X$ be the random variable that denotes the distance from the origin. We first find the cumulative distribution function (CDF) $F_X$, where $F_X(r)$ is the probability that $X \leq r$.

The set of all points with distance to the origin less than or equal to $r$ form a ball of radius $r$. The probability that $X \leq r$ is thus the ratio between the volume of a ball of radius $r$ and that of the unit ball. Since volumes scale appropriately according to dimension, that ratio is $\parens{\frac{r}{1}}^n$. Thus, $F_X(r) = r^n$.

Notice that for a given $r \in (0, 1)$, as $n$ increases, the probability that $X \leq r$ gets smaller and smaller. This implies that as $n$ increases, most of the volume of the ball concentrates near the surface of the ball. For instance, when $n = 500$, the probability that a random point in the unit ball is at least 0.99 units away from the origin is $99.3\%$.

Now we shall use the Inverse Transform Sampling technique that Justin mentioned: If $U$ is a uniform random variable, then $X$ has the same distribution as $F_X^{-1}(U)$. In our case, we have $X = \sqrt[n]{U}$.

\subsection{Python Simulation}

Here's my python code that generates random points in the unit ball using this process.

\begin{verbatim}
    def random_unit_vector(n):
        vec = np.empty(n)
        for i in range(n):
            vec[i] = np.random.normal(0, 1)
        length = np.linalg.norm(vec)
        return vec / length
    
    def normal_sampling(n):
        vec = random_unit_vector(n)
        r = (random.random()) ** (1.0 / n)
        return vec * r
\end{verbatim}

To choose a random unit vector, we first create a vector whose components are sampled from the standard normal distribution. Then we normalize the vector. To find a uniform random point in the ball, we take that random unit vector and scale it by the $n$th root of a uniform random number between 0 and 1.

Even when $n$ is big, the program still finds random points in an $n$-dimensional unit ball very quickly, much more quickly than rejection sampling. Unfortunately, I cannot plot higher-dimensional points in our 3D world, so here's a demonstration of this algorithm in 2 and 3 dimensions.

% Manim animation for 2 and 3 dimensions

\section{Summary}

To recap, we first tried to use rejection sampling to select random points in the unit ball. However, after a simulation and a calculation of the volume of an $n$-dimensional ball, we discovered that the number of trials needed to get a valid point increases faster than exponentially as dimension increases.

We then tried to use spherical coordinates, but decided that it was too difficult to determine the distribution of each variable. Instead, we turned to the method of selecting a random direction and a random distance.

The first task is equivalent to selecting a uniform random point on a sphere. We tried to select points in the box and then normalize them, but that did not result in a uniform distribution. To tackle this problem, we had to think outside the box and use a distribution that is rotationally invariant: the standard multivariate normal distribution.

The second task generalizes easily from the two-dimensional case. We calculated the CDF of the distance from the origin and then used inverse transform sampling.

\section{Reflection}

I personally find the topic of higher-dimensional geometry very intriguing, because many results in higher dimensions are surprising and counterintuitive. I was amazed when I first learned how quickly the volume of the ball decays relative to the box that surrounds it.

% There are many other seemingly strange results which I didn't have time to cover, such as almost all the volume of the high-dimensional ball is concentrated near the equator.

On the other hand, making this video is a rewarding experience for me. Although I learned about this material in class, I had to organize the key points so that I can explain them to you in a relatively short video. In addition, I took up the challenge to teach myself Manim. I was very impressed by the neat style and the cool animations in 3b1b's videos, along with the excellent explanations, so I was determined to learn to use Manim. Despite being only a beginner in this tool, I already feel the power of abstraction compared to simply using PowerPoint.

I hope this was a pleasant journey, and I hope to see you again in next year's Summer of Math Exposition.

\section{References}

\begin{itemize}
    \item Justin's video: \href{https://www.youtube.com/watch?v=4y_nmpv-9lI&list=PLnQX-jgAF5pTkwtUuVpqS5tuWmJ-6ZM-Z&index=6&t=3s}{``The BEST Way to Find a Random Point in a Circle"}
    \item \textit{Vector Calculus, Linear Algebra, and Differential Forms: A Unified Approach} by John Hubbard and Barbara Burke Hubbard
    \item \textit{Foundations of Data Science} by Avrim Blum, John Hopcroft, and Ravindran Kannan
    \item \textit{Introduction to Probability} by David F. Anderson, Timo Sepp\"{a}l\"{a}inen, and Benedek Valk\'{o}
    \item \href{https://scipython.com/blog/visualizing-the-bivariate-gaussian-distribution/}{This website} helped me generate the graph of the PDF of the standard multivariate normal
    \item Manim, numpy, and matplotlib documentations and tutorials
\end{itemize}

\end{document}