Uniform Distributions Part 2: Beware of Transformations

Ok, now getting to the post that I actually meant to do to start the series on uniformity. The previous post was more about the thrilling adventure of bringing together a bunch of disparate math tools to resolve a seemingly simple question about uniform random variables. It was, ironically, not very visual, but it still in the spirit of the blog, to encourage the attitude that experiment, exploration, and trying different approaches to things is as much a part of math as it is for other scientific endeavors. 

Today we bring something much more visual and play around a lot with random number generation. It is a modernization and generalization of this post (it was the featured post on NestedTori for a looong time). In some sense, this might have more fun submission SoME, but given the time constraints (and the admitted surprise of it being earlier this year!), I opted for the adventure story that I had a much more solid narrative for, and something much more self-contained (and in fact, a requirement for the competition). We'll see how I do, but it speaks to the kind of choices and gambles one takes when choosing something to present to a larger event.

But anyway, in this post, we work toward obtaining uniform distributions of spheres and balls (this post will still be primarily setup on the disk in the plane; I'm splitting the already written post up in an effort to cap the length of my posts). We're working up to exploring a somewhat surprising way that uniformity can propagate itself. There has been some online buzz about balls and spheres and their very interesting volume formulas involving one of our frequent friends, the Gamma function, our factorial generalizer. Actually, I'd wanted to do a post about that formula for years, ever since the Pi vs. Tau debate (more specifically: volumes of spheres and balls are where it seems $\pi$ is better), I'm glad it got some visual stuff eventually, but of course, if I wanted to do my own post, maybe it would have more impact if I could rework it in some way in line with the theme of this blog. And that way, as we have been really enjoying lately, is through the lens of probability. However, we will not yet get to the full generalization today, as there is considerable material to get through (considerable FUN material, of course, but still, yours truly does admittedly have a bit of a problem keeping things short).

Summing the Squares: Not So Triangular?

First off, if we ask ourselves the question, what is the distribution of the sum of the squares of two uniform random variables, which is different from just adding them together to form the triangular distribution. This tells us what is the squared distance from the origin in the plane. Two uniform variables in $[-1, 1]$ mean that this radial squared function actually has a very complicated dependence because it has a different value in the corners than on the sides. We'd expect it to be uniform up to the unit disk, touching the boundary, and then it behaves differently when it goes into the corners. To see the real deal, we need to compute the distribution of squares of a uniform $[-1, 1]$ variables, and adding two independent instances of these. This, just as in the un-squared case, uses convolution.

So first let's examine the square of a uniform variable (first on $[0,1]$; the case $[-1,1]$ is a rescaling but causes some complications with having two signs in a square root): the pdf is $dx$ (keeping pdfs multiplied by $dx$ helps me remember how it transforms. One read/understand them more formally as infinitesimals—a "little bit of x"—or differential forms, my favorite mathematical tool in the known universe), but then taking $y = x^2$ means $x = \sqrt{y}$ and $$dx = \frac{1}{2\sqrt{y}} dy,$$ on the interval $[0,1]$ (and zero elsewhere). Write the resulting random variable $Z$ and $\rho_Z(z) = \frac{1}{2\sqrt{z}}$. To get the sum of two independent versions (call it $S$) of this requires us to convolve this function with itself. This turns out to be one tricky problem. First, sums of independent random variables are given by convolution (a fact we have used before):

\[\rho_{S}(\xi) = \int_{\mathbb{R}} \rho_{Z}(\xi - z) \rho_{Z}(z)\; dz\]

Looks straightforward enough to plug in the formula, but the trickiness comes in determining the right bounds of integration. Central to it is to reflect one of them and determine three different stages of overlap. There's 4 stages: no overlap with the reflected-translated function way to the left (so no overlap and so the function is $0$), overlap of the two unbounded sections, overlap of the two bounded sections, and no overlap due to the reflected-translated being far to the right. It's not the most inspiring calculation, so let me just show the result:

\[\rho_Y(y) = \frac{1}{2} \sin^{-1}\left(\frac{2}{y} -1\right)\]

pdf of the sum of two squared uniform random variables

Testing it with random number generators (see accompanying notebook), with a million values of $X$ and $Y$, we similarly bin the squared distance values $S$:

Very cool! You never would have thought of that shape arising from such a simple process. However, if you think about it a bit, this directly reflects that the disk is uniformly covered, but going into the corners is more difficult. The sharp break occurs because the area of the circle hits the boundary square and with distance > 1 cannot be achieved in the middle of the sides, but can be achieved at the corners. 

Now what we do is get rid of the samples $> 1$ which is basically tossing out the fraction $1 - \frac{\pi}{4}$ of the ones that exceed 1. The expected number of trials we need to run to get that is $4/\pi$ (as it's a geometric random variable), so we expect to have $\frac{\pi}{4} \approx 0.7853$ proportion of the generated samples to satisfy that condition. Then now we get to the good part: we can restrict to where $S \leq 1$. Let's just rename it to be $S$ even though technically it is different variable, with a restriction. In our random sample of a million, it did indeed leave us 784560 points (in the notebook version of this, you may get something slightly different; all that is important is it is close to $1-\pi/4$ of a million).

Now if we go back to our random variables $X$ and $Y$ (going from $[-1, 1]$ but whose squares, summed to $S$, gives a distribution exactly like the above by symmetry), also now both restricted so as to never be greater distance than $1$, we scatterplot the random samples corresponding to what produced the above:

Does it look uniform?

A slight lie: we actually showed only 1/100 of the sampled points (about 7845 in this image) because a million is wayy dense for an image like this. Looks very uniform! Again, we should emphasize the fact that although the original random variables for the $x$ and $y$ coordinates were independent uniform random variables their restrictions to the disk, $X$ and $Y$, are NOT independent. It is intuitive to understand at a big-picture, global level: $X$ and $Y$ cannot be independent, because the value of $Y$ is constrained to a very specific range once a value of $X$ is picked. The reason why this can be confusing is because when you look at an integral, the integrand $dx\; dy$ looks like a product. If we explicitly write an indicator function $\chi_D(x, y)$ then it becomes more clear, but often a double integral is written with the domain as a subscript, which hides the fact that there is a nontrivial dependence between the variables. We'll look at this in detail shortly.

Beware of Transformations

Of course, when we routinely dealing with something with circular symmetry, like a disk, often the concept of polar coordinates comes to mind. It would seem intuitive that every radius from $0$ to $1$ is equally possible, and every angle in $[0, 2\pi]$ is, also. But let's take a look at what happens:

Oops! All that uniformity in $\theta$ is scrunched down near the origin, making the sampling very concentrated here. The intuition is that the polar coordinate transformation takes what is evenly spread out in the rectangle $[0,1] \times [0,2\pi]$ (imagine it as a grid of rectangles of equal area), and then warps a bunch of stuff, so that the transformed rectangles have different sizes. This demonstrates exactly why we need Jacobians when we deal with transforming the domains of integration.  To remedy this situation, recall that the area element of the disk is not $\frac{1}{2\pi} dr \; d\theta$ which is effectively what we have here, but rather $\frac{1}{\pi}r \; dr \; d\theta$ instead. The solution to understand things in these terms is to factor it $\left(2r\; dr \right)\left(\frac{1}{2\pi} d\theta\right)$. The CDF for $r$ is simply $r^2$, so we use its inverse n the generation of $r$ from a uniform random variable. A variant of this fix is essentially what we did in our previous post on uniform distribution on the sphere. The moral of the story is, just because we start out with uniform (independed) variables, it doesn't mean we can apply an arbitrary transformation to them and still get something uniform. We may have to tweak one or more of the variables.

Breaking Down Multidimensional Uniformity via Conditioning

We explore another way of generating multidimensional (jointly distributed) random variables, which will be important in the coming posts, and really, the most general way we can understand such variables one dimension at a time. We again begin with an approach that seems natural, but doesn't work. For this we just consider two uniform $[0, 1]$ variables and reject if the sum of their squares $> 1$. This should uniformly fill the upper right quarter disk (it's for convenience so that I don't have to deal with a bunch of plus/minus signs with square roots and state all sorts of caveats with inverse trig functions and calculating the integrals. You already know how much I don't like inverse trig functions).

The reasoning for this sampling method is as follows. If we are already given a uniform sample $y$, then the draw that produced it needs to make it fall within the disk, but with no other restrictions. Namely $x$ would be uniformly be within $[0, \sqrt{1-y^2}]$. Thus a seemingly good strategy would be to get the sample $v$, and synthesize $x$ as the appropriate uniform random variable. The density function for $x$ in this case can be computed by some manipulations of indicator functions:

$$\int_0^1 \frac{1}{\sqrt{1-y^2}} \chi_{[0, \sqrt{1-y^2}]}(x) \; dy$$

The reason for this function is that we are simply using integration to take some weighted averages of sample values continuously indexed by $y$. Namely, given a fixed $y$, we have that $x$ is a some possible value for the random variable, with pdf being the conditional distribution of $x$ given $v$. Then the probability density evaluated at $v$ is the weighting in this weighted average. This is called a mixture distribution. By manipulating the inequality $0 \leq x \leq \sqrt{1-y^2}$ to get it in a form that gives its explicit dependence on the variable we actually want to integrate, $y$, we see it is equivalent to $0 \leq y \leq \sqrt{1-x^2}$, not surprising due to the symmetry of a circle:

$$\chi_{[0, \sqrt{1-y^2}]}(x) = \chi_{[0, \sqrt{1-x^2}]}(y).$$

We need this, because a indicator function on $y$ instead of $x$ enables us to just change the limits of integration. Therefore 

$$p_X(x) = \int_0^{\sqrt{1-x^2}} \frac{1}{\sqrt{1-y^2}} \; dy$$

the integrand of which we recognize as the derivative of $\sin^{-1}(y)$. Thus the density is

$$p_X(x) = \sin^{-1}(\sqrt{1-x^2}) - x \sin^{-1}(0) = \sin^{-1}(\sqrt{1-x^2})$$

on $[0, 1]$. What does this actually look like? (To recap: we select the $y$-coordinate with the usual uniform randomness, and then we select the $x$-coordinate uniformly among all the remaining possible values that are possible once $y$ is fixed):

It totally gets scrunched up top there in the upper left. This is because the uniformity in 2 dimensions means that the frequency of sampling of the horizontal coordinate decreases as the vertical coordinate increases, but this is not accounted for if you just sample the vertical coordinate independently first. This is somewhat subtle and not intuitively obvious; it comes from conflicting definitions of uniformity and equality, and that:

One has to realize uniformity does not persist under transformations and projections.

In order to "fix" the problem of the nonuniformity happening here, we consider the uniform distribution in the disk as having density form $\frac{4}{\pi} dx \; dy$. If the conditioning on a specific value of $y$ gives a uniform distribution in $[0, \sqrt{1-y^2}]$, then we have to split it up between the numerator and denominator: essentially all we're doing is "undoing" the extra scrunching. In some sense, it compromises by thinning out what happens vertically, so that things sample less frequently in the small scrunched up horizontal intervals near the upper left corner.

$$\begin{multline}\frac{4}{\pi} dx\; dy = \frac{4}{\pi} \frac{1}{\sqrt{1-y^2}} \chi_{[0, \sqrt{1-y^2}]}(x) dx\; \sqrt{1-y^2}\; dy \\ = \left( \frac{1}{\sqrt{1-y^2}} \chi_{[0, \sqrt{1-y^2}]}(x) dx\right)\left(\frac{4}{\pi} \sqrt{1-y^2}\;dy\right).\end{multline}$$

The first factor is the conditional distribution $p_{X|Y}(x|y)\; dx$, not $p_X(x)\; dx$. Conditional distributions are how you can factor random variable generation into a sequence of generations each depending on the previous result. To evaluate $p_X(x)\; dx$, one integrates the expression over all $y$. This gives $p_X(x)\; dx = \frac{4}{\pi}\sqrt{1-x^2}\; dx$ which is exactly the same a the distribution of $y$! This should not actually be surprising, given that the circle is symmetric between $x$ and $y$ But it should be noted that using this distribution for $x$ after using it for $y$ will not give correct results, simply because once we've simulated $y$, and we act based on that value, the distribution is by definition conditional. How do we simulate $y$ with the distribution $\frac{4}{\pi}\sqrt{1-y^2}$? The (not so fun) cumulative distribution function corresponding to this is:

$$\frac{2}{\pi} \left(y \sqrt{1-y^2} +\sin^{-1}(y)\right)$$ for $y \in [0, 1]$.

We would need to apply the inverse of this function on a uniform random variable to get what we want. This will take a bit of finessing with Newton's method if you really want to solve it. (At the end of this series, I'll have a Jupyter notebook that will allow you to follow along). We've already described how to get the $x$ coordinate from this, but to recap: we just sample uniformly from $0$ to $\sqrt{1-y^2}$ which is simply generating the uniform from $0$ to $1$ and then rescaling.

