The point of putting data into a NumPy array is to do arithmetic on it without writing a for loop. Element-wise math, aggregations along axes, broadcasting two differently-shaped arrays together, and linear-algebra calls into compiled BLAS routines. This is the core of NumPy. A loop that took 40 milliseconds in pure Python now takes less than a millisecond, and the code that does it is shorter and clearer. This lesson covers the operations that show up in most NumPy code, plus the broadcasting rules required to avoid shape-mismatch surprises.

Arithmetic on a NumPy array applies to every element, in parallel. There's no loop in the code, and no loop in NumPy's Python layer either; the work happens in compiled C.

Each operation produced a new array of the same shape. The original prices is untouched. prices * 2 is shorthand for "multiply every element by 2", and prices * 0.9 applied a 10% discount to every price. The operator on the right is a scalar (a single number), and NumPy broadcast it across the array.

Two arrays of the same shape combine element by element:

prices * quantities multiplied corresponding elements: index 0 with index 0, index 1 with index 1, and so on. The result is a new array of per-row totals. Calling .sum() then collapses it to a single scalar. The whole thing is three lines, one mental step per line, and runs at C speed for a million items.

The four basic arithmetic operators (+, -, *, /) all work this way, as do the rest: ** for exponentiation, // for floor division, % for modulo. Comparison operators (==, !=, <, >, <=, >=) produce a boolean array of the same shape, which can be used as a mask.

The boolean array can drive the next operation: filter, replace, count, or anything else. This composability ("compute a mask, use the mask, compute the next thing") is the rhythm of NumPy code.

Every arithmetic operation on an array allocates a new array for the result. For long chains of operations on huge data, this can be the dominant cost. The in-place forms (arr *= 2, arr += 5) skip the allocation when the original isn't needed.

The augmented operators (+=, -=, *=, /=) modify the array in place instead of allocating a new one. Same semantics, less memory pressure.

prices *= 0.9 rewrote every cell of prices in place. No new array was allocated. For a 1-billion-element array, the difference between in-place and not is the difference between using 8 GB and 16 GB of memory.

In-place arithmetic also enforces dtype compatibility more strictly. prices *= 0.9 works because prices is float64, which can hold the result. If prices were int64, the same expression would raise TypeError because the result needs floats. Out-of-place (prices = prices * 0.9) changes the result's dtype to float64 automatically, which is sometimes desired and sometimes a memory surprise.

NumPy ships with a large library of element-wise functions called universal functions, or ufuncs. np.sin, np.exp, np.log, np.sqrt, np.abs, and many more. They take an array, apply the function to every element, and return a new array of the same shape.

Ufuncs apply a math function to every element, and they are also the foundation for more elaborate operations. The full list is long; the subset that comes up most often:

np.clip handles "cap values between a min and a max" without a loop:

Two-argument ufuncs combine two arrays element-wise. np.maximum(a, b) returns the elementwise maximum, np.minimum the minimum, np.power(a, b) raises a to b. Same shape rules as ordinary arithmetic.

This replaces a Python for loop with min() on each pair, and it runs about 50x faster for arrays large enough to matter.

Aggregations collapse an array down to a smaller one (often a scalar) by combining values. The common ones live as both methods on the array and functions in the top-level np namespace.

sum, mean, std (standard deviation), min, and max are all "reductions" that collapse the whole array into a single number. argmin and argmax return the index of the smallest and largest values, which is often more useful than the values themselves. Knowing that the cheapest product is at index 2 enables a lookup in a parallel array of names.

cumsum is interesting because it doesn't reduce to a scalar. It returns an array of the same length, where each element is the running total up to that point:

After day 0, the running total is 120. After day 1, it's 215. After day 6, it's 900. This is the building block for running revenue dashboards, time-series visualizations, and any report that needs "where am I now" rather than "what did I do today".

The aggregations become more powerful when applied to a multi-dimensional array along a specific axis. Without an axis, they reduce the whole array to a scalar. With axis=0, they collapse the rows, leaving one value per column. With axis=1, they collapse the columns, leaving one value per row.

The rule for axes: axis=k is the axis being collapsed. The remaining axes form the shape of the result. A (7, 3) array collapsed along axis=0 (the rows) leaves a length-3 result (one value per column). Collapsed along axis=1 (the columns) leaves a length-7 result (one value per row).

The same rule scales to higher dimensions. A (4, 7, 3) array (stores by days by categories) summed along axis=0 gives a (7, 3) result (one value per day-and-category, summed across stores). Summed along axis=1 it gives a (4, 3) result. Summed along axes (0, 1) it gives a (3,) result. The shape of the result is whatever was left after removing the collapsed axes.

argmax and argmin follow the same pattern: with no axis, they return a single index into the flattened array. With an axis, they return an array of indices, one per surviving axis position.

The best day for every category is day 5 (Saturday, index 5). The best category every day is category 0 (electronics). The first result has length 3 (one per category), the second has length 7 (one per day). Same rule: the collapsed axis disappears from the result.

Broadcasting is the rule NumPy uses when an operation applies to two arrays of different shapes. Instead of failing or copying data, NumPy stretches the smaller array along the missing dimensions to match the larger one, virtually. The actual memory isn't duplicated; the smaller array is reused at every position along the stretched axis.

The simplest case is scalar broadcasting. prices * 2 doesn't materialize a length-N array of [2, 2, 2, ...]. NumPy multiplies each price by the scalar 2. That's broadcasting in its most basic form.

Broadcasting gets interesting with multi-dimensional arrays. The rule: two arrays are broadcast-compatible if, when their shapes are lined up from the right, every pair of corresponding dimensions is either equal or one of them is 1.

