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.

Can Two Independent Random Variables Sum to Uniform?

As you well know, I'm a big fan of probability and thinking of it in geometric terms. This is because the sample spaces in probability are almost always some sort of high dimensional state space, that describe the possible outcomes (or configurations) of a system. Probability often boils down to measuring what fraction of the overall total space of outcomes, that your particular experiment resulted in. In all of I've read about and studied over the years, understanding the intuition behind probability has been the theme that recurs the most. Today I want to start on an exploration of the notion of uniformity. Now we actually have visited this theme before as well, and the next couple of posts will actually expand upon those, exploring interesting facts about what happens in higher dimensions. There will be plenty of geometry to visualize!

But today I want to start very simple. What can be simpler than a uniform random variable between 0 and 1? It's basically everyone's random number generator: many distributions are simply calculated by combining multiple uniform random variables in different ways, for example, by arithmetic operations (adding, multiplying, raising to powers), applying functions (exponentials, sin, cos, etc), applying an inverse cumulative distribution function, etc. And we will attempt to do something that may look like it breaks all sorts of rules! Break rules? Consummate rule-follower mathematicians BREAKING RULES? What has the world come to? Now I don't encourage literally breaking rules and deliberately getting into trouble. At least, I'm not supposed to. All I am saying to do is to play around a bit and develop intuition, and try to apply things in ways you're not entirely sure of. It's the sense of adventure and some risk-taking, and what the nature of math exploration actually is. At the end of the day, we will need to have firmer justifications for things and we do need to verify. But, thinking of math as living, breathing, experimental science is not what students are used to thinking. Sometimes the intuitions behind things are lost, as well as a sense of adventure. Let me tell you, this was definitely quite the adventure... You will see me get very excited, and I hope it'll rub off on you, too.

Triangular Distribution

Now, adding independent random variables is a very well-known operation. Given $X$ and $Y$ independent uniform random variables (say on $[0, 1]$), if $Z = X + Y$, this is another random variable. What may be surprising at first is that such a sum is actually NOT uniformly distributed. It is easy enough to see that its maximum range is from 0 to 2. Next, we do something like, ask ourselves, when can the two variables sum to something in the range of 1.5 to 2? Once one of the variables is pinned down, this puts a constraint on the possibilities of the other variable. If we get, say, 0.9 for the first, then the second variable can only be in the range of 0.6 and onward, and actually, you would not actually be able to achieve the value 2, since it requires a value of 1.1, which is not an acceptable value. 

In fact, we get something called the triangular distribution. Numbers close to 1 have the highest probability (density) for a sum, while those closer to the extreme ends have lower probabilities, because they rely on each uniform variable to have values clustered very close to 0 or very close to 1. 

So having been sort of surprised about summing two uniform random variables to be non-uniform, one natural question may be to ask, is it actually possible to sum two independent random variables to be uniform? At least, that's the question that occurred to me. If you like, think about it for a bit and pause... Don't spend more than 10 minutes, though.

....

....

Ok. So it turns out that it actually is impossible. There is a very nice, visual proof at Stats StackExchange, which I actually like a lot. It would have made a nice Nested Tori post by itself. Perhaps I'll try to take it apart and explain it someday, but I also think they give a pretty good exposition as it is.

Is That All There Is to it? Should we just give up?

I mean, yes. The theorem says so. At least, we should give up on actually finding those two random variables... And that is well and good. However, should it be the end of the story? What happens if we try anyway? There's many interesting further explorations one can do here to shore up the probabilistic intuition of what is going on. First, knowing that there ARE random variables that can be subdivided in the manner we just described. The exponential (actually, the negative log of a uniform variable), in particular, is the sum of two identically distributed chi-square random variables. The exponential distribution is also implicated in another somewhat related and very interesting theorem about the product of identically distributed uniform random variables, both raised to a power of another such, discussed here by our Summer of Math host.

A second interesting direction is to realize that there are actually ways to sum random variables to uniform. They just cannot be independent. This has come up in some attempts I made to gain intuition on Bose-Einstein statistics, a quantum phenomenon in which indistinguishable particles or states actually behave as if all that binomial stuff collapses, as if, flipping two coins produced three outcomes of equal probability, {heads, heads}, {heads, tails}, and {tails, tails} and somehow that mixed outcome is weighted less. I could do a whole different post about this, as it is also fun, but fundamentally, it boils down to asking this exact same question about summing to uniform. It's just that, in order to solve it, I had to give up independence, in which I would flip different (unfair) coins depending on the outcome. 