But the real takeaway here is basically a granular understanding of a formula of conditional probability that is often committed to memory long before building intuition:

$$\rho_{X,Y}(x,y) = \rho_X(x) \rho_{Y|X}(y|x).$$

The first factor is the marginal pdf, namely, what happens when we forget about $Y$ entirely, and the second is the conditional pdf. This formula always holds, and is a replacement of the more classic formula $\rho_{X,Y}(x,y) = \rho_X(x) \rho_Y(y)$ when $X$ and $Y$ are independent (usually learned much earlier). Anyway, in this formula lies the key to bootstrapping ourselves up to higher dimensions! See you in the next post. (Here's a SNEAK PREVIEW!!)

Would you Like Some Fusilli with your Inverse Trig?

I have a confession to make. I never really liked inverse trigonometric functions. I've had to help a number of students with them over the years, but more or less can relate when they express some distaste for them: it's easy to make mistakes with them, their domains are restricted in what sometimes seems like arbitrary and hard-to-remember ways, sometimes there is more than one valid answer, and conventions differ. Finally, even when you get your hands on a fancy calculator or computer software that supposedly can take care of things for you, it can also give answers that differ from what's expected in class, or otherwise require some interpretation to get right (often finessing with the quadrants, etc.). Actually, dealing with quadrants is cool; that's something akin to the notion of coordinate charts in probably this blog's favorite topic, manifolds. But it all still seems haphazard, and it all adds up to the perception that math is a rigged game in which people enforce rules for seemingly arbitrary reasons just to make you feel bad about yourself.

So I'm going to take you on an adventure involving some inverse trig, and we will get to the bottom of it and understand what is it that makes them so damn hard to deal with. And we'll have a good serving of pasta to go with it. Special edition fusilli. We will also revisit a favorite blog topic: the Riemann surface.

First off, I was very fortunate that my trig class made it a point to (gently) introduce complex numbers and eventually give the big revelation that they're all actually some combinations of complex exponentials. I'm definitely a fan of not making people memorize a lot of trig identities, when conceptually, just one suffices. If there's any lesson one should take from complex numbers besides, that square root of -1 sure is pretty damn useful, it's actually the concept of numbers carrying a generalized sign: neither positive nor negative, but a whole directional space of possibilities. That's personally when the light bulb went off for me.

Despite all of this, and, furthermore, enjoying complex analysis as an undergrad...

I still didn't like inverse trig functions. The problem with it is that there's some rushed discussion on "choosing a principal branch" or "making branch cuts" (I've already given Roger Penrose's opinion on that), they choose it for log and shuffle through how to define it for some functions, and then it's on to the next topic. So the concept never really had the time to "gel". On to grad school, there were many topics that had their origin in figuring out what to do with complex functions, but the overspecialization of topics, plus concerns about choosing a research topic ASAP, made it so that I never really got around to an in-depth study of the real (HAR!) solution to all of this: Riemann Surfaces. Now of course, we've mentioned them on the blog before, but even that was really more of a "beware, functions might not do what you think they will do", rather than actually getting to the bottom of what is really happening.

The Algebra of What's Going On

Inverse sine: blue, Inverse cosine: red.

So let's start off trying to understand the current state of frustration. The "definition" of an inverse trig function is easy enough: it's just whatever angle produces a given ratio of sides. Of course, the problem here is: more than one angle works. The function that produces this is not one-to-one, so therefore there can't be an inverse function. Usually, that's the end of the discussion. When an application comes up to solve for angles, they're usually in some restricted domain in which you can tease out the answer that you really want, by some reasoning about the problem. This is a good habit, but maybe is lost in translation. But instead there are conventions, like inverse cosine always giving an angle between 0 and $\pi$, and inverse sine giving an angle between $-\pi/2$ and $\pi/2$ (see figure above). Don't get me started on the other ones, because I don't even bother (though, you may find a love for the programmer's atan2 after this post). Once you consider the fact the functions are $2\pi$-periodic, it's now all good, right? You can just add a bunch of multiples of $2\pi$ and it's all great! Right? Right??? Oh right. There's like complementary/supplementary angles. Like $\pi$ minus the angle. And stuff. And other stuff. Oof, a headache already, right?

