[Numpy-discussion] "Extended" Outer Product

Gael Varoquaux gael.varoquaux at normalesup.org
Wed Aug 22 03:45:46 EDT 2007

On Tue, Aug 21, 2007 at 02:14:00PM -0700, Timothy Hochberg wrote:
>    I suppose someone should fix that someday. However, I still think
>    vectorize is an attractive nuisance in the sense that someone has a
>    function that they want to apply to an array and they get sucked into
>    throwing vectorize at the problem. More often than not, vectorize makes
>    things slower than they need to be. If you don't care about performance,
>    that's fine, but I live in fear of code like:

>       def f(a, b):
>           return sin(a*b + a**2)
>       f = vectorize(f)

>    The original function f is a perfectly acceptable vectorized function
>    (assuming one uses numpy.sin), but now it's been replaced by a slower
>    version by passing it through vectorize. To be sure, this isn't always the
>    case; in cases where you have to make choices, things get messier. Still,
>    I'm not convinced that vectorize doesn't hurt more than it helps.

I often have code where I am going to loop over a large amount of nested
loops, some thing like:

# A function to return the optical field in each point:

def optical_field( (x, y, z) ):
    loop over an array of laser wave-vector
    return optical field

# Evaluate the optical field on a grid to plot it :

x, y z = mgrid[-10:10, -10:10, -10:10]
field = optical_field( (x, y, z) )

In such a code every single operation could be vectorized, but the
problem is that each function assumes the input array to be of a certain
dimension: I may be using some code like:
r = c_[x, y, z]
cross(r, r_o) 

So implementing loops with arrays is not that convenient, because I have
to add dimensions to my arrays, and to make sure that my inner functions
are robust to these extra dimensions.

Looking at some of my code where I had this kind of problems, I see
functions similar to:

def delta(r, v, k):
    return  dot(r, transpose(k))      
            + Gaussian_beam(r)
            + dot(v, transpose(k))

I am starting to realize that the real problem is that there is no info
of what the expected size for the input and output arguments should be.
Given such info, the function could resize its input and output

Maybe some clever decorators could be written to address this issue,
something like:

@inputsize( (3, -1), (3, -1), (3, -1) )

which would reshape every input positional argument to the shape given in
the list of shapes, and reshape the output argument to the shape of the
first input argument.

As I worked around these problems in my code I cannot say whether these
decorators would get rid of them (I had not had the idea at the time), I
like the idea, and I will try next time I run into these problems.

I just wanted to point out that replacing for loops with arrays was not
always that simple and that using "vectorize" sometimes was a quick and a
dirty way to get things done.


More information about the NumPy-Discussion mailing list