A funny series

One of the reasons I posted a funny number was to tell you the story of a funny series I discovered by chance. It had to do with my fascination for simple constructs like $j^j$ giving unexpected results like ${\rm e}^{-\pi/2}$. I always wondered: are there more numbers connected to powers of $j$ or numbers to the power of $j$?

One day I was using Matlab for a tutorial. Before it started I had the Laptop set up and Matlab running and 10 minutes to fill. Waiting for the students to come I played in Matlab: I made a random complex number, raised $j$ to the power of it, then again $j$ to the power of the last result and so on. In other words, I kept typing x = j^x and looked at x. After a little while I got (rounded to 4 digits) $0.4383 + j 0.3606$ which seemed to be a stationary point since $j^{(0.4383 + j 0.3606)} = 0.4383 + j 0.3606$. I started with another random number and got the same result. This made me think: is this stationary point unique? Do we reach it no matter where we start? If so, does it have an analytical expression? Is there anything else special about it?

I was ripped out of my train of thoughts when the seminar started. I forgot about my little experiment for quite some time until quite recently I found the old script by chance. I really wanted to find out more about that special number.

What is it I am really trying to find? It’s the limit of a recursively defined sequence: it is $\lim_{n\rightarrow \infty} x_{n}$ where

$$x_{n} = j^{x_{n-1}}$$

and $j^x$ refers to the principal value, i.e., $j^x = {\rm e}^{j \cdot \pi/2 \cdot x}$. From this it seems that the sequence can easily be generalized to

$$x_{n} = {\rm e}^{\alpha \cdot x_{n-1}},$$

where $\alpha \in \mathbb{C}$ and the previous series corresponds to the special case $\alpha = j \cdot \pi/2$. This has got to be something veeery fundamental I thought and it does generate the “magic” stationary point I had found earlier.

Unfortunately I didn’t get much further on the analytical side. You can certainly expand $x$ into $y+jz$ so that for $\alpha = j \cdot a$ (purely imaginary), ${\rm e}^{\alpha \cdot x}$ becomes

$${\rm e}^{- a \cdot z + j \cdot a \cdot y} = {\rm e}^{-a \cdot z} \cdot \cos(a \cdot y) + j{\rm e}^{-a \cdot z} \cdot \sin(a \cdot y)$$

but I do not see this going anywhere. A quick google search turned up a lot of results on linear recurrence relations but the one I had is clearly non-linear.

As of now I have no analytical expression for the stationary point for any value of $\alpha$ (trivial solutions like $\alpha=0$ let aside).


So I played a bit more with it in Matlab. It seems that different values of $\alpha$ also give rise to stationary points. For instance, $\alpha=j$ gives $0.5764 + j \cdot 0.3743$. Here is a plot with empirically found fixed points for $\alpha = j \cdot a$, varying $a$ from $0$ to $1.96$ (the series starts to form cycles for $a$ larger than $\approx 1.96125$):


It does not seem to be a circle. This is purely empirical, I have no clue if/how one could find these values analytically (but this would be very nice).

While playing with these numbers I realized that my original observation that one would always converge to one fixed point regardless of the initialization was actually wrong. For some points, the series diverged to infinity. This sparked the question where these points would be located on the complex plane. So I made a second set of experiments: I started the series from initial points $x_0$ chosen in a regular grid in the complex plane, let the series run for a predefined maximum number of iterations and observed what it does: (a) does it diverge to infinity or (b) does it converge to the fixed points I found earlier (i.e., get within a predefined threshold $\delta$ of it)?

I began with $\alpha = j\cdot \pi/2$ (corresponding to the $x_{n+1} = j^x_{n}$ it all started with) and a cartesian grid in ${\rm Re}(x_0)$ and ${\rm Im}(x_0)$ for the initial points. And then, I got this:


The two axis refer to the real part and the imaginary part of our initial point $x_0$. The color indicated the number of iterations until the fixed point is reached up to an accuracy of $10^{-3}$. We notice the periodicity in ${\rm Re}(x_0)$: this is not surprising since the first step in the series gives