To help us understand what is going on, let's first figure out what these inverse functions are in terms of logs. It might be a bit of fun, because I don't think solving for inverse trig in terms of logs is on anyone's radar even after they learn Euler's identity. First consider the cosine:

$$y = \cos(x) = \frac{e^{\mathsf{i}x} + e^{-\mathsf{i}x}}{2}$$.

Now if we multiply through by $2 e^{\mathsf{i}x}$, we get

$$2 y e^{\mathsf{i} x} = e^{2\mathsf{i}x} + 1.$$

Maybe it's not obvious what to do with this, but if we rewrite it like this:

$$e^{2\mathsf{i} x} - 2 y e^{\mathsf{i} x} + 1 = 0$$

Take $u = e^{\mathsf{i}x}$. This gives $u^{2} - 2yu + 1 = 0$. Hopefully this is more familiar. Plugging into the trusty quadratic formula,

$$e^{\mathsf{i} x} = u = \frac{2y \pm \sqrt{4y^{2} - 4}}{2} = y \pm \sqrt{y^{2} - 1}.$$

Then with the logs,

$$x = \frac{1}{\mathsf{i}} \ln\left( y \pm \sqrt{y^{2} - 1}\right) = -\mathsf{i}\ln\left( y \pm \sqrt{y^{2} - 1}\right).$$

