[Numpy-discussion] method to calculate the magnitude squared

[Phillip Feldman]

My apologies for the slow response; I was experiencing some technical
problems with e-mail.

In answer to Antoine's question, my main desire is for a numpy ndarray
method, for the convenience, with a secondary goal being improved
performance.

I have added the function `magsq` to my library, but would like to access
it as a method rather than as a function. I understand that I could create
a class that inherits from NumPy and add a `magsq` method to that class,
but this has a number of disadvantages.

Phillip

[Nathaniel Smith]

Hi Phillip,

My advice would be to stick with the function call. It's consistent with
most other array operations (esp. when you consider that the vast majority
of operations on arrays are functions defined in third party libraries like
yours), and the more things we add to the core array object, the more work
it is for people implementing new array-style containers. I definitely
would not recommend subclassing ndarray for this purpose -- there are all
kinds of subtle problems that you'll run into that mean it's extremely
difficult to do well, and may well be impossible to do perfectly.

Good luck,
-n

Post by Phillip Feldman
My apologies for the slow response; I was experiencing some technical
problems with e-mail.
In answer to Antoine's question, my main desire is for a numpy ndarray
method, for the convenience, with a secondary goal being improved
performance.
I have added the function `magsq` to my library, but would like to access
it as a method rather than as a function. I understand that I could create
a class that inherits from NumPy and add a `magsq` method to that class,
but this has a number of disadvantages.
Phillip
_______________________________________________
NumPy-Discussion mailing list
https://mail.scipy.org/mailman/listinfo/numpy-discussion

[Phillip Feldman]

Hello Nathaniel,

It is hard to say what is normative practice with NumPy, because there are
at least three paradigms:

(1) Some operations are implemented as methods of the `ndarray` class.
`sum` and `mean` are examples.

(2) Some operations are implemented via functions that invoke a private
method of the class. `abs` is an example of this:

In [8]: x= array([1+1J])
In [9]: x.__abs__()
Out[9]: array([ 1.41421356])

(3) Some operations are implemented as functions that operate directly on
the array, e.g., RMS (root-mean-square).

Because calculating the square of the magnitude is such a widely-used
operation, and is often done in a grossly inefficient manner (e.g., by
taking the absolute value, which involves a square-root, and then
squaring), I believe that there is a strong argument for doing either (1)
or (2). I'd prefer (1).

Phillip

Post by Nathaniel Smith
Hi Phillip,
My advice would be to stick with the function call. It's consistent with
most other array operations (esp. when you consider that the vast majority
of operations on arrays are functions defined in third party libraries like
yours), and the more things we add to the core array object, the more work
it is for people implementing new array-style containers. I definitely
would not recommend subclassing ndarray for this purpose -- there are all
kinds of subtle problems that you'll run into that mean it's extremely
difficult to do well, and may well be impossible to do perfectly.
Good luck,
-n

Post by Phillip Feldman
My apologies for the slow response; I was experiencing some technical
problems with e-mail.
In answer to Antoine's question, my main desire is for a numpy ndarray
method, for the convenience, with a secondary goal being improved
performance.
I have added the function `magsq` to my library, but would like to access
it as a method rather than as a function. I understand that I could create
a class that inherits from NumPy and add a `magsq` method to that class,
but this has a number of disadvantages.
Phillip
_______________________________________________
NumPy-Discussion mailing list
https://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
https://mail.scipy.org/mailman/listinfo/numpy-discussion

[Charles R Harris]

On Fri, Oct 9, 2015 at 10:32 PM, Phillip Feldman <

Post by Phillip Feldman
Hello Nathaniel,
It is hard to say what is normative practice with NumPy, because there are
(1) Some operations are implemented as methods of the `ndarray` class.
`sum` and `mean` are examples.
(2) Some operations are implemented via functions that invoke a private
In [8]: x= array([1+1J])
In [9]: x.__abs__()
Out[9]: array([ 1.41421356])
(3) Some operations are implemented as functions that operate directly on
the array, e.g., RMS (root-mean-square).
Because calculating the square of the magnitude is such a widely-used
operation, and is often done in a grossly inefficient manner (e.g., by
taking the absolute value, which involves a square-root, and then
squaring), I believe that there is a strong argument for doing either (1)
or (2). I'd prefer (1).

We tend to avoid adding methods. 2) would be a very easy enhancement, just
a slight modification of sqr.

Chuck

[Marten van Kerkwijk]

Post by Charles R Harris
We tend to avoid adding methods. 2) would be a very easy enhancement,

just a slight modification of sqr.

Did you mean `np.square`? Sadly, that doesn't do the right thing:
`np.square(1+1j)` yields `2j`, while one wants `c*c.conj()` and thus `2`.
Or, for fastest speed, really just `c.real**2 + c.imag**2`.

My guess would be that a new ufunc, say `np.abs2` or `np.modulus2` or so,
would be more appropriate than defining a new method. I'd also be hesitant
to define a new private method -- I like how those usually are just used to
override python basics.

-- Marten

[Charles R Harris]

On Sat, Oct 10, 2015 at 11:14 AM, Marten van Kerkwijk <