$x_1 = {\rm e}^{j \cdot \pi/2 \cdot x_0} = {\rm e}^{-\pi/2 \cdot {\rm Im}(x_0)} \cdot {\rm e}^{j \cdot \pi/2 \cdot {\rm Re}(x_0)}$

and therefore, adding an integer multiple of 4 to ${\rm Re}(x_0)$ gives the same $x_1$. I found the shape looks very interesting and started zooming in.


Further… Frac1_z2

Even further…


This is already a 200x zoom on both axis compared to the original one. Yet more and more details are revealed that look similar to the original ones – we seem to have self-similarity here. I’m not sure if that’s enough to call this a fractal but I’d like to claim so until someone convinces me otherwise.


I have since played with various values of $\alpha$, kept zooming in and looking at the result. So far it seems that the level of details keeps increasing as $\alpha$ approaches $2 \cdot j $. After about $1.961 \cdot j$ it becomes hard to still draw a picture because instead of unique attractors the fractal seems to run into limit cycles consisting of a discrete number of points. There still seems to be a fractal-looking region of divergence but I’m mostly guessing here.

Here is one  for $\alpha = j \cdot 1.96$:


I chose 50 as the maximum number of iterations here, the picture does not change much if I increase it to 1000 or so (it just takes much longer to render).

Here is my best attempt for $\alpha = j \cdot 1.963495$:


I display the plane of initial points in polar instead of cartesian coordinates (for cosmetic reasons only) so this one is periodic in the vertical direction and truncated on the right side horizontally. The maximum number of iterations was set to 2000 here, increasing it further does have an impact on the larger red areas that get filled more and more. I tried to count how long it takes until a limiting cycle is reached but this is based on very rough heuristics and might be wrong. I still enjoy the look of the result.

It’s really a lot of fun to “explore” this little world, especially if you’ve never seen it before and you’re not sure if anyone else ever has. It’s quite possible these fractals are actually well-known, I just don’t know what I should search for. In light of this it would be really great if anyone could give me any hint that brings me closer to any of these open questions:

  • Is this class of series known, does it have a name, is there any published material on it?
  • How can we compute the fixed points of the series for a given value of $\alpha$?
  • Are the resulting images that divide converging from diverging initial points actually a fractal? Does it possibly also have a name, has it been studied before?
  • What happens beyond the “magic” breaking point where the series starts running into limit cycles comprising of a discrete number of points? Does the breaking point have any special meaning? Can the limit cycles be described analytically?

If you have any idea to any of these please leave me a comment!

A funny number

I was giving tutorials for an undergraduate course in Signals and Systems theory. Since students came from quite varying backgrounds I had to start very simple – my first tutorial used to be just on complex numbers. This left some of the students bored, saying they’ve seen it all. To keep them thinking I usually asked them to compute a complex number for me: Let $j$ be the imaginary unit, what is the value of $j^j$?

This question is not as trivial as it may seem at first sight. It requires a generalization of the exponentiation $a^b$ from $a, b \in \mathbb{R}$ to $a, b \in \mathbb{C}$. There are rigorous ways of doing this which I do not want to discuss in detail here. Let us just assume that we found a generalization of $a^b$ to complex numbers which satisfies the laws of powers, in particular the law $a^b = {\rm e}^{b \cdot \ln(a)}$ where $\ln(x)$ is the natural (base-e) logarithm. Then, we can rewrite $j^j$ as $j^j = {\rm e}^{j \cdot \ln(j)}$. What is $\ln(j)$ though? Well, $j$ can be written as ${\rm e}^{j \cdot \pi/2}$, so $\ln(j)$ should be $j \cdot \pi/2$, right? This finally brings us to $j^j = {\rm e}^{j \cdot j \cdot \pi/2} = {\rm e}^{-\pi/2}$.

Some of my students would actually obtain this result. They would usually be surprised to get a real (and sort of strange number) out of such an operation, but always they would be very proud to have the answer. The tech-savvy ones would even check their result with whatever internet-able device they carried and be extra sure to have it right.

