When 1 * x != x

Posted on January 16, 2011 by richierocks in R bloggers, Uncategorized | 0 Comments

[This article was first published on 4D Pie Charts » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Trying to dimly recall things from my maths degree, it seems that in most contexts the whole point of the number one is that it is a multiplicative identity. That is, for any x in your set, 1 * x is equal to x. It turns out that when you move to floating point numbers, in some programming lanugages, this isn’t always true.

In R, try the following

x <- Inf + 0i
1 * x
[1] Inf + NaNi

Something odd is happening, and my first reaction was that this is a bug. In fact, I reported it as such (in a slightly different form) but the terrifyingly wise Brian Ripley came up with the cause.

When you multiply two complex numbers, this is what happens

$(a + bi)(c + di) = (ac - bd) + (ad + bc)i$

In this instance

$(Inf + 0i)(1 + 0i) = (Inf * 1 - 0 * 0) + (Inf * 0 + 0 * 1)i = Inf + NaNi$

So, there is a reasonable argument that R is doing the right thing.

MATLAB works a little differently since all it’s numbers are implicitly complex. The inner workings of MATLAB are somewhat opaque for commercial reasons, but in order for compiled MATLAB code to work, I think it’s reasonable to assume that MATLAB stores its variables in a class that is similar to the MWArray class used by compiled code. As far as I can tell, each double matrix contains two matrices for storing real and imaginary components, and has an IsComplex property that is true whenever the complex part has any nonzero values. If IsComplex is false, only the real matrix is used for arithmetic.

Thus in MATLAB, when you define x = Inf + 0i it immediately resolves itself to simply Inf, and we don’t get the same weirdness.

x = Inf + 0i
1 * x
Inf

Other languages are divided on the subject: python and ruby agree with R but Mathematica (or at least Wolfram Alpha) agrees with MATLAB.

Thinking about this from a purely mathematical sense,

$lim_{n to infty} (n + 0i)(1 + 0i) = lim_{n to infty} (n * 1 - 0 * 0) + (n * 0 + 0 * 1)i = lim_{n to infty} n = Inf$

This concurs with the MATLAB answer and I’m inclined to agree that it makes more sense, but the issue isn’t entirely clear cut. Changing the way complex arithmetic works for a single edge case is likely to introduce more confusion than clarity. It is perhaps sufficient for you to remember that complex infinity is a bit pathological with arithemtic in some languages, and add checks in your code if you think that that will cause a problem.

Tagged: complex numbers, matlab, r