But the third adventure is a real fun one, and that's what the rest of this post is about. Essentially, it had partially been started before learning of that Stats StackExchange proof and I worked it out in different phases. Had I known the StackExchange proof, I might well have never gotten started.

Probability distributions are characterized by their density functions, and we actually very systematically understand what effect does taking a sum of two independent random variables has on their density functions. It's a very well-known operation called convolution.

Convolutions

Given the density functions $\rho_X$ and $\rho_Y$, the density function $\rho_{X+Y}$ is given by the convolution integral

$$\rho_{X+Y}(z) = (\rho_X * \rho_Y)(z) = \int_{-\infty}^\infty \rho_X (x)\rho_Y(z-x) \; dx.$$

The intuition here is using one function to smear out another function by weighted averages. It's somewhat related to the convolutions encountered in machine learning and AI, especially if multidimensional, but convolution is very interesting and important operation in math, useful in all sorts of real analysis contexts. It's also an engineering favorite. This has various nice properties reminiscent of multiplication, like being associative, commutative, and distributive over addition of functions. As for an "identity" ... well... we'll get to that. The reason why convolution is very popular is because of its behavior under the Laplace transform (and Fourier transform as well, but the Laplace transform is what turns out to be the key to really solving this problem). This theorem is the

Convolution Theorem.

Given $f, g$ "nice" functions, the Laplace transform of the convolution, $\mathscr{L}[f*g](s)$, is equal to the product of the Laplace transforms $\mathscr{L}[f](s) \mathscr{L}[g](s)$.

In other words, on the transform side, it goes from being "like" multiplication to just, well, being multiplication. It is for this reason that Laplace transforms are useful in differential equations, and is often quoted as "transforming calculus problems into algebra problems". It should also be noted that convolutions are well-defined for functions that may not necessarily be probability density functions, and the theorem still holds, when the Laplace transforms of such functions exist. Or even, when they don't ... that's a different part of the adventure....

So what is happening in terms of probability density functions (pdfs) is that if one takes Laplace transforms of them, then the Laplace transform (LT) of the pdf of the random variable representing the sum, is the product of the LTs of both pdfs. The LT of uniform with pdf $\chi_{[0,1]}$, the indicator (or characteristic) function that is $1$ on $[0, 1]$ and $0$ elsewhere,

$$\int_{0}^\infty e^{-st} \chi_{[0, 1]} \; dt = \int_{0}^1 e^{-st} \; dt = \frac{1-e^{-s}}{s}.$$

So, now, if we want two independent random variables $\rho_X$ and $\rho_Y$ to represent something that sums to uniform, this means their convolution is the uniform pdf.