However, my (admittedly a bit discouraging) reply would be that the answer is, despite being correct, incomplete. Actually infinitely incomplete. There are (infinitely) more “correct” answers. How so? Well, $j$ can be written as ${\rm e}^{j \cdot \pi/2}$ but it can also be written as ${\rm e}^{j \cdot (\pi/2 + 2\pi)}$ or ${\rm e}^{j \cdot (\pi/2 – 2\pi)}$, and so on. We can add any integer multiple of $2\pi$ to the phase due to the periodicity of complex numbers with respect to their phase.

For the (natural) logarithm of a complex number this means that it is in general ambiguous: $\ln(j) = j \cdot \pi/2 + 2\cdot k \cdot \pi$ for any $k \in \mathbb{Z}$. To avoid the confusion it is common to write ${\rm Ln}(x) = \ln(x) + 2\cdot k \cdot \pi$ where $\ln(x)$ is the (unique) principle value of the logarithm which we obtain by choosing the principle phase of $x$.

In that sense, the full answer would be $j^j = {\rm e}^{-\pi/2} \cdot {\rm e}^{-2 \cdot k \cdot \pi}$, where $k=0$ corresponds to the principle value of $j^j$ given by ${\rm e}^{-\pi/2}$.

Infinitely many solutions, all of them are real, and all of them connect ${\rm e}$ and $\pi$. Pretty cool, right?

Long time no post…

… I apologize. But I have a good excuse: I was really really busy trying to finish my thesis. That’s right folks, this guy has submitted his dissertation recently and the last weeks towards this goal were quite … intense.

I’m not entirely happy with the result but that’s normal I think, noone ever is. Still I’m happy I can again focus on other things in life.

I will make the thesis publicly available as soon as I am officially allowed to do so, likely some time in October. I’ve decided to publish it under a Creative Commons BY-NC-ND license, so you will be free to copy/print/read it as often as you like. 😉


Until then, stay tuned and keep your fingers crossed for the defense that I still have to go through. 😉

Fun with statistics – transformations of random variables part 2

I recently posted on how to find the distribution of functions of random variables, i.e., the distribution of $Y=g(X)$, where $X$ is a random variable with known distribution and $y=g(x)$ is some function.


As a natural extension of this concept we may ask ourselves what happens if we have two random variables involved. Let us start with one function of two random variables, i.e., given $X$ and $Y$ and knowing their joint PDF $f_{X,Y}(x,y)$ or their joint CDF $F_{X,Y}(x,y) = {\rm Pr}[X \leq x, Y \leq y]$ we would like to calculate the distribution of $Z = g(X,Y)$ where $z=g(x,y)$ is a function with two arguments, e.g., $z=x+y$.

Again, there are multiple ways of addressing this problem. A natural way would be to calculate the CDF of $Z$ directly, i.e., $F_Z(z) = {\rm Pr}[Z \leq z] = {\rm Pr}[g(X,Y) \leq z]$. In other words, we need to compute the probability of the event that relates to all realization of $X$ and $Y$ which satisfy $g(X,Y) \leq z$. This is easily done by integrating the joint PDF $f_{X,Y}(x,y)$ over all points in the set ${\mathcal D}_z$ which contains all points $(x,y)$ for which $g(x,y) \leq z$. Written out, we have

$$F_Z(z) = {\rm Pr}[Z \leq z] = {\rm Pr}[g(X,Y) \leq z] = \iint_{{\mathcal D}_z} f_{X,Y}(x,y) {\rm d}x {\rm d} y$$

Whether or not this approach is easy to follow depends on two things: (1) how easy it is to parametrize the set ${\mathcal D}_z$ and (2) how easy it is to integrate the joint PDF over ${\mathcal D}_z$.

Let us make an example considering the simple function $g(x,y) = x+y$. Then ${\mathcal D}_z$ contains all points $(x,y)$ for which $x+y \leq z$, i.e., $y\leq z-x$ or $x \leq z-y$. Geometrically, this is the set of points that are on the lower-left of a line with slope -1 and offset $z$, i.e., a line passing through $(z,0)$ and $(0,z)$. The integral over this set is relatively simple, as we can directly write it as