First off, whoa. We have inverse cosine in terms of logs and algebraic operations. Maybe it's old hat to the math majors, but we should remember complex numbers are often just introduced as "Oh you can't take the square root of -1? Can't stop me! 🤪🤪". It's a big unifying concept.

... but I hope you also see some, um, issues with this. First of all, that nasty $\pm$. Which one is it??? If you want purely real numbers... well... some bad news there. If $y$ is positive, then since $\sqrt{y^2-1}$ is always of smaller magnitude than $y$, both the plus and the minus give you something legal to take the log of. But then the $-\mathsf{i}$ stops you afterward. Ok, maybe we want to make the log be purely imaginary, so that the $-\mathsf{i}$ will cancel it. So we are still made to venture out into the nuances of how logs and complex numbers work. Another fact that is not obvious to start (and actually, this will be the thing that makes our pasta more interesting): it turns out that for any complex numbers $y$, the two numbers $y + \sqrt{y^{2} - 1}$ and $y - \sqrt{y^{2} - 1}$ are reciprocal (or: since there are two square roots for every complex number, this says the two possible values gotten by the square root, are reciprocal). This is easy to verify: $(y + \sqrt{y^{2} - 1})(y - \sqrt{y^{2} - 1}) = y^2 - (y^2 -1) = 1$.

Reciprocal numbers pass through the logarithm to become a minus sign, so we can (rather surprisingly) rewrite it as

$$x = \pm \mathsf{i} \ln\left(y + \sqrt{y^{2}-1}\right).$$

(Surprising, because usually it is NOT legal to take out plus/minus signs through a log like that). Now if in our classic situation we have $-1 \leq y \leq 1$, then $y^2 - 1$ is going to be negative, and thus have a complex square root. $y + \sqrt{y^2-1} = y +\mathsf{i}\sqrt{1-y^2}$. We should note that if you take the complex modulus of that, you get $y^2 + (1-y^2) = 1$. What happens when you take a log of something of complex modulus 1? The Pythagorean Trig identity and Euler's identity, the only two you need, show that you get something purely imaginary. Which, when combined with the outside factor of $\pm\mathsf{i}$, gets you two real solutions of opposite sign. If you then think about it some, it helps to recall that cosine is an even function, i.e. it gives the same result when switching the sign of its argument. So it makes sense you can have oppositely signed results for the inverse.

Finally, now the $2\pi$-periodicity of trig functions can be brought in, since Euler's identity is valid for angles that keep wrapping 'round and 'round: you get a bunch of results separated by multiples of $2\pi$ and also the result of the opposite sign separated by multiples of $2\pi$.

That's all well and good. It was a bunch of algebra and symbol wrangling. The whole point of this blog is, what the hell does this actually look like? This requires a bit more finessing, but it shows up on the teaser title image: note that it's a surface with ramps moving up and down both in a counterclockwise and a clockwise direction. This is in contrast to the usual depictions of fusilli pasta, which only has one ramp spiraling up.

The Full Complex Definition

To understand how to get a full-blown surface from all of this, we will have to stop confining ourselves to real numbers, or simple images of the real number line (i.e. $1$-dimensional subsets), such as the unit circle. In complex analysis, we define for all complex $z$,

$$\cos(z) = \frac{e^{\mathsf{i} z} + e^{-\mathsf{i} z}}{2}.$$

where you can use the complex exponential. We don't quite have the space to define that here, but one quick way is to use power series. What's important is to realize it still satisfies the same laws of exponents we know and love: $e^{z + w} = e^{z} e^{w}$, etc. Then when writing down the formula for the inverse of $w = \cos(z)$, we get

$$z = \mp\mathsf{i} \ln\left(w+ \sqrt{w^2 - 1}\right) = \mp\mathsf{i} \ln\left(w+ \mathsf{i}\sqrt{1-w^2}\right).$$

Now the question of which square root and which logarithm becomes an issue. It's easy to say, I'm going to ask you to ignore all that we drilled into your head about functions, the vertical line test, the horizontal line test, etc. Or at least, temporarily suspend it... We'll answer questions about how to precisely, unabiguously choose a different branch of a square root in programming in a bit, but for now: the most common complex square root that is implemented is, effectively, taking the polar form of the complex number, with its angle $\pi < \theta \leq \pi$ (yes, the less than or equal on the right, unless you're in Apple's grapher, for which it is $-\pi \leq \theta < \pi$, which made me have to jump through additional hoops: so please beware of this if you're going to go off on some explorations on your own. I told you that conventions differ!). The picture of the square root is to take the whole complex plane, and map it to the right half-plane, including the upper segment of the boundary, but not the lower segment. Introductory complex analysis texts pay no heed to this extra boundary happening, because the theory gives preference to open sets.

