NumPy Fundamentals

The ndarray: what makes it different from a list

In M47 you represented vectors as Python lists, and saw that a pure-Python matrix multiply spends nearly all its time on interpreter overhead. A list stores pointers to separate Python objects scattered across memory, and each object carries its own type, reference count, and value. That flexibility (any mix of types, any size) is exactly what makes lists slow for arithmetic.

A NumPy ndarray ("n-dimensional array") makes three promises instead:

• Homogeneous: every element has the same dtype (data type), such as float64 (8-byte float) or int32 (4-byte integer).
• Contiguous: the raw values sit side by side in one block of memory (the buffer), with no per-element Python objects.
• Typed metadata: a shape (how many elements along each dimension) and strides (how many bytes to step to reach the next element along each dimension) describe how to interpret that block.

Because the type is known once for the whole array and the values are contiguous, NumPy can run its loops in compiled C, use SIMD instructions, and stay friendly to the CPU cache. The price is that every element must be the same type and the size is fixed at creation.

NumPy upcasts to a single type that can hold everything: integers plus a float become floats; anything with a string becomes strings (and arithmetic on that array then fails). Mixed data is a sign you want a pandas DataFrame (M51), which stores one array per column.

Remember: ndarray = one dtype + one buffer + shape + strides. Operators work element-wise.

Creating arrays

You rarely type arrays out by hand. These constructors cover most needs:

Remember: shapes are tuples; linspace for float ranges; default_rng(seed) for random numbers.

dtypes, overflow, and precision

A dtype fixes how many bytes each element uses and how those bits are read. Common ones:

Python's own int never overflows: it grows as needed. NumPy integers are fixed-width hardware integers, so they wrap around silently when an array operation exceeds their range, exactly like a car odometer rolling over. Floats don't wrap; they lose precision, and they overflow to inf.

Remember: NumPy ints wrap silently on overflow; floats round. Pick the dtype on purpose and convert with astype before arithmetic that might exceed the range.

Shape, reshape, and axis

The shape is a tuple with one entry per dimension: (5,) is a 1-D array of 5 elements, (3, 4) is 3 rows by 4 columns, (10, 28, 28) could be ten 28 × 28 grayscale images. ndim is the length of the shape. reshape reinterprets the same elements in a new shape, as long as the total count is unchanged; pass -1 for one dimension and NumPy computes it.

Many functions take an axis argument. The rule that never fails: axis=k is the dimension that disappears. For a 2-D array, sum(axis=0) collapses the rows (axis 0), going down each column and leaving one result per column; sum(axis=1) collapses the columns, going across each row and leaving one result per row.

Remember: axis=k is the dimension that vanishes. keepdims=True keeps it as size 1, which is handy for broadcasting the result back.

Indexing: basic, slicing, boolean masks, fancy

NumPy indexing extends list indexing to many dimensions, separated by commas inside one pair of brackets: m[i, j] rather than m[i][j].

• Basic: m[1, 2] is one element; m[1] (or m[1, :]) is row 1; m[:, 2] is column 2; negative indices count from the end.
• Slicing: start:stop:step in each dimension, e.g. m[::2, 1:]. Slices return views (next section).
• Boolean mask: an array of True/False with the same shape selects the elements where it is True: a[a > 0]. Combine conditions with &, |, ~, and parentheses (not and/or).
• Fancy (integer-array) indexing: a list or array of indices picks those positions, in that order, possibly repeated: a[[3, 0, 0]].

Remember: m[i, j], slices per axis, masks for "where condition", integer arrays for "these positions".

Views, copies, and strides

A view is a new array object that shares the same buffer as the original, with different shape or strides. Basic slicing, reshape (when possible), .T, and ravel (when possible) return views: no data is copied, so they are instant even for huge arrays, but writing to a view changes the original. Boolean and fancy indexing always return copies, because the selected elements generally can't be described by a regular stride pattern.

This is the opposite of lists, where lst[1:3] makes a copy. It surprises everyone once.

Strides explained

Strides are the byte steps along each axis. For a (3, 4) array of 8-byte floats in row-major order, moving one column right is 8 bytes and moving one row down is 4 × 8 = 32 bytes: strides (32, 8). The transpose just swaps them to (8, 32), which is why .T is free. Taking every other column, m[:, ::2], doubles the column stride to 16. The element at [i, j] lives at byte offset i·stride₀ + j·stride₁: exactly the row-major formula from M47, generalized.

Remember: slices and transposes are views (shared, free); masks and fancy indexing are copies. Strides say how to walk the buffer.

Vectorization, ufuncs, and aggregation

Vectorization means expressing a computation on whole arrays instead of writing a Python loop over elements. The loop still happens, but inside compiled code. NumPy's element-wise functions are called ufuncs (universal functions): arithmetic operators, np.sqrt, np.exp, np.log, np.sin, np.maximum, comparisons, and many more. They work on any shape and return an array of the same shape.

Aggregations (reductions) collapse values: sum, mean, std, min, max, argmin, argmax, any, all, all with an optional axis. NaN-aware versions (np.nanmean, np.nansum) skip missing values.

Remember: write the formula once for the whole array. If you find yourself writing for i in range(len(a)) over NumPy data, look for a ufunc, a reduction, or broadcasting.

Broadcasting, step by step

