Subtracting a shape-(3,) mean from a shape-(5, 3) matrix should raise an error. The shapes do not match. NumPy does it anyway, without complaining, without copying data. If you do not know why, you will write code that works by accident until the shapes change and it suddenly fails in production with no obvious cause.
The answer is broadcasting: NumPy’s rule for combining arrays of different shapes by stretching size-1 dimensions without copying data. Once you understand that rule, the gap between code that works and code that works quickly closes fast.
The key concept is broadcasting. Once you understand it, the rest follows: memory layout, vectorisation, linear algebra, saving and loading arrays, and the gotchas that bite you when shapes are almost right but not quite. Chapter 6 built the array fundamentals; this chapter shows what those arrays can do at speed. By the end you will have everything you need for Chapter 14: matplotlib visualisations of the arrays you have been building.
Callout markers: book cover page.
Re-run this cell to restore the imports and student data matrix from Chapter 6.
import numpy as np
# Student data matrix: 5 students x 3 measurements (study_hours, attendance_pct, prior_gpa)
student_data = np.array(
[
[12.0, 85.0, 3.1],
[5.0, 60.0, 2.4],
[18.0, 95.0, 3.8],
[9.0, 70.0, 2.9],
[22.0, 98.0, 3.9],
]
)
column_names = ["study_hours", "attendance_pct", "prior_gpa"]
print(f"student_data.shape={student_data.shape} student_data.dtype={student_data.dtype}")student_data.shape=(5, 3) student_data.dtype=float64Broadcasting is the rule NumPy uses to apply an operation between two arrays of different shapes, by virtually “stretching” the smaller one, without actually copying any data. It is what lets you write student_data - student_data.mean(axis=0) instead of a loop over rows.
Key Concept: the broadcasting rule
Compare shapes from the right-hand side. Two dimensions are compatible when they are equal, or when one of them is 1 (it gets stretched to match). Missing leading dimensions are treated as 1.
(5, 3) and (3,) → treat (3,) as (1, 3) → stretch to (5, 3). (valid)
(5, 3) and (5,) → treat (5,) as (1, 5) → 3 != 5 and neither is 1. (error)
The code below normalises each column of the student data matrix. (5, 3) - (3,) broadcasts without a loop: NumPy treats the (3,) mean as (1, 3) and stretches it across all five rows.
student_data = np.array(
[
[12.0, 85.0, 3.1],
[5.0, 60.0, 2.4],
[18.0, 95.0, 3.8],
[9.0, 70.0, 2.9],
[22.0, 98.0, 3.9],
]
)
col_mean = student_data.mean(axis=0) # shape (3,): one mean per column
col_std = student_data.std(axis=0) # shape (3,)
print(f"student_data.shape : {student_data.shape}")
print(f"col_mean.shape : {col_mean.shape}")
# (5, 3) - (3,) broadcasts the mean across every row: no loop needed
student_data_z = (student_data - col_mean) / col_std
print(f"normalised:\n{student_data_z}")
print(f"new per-column mean (~0): {student_data_z.mean(axis=0).round(6)}")student_data.shape : (5, 3)
col_mean.shape : (3,)
normalised:
[[-0.196326 0.23338089 -0.21293203]
[-1.34156098 -1.48265509 -1.45503555]
[ 0.78530399 0.91979529 1.02917149]
[-0.68714099 -0.7962407 -0.56781875]
[ 1.43972398 1.1257196 1.20661485]]
new per-column mean (~0): [ 0. 0. -0.]The rule is: align shapes from the right. Dimensions of size 1 stretch to match. Everything else must match exactly or NumPy raises an error.
Common Mistake: broadcasting a (5,) array against a (5, 3) matrix fails for a subtle reason
If you compute per_row_mean = student_data.mean(axis=1) (shape (5,)) and try student_data - per_row_mean, NumPy raises ValueError: operands could not be broadcast together. It is comparing the trailing dimensions 3 vs 5, not what you intended. Fix it by giving the per-row result an explicit column shape with keepdims=True: student_data.mean(axis=1, keepdims=True) has shape (5, 1), which broadcasts correctly against (5, 3).
row_mean = student_data.mean(axis=1, keepdims=True) # shape (5, 1), NOT (5,)
print(f"row_mean.shape : {row_mean.shape}")
# (5, 3) - (5, 1) broadcasts the per-row mean across every column
centered_per_row = student_data - row_mean
print(f"centered_per_row:\n{centered_per_row}")row_mean.shape : (5, 1)
centered_per_row:
[[-21.36666667 51.63333333 -30.26666667]
[-17.46666667 37.53333333 -20.06666667]
[-20.93333333 56.06666667 -35.13333333]
[-18.3 42.7 -24.4 ]
[-19.3 56.7 -37.4 ]]Every NumPy array is, at bottom, a 1D block of bytes. arr.strides exposes the step sizes that let NumPy navigate a multi-dimensional shape through that flat block. arr.itemsize is the number of bytes in one element.
Key Concept: strides map a multi-dimensional shape onto flat memory
For a float64 array of shape (5, 3): itemsize is 8 bytes, strides are (24, 8). Moving one column right costs 8 bytes (= itemsize). Moving one row down costs 24 bytes (= 3 columns × 8 bytes). Basic slicing adjusts strides without copying data: that is why reshape() and [::-1] return views instantly.
student_data = np.array(
[
[12.0, 85.0, 3.1],
[5.0, 60.0, 2.4],
[18.0, 95.0, 3.8],
[9.0, 70.0, 2.9],
[22.0, 98.0, 3.9],
]
)
print(f"shape : {student_data.shape}")
print(f"itemsize : {student_data.itemsize} bytes per element") # 8 for float64
print(f"strides : {student_data.strides}") # (bytes/row, bytes/col) = (24, 8)
# A reversed view uses a negative stride: no data is copied
rev = student_data[::-1]
print(f"reversed strides : {rev.strides}") # (-24, 8)
print(f"shares memory original : {np.shares_memory(rev, student_data)}")shape : (5, 3)
itemsize : 8 bytes per element
strides : (24, 8)
reversed strides : (-24, 8)
shares memory original : True Activity 1: min-max scale every column
student_data to the [0, 1] range using the formula (student_data - student_data.min(axis=0)) / (student_data.max(axis=0) - student_data.min(axis=0)). Confirm the result’s per-column min is 0 and max is 1.
student_data_scaled.min(axis=0) # -> array([0., 0., 0.]) student_data_scaled.max(axis=0) # -> array([1., 1., 1.])
student_data = np.array(
[
[12.0, 85.0, 3.1],
[5.0, 60.0, 2.4],
[18.0, 95.0, 3.8],
[9.0, 70.0, 2.9],
[22.0, 98.0, 3.9],
]
)
student_data_scaled = ... # TODO: min-max scale every column to [0, 1]
# print(f"min per column: {student_data_scaled.min(axis=0)}")
# print(f"max per column: {student_data_scaled.max(axis=0)}")Now that you can express normalisation as (student_data - mean) / std with no loop, it’s worth seeing why that matters. Every NumPy operation you have used so far runs as a single compiled loop over contiguous memory; a Python for loop pays the cost of the interpreter on every single element.
Key Concept: replace loops with vectorised expressions, 10-100x faster
A Python loop over one million floats takes ~500ms. The same operation as (arr - arr.mean()) / arr.std() takes ~5ms. NumPy calls compiled C code: the Python interpreter never enters the inner loop. The speedup grows with array size.
import time
big = np.random.default_rng(0).normal(size=1_000_000)
def zscore_loop(values: np.ndarray) -> list[float]:
mean = sum(values) / len(values)
variance = sum((v - mean) ** 2 for v in values) / len(values)
std = variance**0.5
return [(v - mean) / std for v in values]
start = time.perf_counter()
_ = zscore_loop(big)
loop_time = time.perf_counter() - start
start = time.perf_counter()
_ = (big - big.mean()) / big.std()
vector_time = time.perf_counter() - start
print(f"Python loop time : {loop_time:.4f}s")
print(f"Vectorised time : {vector_time:.4f}s")
print(f"Speedup : {loop_time / vector_time:,.0f}x")Python loop time : 0.4365s
Vectorised time : 0.0144s
Speedup : 30x Pro Tip: if you’re writing a for loop over an array, stop and look for a vectorised way
Almost every elementwise transformation, filter, or aggregation you would write as a Python loop already has a NumPy equivalent: arithmetic operators, np.where, boolean masks, axis aggregations. Reach for those first: a hand-written loop over a large array is one of the most common NumPy performance bugs, and it’s usually a 10-100x slowdown for no benefit.
Activity 2: benchmark loop vs vectorised z-score
Measure the speedup of vectorisation.
1. Create arr = np.random.default_rng(0).normal(size=100_000).
2. Implement a Python loop z-score (iterate element by element).
3. Implement the vectorised version: (arr - arr.mean()) / arr.std().
4. Time both with import time and print the speedup ratio.
Expected: vectorised is 50-200x faster.
A linear model’s prediction is a dot product: multiply each measurement by a learned weight, sum the results, add a bias. @ (matrix multiplication) computes this for every student row at once, no loop required.
Key Concept: @ is matrix multiplication; use it, never loop
student_data @ w computes every row’s dot product with w in one call. Manually looping for row in student_data: sum(row * w) is correct but 100x slower. Every linear model, transformer attention layer, and PCA uses @ internally.
student_data = np.array(
[
[12.0, 85.0, 3.1],
[5.0, 60.0, 2.4],
[18.0, 95.0, 3.8],
[9.0, 70.0, 2.9],
[22.0, 98.0, 3.9],
]
) # shape (5 students, 3 measurements)
# Suppose a (already-fitted) linear model has these learned weights and bias
weights = np.array([1.5, 0.3, 8.0]) # one weight per column, shape (3,)
bias = 10.0
# student_data @ weights: (5, 3) @ (3,) -> (5,): one prediction per student
predicted_scores = student_data @ weights + bias
print(f"predicted_scores: {predicted_scores.round(1)}")predicted_scores: [ 78.3 54.7 95.9 67.7 103.6]@ on a matrix and a vector is shorthand for: for each row, multiply element-wise by weights and sum: exactly (student_data * weights).sum(axis=1). Verify the two are equivalent, then measure prediction error against the true scores with np.linalg.norm (the Euclidean / RMS-style distance):
# @ is equivalent to elementwise multiply + sum along axis=1
manual = (student_data * weights).sum(axis=1) + bias
print(f"@ matches manual sum: {np.allclose(predicted_scores, manual)}")
actual_scores = np.array([88.0, 65.0, 95.0, 78.0, 99.0])
errors = predicted_scores - actual_scores
rmse = np.linalg.norm(errors) / np.sqrt(len(errors))
print(f"errors : {errors.round(1)}")
print(f"RMSE : {rmse:.2f}")@ matches manual sum: True
errors : [ -9.7 -10.3 0.9 -10.3 4.6]
RMSE : 8.10 Common Mistake: comparing floats with ==
np.allclose(a, b) was used above instead of (a == b).all() on purpose: floating-point arithmetic accumulates tiny rounding errors, so two mathematically-equal results can differ in their last bit. Always compare floats with a tolerance: np.allclose, or abs(a - b) < 1e-9, never with exact ==.
Activity 3: linear prediction with @
Use matrix multiplication to compute predictions.
Using student_data from the setup section:
1. Define weights = np.array([0.3, 0.4, 1.5]) and bias = -2.0.
2. Compute predicted scores as preds = student_data @ weights + bias.
3. Confirm preds.shape == (5,).
4. Print the index of the highest-scoring student using np.argmax(preds).
A typical pipeline computes a student data matrix once and reuses it across many later steps (training, evaluation, serving). .npy stores a single array in NumPy’s own binary format: far smaller and faster to read than CSV for numeric data, and it preserves dtype and shape exactly. .npz bundles several named arrays together.
Key Concept: use .npy for one array, .npz for multiple
np.save(‘student_data.npy’, student_data) and np.load(‘student_data.npy’) preserve dtype and shape exactly. np.savez(‘data.npz’, features=student_data, target=y) bundles multiple arrays in one file. Both load in microseconds: much faster than re-parsing a CSV.
from pathlib import Path
tmp_dir = Path("tmp_numpy_activity")
tmp_dir.mkdir(exist_ok=True)
student_data = np.array([[12.0, 85.0, 3.1], [5.0, 60.0, 2.4], [18.0, 95.0, 3.8]])
y = np.array([88.0, 65.0, 95.0])
# Save a single array
np.save(tmp_dir / "student_data.npy", student_data)
student_data_loaded = np.load(tmp_dir / "student_data.npy")
print(f"round-trip equal: {np.array_equal(student_data, student_data_loaded)}")
# Save several named arrays together in one file
np.savez(tmp_dir / "dataset.npz", features=student_data, target=y)
bundle = np.load(tmp_dir / "dataset.npz")
print(f"keys : {list(bundle.keys())}")
print(f"target : {bundle['target']}")round-trip equal: True
keys : ['features', 'target']
target : [88. 65. 95.]import shutil
shutil.rmtree(tmp_dir)
print(f"cleaned up: {tmp_dir.exists()}")cleaned up: FalseLike the Python gotchas in Chapter 4, none of these raise an exception. They silently produce a wrong (or surprising) result. Recognise them now so you do not lose hours to them later.
Key Concept: three gotchas that never raise an exception
1. View vs copy: sub = student_data[:2] shares memory; mutating sub mutates student_data. Fix: student_data[:2].copy().
2. Integer overflow: np.array([120, 10], dtype=np.int8) wraps silently. Fix: check dtype before arithmetic.
3. Shape mismatch in broadcasting: (5, 3) - (3, 5) may broadcast unexpectedly. Always check shapes before subtracting.
# GOTCHA 1: basic slicing returns a VIEW, not a copy
scores = np.array([62, 78, 85, 91, 55])
top_three = scores[:3]
top_three[0] = 0 # mutating the "slice"...
print(f"top_three : {top_three}")
print(f"scores : {scores}") # ...changed the original too!
# Fix: copy explicitly when you need an independent array
safe_copy = scores[:3].copy()top_three : [ 0 78 85]
scores : [ 0 78 85 91 55]# GOTCHA 2: integer dtype truncates on division-like ops, and can overflow
small = np.array([120, 10], dtype=np.int8)
print(f"int8 + 50 : {small + 50}") # 120 + 50 = 170, overflows int8's max of 127!
# Fix: use a wide-enough dtype, or let NumPy infer (default is int64/float64)
safe = small.astype(np.int64) + 50
print(f"int64 + 50: {safe}")int8 + 50 : [-86 60]
int64 + 50: [170 60]# GOTCHA 3: {} is a dict, not a set: same trap as in plain Python (Chapter 4)
empty_dict = {}
empty_set = set()
print(f"type({{}}) : {type(empty_dict)}")
print(f"type(set()) : {type(empty_set)}")type({}) : <class 'dict'>
type(set()) : <class 'set'># GOTCHA 4: comparing floats with == (see Sec. 9 for the fix: np.allclose)
a = np.array([0.1 + 0.2])
print(f"0.1 + 0.2 == 0.3 : {a == 0.3}") # False: 0.30000000000000004 != 0.3
print(f"np.allclose : {np.allclose(a, 0.3)}") # True: tolerant comparison0.1 + 0.2 == 0.3 : [False]
np.allclose : TrueApply everything from this notebook together. Each exercise is self-contained.
Activity 4: vectorised anomaly detector
2 standard deviations from the overall mean using a single boolean mask.
readings = [36.5, 36.7, 36.8, 36.6, 36.9, 39.5, 36.7, 36.8] # Expected: reading 39.5 (index 5) flagged as anomaly
Hint: z = (readings - readings.mean()) / readings.std(), then mask np.abs(z) > 2.
readings = np.array([36.5, 36.7, 36.8, 36.6, 36.9, 39.5, 36.7, 36.8])
z_scores = ... # TODO
anomaly_mask = ... # TODO
if z_scores is not ...:
print(f"z_scores : {z_scores.round(2)}")
print(f"anomaly_mask : {anomaly_mask}")
print(f"anomalies : {readings[anomaly_mask]}")
else:
print("(implement z_scores and anomaly_mask above to see output)")(implement z_scores and anomaly_mask above to see output)NumPy 2.0 (released June 2024) is the first major version bump in almost two decades. Two changes matter most for day-to-day data science work:
The old Python-builtin aliases (np.int, np.float, np.bool, np.complex, np.object, np.str) were deprecated for years and are now fully removed. They were just aliases to the Python built-ins anyway, so the fix is mechanical:
| Old (removed) | Replacement |
|---|---|
np.bool |
np.bool_ or Python bool |
np.int |
np.intp or Python int |
np.float |
np.float64 or Python float |
np.complex |
np.complex128 or Python complex |
np.object |
np.object_ or Python object |
np.str |
np.str_ or Python str |
StringDType for proper string arraysNumPy 2.0 introduced np.dtypes.StringDType(), a real variable-length string dtype backed by UTF-8 memory. The old np.str_ stored fixed-width UCS-4 strings (one array dtype for every character count), StringDType stores arbitrary-length strings efficiently.
Pro Tip: use np.float64 not np.float in new code
If you see a module ‘numpy’ has no attribute ‘float’ error, the codebase was written against NumPy 1.x. A global search-and-replace of np.float → np.float64 (and so on for the others in the table) is the complete fix.
import numpy as np
# Old aliases are removed in NumPy 2.0: this would raise AttributeError:
# arr = np.array([1, 2, 3], dtype=np.float) # removed in NumPy 2.0
# Use the explicit dtype names instead:
arr = np.array([1, 2, 3], dtype=np.float64) # correct
print(f"dtype: {arr.dtype}")
# StringDType: variable-length strings, efficient UTF-8 storage (NumPy 2.0+)
names = np.array(["Alice", "Bob", "Charlie"], dtype=np.dtypes.StringDType())
print(f"names : {names}")
print(f"dtype : {names.dtype}")
# The old fixed-width str_ is still available but StringDType is the modern choice
old_style = np.array(["Alice", "Bob", "Charlie"]) # infers str_ (fixed width)
print(f"old dtype: {old_style.dtype}") # <U7 means Unicode, 7 chars widedtype: float64
names : ['Alice' 'Bob' 'Charlie']
dtype : StringDType()
old dtype: <U7np.strings: vectorised string operationsNumPy 2.0 ships a dedicated np.strings module for vectorised string operations on np.ndarray values backed by StringDType. Before this module existed, string operations on object arrays required Python-level loops or np.vectorize, both of which are slow. np.strings routes each operation through a compiled C loop, making it roughly 10x faster on large arrays than the equivalent Python loop.
The module mirrors the Python str method set: np.strings.upper(), np.strings.lower(), np.strings.split(), np.strings.startswith(), np.strings.replace(), and others all work element-wise, broadcasting across the entire array in one call.
import numpy as np
# Create a StringDType array of course codes
course_arr = np.array(["DATA101", "MATH201", "STATS301", "ML401"], dtype=np.dtypes.StringDType())
# np.strings.upper() broadcasts across the whole array in one compiled call
upper = np.strings.upper(course_arr)
print("upper:", upper)
# np.strings.startswith() returns a boolean array, usable as a mask
prefix_mask = np.strings.startswith(course_arr, "DATA")
print("DATA prefix mask:", prefix_mask)
print("courses starting with 'DATA':", course_arr[prefix_mask])upper: ['DATA101' 'MATH201' 'STATS301' 'ML401']
DATA prefix mask: [ True False False False]
courses starting with 'DATA': ['DATA101'] Pro Tip: use np.dtypes.StringDType() for string arrays that need fast operations
The default np.array([“a”, “b”]) stores strings as object dtype: each element is a Python object pointer, so any operation loops at the Python level. np.dtypes.StringDType() stores strings as contiguous UTF-8 bytes, avoids the object-dtype overhead, and unlocks the entire np.strings.* API for compiled, vectorised string work.
| Resource | Why it matters |
|---|---|
| Harris, C.R. et al. (2020). Array programming with NumPy. Nature 585, 357–362. | The primary citation for NumPy; the paper explains the design decisions behind broadcasting and ufuncs |
| VanderPlas, J. (2016). Python Data Science Handbook, Ch. 2. O’Reilly. | Free online: the most readable treatment of fancy indexing, broadcasting, and structured arrays |
| NumPy documentation: Broadcasting | Official broadcasting rules with diagrams; bookmark for the next time the shapes don’t align |
| NumPy documentation: Indexing | Covers basic, advanced, and boolean indexing in one place |
| Concept | Key rule |
|---|---|
ndarray |
One dtype, contiguous memory, fast because it skips the Python interpreter |
| Creation | np.array, arange, linspace, zeros/ones; np.random.default_rng(seed) for reproducible random data |
shape/dtype |
Always check before trusting a result; mixed-type input silently upcasts |
reshape |
Returns a view, same data, same size, different shape |
flatten |
Always returns a copy |
| Slicing | Basic slices (student_data[:, 0]) are views; fancy/boolean indexing always copies |
| Boolean masks | & / \| (not and/or) on arrays, each side parenthesised |
np.where / np.select |
Vectorised if/else and multi-branch labelling |
axis=0 vs axis=1 |
0 collapses rows -> one value per column; 1 collapses columns -> one value per row |
| Broadcasting | Compare shapes right-to-left; dims match if equal or one is 1; use keepdims=True to broadcast a per-row stat |
| Vectorisation | A NumPy expression beats a Python loop by 10-100x, look for one before writing for |
@ / np.dot |
Matrix multiplication: student_data @ weights predicts every row in one call |
np.allclose |
Always compare floats with a tolerance, never == |
.npy / .npz |
Compact, dtype/shape-preserving array storage between pipeline stages |
Next: Chapter 14 (14-matplotlib.ipynb), covering how to visualise arrays and DataFrames with matplotlib.