$$\displaystyle F_Z(z) = \int_{-\infty}^{+\infty} \int_{-\infty}^{z-y} f_{X,Y}(x,y) {\rm d}x {\rm d} y = \int_{-\infty}^{+\infty} \int_{-\infty}^{z-x} f_{X,Y}(x,y) {\rm d}y {\rm d} x$$.

Another example is $g(x,y) = \max(x,y)$. Since $\max(x,y) \leq z \Leftrightarrow ((x \leq  z) \;\mbox{and} \; (y \leq z))$ we can argue

$$F_Z(z) = {\rm Pr}[\max(X,Y) \leq z] = \int_{-\infty}^z \int_{-\infty}^z f_{X,Y}(x,y) {\rm d}x {\rm d} y$$.

Geometrically, ${\mathcal D}_z$ contains all points on the “lower left” of the point $(z,z)$, i.e., the intersection of the half-planes below $y=z$ and left of $x=z$.


The second extension is to consider two functions of two random variables. Say we are given the distribution of $X$ and $Y$ via their joint PDF, we would like to find the joint PDF of $Z=g(X,Y)$ and $W=h(X,Y)$. There is a closed-form expresion for it as a direct extension of the closed-form expression for the PDF of one function of one random variable. It reads as

$$f_{Z,W}(z,w) = \sum_{i=1}^N \frac{1}{|{\rm det} \ma{J}(x_i,y_i)|} f_{X,Y}(x_i,y_i)$$,

where $(x_i,y_i)$ are all solutions to the system of equations $z=g(x,y)$, $w=h(x,y)$ in $x$ and $y$. Here, $\ma{J}$ is the Jacobian matrix given by

$$ \ma{J} = \left[ \begin{array}{cc}
\frac{\partial g}{x} & \frac{\partial g}{y} \\
\frac{\partial h}{x} & \frac{\partial h}{y}

Moreover, the term ${\rm det} \ma{J}(x_i,y_i)$ means that we first compute the determinant of the Jacobian matrix (in terms of $x$ and $y$) and then insert $x_i(z,w)$ and $y_i(z,w)$.


Example? How about the joint distribution of $X+Y$ and $X-Y$? In this case, solving for $z=x+y$ and $w=x-y$ for $x$ and $y$ is simple, we have one solution given by $x_1 = (z+w)/2$ and $y_1 = (z-w)/2$. The Jacobian matrix is given by

$$ \ma{J} = \left[ \begin{array}{cc}
1& 1 \\ 1 & -1
and hence its determinant is $-2$ everywhere. This gives the solution for $f_{Z,W}(z,w)$ in the form

$f_{Z,W}(z,w) = \frac{1}{2} f_{X,Y}((z+w)/2,(z-w)/2)$.

As in the 1-D case, this direct solution depends heavily on our ability to solve the given functions for $x$ and $y$, which may be tedious for complicated functions.


Interestingly, the first case where we considered one function of one random variable can be solved also via this approach, simply by creating another “auxiliary” variable, and then marginalizing over it. So once we have $Z=g(X,Y)$ we make up another $W=h(X,Y)$, choosing it such that the remaining calculations are simple. For instance, for $g(x,y) = x+y$ we may choose $h(x,y) = y$. Then, the Jacobian matrix becomes

$$ \ma{J} = \left[ \begin{array}{cc}
1& 1 \\ 0 & 1

with determinant one. Moreover, we have $x_1 = z-w$ and $y_1 = w$. Therefore, we get

$$f_{Z,W}(z,w) = f_{X,Y}(z-w,w)$$.

The final step is marginalizing out the auxiliary $W$ which gives

$$f_Z(z) = \int_{-\infty}^{+\infty} f_{X,Y}(z-w,w) {\rm d}w.$$

Looks much like a convolution integral, doesn’t it? In fact, if $X$ and $Y$ are statistically independent, we can write $f_{X,Y}(x,y) = f_X(x) \cdot f_Y(y)$ and hence we obtain

$$f_Z(z) = \int_{-\infty}^{+\infty} f_{X}(z-w)\cdot f_Y(w) {\rm d}w = f_X(x) * f_Y(y),$$

where $*$ denotes convolution. This shows very easily that the PDF of the sum of two random variables is the convolution of their PDFs, if they are statistically independent.