Broadcasting is how NumPy combines arrays of different shapes without copying data. It is why a * 2 works (the scalar "stretches" to a's shape) and why you can subtract a row of column means from every row of a matrix in one expression. The rules, applied to the two shapes:

1. Align from the right. Write the shapes one above the other, right-aligned. If one has fewer dimensions, pad it on the left with 1s.
2. Compare each dimension. Two sizes are compatible if they are equal, or if one of them is 1.
3. Stretch the 1s. A size-1 dimension is (virtually) repeated to match the other. The result shape is the maximum in each position. If any pair is incompatible, NumPy raises ValueError: operands could not be broadcast together.

Remember: right-align, equal-or-1, stretch the 1s. Use [:, None] or keepdims=True to turn a per-row result into a column that broadcasts back.

The everyday toolbox: @, where, clip, cumsum, sort, unique

A handful of functions show up in nearly every NumPy program:

• A @ B (or np.matmul): matrix multiplication, exactly the "rows dot columns" product from M47, handed to an optimized BLAS library. On 1-D arrays, v @ w is the dot product. Don't confuse it with A * B, which is element-wise.
• np.where(cond, x, y): element-wise "if cond then x else y". With one argument it returns the indices where cond is true.
• np.clip(a, lo, hi): limit values to a range.
• np.cumsum, np.cumprod, np.diff: running totals and differences (each undoes the other, up to the first element).
• np.sort returns a sorted copy (a.sort() sorts in place); np.argsort returns the indices that would sort, so you can reorder another array by this one.
• np.unique(a, return_counts=True): distinct values (sorted) and how often each occurs.

Remember: @ is matrix product, * is element-wise. argsort gives positions, which you then use for fancy indexing.

Random generators, seeding, and saving arrays

Random numbers the modern way

Create a Generator with rng = np.random.default_rng(seed) and call methods on it: rng.random(n) (uniform in [0, 1)), rng.integers(lo, hi, size), rng.normal(mean, sd, size), rng.choice(a, size, replace=False), rng.permutation(a), rng.shuffle(a). A seed makes the sequence reproducible: the same seed gives the same numbers, which you need for tests, debugging, and reproducible science. Pass the generator into functions instead of using global state. The older np.random.seed() / np.random.rand() API uses one hidden global generator; avoid it in new code.

Saving and loading

np.save("a.npy", a) writes one array in NumPy's binary format (dtype and shape included, exact values, fast). np.savez("data.npz", x=x, y=y) bundles several named arrays (savez_compressed zips them). np.load reads both. For human-readable text, there's np.savetxt/np.loadtxt, which is slower and can lose precision.

Remember: one default_rng(seed) per program (or test), passed around explicitly; .npy/.npz for exact, fast storage.

Performance and common pitfalls

Remember: preallocate or build lists then convert; check shapes when results look "too big"; compare floats with tolerance.

Practice

Use NumPy and no Python loops over elements unless the exercise allows it.

import numpy as np

def checkerboard(n):
    i = np.arange(n)
    return (i[:, None] + i[None, :]) % 2

def normalize_rows(X):
    X = np.asarray(X, dtype=float)
    return X / X.sum(axis=1, keepdims=True)

import numpy as np

def zscore(X):
    X = np.asarray(X, dtype=float)
    mu = X.mean(axis=0)
    sd = X.std(axis=0)
    sd = np.where(sd == 0, 1.0, sd)
    return (X - mu) / sd

import numpy as np

def one_hot(labels, n_classes):
    labels = np.asarray(labels)
    out = np.zeros((labels.size, n_classes), dtype=int)
    out[np.arange(labels.size), labels] = 1
    return out

def top_k(scores, k):
    return np.argsort(scores)[::-1][:k]

import numpy as np

def pairwise_distances(A, B):
    diff = A[:, None, :] - B[None, :, :]
    return np.sqrt((diff ** 2).sum(axis=-1))

def nearest(A, B):
    return pairwise_distances(A, B).argmin(axis=1)

import numpy as np

def moving_average(x, w):
    x = np.asarray(x, dtype=float)
    if w < 1 or w > x.size:
        raise ValueError(f"window must be in 1..{x.size}, got {w}")
    c = np.concatenate(([0.0], np.cumsum(x)))
    return (c[w:] - c[:-w]) / w

Mini-project: image manipulation with plain arrays

A grayscale image is just a 2-D array of brightness values; a color image is a 3-D array of shape (height, width, 3) holding red, green, and blue. Everything you learned applies: slicing flips and crops, broadcasting weights channels, and a blur is a convolution, where each output pixel is the average of a small neighborhood. The cell below builds a synthetic picture so it runs offline, then processes it.

pixi init imgtools && cd imgtools
pixi add python=3.14 numpy matplotlib pytest
pixi run python demo.py photo.png      # writes compare.png
pixi run pytest -q

Quiz

Recap

• An ndarray is a homogeneous, contiguous buffer plus dtype, shape, and strides; operators are element-wise.
• Create with array, zeros, ones, arange, linspace, eye, and default_rng(seed).
• Fixed-width integers wrap on overflow; choose dtypes deliberately and astype before risky math.
• axis=k is the dimension that disappears; keepdims keeps it for broadcasting.
• Slices and transposes are views; masks and fancy indexing copy. np.shares_memory tells you.
• Broadcasting: right-align shapes, sizes must be equal or 1, the 1s stretch.
• Vectorize with ufuncs and reductions; avoid Python loops over elements and np.append in loops.

Going further

• NumPy: the absolute basics for beginners
• NumPy user guide: broadcasting
• NumPy user guide: indexing on ndarrays (views vs copies)
• Random Generator reference
• Next: M49 · Linear Algebra with NumPy & SciPy