A concrete example. Subtract the per-category average from every cell in a (7, 3) revenue grid:

revenue has shape (7, 3). category_avg has shape (3,). NumPy lines them up from the right: 3 matches 3. There's no second dimension on category_avg, so NumPy treats it as if it were (1, 3), and broadcasts it down all 7 rows. The result is a (7, 3) array where each row has had the category averages subtracted. No loops, no for over the seven days.

The same trick works for "subtract the per-day average from every column". The average needs an extra dimension first, because the shape mismatch must align correctly:

The (7, 1) shape lines up with (7, 3) because both have 7 on the left and the column on the right is 1, which broadcasts. Without the reshape, aligning (7,) against (7, 3) puts the 7 in the wrong slot (the column slot) and fails with a ValueError.

The rule visualized:

The whole rule: line shapes up from the right, walk each pair, accept if they're equal or one is 1, fail otherwise. Missing dimensions on the left of a shape are treated as 1. A (3,) array is equivalent to (1, 3) for broadcasting purposes, which is why it lined up with (7, 3) on the right.

When broadcasting fails, the error is explicit:

b has shape (2,), which aligns to the right of a's shape: (2, 3) and (2,) become (2, 3) and (_, 2). The 3 and the 2 don't match, neither is 1, so broadcasting fails. To make it work, b would need shape (3,) (broadcast across rows) or (2, 1) (broadcast across columns).

The main reason to learn NumPy is the speedup over loops. Concrete numbers help. Compute the sum of squares of a million numbers:

The exact ratio depends on the CPU and Python version, but 30x-100x is typical for arithmetic-heavy work. Two factors drive the speedup. First, the inner loop runs in compiled C instead of the Python interpreter, which strips out per-iteration overhead (no bytecode dispatch, no reference counting, no boxed object lookup). Second, the data is laid out contiguously, so the CPU's cache and prefetcher work efficiently.

A NumPy loop of the form for x in arr: do_something(x) is usually no faster than a list loop and sometimes slower, because of the overhead of extracting each element from the array. Look for a vectorized equivalent before writing such a loop.

Not every operation vectorizes neatly. When per-element logic needs branches, exception handling, or external state, a loop is the right answer. The goal isn't to eliminate loops everywhere, it's to use vectorized ops when the work is "do the same arithmetic to every element".

np.sort returns a sorted copy of an array. arr.sort() sorts in place and returns None (same pattern as Python's sorted and list.sort).

For a 2D array, np.sort sorts along the last axis by default. Passing axis=0 sorts each column independently, and axis=None flattens then sorts.

The first call sorts each row from low to high. The second sorts each column from low to high. Both keep the array's shape.

np.argsort returns the indices that would sort the array, which is useful for "sort a parallel array by this one":

argsort returned [2, 1, 4, 0, 3], meaning "to sort prices, take index 2 first, then 1, then 4, then 0, then 3". Applying that same index to names gives the names sorted by price. This is a common NumPy pattern: sort one array, carry parallel arrays along with it.

NumPy has a full linear-algebra module under np.linalg, plus a handful of operations available at the top level. The two most common operations are matrix multiplication and solving linear systems.

@ is the matrix multiplication operator. It works on 2D arrays as matrix multiplication, on 1D arrays as a dot product, and on mixed cases for matrix-vector products.

quantities @ prices did the matrix-vector product: for each row of quantities, it dotted with prices, producing one total per customer. Compared with (quantities * prices).sum(axis=1), both work, but @ is more idiomatic for matrix-style operations and dispatches to optimized BLAS routines for large matrices.

np.dot is the older form; it does the same job and is still common in legacy code:

Both return 1*4 + 2*5 + 3*6 = 32. The @ operator was added in Python 3.5 and is the preferred form today.

For inverting matrices or solving linear systems, np.linalg has the heavy machinery. The classic case: given A @ x = b, with A and b known, solve for x.

So x = 2.2, y = 1.2. Check: 2*2.2 + 3*1.2 = 4.4 + 3.6 = 8. Correct.

This could also be written as np.linalg.inv(A) @ b, which would explicitly compute the inverse of A and multiply by b. Mathematically equivalent, numerically worse. solve is more accurate and faster because it uses a specialized factorization rather than explicitly inverting.

np.linalg.inv is O(n^3) for an n-by-n matrix and is numerically unstable for ill-conditioned matrices. Use np.linalg.solve(A, b) instead of np.linalg.inv(A) @ b whenever possible. The two compute the same result mathematically; solve computes it better.

Other useful linear algebra routines: np.linalg.det(A) for the determinant, np.linalg.eig(A) for eigenvalues and eigenvectors, np.linalg.norm(v) for vector norms, np.linalg.svd(A) for singular value decomposition. Most show up rarely in everyday work, but solve and @ are common.

The four operations come together in most NumPy tasks. Load the data, reshape if needed, do the math, summarize.

(The exact numbers depend on the seed and platform, but the structure is the same.)

Every step is a one-liner. sum(axis=0) collapses days, leaving per-store totals. mean(axis=1) collapses stores, leaving per-day averages. argmax() on the whole array gives a flat index; np.unravel_index converts it back into row-and-column coordinates. The final share = per_store / per_store.sum() is a one-line normalization. None of this needs a for loop. None of it builds an intermediate list. And it runs at C speed on data 1000x larger than the example.

This is the rhythm of practical NumPy code. Thinking in shapes and axes rather than loops, the same patterns repeat: filter with a mask, transform with arithmetic, aggregate along an axis, sort or argsort to carry along parallel data. The exact set of operations is small; the combinations cover most of what data work requires.