Home > MATLAB, R > When 1 * x != x

When 1 * x != x

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.

About these ads
  1. 17th January, 2011 at 6:38 am

    A more rigorous mathematical proof would be to look at the norm. There is only one “infinity” in the complex plane. And R knows it:

    > is.finite(Inf+3i)
    [1] FALSE

    even though its maths are wrong:

    > Inf*1i
    [1] NaN+Infi

    (same reason as yours!)

  2. 17th January, 2011 at 19:29 pm

    Another way of understanding the two different results comes from thinking about limits. I’ve discussed your example in that context in this post.

  1. 17th January, 2011 at 22:57 pm

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 229 other followers

%d bloggers like this: