Search code examples
pythonnumpycomplex-numbers

Inconsistent numpy complex multiplication results


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?


Solution

  • 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.