Exercises — Vectorization — avoiding Python loops, speed comparison

B is vectorized. In B, data * 2 hands the entire array to a single compiled C loop — the multiplication happens inside NumPy, not in the Python interpreter. A loops in Python: the for runs one interpreter step per element, and each step unboxes a np.float64 object, checks its type, and appends to a Python list. That is exactly the "reading the manual once per brick" pattern the parent note warns about.

• (a) Yes → data.sum(); a pure reduction over a numeric array.
• (b) No → file I/O and string parsing are not uniform numeric array ops; there is no contiguous float64 block to sweep.
• (c) Yes → data**2 + 1, two element-wise C loops. Rule of thumb: vectorization pays off when you have the same numeric operation applied to every element of a big packed array.

Way 1 — boolean mask (in place):

a[a < 0] = 0

a < 0 builds a boolean mask array in C (one True/False per element); indexing with it writes 0 only where the mask is True. Trade-off to notice: the mask a < 0 is itself a full-length temporary array — one boolean per element, so it costs $N$ extra elements of memory and one extra sweep of memory traffic to build and read. It is cheap (a byte, sometimes a bit, per element rather than 8 bytes) but not free, and on bandwidth-bound work that extra pass matters. That is why fusing the comparison into the write (as a[a < 0] = 0 does) is preferred over materialising the mask, storing it in a variable, and reusing it repeatedly. Way 2 — np.where (new array):

a = np.where(a < 0, 0, a)

Reads "where the condition is true pick 0, else keep a", selecting per element inside C. Note this also builds the mask internally and allocates a new output array, so it moves more memory than the in-place mask version. Both remove the Python-level if, which was the expensive part.

M = row[:, None] + col[None, :]   # shape (4, 3)

row[:, None] reshapes row to a column of shape $(4,1)$; col[None, :] reshapes col to a row of shape $(1,3)$. NumPy broadcasts them: the column is copied across 3 columns, the row across 4 rows, and they add element-wise in C — a full nested loop's worth of work with zero Python iteration.

• $T_{\text{py}}=N{\tau }_{\text{py}}=10^6\times 100\text{ ns}=10^8\text{ ns}=0.1\text{ s}.$
• $T_{\text{vec}}=c_{\text{setup}}+N{\tau }_{\text{vec}}=1000+10^6\times 1=1,001,000\text{ ns}\approx 1.001\text{ ms}.$
• $\text{Speedup}=\frac{10^8}{1,001,000}\approx 99.9\times .$ The setup cost is negligible here, so we sit essentially at the ${\tau }_{\text{py}}/{\tau }_{\text{vec}}=100\times$ ceiling.

Set $T_{\text{py}}=T_{\text{vec}}$: $N{\tau }_{\text{py}}=c_{\text{setup}}+N{\tau }_{\text{vec}}.$ WHAT we did: equated the two lines of the figure. WHY: the crossing point is exactly where neither method wins. Solve for $N$: $N({\tau }_{\text{py}}-{\tau }_{\text{vec}})=c_{\text{setup}}\Rightarrow N=\frac{c_{\text{setup}}}{{\tau }_{\text{py}}-{\tau }_{\text{vec}}}=\frac{1000}{100-1}=\frac{1000}{99}\approx \mathrm{10.1.}$ So for arrays smaller than about 10 elements the loop is actually faster — precisely the "don't vectorize tiny arrays" warning, now derived rather than asserted.

np.append loop: total copies $\approx {\sum }_{k=1}^{N}k=\frac{N(N+1)}{2}$. For $N=10^4$ that is $\frac{10^4\times 10001}{2}=50,005,000\approx 5\times 10^7\text{ element-copies.}$ This is the hallmark $O(N^2)$ behaviour (Big-O Complexity). List then np.array: the list append is $O(1)$ amortized and the final conversion copies each element once → $\approx 10^4=N$ copies, i.e. $O(N)$. Ratio $\approx \frac{5\times 10^7}{10^4}=5000\times$ less copying. The "obvious" array-growing loop is thousands of times more work.