But that means ... their Laplace transforms multiply to $ \frac{1-e^{-s}}{s}$! (that's me exclaiming, not taking the factorial). But if they are identically distributed... This means the our purported density function for this is $$\mathscr{L}[\rho_X] = \mathscr{L}[\rho_Y] = \sqrt{\frac{1-e^{-s}}{s}}.$$

Great!!! All we have to do, then, is find the inverse Laplace Transform! No problem at all!!

Except... well... ::opens up a reference book:: looks at the Bromwich Integral... ok where am I supposed to put the contours? What are the singularities? What about that branch cut? Didn't we spend a lot of time last year and years before that hating on branch cuts? ARRRRGH WHAT?? :: slinks away in defeat... ::

Highly Suspect Distributional Calculus

So back to the drawing board. Let's pick apart what that Laplace transform of a uniform pdf really means. It is $\chi_{[0,1]}$, with transform $\frac{1 - e^{-s}}{s}$. As Laplace transform fans know, multiplication by $s$ is differentiation, and division by $s$ is integration starting at some point. Now let's forget that denominator for a moment, and look at pieces separately.

Laplace (and Fourier) transform fans are may be familiar with some (Schwartzian) distributions, or generalized functions (we'll years before that hating on branch cuts the latter term since the word "distribution" is overused) namely the Dirac delta (impulse "function") and its derivatives. I actually first learned about that function in an engineering class, not in a math class. And the dang thing is sure useful. The important things to know about this "function" (As a mathematician, I'm obligated to put function in quotes):

1. $\int_{-\infty}^{\infty} \delta(t) \; dt = 1$.
2. $\delta * f = f*\delta = f$ (it is the identity for convolution).
3. The Laplace transform of $\delta$ is $1$ (actually this follows from the above by the Convolution theorem).
4. $\delta$ is the (suitably generalized) derivative of the unit step function (Heaviside function) $\chi_{[0, \infty)}$, the function that is $1$ on the positive real line, and $0$ elsewhere.

It can be visualized as having a graph with infinitely skinny, infinitely high peak at $0$, a value of $0$ elsewhere, yet has a total area under it equal to $1$ (no actual function on the real line can do this: being nonzero at only one point always means it integrates to $0$). When the derivative of something is a $\delta$, that represents a jump discontinuity of size $1$. Yes, I know some of those listed properties do need some more justification, since $\delta$ is not in fact actually a function, but rest assured, there are firm foundations for this.

So if we just look at $1 - e^{-s}$, this corresponds to the Dirac delta "function", minus the Dirac delta "function" shifted to the right by $1$. This is also intuitive: the uniform pdf $\chi_{[0,1]}$  has two jumps, the value going from $0$ to $1$ at the point $0$, and from value $1$ back to $0$, at the point $1$. In general, $e^{-as}\mathscr{L}[f](s)$ the transform of the function $f(x-a)$, another useful property (the translation property)... we'll be using this a lot, soon.

Now if we take the square root:

$$\sqrt{1-e^{-s}}$$

then corresponds to the convolutional square root of that difference of deltas. What can we do with that? There is another thing we know: we know how to expand $\sqrt{1 - e^{-s}}$ in a power series (binomial series) (but crucially, and this is what took some time for me to see, in powers of $e^{-s}$, not $s$). For $\operatorname{Re}(s) > 0$, the series converges uniformly, so there are a number of good analytic operations we can do. Namely, the series is

$$\sqrt{1 - e^{-s}} = (1-e^{-s})^{1/2} = \sum_{n=0}^\infty {1/2 \choose n}(-1)^n e^{-ns}$$

The suspicious looking binomial coefficent ${ 1/2 \choose n}$ is defined using the Gamma function, the usual function that generalizes factorial to non-integers. The notation really is the generalization of binomial coefficient to the non-integer case and is written that way precisely because it generalizes a binomial expresion like $(1+y)^n$ to a non-integer power. It's just that the coefficients, once nonzero, never become zero again, so expanding it out makes it an infinite series. We find that we actually can massage the coefficient a bit to express it in terms of the more normal (integer) binomial coefficients:

$${1/2 \choose n} = \frac{1}{2^{2n-1} n} (-1)^{n+1} { 2n-2 \choose n-1}$$ (and conventionally taking it to be $1$ at $0$). Then $e^{-ns}$ is the transform of $\delta$ translated by $n$ (it's using the property that its LT is 1, and the translation property). Taking this inverse transform we have:

$$\mathscr{L}^{-1}\left[\sqrt{1-e^{-s}}\right] = \sum_{n=0}^\infty { 1/2 \choose n}(-1)^n\delta_n = \delta_0 - \sum_{n=1}\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \delta_n =: \rho(t).$$

The latter sum rather nicely gets rid of the signs, so that only the first term has a plus sign and the rest of the terms are subtracted from it.

A more "illegal" way of seeing this is to consider $\delta_0$ the identity operator as a convolution, and $\delta_0 - \delta_1$ as being a "perturbation" of the identity, and we formally take the square root of that difference using Binomial series:

$$\sqrt{\delta_0 - \delta_1} = \delta_0 - \frac{1}{2} \delta_1 - \frac{1}{8} \delta_1^2 - \dots = \sum_{n=0}^\infty { 1/2 \choose n}(-1)^n\delta_1^n$$

and using the fact that $\delta_1^n = \delta_n$ (because convolution by $\delta_1$ multiple times will translate by 1 over and over again). This is the best of Euler and he would do it with abandon, without regard to convergence. The advantage of the Laplace transform method is that we have all the series converging and we know that stuff like term-by-term integration is guaranteed to work in such circumstances.

In either case, this is indeed a generalized function in the same sense as $\delta$ is, and in fact, a tempered one: given a function $\phi$ the very nice Schwartz class of smooth functions vanishing at infinity after multiplication by polynomial of any degree, when $\rho$ operates on it, the resultant infinite series of translations of $\phi$ must converge, since each discrete value will be dominated by any negative $n$th power. The situation is even nicer with smooth functions of compact support (vanishing outside a finite interval): $\rho$ will always be a finite sum at any evaluation point. Suffice to say, the convergence is good enough.

Next We Will Need: The Square Root of .... Integral?

Ah. But we're missing something. We said to ignore that denominator $s$, so we're still not done yet. It's time to bring it back. Now, just like multiplying by $s$ corresponds to differentiation, dividing by $s$ corresponds to integration (from 0). Since the convergence is nice, dividing by $s$ we can just integrate that crazy sum above term-by-term, and be done..........

Oh wait. We divide by the square root of $s$, not $s$. Oops. What does dividing by the square root of $s$ mean? Dividing by $s^2$ is integrating twice. So dividing by the square root of $s$ is, of course, just integrating half a time. WHAT??? Yes the square root of integral!! (a double exclamation. Not the double factorial of integral, although I wonder if it's possible to..........)

One additional way to motivate things is that differentiation plays very well with convolution. Namely $\partial^{i+j} ( f * g) = (\partial^{i+j} f) * g = f*(\partial^{i+j} g) = \partial^i f * \partial^j g$ in any combination that makes sense. You can split the differentiations up to happen on each component that can take that many derivatives, you can put them all on one factor, or the other factor, you can re-split it up in different ways, etc. Similarly, if we think of integration from $0$ to $t$ as a negative order differential, it all continues to work out, too.

So here if we define $\rho$ to be the binomial-series sum of deltas above, we have $\rho * \rho = \delta_0 - \delta_1$ by construction. If we have to do it on a test function, fine. The second convolution will stagger the additional terms, so they all cancel exactly outside of $[0, 1]$. Then the function we seek, $\tilde \rho$ satisfies $\tilde \rho * \tilde \rho = \chi_{[0,1]} = \int_0^x \delta(t) \; dt$.

Thus $\tilde \rho * \tilde \rho = \int_0^x \rho * \rho \; dt = \partial^{-1}(\rho * \rho)$, that is, it is the integral of the whole thing. For our next illegal trick, we're going to split that integration in exactly half, and put it on each individual $\rho$:

$$\tilde \rho * \tilde \rho = \partial^{-1}(\rho * \rho) = \partial^{-1/2} \rho * \partial^{-1/2} \rho.$$

Total nonsense?? Maybe. I have not actually justified for generalized functions in general, but we'll see in a short bit how Laplace transforms take care of all of it, at least with the exponent of ½. For functions vanishing at $0$, it sure is convenient (vanishing at $0$ makes it work well with Laplace transforms, since we often need to subtract derivatives of initial conditions). It turns out the iterated integral operation has, what is called, a kernel representation: you can represent iterated integrals as one single integral. The "kernel" in kernel representation is analogous to continuous version of a matrix, and the integral against a function is like a continuous sum version of matrix multiplication.

The kernel representation of iterated integration is:

$$\int_0^t \int_0^{t_{n-1}} \dots \int_0^{t_0} f(t_0) \; dt_0 \; dt_1 \; \dots \; dt_{n-1} = \frac{1}{(n-1)!}\int_{0}^t f(\tau)(t-\tau)^{n-1}\; d\tau.$$

The proof is actually very pretty straightforward: interchanging the order of integration and proof by induction. When the innermost variables are switched, the inside is constant with respect to the switched variable, so integrates to a difference. Stacking them up will give you the powers, so now it's only one integral.

But the amazing thing is... that expression is defined for all real numbers $\nu$ (we can rely on our trusty old friend, the Gamma function, for the factorial). We DEFINE the fractional integral

$$I^{\nu}(f)(t) = \partial^{-\nu}(f)(t) = \frac{1}{\Gamma(\nu)}\int_0^t (t - \tau)^{\nu - 1}f(\tau) \; d\tau.$$

The important thing to know here is $I^{\nu + \mu} = I^\nu I^\mu$. This can be shown by composing two integrals and using beta functions (of Bayesian statistics fame) to rewrite the product $\Gamma(\alpha) \Gamma(\beta)$. It will work for non-integer $\nu$ and $\mu$. It's a bit of a mess, so we won't reproduce all of it here. Anyway, this makes integration operations of all orders a form of convolution—the telltale sign being, the $\tau$ and the $t-\tau$, though there is some funny business with the limits of integration. Also, there's trouble if $\nu$ is a nonpositive integer (because the gamma function is infinite in such places), with the function $\frac{1}{\Gamma(\nu)}t^{\nu-1} \chi_{[0, \infty)}(t)$. Thus if we want to write it for $\nu = -\frac{1}{2}$, we have (recalling $\Gamma(1/2) = \sqrt{\pi}$, another amazing fact, and another place where $\pi$ has a habit of just inviting itself to random parties):

$$I^{1/2}(f)(t) = \frac{1}{\sqrt{\pi t}} \chi_{[0, \infty)}(t) * f(t).$$

We can also justify this a bit more by showing that the Laplace transform of $\frac{1}{\sqrt{\pi t}}\chi_{[0, \infty)}(t)$ is in fact $\frac{1}{\sqrt{s}}$ so that it just becomes $\frac{1}{\sqrt{s}}$ times the original series with $\sqrt{1-e^{-s}}$. Then this means the original transform is just a series of translated $\frac{1}{\sqrt{\pi t}}\chi_{[0, \infty)}(t)$'s. But STILL, ISN'T IT TOTALLY AWESOME THAT WE TOOK THE SQUARE ROOT OF INTEGRAL???

AND FINALLY...

So what's left to do? We need to take the ½th integral of that sum of delta's above. So this is just convolving the function $$\frac{1}{\sqrt{\pi t}} \chi_{[0, \infty)}(t)$$ with the majestic

$$\delta_0 - \sum_{n=1}\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \delta_n$$

Convolutions with deltas are supposed to be trivial, because they're the identity op, remember? We have some shifted deltas, however. Convolving with it gives a shifted function, namely, we shift that $\frac{1}{\sqrt{\pi t}} \chi_{[0,\infty)}(t)$ over by $n$ when we apply $\delta_n$. The result is therefore

$$\tilde \rho(t) = \frac{1}{\sqrt{\pi t} } \chi_{[0, \infty)}(t)* \rho(t) = \frac{1}{\sqrt{\pi t} }\chi_{[0, \infty)}(t)- \sum_{n=1}^\infty\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \frac{1}{\sqrt{\pi (t-n)} }\chi_{[0, \infty)}(t-n)$$

Those indicator functions are bloody annoying looking. I kept them in specifically because it allows you to keep track of the bounds of integration better. But there's a nice reformulation, since here, the action of these indicator functions is merely to cut off the sums depending on whether we have encountered the relevant delta term. In fact, its effect is to finitize the sum (which turns out to allow us to dodge some convergence questions): at any given $t$, only the terms up to $\lfloor t \rfloor$ have been added! We thus rewrite:

$$\tilde \rho(t) = \frac{1}{\sqrt{\pi t} } - \sum_{n=1}^{\lfloor t \rfloor}\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \frac{1}{\sqrt{\pi (t-n)} }$$

for $t > 0$ and $0$ otherwise. This includes a value of $\pm \infty$ at every integer. To avoid that, we can use $\lceil t \rceil -1$ instead of $\lfloor t \rfloor$, but infinite values on a discrete set of measure zero (set with no length) don't affect the result, especially when it comes to convolution, so we'll not bother.

Of course, we should actually probably take a look at what we just constructed. This is the graph of $\tilde{\rho}$:

You can very plainly see it has negative values, which prevents it from being a pdf. We've shown that a convolutional square root does in fact exist, but it is not a pdf. Because the uniform pdf has a downward jump, no fix can be made to eliminate the negative values (since a downward jump requires convolution by something with vertical asymptotes, whose values approach $-\infty$). And that, my friends, is my very satisfying negative answer to the this post's title! The the number of different tools and areas of math that this touches upon, not to mention, having just a random factor of $\pi$ showing up where there are no circles, is a real delight!

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.

Interstellar Lily Pads

 

I'll someday explain what this is, but for now just enjoy! (Well, this is ultimately a visualization blog. The math explanations are bonus.)

Some Loxodromes

Loxodromes are generated by curves of constant bearing on a sphere (always making a constant angle with circles of latitude). If you go in a direction and proceed forever in something that's NOT one of the 4 cardinal directions, you'll get something very interesting. It's not really obvious that you actually will spiral forever around a pole. These have the simplest expression in isothermal coordinates we talked about in the last post: they can be represented simply as straight lines in the $uv$ plane (well, if you allow $v$ to range over all of $\mathbb R$ and wrap around; if you insist on $v$ only being in an interval of length $2\pi$, then after the line goes off the top edge, you make it jump back to the bottom edge so you get something that looks like this: $/////$. So the full parametrization would be: $(\operatorname{sech}u \cos(au+b), \operatorname{sech}u \sin(au+b), \tanh u)$. It is easy to intuitively see why the angle between these curves and the $\theta$ curves are constant: the $uv$-plane is first mapped (conformally) to the whole $z$-plane by $(e^u\cos v, e^u\sin v)$, and then the (conformal) stereographic projection wraps it all up. The same cannot be said of the $\varphi\theta$ plane (to confirm, of course, you need to mess with metric coefficients in that last post).

The Uniform-Distribution-on-the-Sphere's Nemesis: Isothermal Coordinates

Inspired by the last post where the conformal version of the polar grid is $(x,y) = (e^u \cos v, e^u \sin v)$ i.e. the radius is now exponentiated (namely what looks like a grid of squares in the domain gets mapped to what looks like squares in the range: the concentric circles bunch up much more densely close to the origin, to correspond to radial lines converging at the origin), I set out to do it on the sphere as well (that conformal polar grid is simply the image of the complex plane under the complex exponential $e^z$). This is in some sense the polar (hah) opposite of something that would produce a uniform distribution of points on a sphere, because this map would assign each square an equal chance (so of course, the closer you get to the poles, the squares are smaller and more of them encompass the same space in the image, so they'll be much more densely packed).

The quickest way to see this is to note the stereographic projection is a conformal map, and has the formula

\[(x,y) \mapsto \left(\frac{2x}{x^2+y^2+1}, \frac{2y}{x^2+y^2+1}, \frac{x^2+y^2-1}{x^2+y^2+1}\right)\]

and with a substitution of the complex exponential, you get \[(u,v) \mapsto \left(\frac{2e^u \cos v}{e^{2u}+1}, \frac{2e^u \sin v}{e^{2u}+1}, \frac{e^{2u}-1}{e^{2u}+1}\right) = \left(\frac{2}{e^u+e^{-u}}\cos v, \frac{2}{e^u+e^{-u}}\sin v, \frac{e^u-e^{-u}}{e^u+e^{-u}}\right)\]

I rewrote that in that form for fans of hyperbolic functions, which gives the map

\[ (u,v) \mapsto(\operatorname{sech} u \cos v, \operatorname{sech} u \sin v, \tanh u)\] (hyperbolic secant seems to be one of the least used ones; $\LaTeX$ doesn't recognize it!). But the upshot is, this formula looks very similar to the usual geographic coordinate system (just swapping out some transcendental functions).

Another version with finer grid lines. Notice how they cluster at the top.

Another way of looking at this is via manipulation of the metric coefficients directly, which raises some interesting points on the nature of the Laplace equation. If we consider the standard geographic parametrization$(x,y,z) = (\sin\varphi\cos\theta, \sin\varphi\sin\theta, \cos\varphi)$ with Euclidean metric $g = ds^2 = d\varphi^2 + \sin^2 \varphi\; d\theta^2$. To say that the parametrization is conformal (grid squares to grid squares) is to say that there are some parameters $u$ and $v$ and a function $\sigma$ such that $g = \sigma^2(du^2 + dv^2)$ (forcing the matrix to not only be diagonal, but a multiple of the identity). If we set things up so that $du = -\csc\varphi d\varphi$ and $dv = d\theta$, then $\sin^2\varphi (du^2 + dv^2) = d\varphi^2 + \sin^2\varphi\; d\theta^2)$. Consulting trusty integral tables, Wolfram Alpha'ing it, or whatever, we find that this differential equation gives $u = \ln(\csc \varphi + \cot \varphi) + C$. One can make it this derivation a little less "out of the blue" situation in a very interesting derivation of Laplace's equation using certain complex analysis hackery (George Springer's Introduction to Riemann Surfaces, section 1-3). We'll derive this another day, but the summary is, we end up deriving that $u$ has to satisfy the spherical Laplace equation \[\frac{1}{\sin\varphi} \frac{\partial}{\partial \varphi}\left(\sin \varphi \frac{\partial u}{\partial \varphi}\right) + \frac{1}{\sin \varphi} \frac{\partial^2 u}{\partial \theta^2} = 0,\] or really, it only has to satisfy this on an open subset of the sphere, in this case, a sphere with two points deleted. This fact often gets glossed over by people talking about geometry; we know that on a sphere, the only harmonic functions that are actually defined on the sphere everywhere are constant functions—so we can only find "interesting" harmonic functions on a sphere if they blow up at certain points. Using the very reasonable assumption that such a function $u$ should be independent of $\theta$, we derive $\frac{d}{d \varphi} \left(\sin\varphi \frac{d u}{d \varphi}\right) = 0$, making $\sin\varphi \frac{d u}{d \varphi} = C_1$ a constant, and $\frac{du}{d\varphi} = C_1 \csc \varphi$ to finally get $u = C_1 \ln(\csc \varphi + \cot \varphi) + C_2$, as before (just with two degrees of freedom, namely $C_1$ controls what the initial spacing of the grid lines, and $C_2$ controls what the value is at the equator). This function blows up at the poles, a fact which manifests itself visually by the grid curves $u =$ constant become more closely spaced together at the poles (and become infinitely dense there, meaning, this is only a valid coordinate system in the sphere outside the poles). We just set $C_1 = 1$ and $C_2 = 0$, which assigns the equator to $u=0$.

Finally, we have to do some algebra, though, to get to where we want: trying to find the inverse of $u = \ln(\csc \varphi + \cot \varphi)$. This says $e^u = \csc\varphi+\cot\varphi = \frac{1+\cos\varphi}{\sin\varphi}$. This is just $\frac{1+z}{r}$ where $r$ is the cylindrical coordinate $\sqrt{x^2+y^2}$. But if we're on the sphere, $z^2 + r^2 =1$, so this gives us \[e^{2u} = \frac{(1+z)^2}{1-z^2} = \frac{1+z}{1-z}.\] Rearranging we get $e^{2u}-e^{2u}z = 1+z$, or $z(1+e^{2u}) = e^{2u}-1$. This gives \[z=\frac{e^{2u}-1}{e^{2u}+1} = \frac{e^{u}-e^{-u}}{e^{u}+e^{-u}} = \tanh u.\]. To get the $x$ and $y$, we note that $\sin\varphi = r = e^{-u}(1+z) = e^{-u}(1+\tanh u)$. Finally,

\[ 1+ \tanh u = \frac{\cosh u + \sinh u}{\cosh u} = \frac{e^u}{\cosh u}, \] which, multiplying by the $e^{-u}$ gives $\frac{1}{\cosh u}$ or $\operatorname{sech} u$.

The coordinates we derived are called isothermal coordinates, and knowing what a big fan I am of stat mech, I'd like to fully understand the origin of the term someday. Finally, one might ask if there are coordinate systems (conformal or otherwise) that will cover up all the sphere without weird coordinate singularities like that. The answer is no, and the short explanation of why is because a smooth vector field on the sphere must vanish somewhere (the hairy ball theorem), and since coordinates always set up vector fields (take the vector in the direction of a coordinate curve) and aren't supposed to collapse at good points. This is a long topic for another day (we'll show some pretty cool coordinate systems on the sphere, though!)

Happy Halloween from Nested Tori!

Spooky greetings. A quick post with more complex analysis for you. Though maybe not as visually interesting as fractals, but getting back to some basics and/or nuance regarding complex functions. For one thing, I just never had a chance to really fit in a good study of Riemann surfaces (interest in them was partially stirred on by all that recent cubic hackery). Also, a lot of stuff involving branch cuts and branch points are something quickly gone over in standard complex analysis and never really explored in depth, other than to say "we take the principal branch of the logarithm/root/power" function to be... and you don't get into the nuance with the details of differences between different branches of functions. This example here is exploring the mapping $z \mapsto \sqrt{z^2 + 1}$. This is already interesting from the Riemann surface point of view, because neither the function nor its inverse is definable in its most interesting form by mapping (subsets of) the complex plane alone. Formally, to work on $\mathbb{C}$, we have to define things by taking a continuous choice of square root on $\mathbb{C}$ minus some ray from $0$ to $\infty$. The choice here is between a square root and its opposite (which is also a square root of the same complex number!) is not natural in the complex plane (unlike in the real numbers where you can always take the positive square root of a positive number, which is in some sense, a defining characteristic of real numbers, cementing why it is the unique complete ordered field), and being forced to choose forces a discontinuity somewhere. Riemann surface theory is, in some sense, about never having to make the choice. To quote Sir Roger Penrose (congrats on that Nobel Prize!): 

In particular, the domain of the logarithm function would be ‘cut’ in some arbitrary way, by a line out from the origin to infinity. To my way of thinking, this was a brutal mutilation of a sublime mathematical structure. (The Road to Reality, p. 136)

A couple of references (e.g. Penrose's student Tristan Needham, in his otherwise fantastic Visual Complex Analysis) suggest the most common branch cut for this is function is two horizontal lines to the left of $\mathsf{i}$ and $-\mathsf{i}$. I wrangled with this a bit, but I believe a subtlety is left out. To study $\sqrt{z^2 + 1}$ as truly the composition of a square root function, and the mapping $z^2 + 1$, your branch cuts will have to be the corresponding inverse image of the branch cut of the square root, under $z^2 + 1$. Using the negative real axis, this is the ray above $\mathsf{i}$ and the ray below $-\mathsf{i}$. Two rays to the left of $\pm \mathsf{i}$ cannot be realized as such an inverse image, for any ray to infinity (the image under $z^2 + 1$ is a whole parabola to infinity). Of course, the essential definition for complex analysis is that you define the function in this set by taking paths from an origin point and analytically continuing, and entrusting your result to algebraic topology. But this is not very effective, computationally. This is basically saying that for some selections of branch cuts, the function composition $\sqrt{z^2 + 1}$ is a lie: you cannot get it by assembling square root, squaring, and plus one. What turns out to work, in terms of assembling functions rather than continuing along curves, is not $\sqrt{z^2 + 1}$, but $\sqrt{z+\mathsf{i}} \sqrt{z-\mathsf{i}}$, that is, if you take the branch of square root, along the negative axis, subtract and add $\mathsf{i}$, take the square root of each, and then multiply, you do get something that reproduces what you theoretically get by paths. One might wonder what is the difference between $\sqrt{z^2 + 1} = \sqrt{(z+\mathsf{i})(z-\mathsf{i})}$ and $\sqrt{z+\mathsf{i}}\sqrt{z- \mathsf{i}}$; these are not the same because $\sqrt{ab} \neq \sqrt{a}\sqrt{b}$ in the complex plane unless you regard it as a multivalued set equality.

The best rendering of this function is, I've discovered, is taking the branch $-\mathsf{i}$ to $\mathsf{i}$ along the imaginary axis—something that doesn't even go out to infinity. To his credit, Needham does mention that connecting the two branch points in a finite part of the plane is possible, and it just has to be done in a manner that disables you from being able to complete a turn around the two branch points $\pm \mathsf{i}$. But there's no explicit computation that develops a visceral understanding of this fact (and so I made this post).

The title picture is the image of a large disk (parametrized by the usual polar coordinates) under the mapping $\sqrt{z-\mathsf{i}} \sqrt{z+ \mathsf{i}}$. The cut (deleted out of the disk) is taken on the interval $-\mathsf{i}$ to $\mathsf{i}$ ($y=-1$ to $y=1$ in the picture). Computationally this means the square roots were taken with the angles in the interval $-\pi/2$ to $3\pi/2$, which required me to define the square root using some hackery of the arctan function. I leave you with the conformal version of this (where the concentric ring radii increase exponentially rather than linearly; the linear increase looks better for looking like a cobweb, but this version looks better in preserving the proportions of the grid rectangles formed by the intersecting lines).

Or, take a look at it in Desmos: