Inconsistent numpy complex multiplication results

Question

Consider the following Python code that multiplies two complex numbers:

import numpy as np
a = np.matrix('28534314.10478439+28534314.10478436j').astype(np.complex128)
b = np.matrix('-1.39818115e+09+1.39818115e+09j').astype(np.complex128)

#Verify values
print(a)
print(b)

c=np.dot(a.getT(),b)

#Verify product
print(c)

Now the product should be -7.979228021897728000e+16 + 48j which is correct when I run on Spyder. However, if I receive the values a and b from a sender to a receiver via MPI on an MPI4py program (I verify that they have been received correctly) the product is wrong and specifically -7.97922801e+16+28534416.j. In both cases I am using numpy 1.14.3 and Python 2.7.14. The only difference in the latter case is that prior to receiving the values I initialize the matrices with:

a = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
b = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)

and then the function MPI::Comm::Irecv() gives them the correct values.

What is going wrong in the latter case if the a and b are correct but c is wrong? Is numpy arbitrarily setting the imaginary part since it's quite smaller than the real part of the product?

Would recommend changing from using `np.matrix()` to `np.array()`s first off and see if the behavior continues. The short of it: SciPy recommends using `ndarray`s, the long of it: https://stackoverflow.com/q/4151128/5087436. Also: "Is numpy arbitrarily setting the imaginary part since it's quite smaller than the real part of the product?" no. `complex128` in numpy are just stored as two separate 64-bit floats. — alkasm, May 22 '18 at 03:24
Thanks but it didn't work. I replaced `np.matrix()` with `np.array()` and `getT()` with `transpose()` and it is still giving the same result in the MPI case. I even removed the `np.empty_like()` to see if this changes anything but no. — mgus, May 22 '18 at 04:00
@AlexanderReynolds Any other suggestions? Please also note that the machine is the same in both cases. I am thinking of using sympy but as I remember from the past this was a bit slow for large matrix multiplications. — mgus, May 22 '18 at 04:30
I don't understand the use of `np.empty_like()` in this example, why not just use `np.empty((1, 1), dtype=np.complex128)`? Either way, I'm going to guess that the junk values that get inserted into `np.empty()` are actually what is multiplying your number to give you not what you expect. Try initializing the array as zeros instead with `np.zeros()` and see if you end up with zeros. Why initialize at all anyways? — alkasm, May 22 '18 at 04:31
@AlexanderReynolds The reason for initialization is that the received operands are n-by-n matrices in general and not just 1-by-1. Thus, the MPI needs to know that the receiver buffer is `complex128` and also preallocated. By the way, I also tried all of your suggestions but the result didn't change. — mgus, May 22 '18 at 04:41

hpaulj · Accepted Answer · 2018-05-22T05:19:53.343

First, this doesn't address the mp stuff, but since it was raised in comments:

np.matrix can takes a string argument, and produce a numeric matrix from it. Also notice that the shape is (1,1)

In [145]: a = np.matrix('28534314.10478439+28534314.10478436j')
In [146]: a
Out[146]: matrix([[28534314.10478439+28534314.10478436j]])
In [147]: a.dtype
Out[147]: dtype('complex128')

String input to np.array produces a string:

In [148]: a = np.array('28534314.10478439+28534314.10478436j')
In [149]: a
Out[149]: array('28534314.10478439+28534314.10478436j', dtype='<U36')

But omit the quotes and we get the complex array, with shape () (0d):

In [151]: a = np.array(28534314.10478439+28534314.10478436j)
In [152]: a
Out[152]: array(28534314.10478439+28534314.10478436j)
In [153]: a.dtype
Out[153]: dtype('complex128')

And the product of these values:

In [154]: b = np.array(-1.39818115e+09+1.39818115e+09j)
In [155]: a*b       # a.dot(b) same thing
Out[155]: (-7.979228021897728e+16+48j)

Without using mp, I assume the initialization and setting is something like this:

In [179]: x=np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
In [180]: x[:]=a
In [181]: x
Out[181]: matrix([[28534314.10478439+28534314.10478436j]])
In [182]: y=np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
In [183]: y[:]=b
In [184]: y
Out[184]: matrix([[-1.39818115e+09+1.39818115e+09j]])
In [185]: x*y
Out[185]: matrix([[-7.97922802e+16+48.j]])

It may be worth trying np.zeros_like instead of np.empty_like. That will ensure the imaginary part is 0, instead of something random. Then if the mp process is just setting the real part, you should get something different.

Inconsistent numpy complex multiplication results

1 Answers1