conj and element-by-element multiplication
This is really more of a NumPy question, but I was calculating coherence and I need to multiply a vector by its own complex conjugate element-by-element. This should return a vector of real numbers in my mind, but instead returns a vector with imaginary parts identically equal to zero. Now maybe I have been spoiled by Matlab, but I ended up having to test that the imaginary parts really are 0 and then just kept the real part. Should this happen automatically, or am I just being a lazier whiner? Ryan
Ryan Krauss wrote:
This is really more of a NumPy question, but I was calculating coherence and I need to multiply a vector by its own complex conjugate element-by-element. This should return a vector of real numbers in my mind, but instead returns a vector with imaginary parts identically equal to zero.
Are you absolutely sure? I do precisely the same calculations, often, and the imaginary parts are not zero; they are small, well within roundoff error. I usually just take scipy.real() of the result when I'm done.
So far, it seems to be true. I calculated the coherence of about 10 sets of data so far and included the following code to check this: cohnum=multiply(Gxyave,conj(Gxyave)) cohden=multiply(Gxxave,Gyyave) curavebode.coh=divide(cohnum,cohden) print('max(abs(imag(coh)))='+str(max(abs(imag(curavebode.coh))))) print(str(max(abs(imag(curavebode.coh)))==0.0)) and I am getting 0.0 as the max(abs(imag(coherence))) and True for the boolean test. Ryan Stephen Walton wrote:
Ryan Krauss wrote:
This is really more of a NumPy question, but I was calculating coherence and I need to multiply a vector by its own complex conjugate element-by-element. This should return a vector of real numbers in my mind, but instead returns a vector with imaginary parts identically equal to zero.
Are you absolutely sure? I do precisely the same calculations, often, and the imaginary parts are not zero; they are small, well within roundoff error. I usually just take scipy.real() of the result when I'm done.
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Ryan Krauss <ryanfedora@comcast.net> writes:
This is really more of a NumPy question, but I was calculating coherence and I need to multiply a vector by its own complex conjugate element-by-element. This should return a vector of real numbers in my mind, but instead returns a vector with imaginary parts identically equal to zero. Now maybe I have been spoiled by Matlab, but I ended up having to test that the imaginary parts really are 0 and then just kept the real part.
There's no downcasting; if you want a vector of real numbers, you have to ask for it. The multiply doesn't actually know that the two things it's multiplying are conjugates of each other, and checking for that would slow down the general case. In this case, you're probably better off doing the norm explicitly. Using scipy, it'd look like this: def norm2(x): """compute |x|^2 = x*conjugate(x)""" if scipy.iscomplexobj(x): return x.real**2 + x.imag**2 else: return x**2 which gives you a real vector back. This is also theoretically twice as fast, as you don't do the imaginary part (which comes to zero). -- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
Ryan Krauss <ryanfedora@comcast.net> writes:
Thanks David, that seems like the right answer. But to get the speed benefit you mention, I need to use multiply(x.real,x.real) instead of x.real**2 - same for imag. There is about a factor of 10 time difference between the two. So my norm2 function now looks like this:
Are you on Windows? The pow() function there, well, sucks big time. On other platforms the two are pretty close in speed. I'm working on a better version of the power ufunc for Numeric, so hopefully in a later version this will be fixed (or rather, worked around). -- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
Ryan Krauss <ryanfedora@comcast.net> writes:
No, I am running Linux.
Odd; I do too, and don't see this. glibc's pow doesn't have the same problem as Window's does. x**2 is at most 50% slower, not 10x. Which version of Numeric and python are you using, and which Linux distribution? (Numeric's version you can get from Numeric.__version__) Ideally, x**2 should run in the same time as x*x You're using scipy's multiply, which is the same (I think) as the multiply in Numeric 24.0b2, which are different from 23.8 and earlier (which does a time-wasting postcheck of the array for infinite values).
David M. Cooke wrote:
Ryan Krauss <ryanfedora@comcast.net> writes:
Thanks David, that seems like the right answer. But to get the speed benefit you mention, I need to use multiply(x.real,x.real) instead of x.real**2 - same for imag. There is about a factor of 10 time difference between the two. So my norm2 function now looks like this:
Are you on Windows? The pow() function there, well, sucks big time. On other platforms the two are pretty close in speed.
I'm working on a better version of the power ufunc for Numeric, so hopefully in a later version this will be fixed (or rather, worked around).
-- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
On May 20, 2005, at 18:35, Ryan Krauss wrote:
I am running Numeric 24.02b on python 2.3.4. I am using Fedora Core 3 and kde.
I have modified the norm2 function to include printing both times:
def norm2(x): """compute |x|^2 = x*conjugate(x)""" if iscomplexobj(x): t1=time.time() mat1=x.real**2 + x.imag**2 t2=time.time() mat2=multiply(x.real,x.real) + multiply(x.imag,x.imag) t3=time.time() print('pow time='+str(t2-t1)) print('multiply time='+str(t3-t2)) return mat2 else: return multiply(x,x)
and get these times printed to the screen: pow time=0.00401091575623 multiply time=0.000787019729614 pow time=0.00442790985107 multiply time=0.000679016113281 pow time=0.0130319595337 multiply time=0.000170946121216 pow time=0.00465989112854 multiply time=0.000912189483643 pow time=0.00434708595276 multiply time=0.000770807266235 pow time=0.00153613090515 multiply time=0.000103950500488
Most of the calls are with 4200x3 matrices. The last one is with a 4200x1 vector (I think).
Curious. One more question: are these Complex32 matrices? I.e., are they matrices using C floats (single precision), as opposed to C doubles (double precision)? The reason is that if x is a vector of floats, x**2 upcasts x to doubles first (b/c 2 is an int -- ints are cast to doubles), so that kills the speed, and I can get timings like this. This isn't a good behaviour, and I'm going to see what I can do to fix this for the next version of Numeric.
David M. Cooke wrote:
No, I am running Linux. Odd; I do too, and don't see this. glibc's pow doesn't have the same problem as Window's does. x**2 is at most 50% slower, not 10x. Which version of Numeric and python are you using, and which Linux distribution? (Numeric's version you can get from Numeric.__version__) Ideally, x**2 should run in the same time as x*x You're using scipy's multiply, which is the same (I think) as
Ryan Krauss <ryanfedora@comcast.net> writes: the multiply in Numeric 24.0b2, which are different from 23.8 and earlier (which does a time-wasting postcheck of the array for infinite values).
David M. Cooke wrote: Ryan Krauss <ryanfedora@comcast.net> writes: Thanks David, that seems like the right answer. But to get the speed benefit you mention, I need to use multiply (x.real,x.real) instead of x.real**2 - same for imag. There is about a factor of 10 time difference between the two. So my norm2 function now looks like this: Are you on Windows? The pow() function there, well, sucks big time. On other platforms the two are pretty close in speed. I'm working on a better version of the power ufunc for Numeric, so hopefully in a later version this will be fixed (or rather, worked around).
-- |>|\/|< /------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
Hi David, You can prevent the upcasts using the .savespace() method. For example: In [8]:a = arange(10.0, typecode='F') In [9]:a**2 Out[9]: array([ 0.+0.j, 1.+0.j, 4.+0.j, 9.+0.j, 16.+0.j, 25.+0.j, 36.+0.j, 49.+0.j, 64.+0.j, 81.+0.j]) In [10]:a.savespace() In [11]:a**2 Out[11]: array([ 0.+0.j, 1.+0.j, 4.+0.j, 9.+0.j, 16.+0.j, 25.+0.j, 36.+0.j, 49.+0.j, 64.+0.j, 81.+0.j],'F') This is _very_ useful when dealing with very large arrays... Scott
Curious. One more question: are these Complex32 matrices? I.e., are they matrices using C floats (single precision), as opposed to C doubles (double precision)? The reason is that if x is a vector of floats, x**2 upcasts x to doubles first (b/c 2 is an int -- ints are cast to doubles), so that kills the speed, and I can get timings like this. This isn't a good behaviour, and I'm going to see what I can do to fix this for the next version of Numeric.
-- Scott M. Ransom Address: NRAO Phone: (434) 296-0320 520 Edgemont Rd. email: sransom@nrao.edu Charlottesville, VA 22903 USA GPG Fingerprint: 06A9 9553 78BE 16DB 407B FFCA 9BFA B6FF FFD3 2989
So I am not sure about the single or double precision complex question. I changed my code to print out the types and some more diagnostic information. They all seem to be doubles. Most of the time multiply is 5-7 times faster for 4200x3 matrices and 10-12 times faster for 4200x1 vectors. Sometimes they are strangely nearly the same speed, but this just happens once in a while with no reason - I suspect other system processes?. Take a look at the code and the output: def norm2(x): """compute |x|^2 = x*conjugate(x)""" if iscomplexobj(x): t1=time.time() mat1=x.real**2 + x.imag**2 t2=time.time() mat2=multiply(x.real,x.real) + multiply(x.imag,x.imag) t3=time.time() print('---------------------------') print('pow time='+str(t2-t1)) print('multiply time='+str(t3-t2)) print('pow time/multiply time='+str((t2-t1)/(t3-t2))) print('shape(x)='+str(shape(x))) print('x.typecode='+str(x.typecode())) print('mat1.typecode='+str(mat1.typecode())) print('mat2.typecode='+str(mat2.typecode())) if len(shape(x))==1: print('type(x[0,0])='+str(type(x[0]))) print('type(mat1[0,0])='+str(type(mat1[0]))) print('type(mat2[0,0])='+str(type(mat2[0]))) else: print('type(x[0,0])='+str(type(x[0,0]))) print('type(mat1[0,0])='+str(type(mat1[0,0]))) print('type(mat2[0,0])='+str(type(mat2[0,0]))) return mat2 else: return multiply(x,x) and I now get this output: --------------------------- pow time=0.00702810287476 multiply time=0.00100302696228 pow time/multiply time=7.00689327312 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00710797309875 multiply time=0.00673699378967 pow time/multiply time=1.05506600134 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.0021538734436 multiply time=0.000178098678589 pow time/multiply time=12.093708166 shape(x)=(4200,) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00687909126282 multiply time=0.0013599395752 pow time/multiply time=5.05838008415 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00610399246216 multiply time=0.00104904174805 pow time/multiply time=5.81863636364 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00214910507202 multiply time=0.000160932540894 pow time/multiply time=13.3540740741 shape(x)=(4200,) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> curamp=0.5 --------------------------- pow time=0.00592398643494 multiply time=0.0013701915741 pow time/multiply time=4.32347311641 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00587201118469 multiply time=0.000941038131714 pow time/multiply time=6.23992906005 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00174999237061 multiply time=0.0001380443573 pow time/multiply time=12.677029361 shape(x)=(4200,) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00701999664307 multiply time=0.00106501579285 pow time/multiply time=6.59144839937 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00612807273865 multiply time=0.000999927520752 pow time/multiply time=6.12851692895 shape(x)=(4200, 3) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> --------------------------- pow time=0.00207996368408 multiply time=0.0001540184021 pow time/multiply time=13.5046439628 shape(x)=(4200,) x.typecode=D mat1.typecode=d mat2.typecode=d type(x[0,0])=<type 'complex'> type(mat1[0,0])=<type 'float'> type(mat2[0,0])=<type 'float'> curamp=0.75 --------------------------- David M. Cooke wrote:
On May 20, 2005, at 18:35, Ryan Krauss wrote:
I am running Numeric 24.02b on python 2.3.4. I am using Fedora Core 3 and kde.
I have modified the norm2 function to include printing both times:
def norm2(x): """compute |x|^2 = x*conjugate(x)""" if iscomplexobj(x): t1=time.time() mat1=x.real**2 + x.imag**2 t2=time.time() mat2=multiply(x.real,x.real) + multiply(x.imag,x.imag) t3=time.time() print('pow time='+str(t2-t1)) print('multiply time='+str(t3-t2)) return mat2 else: return multiply(x,x)
and get these times printed to the screen: pow time=0.00401091575623 multiply time=0.000787019729614 pow time=0.00442790985107 multiply time=0.000679016113281 pow time=0.0130319595337 multiply time=0.000170946121216 pow time=0.00465989112854 multiply time=0.000912189483643 pow time=0.00434708595276 multiply time=0.000770807266235 pow time=0.00153613090515 multiply time=0.000103950500488
Most of the calls are with 4200x3 matrices. The last one is with a 4200x1 vector (I think).
Curious. One more question: are these Complex32 matrices? I.e., are they matrices using C floats (single precision), as opposed to C doubles (double precision)? The reason is that if x is a vector of floats, x**2 upcasts x to doubles first (b/c 2 is an int -- ints are cast to doubles), so that kills the speed, and I can get timings like this. This isn't a good behaviour, and I'm going to see what I can do to fix this for the next version of Numeric.
David M. Cooke wrote:
Ryan Krauss <ryanfedora@comcast.net> writes:
No, I am running Linux.
Odd; I do too, and don't see this. glibc's pow doesn't have the same problem as Window's does. x**2 is at most 50% slower, not 10x. Which version of Numeric and python are you using, and which Linux distribution? (Numeric's version you can get from Numeric.__version__) Ideally, x**2 should run in the same time as x*x You're using scipy's multiply, which is the same (I think) as the multiply in Numeric 24.0b2, which are different from 23.8 and earlier (which does a time-wasting postcheck of the array for infinite values).
David M. Cooke wrote: Ryan Krauss <ryanfedora@comcast.net> writes: Thanks David, that seems like the right answer. But to get the speed benefit you mention, I need to use multiply (x.real,x.real) instead of x.real**2 - same for imag. There is about a factor of 10 time difference between the two. So my norm2 function now looks like this: Are you on Windows? The pow() function there, well, sucks big time. On other platforms the two are pretty close in speed. I'm working on a better version of the power ufunc for Numeric, so hopefully in a later version this will be fixed (or rather, worked around).
participants (5)
-
cookedm＠physics.mcmaster.ca -
David M. Cooke
-
Ryan Krauss -
Scott Ransom -
Stephen Walton