Post by Marten van Kerkwijk

Post by Charles R Harris
We tend to avoid adding methods. 2) would be a very easy enhancement,

just a slight modification of sqr.
`np.square(1+1j)` yields `2j`, while one wants `c*c.conj()` and thus `2`.
Or, for fastest speed, really just `c.real**2 + c.imag**2`.

Yes, I meant the new function could made by reusing the square code with
slight modifications.

Post by Marten van Kerkwijk
My guess would be that a new ufunc, say `np.abs2` or `np.modulus2` or so,
would be more appropriate than defining a new method. I'd also be hesitant
to define a new private method -- I like how those usually are just used to
override python basics.

Julia uses abs2.

Chuck

[Nathaniel Smith]

Post by Charles R Harris
On Sat, Oct 10, 2015 at 11:14 AM, Marten van Kerkwijk <

Post by Marten van Kerkwijk

Post by Charles R Harris
We tend to avoid adding methods. 2) would be a very easy enhancement,

just a slight modification of sqr.
`np.square(1+1j)` yields `2j`, while one wants `c*c.conj()` and thus `2`.
Or, for fastest speed, really just `c.real**2 + c.imag**2`.

Post by Charles R Harris
Yes, I meant the new function could made by reusing the square code with

slight modifications.

Post by Charles R Harris

Post by Marten van Kerkwijk
My guess would be that a new ufunc, say `np.abs2` or `np.modulus2` or

so, would be more appropriate than defining a new method. I'd also be
hesitant to define a new private method -- I like how those usually are
just used to override python basics.

Post by Charles R Harris
Julia uses abs2.

I don't have an opinion on whether abs2 is important enough to bother with
(I don't work much with complex numbers myself, nor have I run any
benchmarks), but I agree that if we do want it then adding it as a regular
ufunc would definitely be the right approach.

-n

[Phillip Feldman]

The ufunc approach makes sense.

Something like abs2 is essential for anyone who does signal processing
simulations using NumPy.

Phillip

Post by Marten van Kerkwijk

Post by Charles R Harris
On Sat, Oct 10, 2015 at 11:14 AM, Marten van Kerkwijk <

Post by Marten van Kerkwijk

Post by Charles R Harris
We tend to avoid adding methods. 2) would be a very easy enhancement,

just a slight modification of sqr.
`np.square(1+1j)` yields `2j`, while one wants `c*c.conj()` and thus `2`.
Or, for fastest speed, really just `c.real**2 + c.imag**2`.

Post by Charles R Harris
Yes, I meant the new function could made by reusing the square code with

slight modifications.

Post by Charles R Harris

Post by Marten van Kerkwijk
My guess would be that a new ufunc, say `np.abs2` or `np.modulus2` or

so, would be more appropriate than defining a new method. I'd also be
hesitant to define a new private method -- I like how those usually are
just used to override python basics.

Post by Charles R Harris
Julia uses abs2.

I don't have an opinion on whether abs2 is important enough to bother with
(I don't work much with complex numbers myself, nor have I run any
benchmarks), but I agree that if we do want it then adding it as a regular
ufunc would definitely be the right approach.
-n
_______________________________________________
NumPy-Discussion mailing list
https://mail.scipy.org/mailman/listinfo/numpy-discussion

[Nils Becker]

Hey,

I use complex numbers a lot and obviously need the modulus a lot. However,
I am not sure if we need a special function for _performance_ reasons.

At 10:01 AM 9/20/2015, you wrote:

It is, but since that involves taking sqrt, it is *much* slower. Even now,
```
In [32]: r = np.arange(10000)*(1+1j)

In [33]: %timeit np.abs(r)**2
1000 loops, best of 3: 213 µs per loop

In [34]: %timeit r.real**2 + r.imag**2
10000 loops, best of 3: 47.5 µs per loop

This benchmark is not quite fair as the first example needs a python
function call and the second doesn't. If you benchmark a modulus function
against np.abs(x)**2 the performance gain is ca. 30% on my machine. This
means that for such a basic operation most of the time is spent in the
function call.
In my opinion if you want to have speed you write the modulus explicitly in
your expression (3-4x speedup on my machine). If you don't need speed you
can afford the function call (be it to abs2 or to abs).

By not providing abs2 in numpy, however, people do not loose out on a lot
of performance...

There may be reasons to provide abs2 related to accuracy. If people (for
not knowing it better) use np.abs(x)**2 they lose significant digits I
think (may be wrong on that...). I did not look into it, though.

Cheers
Nils

[Marten van Kerkwijk]

Hi Nils,

I think performance will actually be better than I indicated, especially
for larger arrays, since `r.real**2 + r.imag**2` makes a quite unnecessary
intermediate arrays. With a `ufunc`, this can be done much faster. Indeed,
it should be no slower than `np.square` (which does more operations):

%timeit b = np.square(a)
100000 loops, best of 3: 16.6 µs per loop

[This is on same laptop as the timings above.]

Chuck: agreed that a future np.abs2 could just reuse the internal ufunc
loops for np.square except for the complex case.