Memory bandwidth (Memory Layout and Cache). For $N=10^7$ float64 arrays, a + b must read two arrays and write one → $3\times 10^7\times 8=2.4\times 10^8$ bytes moved. At, say, ~15 GB/s effective bandwidth that is $\approx 16$ ms — matching the observation. In the time model this means ${\tau }_{\text{vec}}$ is no longer set by the add instruction but by the time to fetch/store one element from memory. The CPU waits on data, not on math — so making the arithmetic cheaper wouldn't help. This is the ceiling the parent note calls "memory-bandwidth bound".

Prediction: essentially no speedup (~1×). np.vectorize is a convenience wrapper that still calls your Python function once per element — the loop remains in the interpreter. So ${\tau }_{\text{per-element}}$ stays at the slow Python value ${\tau }_{\text{py}}$; only $c_{\text{setup}}$ changes slightly. $\text{Speedup}=\frac{N{\tau }_{\text{py}}}{c_{\text{setup}}+N{\tau }_{\text{py}}}\underset{N\to \infty }{\to }1.$ Fix: rewrite the function using ufuncs (so the per-element cost drops from ${\tau }_{\text{py}}$ to ${\tau }_{\text{vec}}$), or compile it with Numba and Cython.

import numpy as np, timeit
for N in [1, 3, 10, 30, 100, 1000, 10000]:
    a = np.random.rand(N)
    t_loop = timeit.timeit(lambda: [x*2 for x in a], number=1000)
    t_vec  = timeit.timeit(lambda: a*2,             number=1000)
    print(N, t_loop/1000, t_vec/1000, t_loop/t_vec)

Expected curve: the speedup column starts below 1 for tiny $N$ (loop wins because $c_{\text{setup}}$ dominates the vectorized denominator), crosses 1 near the break-even $N$ from Exercise 3.1, then rises toward the ratio ${\tau }_{\text{py}}/{\tau }_{\text{vec}}$ and finally flattens at large $N$ when memory bandwidth caps ${\tau }_{\text{vec}}$. Plotting speedup vs N on a log-$N$ axis reproduces the crossing in the page figure.

Design A: $T_A=c_{\text{setup}}+N{\tau }_{\text{vec}}=1000+10^6\times 1=1,001,000\text{ ns}\approx 1.001\text{ ms}.$ Design B: each call handles $10^3$ elements: $t=1000+10^3\times 1=2000\text{ ns}$. With 1000 calls: $T_B=1000\times 2000=2,000,000\text{ ns}=2.0\text{ ms}.$ A is faster, by $T_B/T_A=2,000,000/1,001,000\approx 2.0\times .$ Lesson: fragmenting a vectorized job re-pays $c_{\text{setup}}$ on every chunk. Keep arrays whole and calls few — chunk only when memory forces you to.

Straightforward one-liner:

total = (x**2 + 3*y).sum()

Temporaries of length $N$: x**2 (1), 3*y (2), and their sum x**2 + 3*y (3) → 3 temporaries before the reduction collapses it to a scalar. Each temporary is $N\times 8=10^7\times 8=8\times 10^7$ bytes $=80$ MB. Three of them $=240$ MB of transient allocation. Lower-memory version using fused reductions:

total = np.dot(x, x) + 3.0 * y.sum()     # (sum of x_i^2) + 3*(sum y_i)

Here $\sum x_i^2=x\cdot x$ and $\sum (3y_i)=3\sum y_i$ — algebra lets each reduction stream through without building a full-length temporary. Temporaries: 0 full-length arrays. Bytes saved: $\approx 240$ MB of temporary array memory (3 × 80 MB) eliminated, and far less traffic across the memory bus — directly attacking the bandwidth bottleneck from Exercise 3.3.