[Numpy-discussion] Apparently non-deterministic behaviour of complex array multiplication
Karl Kappler
magnetotellurics at gmail.com
Wed Nov 30 20:44:41 EST 2011
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?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20111130/154d3cad/attachment.html>
More information about the NumPy-Discussion
mailing list