[Numpy-discussion] Apparently non-deterministic behaviour of complex array multiplication

kneil magnetotellurics at gmail.com
Thu Dec 1 15:17:19 EST 2011


Hi Oliver, indeed that was a typo, I should have used cut and paste.   I was
using .transpose()

Olivier Delalleau-2 wrote:
> 
> I guess it's just a typo on your part, but just to make sure, you are
> using
> .transpose(), not .transpose, correct?
> 
> -=- Olivier
> 
> 2011/11/30 Karl Kappler <magnetotellurics at gmail.com>
> 
>> Hello,
>> I am somewhat new to scipy/numpy so please point me in the right
>> direction
>> if I am posting to an incorrect forum.
>>
>> The experience which has prompted my post is the following:
>> I have a numpy array Y where the elements of Y are
>> type(Y[0,0])
>> Out[709]: <type 'numpy.complex128'>
>>
>> The absolute values of the real and complex values do not far exceed say
>> 1e-10.  The shape of Y is (24, 49218).
>> When I perform the operation: C = dot(Y,Y.conj().transpose), i.e. I form
>> the covariance matrix by multiplying T by its conjugate transpose, I
>> sometimes get NaN in the array C.
>>
>> I can imagine some reasons why this may happen, but what is truly
>> puzzling
>> to me is that I will be working in ipython and will execute for example:
>> find(isnan(C)) and will be returned an list of elements of C which are
>> NaN,
>> fine, but then I recalculate C, and repeat the find(isnan(C)) command and
>> I get a different answer.
>>
>> I type:
>> find(isnan(dot(Y,Y.conj().transpose)))
>> and an empty array is returned.  Repeated calls of the same command
>> however result in a non-empty array.  In fact, the sequence of arrays
>> returned from several consecutive calls varies. Sometimes there are tens
>> of
>> NaNs, sometimes none.
>>
>> I have been performing a collection of experiments for some hours and
>> cannot get to the bottom of this;
>> Some things I have tried:
>> 1. Cast the array Y as a matrix X and calculate X*X.H --- in this case i
>> get the same result in that sometimes I have NaN and sometimes I do not.
>> 2. set A=X.H and calculate X*A --- same results*
>> 3. set b=A.copy() and calc X*b --- same results*.
>> 4. find(isnan(real(X*X.H))) --- same results*
>> 5. find(isnan(real(X)*real(X.H))) - no NaN appear
>>
>> *N.B. "Same results" does not mean that the same indices were going NaN,
>> simply that I was getting back a different result if I ran the command
>> say
>> a dozen times.
>>
>> So it would seem that it has something to do with the complex
>> multiplication.   I am wondering if there is too much dynamic range being
>> used in the calculation?  It absolutely amazes me that I can perform the
>> same complex-arithmetic operation sitting at the command line and obtain
>> different results each time.  In one case I ran a for loop where I
>> performed the multiplication 1000 times and found that 694 trials had no
>> NaN and 306 trials had NaN.
>>
>> Saving X to file and then reloading it in a new ipython interpreter
>> typically resulted in no NaN.
>>
>> For a fixed interpretter and instance of X or Y, the indices which go NaN
>> (when they do) sometimes repeat many times and sometimes they vary
>> apparently at random.
>>
>> Also note that I have had a similar problem with much smaller arrays, say
>> 24 x 3076
>>
>> I have also tried 'upping' the numpy array to complex256, I have like
>> 12GB
>> of RAM...
>>
>> This happens both in ipython and when I call my function from the command
>> line.
>>
>> Does this sound familiar to anyone?  Is my machine possessed?
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
> 
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> 

-- 
View this message in context: http://old.nabble.com/Apparently-non-deterministic-behaviour-of-complex-array-multiplication-tp32893004p32898294.html
Sent from the Numpy-discussion mailing list archive at Nabble.com.




More information about the NumPy-Discussion mailing list