# Is there a way to do array broadcasting faster, or without using loops?

I would like to optimize the following code. In real applications, the number of snapshots can be 10000, and for each of them, the lengths can be length_list = np.arange(0, 13, 0.01). So, it takes ages for it to calculate the properties of matter that I want.

import numpy as np

list_shot = list(np.arange(0, 100, 50))
list_lens = list(np.arange(10))
list_total = []

for shot in list_shot:
print("-", shot)
list_out = []

for lens in list_lens:
val = shot * lens
list_out.append(val)
print(list_out)

list_total.append(list_out)


Basically, the above code is just a simple representation of a code that I need running much faster than it's presently running.

• Why does this need to go faster? It looks like it shouldn't take that long anyways, other than maybe you keep appending. You could define the list at the start since you should know the length of it I think. Jan 4 at 19:36
• @TristanMaxson Actually, you can see in the code below the sample one, the number of snapshots can be 10000 and for each of them, the lengths can be length_list = np.arange(0, 13, 0.01). So, it takes ages for it to calculate the desired properties. Jan 4 at 19:49
• Have you strong evidence by a measurement (e.g., by a line profiler) that this snippet is the bottle neck of your computation where most time is spent (single operation step times application of this step)? How much time in this snippet is spent vs. the overall computation? How many times faster (over all steps) realistically is fast enough for you? Per your profile you know about sibling site codereview.se, it might be a better venue. Jan 4 at 21:21
• I think this should be closed and moved to a different SE after reformulating it. You need to properly profile this code and provide a test dataset etc. We cannot improve this code in its current form (and this isn't a coding SE). Jan 4 at 21:36
• Cross posted on StackOverflow: stackoverflow.com/q/70584020/1271772. Generally, it's discouraged to post a question on multiple sites without linking between the different versions. Otherwise, it's difficult to know if a question has already been answer somewhere else.
– Tyberius
Jan 5 at 13:38

What you're trying to do can actually be done very quickly, without any loops, and can even be ported to make valuable use of a GPU.

Here's a quick example that you can test yourself in octave-online:

shots=0:50:250
lens = 0:10
bsxfun(@times,shots',lens)


The outputs are as follows:

shots=[0    50   100   150   200   250]
lens =[0    1    2    3    4    5    6    7    8    9   10]
bsxfun(@times,shots',lens)
ans =

0      0      0      0      0      0      0      0      0      0      0
0     50    100    150    200    250    300    350    400    450    500
0    100    200    300    400    500    600    700    800    900   1000
0    150    300    450    600    750    900   1050   1200   1350   1500
0    200    400    600    800   1000   1200   1400   1600   1800   2000
0    250    500    750   1000   1250   1500   1750   2000   2250   2500


For shots up to 10,000 in steps of 50, the calculation takes 5.6e-5 seconds on a CPU:

tic;shots=0:50:10000;lens=0:10;bsxfun(@times,shots',lens);toc;
Elapsed time is 5.79357e-05 seconds.


With an extremely un-optimized FOR LOOP in octave-online, I repeated the calculation 100,000 times and the total CPU time was 1.97 seconds:

tic;for i=1:1e5;shots=0:50:10000;lens=0:10;bsxfun(@times,shots',lens);end;toc;
Elapsed time is 1.96675 seconds.


Believe it or not, this single function called bsxfun which stands for "binary singleton expansion function" was the basis for porting my code for calculating Feynman integrals for open quantum systems numerically from CPUs onto GPUs (though I used bsxfun for the CPU version too since this saves enormous time compared to using loops or other functions), which lead to huge speed-ups (more and more as the arrays got bigger). A paper about this: "Dattani (2013) FeynDyn: A MATLAB program for fast numerical Feynman integral calculations for open quantum system dynamics on GPUs" explains the story.

If you need to run this on a GPU you can ask a separate question about that and I'll answer that (it's been a while since I've done development for GPUs).

# In Python

In Python the task accomplished by bsxfun is known as "broadcasting" and is implemented in numpy.broadcast so that when you multiply two arrays, you automatically get nearly exactly the result that you desire here:

shots = np.array([1.0, 2.0, 3.0])
lens = np.array([2.0, 2.0, 2.0])
shots * lens
array([ 2.,  4.,  6.])


The numpy.broadcast documentation will then tell you how to generalize to work with more elements in the way that bsxfun does in MATLAB/Octave.

• These days MATLAB/Octave also natively supports broadcast operations: Try it online! Jan 5 at 11:06
• @Sanchises I was going to talk about that but forgot! That example you gave, would be excellent as another supplementary answer to this question, just like antimon's improvement to vtan707's answer! Can you write an answer including timing? Jan 5 at 16:48
• Did somebody say something about numpy on GPU? How amenable is this 2D Frenkel–Kontorova-like energy minimization problem in Python to the use of a modest PC + GPU? (Heavy reliance on indexing) I've been asked to profile it but don't know how to choose which profiler to use.
– uhoh
Jan 6 at 0:51
• @uhoh any profiler would be good enough! I'd suggest you to ask a question here with the high-performance-computing and python tags to find out which profiler to use, but currently there's a bit of a debate on our site meta about whether or not to allow such questions (myself and all 3 diamond moderators have been okay with such questions so far, but a user complained). Jan 6 at 0:59
• @NikeDattani Oh I see! Yes great idea. I can try to tailor the question to avoid some concerns in meta to the extent that I can.
– uhoh
Jan 6 at 1:08

Assuming that your sample code is really what you are looking for (which I doubt judging by the full code example you had posted originally), the answer by @vtan707 would work, but is very awkward because of its requirement for specifically shaped 2D arrays. Having to do nested calls of np.array on a single-element list of a function output kinda gives it away.

In fact, your example code is essentially a verbose version of the outer vector product where you multiply each element of one vector with each element of another one. You could simplify it to:

shot = np.arange(...)
lens = np.arange(...)

total = np.outer(shot, lens)


As per @NikeDattani's request in the comments, here are some ipython timings with their larger arrays. The first one is using np.outer, the second one reshapes the arrays into row and column vectors and lets numpy broadcast the result when using the * operator:

shots = np.arange(0,11) # 0 to 10 including
lens = np.arange(0,10001,50) # 0 to 10000 including

%timeit -n 100000 np.outer(shots, lens)
4.33 µs ± 17.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit -n 100000 shots[:,None] * lens[None,:]
3.38 µs ± 19.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


(I've used the slightly obscure indexing using None to create new dimensions here. For more info, see Dimensional Indexing Tools in the numpy docs.)

Because numpy is heavily optimized behind the scenes, there is no big performance difference between np.outer and the reshaped example -- presumably because it falls back to the exact same operations anyway. The slight overhead we see might come from function calls.

Interesting side note: The processing times increase slightly for both strategies (to about 5 to 6 µs) if I reverse the two arrays.

(Thanks also to @wflynny for doing some more timings in their answer.)

• I'm curious if you could compare the timing of this for 100,000 iterations on the larger arrays in my timed example? I know this will be comparing apples and oranges but I'm curious about the results! Jan 5 at 6:33
• Thanks for the timings! How many shots and lens do you have? I just see shot = np.arange(...) and same for lens. Jan 5 at 19:33
• Ugh, forgot to add that detail and then also realized I got your example wrong. Let me correct that. Jan 5 at 19:36

You can use dot multiplication to achieve this:

import numpy as np

shot = np.array([np.arange(0, 100, 50)]).T
lens = np.array([np.arange(10)])

total = np.dot(shot,lens)
print(total)


This gives you the below results:

[[  0   0   0   0   0   0   0   0   0   0]
[  0  50 100 150 200 250 300 350 400 450]]

• This works, but relies on some awkward array shaping. This is the exact use case that the np.outer function is for. Just use the 1D arrays from arange. Jan 5 at 0:45
• @Antimon I've mentioned np.broadcast, this answer uses np.dot and you seem to have a third (perhaps better) suggestion of using np.outer: Could you write this as a third answer? I think it would be very valuable to have that here as an answer! Jan 5 at 1:20
• @NikeDattani Sure, can do... although I'm still not sure if we are actually solving OP's problem here (see my comment on the original post). Jan 5 at 1:39

Just adding some timings for @Antimon's answer. Using numpy.outer is definitely the way to go IMO.

def list_version(shots, lens):
list_total = []
for shot in list_shot:
list_out = []

for lens in list_lens:
val = shot * lens
list_out.append(val)

list_total.append(list_out)
return list_total

def numpy_version1(shots, lens):
_shots = np.array(shots)
_lens = np.array(lens)
return np.outer(_shots, _lens)

def numpy_version2(shots, lens):
return np.outer(shots, lens)


Depending on the input type (meaning the original arrays are numpy arrays and don't need to be converted to lists), numpy.outer can give a substantial (300 fold) increase in performance:

arr_shot = np.arange(0, 10000)
arr_lens = np.linspace(0, 13, 0.01)
list_shot = list(arr_shot)
list_lens = list(arr_lens)

%timeit list_version(list_shot, list_lens)
21.5 s ± 24.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit numpy_version1(list_shot, list_lens)
73 ms ± 75.2 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit numpy_version2(arr_shot, arr_lens)
72.2 ms ± 46.6 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Method Time (ms)
Original 21500
Native Numpy 72.2
Numpy convert from list 73

In response to @Nike's comment below, I timed the numpy version over several scales of input array sizes:

results = []
for m in np.logspace(1, 4, 10, base=10, dtype=int):
for n in np.logspace(1, 4, 10, base=10, dtype=int):
x = np.random.rand(m)
y = np.random.rand(n)
t = timeit.timeit("numpy_version2(x, y)", number=10, globals=globals())
results.append((m, n, t))
data = pd.DataFrame(results, columns=["size_shots", "size_lens", "time_numpy"])
data["size_total"] = data.size_shots * data.size_lens


This should produce 2 measurements of 10 iterations at every combination of array sizes (whose log10s are equally spaced) from 10 to 10,000.

TLDR: OP's arrays are of the order of 10^5 and 10^3 elements which should take <10s (10^8 output elements takes <10s above).

• +1 and welcome to our new community! We hope to see more of you in the future, if it interests you enough! Since your final time was still 0.0722 seconds, which could be on the order of noise (in some cases), I wonder if you could time it for 100,000 repetitions like I did in my answer? I'm a bit surprised by the speed since I got: tic;for i=1:1e5;shots=0:50:10000;lens=0:10;bsxfun(@times,shots',lens);end;toc; Elapsed time is 1.96675 seconds. Maybe it's because you did arr_shot = np.arange(0, 10000); arr_lens = np.linspace(0, 13, 0.01) whereas I had only 200 different shots but 10 lens. Jan 5 at 19:31
• @NikeDattani, thanks! I found my way here through a comment on an answer I made elsewhere ~8 (wow!) years ago. But I added some additional timings of the numpy output and edited the comment! Jan 6 at 15:10