The strategy for choosing the multiple values is: take a principal branch, and then adjust based on two parameters: one, the sign of the square root, and two, which $2\pi$-period it comes from. We'll call it its ladder position. These parameters actually have a group structure which definitely surprised me the first time I saw it. We'll get to it. For now, we take the principal branch of inverse cosine to be:

$$\cos^{-1}(z) = -\mathsf{i} \ln\left(z + \mathsf{i} \sqrt{1-z^2}\right),$$

where we use the standard logarithm and square root that do the funny stuff on the negative axis. We derive this next (this choice is what reduces to that $0$ to $\pi$ range originally given, when we restrict to real values between $-1$ and $1$). The actual place where one needs to worry about a discontinuity in the complex plane (which is what branch cut means), in this case is $(-\infty, -1] \cup [1, \infty)$. 

The Graph

The classic visualization is $y = f(x)$ meaning the set of points $(x, f(x))$ in the plane, where $x$ lies in some portion of $\mathbb R$. This we've talked about many times before. It is the graph parametrization.

For complex-valued functions, ideally we would be able to visualize $(z, f(z))$ as something in 4 dimensions, two for the domain and two for the range. And visualizing 4-dimensional things is a favorite thing for mathematicians to try to do in various ways. A common way to do it for Riemann surfaces is by simply taking real and imaginary parts: $(z, \operatorname {Re}(f(z))$ and $(z, \operatorname {Im}(f(z))$ as separate 3D graphs. What's nice is that this tells you one method of visualizing the inverse function: $(f(z), z)$, or for a pair of 3D functions, in terms of parametrizations as $(\operatorname {Re}(f(z)), \operatorname{Im}(f(z)), \operatorname{Re}(z))$ and $(\operatorname {Re}(f(z)), \operatorname{Im}(f(z)), \operatorname{Im}(z))$. All of this is with $z = x+\mathsf{i}y$. The former is our title image with $-3 \leq x \leq 3$ and $-3 \leq y \leq 3$.

This is a similar graph of its imaginary part:

which rather surprisingly does not spiral around, namely, it has no more than two values per input value $z$. This comes from the fact that the real part of the complex logarithm comes from the complex number's radius. Moving things in the forward direction, though, namely the approach $(z, f(z))$ with the two coordinates of $z$ forming the $xy$ part of the parametrization, is a bit trickier, because we need to confront head on the multivaluedness of the inverse. But it does have an advantage that we can more finely control the domain coordinates. For that, we now talk about...

Learning to Live with the Branch Cuts We've Got

And just so how do we exactly how to deal with branch cuts? How do we deal systematically with them? One way was described here. Basically, the key is to consider the "problematic" part of the standard functions (both along the negative real axis) and work backward, computing what the problematic part looks like when mapped this way and that. The classical way of dealing with it is to consider continuous values of the functions along curves. But standard functions provided by computers don't take continuous curves as a parameter. However, for certain parametrizations of our surfaces, such as cylindrical coordinates, coordinate curves like the polar coordinate $\theta$ will cross the cuts at very predictable values of the coordinate. Then the key to evaluating the function is to strategically switch the sign of the square root, and the ladder position, based on the coordinate, in order to maintain the continuity. For parametrizing in terms of ladder steps going around multiple times, we will parametrize as follows: taking $z = x+iy = u e^{2\pi t}$, we take 

$$\begin{pmatrix}u \cos(2\pi t) \\ u \sin(2 \pi t) \\ \operatorname{Re}\left(-\mathsf{i} \ln \left(z + \mathsf{i} \sqrt{1-z^2}\right)\right) \end{pmatrix}$$

Now to deal with the multivaluedness of the function, we have that the crossing happens when $u \geq 1$ and $t$ at every integer and half-integer. The half-integer induces both a sign switch in the square root, and a bump along the ladder of the log, and at each integer, only the sign of the square root needs to switch in order to maintain continuity.

This changes the formula to something that's definitely less pretty... BUT it'll be much more concrete, in terms of being able to use common, readily available complex logs and square roots on your system without further custom hacks:

$$-\mathsf{i} \ln \left(z + \mathsf{i} \sqrt{1-z^2}\right)$$

is realized as

$$(-1)^{k+1}\mathsf{i} \ln \left(z + (-1)^{\lceil 2t\rceil }\mathsf{i} \sqrt{1-z^2}\right) + 2\pi\lceil t - \tfrac{1}{2}\rceil.$$

where $k$ is $0$ or $1$ to choose which ramp to go on, and the ceiling functions $2\pi\lceil t - \tfrac{1}{2}\rceil$ round you up to the next half-integer needed to make the log continuous after a full trip around, and there's an inner sign switch of the square root that keeps track both crossings.

There's actually so much more to say here about getting down into the weeds and nitty gritty to really practice wrangling with the branch cuts, but this can be an entire post in itself! I know it can be unsatisfying! As the YouTubers would say, let me know in the comments. A hidden gem in this is the intertwined nature of how the sign interacts with the log's ladder step. I originally thought the two were independent. They are not, and it turns out that it's the action of an infinite dihedral group. This would be material for yet another post! To me, this is the ultimate explanation of the double ramp, the true way I learned to love the inverse trig functions ... but hopefully the awesomeness hits sooner than having to do abstract algebra!

The result of all of this:

Notice how the coordinate curves are more uniform.

Sneak Peek

A fusilli with ramps going in two opposite directions is all well and good. But what about three directions in a 3-fold symmetry? Here's a sneak peek of that. And that's another adventure, for another post.

There's lots of directions we can go from here, and it'll be a big subject for the blog because it intertwines so many subjects! This 3-fold symmetry version will have us getting into the nitty gritty of defining branches, and a review of the cubic formula (also talked about here previously). We could talk about the algebraic topology aspect, where taking various paths around the branch points forms a group. And some unexpected cool results surrounding that.

The Configuration Space of All Lines In the Plane

Time for another classic post from the old site, while I'm at it… As we have seen in a bunch of previous posts, the notion of configuration space has held a prominent place in my mathematical explorations, because, I can't emphasize enough, it's not just the geometry of things that you can directly see that are important; it's the use of geometric methods to model many things that is.

Consider all possible straight lines in the plane (by straight line I mean those that are infinitely long). This collection is a configuration space—a catalogue of these lines. What does that mean? Of course, in short words, it's a manifold, and we've talked about many of those before. But it always helps to examine examples in detail to help develop familiarity. Being a manifold means we can find local parametrizations of our collection by $n$-tuples of real numbers. How do we do this for lines in the plane?

Finding Some Coordinates: The Slope-Intercept Equation

First off, almost anyone reading this, regardless of mathematical background, probably was drilled at one time or another on the following famous equation, called the slope-intercept equation:
\[
y = mx + b.
\] The variable m if you recall is the slope of a line, and b is the y-intercept which means where the line crosses the y-axis. However few people realize that what they're doing, really, is indexing every non-vertical line in the plane with some point in another plane (called, say, the mb-plane rather than the old traditional xy-plane). That is, we have a correspondence between points in one abstract mb-plane to the lines in the xy-plane, where m is the slope, and b the y-intercept. It doesn't get all lines, though—vertical lines have infinite slope, and you can't have a coordinate on any plane with an infinite value. That is, you have charted out the space of all possible nonvertical lines with points in a different plane.

Of course how would we catch the vertical lines? We chart it out using different coordinates. There's the possibility of using inverse slope and x-intercept, which merely means we write the equation using x as a function of y. In other words, all non-horizontal lines are given by:
\[
x=ky + c.
\] It's essentially obtained by reflecting everything about the diagonal line x = y and finding the usual slope and intercept. In other words we can chart out all non-horizontal lines on this new kc-plane. So you could represent the collection of all lines in the plane by having two charts... an atlas or catalogue of all lines, with these two sheets, the mb-plane and the kc-plane.

However, note that most lines in the plane—ones that are neither vertical nor horizontal, can be represented in both forms. That is, they correspond to points on both sheets. So, we would say in our usual terminology that these two sheets, or charts, determine a manifold of all possible lines; we just need to check that the overlap is smooth. But let's stop to think about what that means. We can think of it as trying to glue together both sheets to form just a single catalogue. The method of gluing is that we glue together the points that represent the same line in the plane. There is a nice, exact formula for the gluing, which is very simply determined: $(m,b)$ and $(k,c)$ represent the same line if and only if $y = mx + b$ and $x = ky + c$ are equations of the same line. All we have to do to convert from one to the other is solve for x in terms of y, that is, invert the functions. So we solve
\[
y = mx + b \iff y-b = mx \iff y/m - b/m = x \iff x = (1/m) y - b/m
\] that is, if $k = 1/m$ and $c = -b/m$, then $(m,b)$ and $(k,c)$ represent the same line in the two sheets. So to "glue" the sheets together, we glue every $(m,b)$ in the mb-plane to $(1/m,-m/b)$ in the kc-plane. Obviously the sheets need to be made of some very stretchable material, because it is going to be awfully hard to glue points together. Actually it's pretty hard to physically do this so don't try this at home, but just try to imagine it (don't you just love thought experiments?). For example, points in the mb-plane $(1,1)$, $(2,2)$, $(3,3)$, and $(4,4)$ get glued to the corresponding points $(1,-1)$, $(1/2,-1)$, $(1/3,-1)$, and $(1/4,-1)$. You glue them in a very weird way, but if you suppose for a moment, that you allow all sorts of moving, rotating, shrinking, stretching in this process (topological deformations), but without tearing, creasing, or collapsing, you can preseve the "shape" of this space, and yet make it look like something more familiar. This would be our new catalogue of lines. In addition, the catalogue has an additional property: nearby points in the catalogue translate to very similar-looking lines.

So, What Is It?

One should wonder what kind of overall "shape" our nice spiffy catalogue has, after gluing together the two possible charts we've made for it. As it turns out, its shape is the Möbius strip! That's right, the classic one-sided surface (without, as it turns out, its circle-boundary).

That is to say, if you give me a point on the Möbius strip, it specifies one and only one line in the plane. One would not, initially, be able see why non-orientability enters the picture. But a little interpretation is in order. First, if we take a particular line and rotate it through 180 degrees, we get the same line back. Everything in between gives every possible (finite) slope. It so happens that as far as slopes of lines is concerned, $\infty=-\infty$, and if you go "past" this single projective infinity, as they call it, you go to negative slopes. In other words, if you start on a journey on rotating a line through 180 degrees, from vertical back to vertical, you come back to the same line, except with orientation reversed (because what started out as pointing up now points down).

If you fix an origin and declare that it correspond to a certain special line in the plane, and then select a "core circle" for the Möbius strip, then as you travel around this circle, the distance traveled represents rotation angle for this special line. Traveling from the origin along the core circle and making one full loop should correspond to rotating the special line by 180 degrees. If you instead move up or down on the core circle, you instead end up sliding the line along a perpendicular direction, without changing its angle. So moving up and down the strip corresponds to parallel sliding of lines, and moving around the strip along a circle corresponds to rotating a line.

The Derivation

The specific formula we use is \[
F(m,b) = \left(\cot^{-1} m, \frac{b}{\sqrt{m^2 + 1}}\right),
\] which sends the line to the angle it makes with the $y$-axis, and its signed perpendicular distance to the origin (the sign is determined by $b$). For the other chart, \[ F(k,c) = \begin{cases}\left( \cot^{-1} \left( \frac{1}{k}\right), \dfrac{c}{\operatorname{sgn}(-k)\sqrt{k^2+1}}\right) \quad & \text{ if } k \neq 0 \\
(0,-c) & \quad \text{ if } k = 0
\end{cases}
\] We'll show how we got this in an update to this post, or perhaps a "Part 2." This can be readily checked by simply substituting the transition charts for $(m,b)$ and $(k,c)$. However, the extra case for $k=0$ here is simply gotten by taking the limit as $k$ goes to zero from above in the other case. What proves that it is a Möbius strip is that, if we take the limit as $k$ goes to zero from below, it will approach $(\pi,c)$ instead of $(0,-c)$. This would make it discontinuous, unless we decide to identify $(\pi,c)$ with $(0,-c)$: the $c$ going to $-c$ means we take the strip at $\pi$ and flip it around to glue it to the strip at $0$ (see this post for another example of defining a Möbius strip this way). Technically, we need an infinitely wide Möbius strip for this, but we can always scrunch it down into a finite-width strip without its boundary circle (using something like arctangent). It's just that the closer you get to the edge, the quicker things go off to infinity.

The animation is an example "path" through "line space," The blue dot travels around the white circle, and the line in the plane that corresponds to it is the blue line. The red line is a reference line perpendicular to the blue one, and always passes through the origin. The distance the blue dot from the core circle (in turquoise) indicates how far from the origin that the blue and red lines intersect. Because of the "scrunching down," though, the closer to the edge of the strip we get, the more dramatic the change in distance the blue line is from the origin. Here it is again in the plane with the Möbius strip as a picture-in-picture reference:

The more general object here we are describing is closely related to the notion of Grassmannian manifold, which are all $k$-dimensional subspaces in an $n$-dimensional vector space (the only difference is that Grassmannians only consider spaces through the origin).

Clifford Circles

30 Clifford circles with with $\varphi = \pi/8$ and $\theta_0$ ranging from $0$ to $2\pi$

A discussion about Clifford tori would never be complete without a corresponding discussion about Clifford circles. These were featured as the logo of the UCSD math site for many years (not the case anymore, though, but I saved a screenshot!):

Just as the Clifford tori foliate $S^3$, the Clifford circles foliate each of the Clifford tori. As part of ongoing efforts to revisit algebraic topology, this example is one of the best to explore. We begin with classical mapping, called the Hopf fibration, which maps $S^3$ to $S^2$ by

\[
p\begin{pmatrix}x_1\\x_2\\x_3\\ x_4\end{pmatrix} =\begin{pmatrix}2(x_1 x_3 + x_2 x_4)\\ 2(x_2 x_3 -x_1x_4)\\ x_1^2 + x_2^2 -x_3^2 -x_4^2\end{pmatrix}.
\]
In fancy-schmancy homotopy-theory speak, $p$ is a generator of $\pi_3(S^2)$. Considering our previous parametrization $F$ from last time:
\[
F\begin{pmatrix}\varphi \\ \alpha \\ \beta\end{pmatrix} = \begin{pmatrix} \cos \alpha \cos \varphi \\ \sin\alpha\cos\varphi \\ \cos\beta\sin\varphi \\ \sin\beta\sin\varphi \end{pmatrix},
\]
we recall that the last coordinate of $p$ is simply $A^2 - B^2$, or $\cos^2(\varphi) - \sin^2(\varphi)$. Although I tell my students to not bother memorizing trigonometric identities, we derived that from a $\sin^2 \alpha + \cos^2 \alpha = 1$ and $\sin^2 \beta + \cos^2 \beta = 1$. We can further simplify that last coordinate to $\cos(2\varphi)$.

The key property that we want to demonstrate is that each point of $S^2$ corresponds to a whole circle in $S^3$ (for any $\boldsymbol \xi$ in $S^3$, the circle in question is the inverse image $p^{-1}(\boldsymbol\xi)$), and that each of these circles corresponding to distinct $\boldsymbol \xi$, though disjoint, are nevertheless linked together.

To do this, we visit the first two coordinates:
\[
2(x_1x_3 + x_2x_4) = 2(\cos\alpha \cos\varphi \cos\beta \sin\varphi+ \sin\alpha\cos\varphi \sin\beta\sin\varphi)
\] \[= 2\sin\varphi\cos\varphi(\cos\alpha\cos\beta +\sin\alpha\sin\beta) = \sin(2\varphi) \cos(\alpha-\beta),\] where the last equation is gotten either by trolling the back of a calculus book for some trig identities, or using complex numbers. Similarly,
\[
2(x_2x_3 - x_1x_4) = 2(\sin\alpha \cos\varphi \cos\beta \sin\varphi- \cos\alpha\cos\varphi \sin\beta\sin\varphi)
\] \[= 2\sin\varphi\cos\varphi(\sin\alpha\cos\beta -\cos\alpha\sin\beta) = \sin(2\varphi) \sin(\alpha-\beta).\]
All together, we have \[
p \circ F\begin{pmatrix}\alpha\\ \beta\\ \varphi\end{pmatrix} = \begin{pmatrix}\sin(2\varphi)\sin(\alpha-\beta) \\ \sin(2\varphi)\cos(\alpha-\beta) \\ \cos(2\varphi)\end{pmatrix}.
\]
Letting $\theta = \alpha-\beta$, we see that this looks almost like our standard parametrization of a $3$-sphere, except that the polar angle $\varphi$ is off by a factor of $2$. No matter, it's still a sphere; we just have to take care to remember that the range of this $\varphi$ is $[0,\pi/2]$ rather than $[0,\pi]$. This confirms, incidentally, that $p$ really maps onto the sphere, rather than merely mapping into some amorphous blob in $\mathbb{R}^3$, as one can only assume at first, because the destination of the map $p$ has $3$ coordinates. The important thing to realize is that given any $\boldsymbol \xi$ in $S^2$, there is a unique $\varphi_0$ in $[0,\pi/2]$ and $\theta_0$ in $[0,2\pi]$ that correspond, under the usual parametrization, to $\boldsymbol \xi$. So, to calculate $p^{-1}(\boldsymbol \xi)$, we have to see how much we can change $\alpha$, $\beta$, and $\varphi$ in order to always get $(\varphi_0,\theta_0)$. By our above computations, this is easy: $\alpha$ and $\beta$ must satisfy $\alpha - \beta = \theta_0$, and $\varphi = \varphi_0$ is already completely determined. So our only degree of freedom here is $\alpha-\beta$: given some point in the fiber $p^{-1}(\boldsymbol \xi)$, if we add the same thing to both $\alpha$ and $\beta$, we will stay in the fiber. This means the fiber has a parametrization that looks like
\[
\beta \mapsto F\begin{pmatrix}\theta_0 + \beta \\ \beta \\ \varphi_0\end{pmatrix} =  \begin{pmatrix} \cos (\theta_0 + \beta) \cos \varphi \\ \sin(\theta_0 + \beta)\cos\varphi \\ \cos\beta\sin\varphi \\ \sin\beta\sin\varphi \end{pmatrix}.
\]
We finally finish things off by composing with the stereographic projection $P$ as before:
\[
\beta \mapsto \frac{1}{1-\sin\beta\sin\varphi } \begin{pmatrix}\cos (\theta_0 + \beta) \cos \varphi \\ \sin(\theta_0 + \beta)\cos\varphi \\ \cos\beta\sin\varphi \end{pmatrix}.
\]
Plotting this out, this gives us the nice circles shown at the start of the post. We'll continue to explore the properties of $p$ and its visualizations as we move along in algebraic topology.

Some Cool Views of Projective Space

While finally revisiting one of the geometry books on my shelf, Glen Bredon's Topology and Geometry, I encountered an exercise about showing that the projective space is homeomorphic to the mapping cone of a map that doubles a circle on itself (the complex squaring map $z \mapsto z^2$). The mapping cone has a nice visualization, first as a mapping cylinder, which takes a space $X$ and crosses it with the interval $I$ to form $X \times I$ (thus forming a "cylinder"), and then glues the bottom of it to another space $Y$ using a given continuous map $f : X \to Y$. Finally, to make the cone, it collapses the top to a single point. Of course, this can be visualized as deforming the bottom part of $X \times I$ through whatever contortion $f$ does, which might include self-intersection (and of course, it could be more gradual). So I used a good old friend, parametrizations, to help set up an explicit example. Take a look!

An immersion of projective space into $\mathbb R^3$. Shown as a mesh to make the self-intersecting portion visible. Looks a little like a molar, although one would hope I take good enough care of my teeth to not have that many holes in it…

A cutaway view, now as a more solid surface, basically illustrating it now as a mapping cylinder (it is homeomorphic to projective space minus a disk, which is a Möbius strip). Anyone know a good glassblower so we can make vases that look like this?

A view from the open top, allowing us to see the self-intersecting part of the surface from above

The equation, in "cylindrical coordinates", is $r = (2z+\cos\theta) \sqrt{1-z^2}$ for $0\leq z \leq 1$ (for the closed surface), $0 \leq z \leq 0.97$ (for the cutaway), and $0 \leq \theta \leq$ (what else?) $2\pi$. I say it in scare quotes because technically, it allows negative values of $r$. For fun, though (and to make it totally legit, even if you have qualms about negative radii), we rewrite it as a (Cartesian) parametrization (by substituting for $r$):
\[
\begin{pmatrix} x \\ y \\ z\end{pmatrix} = \begin{pmatrix} (2u + \cos v)\sqrt{1-u^2} \cos v\\ (2u + \cos v) \sqrt{1-u^2}\sin v \\ u\end{pmatrix}
\]
with $0\leq u \leq 1$ or $0.97$, and $0 \leq v \leq 2\pi$.

The motivation for this is that the equations $r = a + \cos\theta$, as $a$ varies, goes from a loop traversed once to a loop that folds over itself exactly in a $2:1$ manner:

but now visualized stacked on top of one another. Then to collapse the top, we shrink the diameter by a function that vanishes at $1$ and approaches $0$ at infinite slope to make it smooth (topologically speaking, it is valid to let it collapse to a corner point, though). Flowers, however, are not included (you'll have to fiddle with things like $r = \cos(k\theta)$ if you want that…)

The Configuration Space of a Double Pendulum

One of the classic ways of understand the torus is by anything that involves two circles, such as a system with two independent angular variables—roughly speaking, it is because the torus has two families of circles, one surrounding the tube, and another surrounding the hole. The classic example of such a system is a double pendulum, which is one pendulum hung off the pendulum bob of another.

Double pendulum: the path is traced by the second bob and fills a path within the annulus. If we imagine that the second circle goes into and out of the page instead, we can see it as a path on a torus.

Of course, it is not immediately apparent why one would find it advantageous to consider a whole torus, as opposed to two distinct circles, other than for cool visualization purposes, so we should explain this. First, many interesting possibilities arise when the two angular variables are not independent: they bear some complicated relationship imposed by, say, physics, in this case. Not all the possibilities are there, at least, in a short amount of time, and it merely traces some path through this grand torus. The dynamics of the interaction is something that makes the system transcend merely "having two angles".

This notion of "where all the variables live together" is historically what led to the development of the theory of manifolds. Think again to parametrizations and equations, where variables that we might only know, initially, lie in some Euclidean space, either get warped into more complicated spaces, or, when constrained, lead to funny subsets where, and only where, the constraint is satisfied. These state spaces (or configuration spaces) are everywhere, at least implicitly in any kind of modeling of real-world objects and their dynamics.

This visualization was coded in MATLAB, using a symplectic integrator (in order more accurately reflect the energy-conservation). We delve into the calculational specifics of this example in a future post.

Here's a still photo for the featured post image: