NumPy
NumPy lets you do one operation to millions of numbers at once instead of one at a time. It is the foundation every AI library in Python is built on.
- 14 min read
- 3 reading levels
- Published
Read these first
On this page 7
One lesson, three depths. Pick the one that fits you today — you can switch any time.
Beginner — No maths. Plain English.
NumPy lets you do one operation to a whole pile of numbers at once. Not one by one — all of them, together.
Think about an ice cube tray. You do not fill each little square with a spoon, one after another. You hold the tray under the tap and pour once. Every square fills at the same time, and you are done in seconds.
Plain Python fills the squares with a spoon. NumPy pours.
That difference sounds small. It is the difference between an AI model that trains overnight and one that would take three years.
Why you should care
AI is mostly arithmetic done to enormous piles of numbers. One photo on your phone is roughly twelve million brightness values. One sentence going into ChatGPT becomes thousands of values. A model looks at millions of examples.
If each value had to be handled separately by Python, nothing would finish. NumPy is what makes the pile handleable. Every library you will meet later — pandas, scikit-learn, PyTorch — stands on it.
Why it exists
Python was built to be readable, not fast. Handling one value at a time in Python costs a lot of hidden bookkeeping. Doing that a million times is painfully slow.
Meanwhile, the languages that scientists used before Python were fast but unfriendly. Nobody wanted to give up readable code, and nobody could accept the slowness.
NumPy is the peace treaty. You keep writing friendly Python. Behind the scenes, the heavy work goes to fast compiled code. That code chews through the whole pile in one go.
How it works
A normal Python list keeps each value in its own little parcel, scattered around memory. NumPy keeps them shoulder to shoulder in one long tidy row.
Python list → value ... elsewhere ... value ... elsewhere ... value
(scattered, each one wrapped in its own parcel)
NumPy array → value value value value value value value value
(one unbroken row, all the same kind of thing)Because the values are neat and identical, the computer can work on many of them in a single sweep:
Plain Python: take one value → change it → store it → repeat a million times
NumPy: say what you want once → the whole row changes togetherThat "say it once" style has a name: vectorised. It means the operation applies to the whole collection, not to one element at a time.
A real example you have seen
When you slide the brightness control on a photo in your phone's gallery, the picture changes instantly. Behind that slider, millions of brightness values are being raised at the same moment.
Nothing walks through them one at a time. The whole grid is lifted in one sweep. It is the same as tilting a full tray, rather than moving each cube.
An honest word
The idea of "do it to everything at once" takes a while to sink in. Most people write a loop first, then learn to remove it. That is a normal path, not a failure.
If you find yourself writing a loop over a NumPy array, pause and ask whether one line could replace it. The answer is often yes, and the answer being no is also fine.
Remember this
- NumPy stores many values of one type in a tidy block, so the computer sweeps through them fast.
- You describe the operation once, and it applies to everything.
- Every serious AI library in Python is built on top of this idea.
What to learn next
- Pandas — tables with names and labels, built on NumPy.
- Matplotlib — turning those numbers into pictures.
- Vectors and matrices — what those rows and grids mean.
Developer — Code and libraries.
Setup
pip install numpyEverything below runs on a laptop CPU in under a second. No GPU, no download.
Your first array
import numpy as np
# Five temperature readings, one per minute, in degrees Celsius.
temps = np.array([31.5, 32.0, 34.0, 33.5, 30.0])
print("shape:", temps.shape, " dtype:", temps.dtype)
print("mean :", temps.mean(), " max:", temps.max())
# Vectorised: one line touches all five values, with no Python loop.
fahrenheit = temps * 9 / 5 + 32
print("degF :", fahrenheit)
# A mask is an array of True/False, one entry per element.
hot = temps > 32
print("mask :", hot)
print("hot :", temps[hot])
print("count:", hot.sum())shape: (5,) dtype: float64 mean : 32.2 max: 34.0 degF : [88.7 89.6 93.2 92.3 86. ] mask : [False False True True False] hot : [34. 33.5] count: 2
Two dimensions, axes, and broadcasting
import numpy as np
# Three students, four subjects. Rows = students, columns = subjects.
marks = np.array([
[88.0, 71.0, 95.0, 60.0],
[79.0, 92.0, 68.0, 84.0],
[55.0, 64.0, 72.0, 91.0],
])
print("shape:", marks.shape)
print("per-student mean:", marks.mean(axis=1))
print("per-subject mean:", marks.mean(axis=0))
# Broadcasting: one number per subject is stretched across all three rows.
bonus = np.array([2.0, 0.0, 5.0, 0.0])
print("after bonus:\n", marks + bonus)
print("reshaped:\n", marks.reshape(2, 6))
print("transposed shape:", marks.T.shape)shape: (3, 4) per-student mean: [78.5 80.75 70.5 ] per-subject mean: [74. 75.66666667 78.33333333 78.33333333] after bonus: [[ 90. 71. 100. 60.] [ 81. 92. 73. 84.] [ 57. 64. 77. 91.]] reshaped: [[88. 71. 95. 60. 79. 92.] [68. 84. 55. 64. 72. 91.]] transposed shape: (4, 3)
Line-by-line walkthrough
shape is a tuple describing the size along each dimension. (3, 4) means three rows and four columns. Almost every NumPy error you will hit is a shape you did not expect, so print shapes often.
dtype is the single type every element shares — here float64, a 64-bit decimal number. A NumPy array cannot hold a mix of text and numbers the way a list can. That restriction is what buys the speed.
axis answers "which direction do I collapse?" axis=0 walks down the rows and gives you one value per column. axis=1 walks across the columns and gives you one value per row. This is confusing for almost everyone the first week. A trick that helps: axis=n is the dimension that disappears from the shape.
Broadcasting is NumPy stretching a smaller array to match a larger one without copying it. marks is (3, 4) and bonus is (4,). NumPy lines them up from the right, sees the four columns match, and reuses bonus for every row.
reshape(2, 6) rearranges the same twelve values into a different grid. The total count must match, and the values are read row by row.
marks.T is the transpose, rows and columns swapped. It does not copy anything; it is a different view of the same memory.
Speed, measured
import time
import numpy as np
n = 5_000_000
py_list = list(range(n))
np_arr = np.arange(n, dtype=np.float64)
t0 = time.perf_counter()
squared_list = [x * x for x in py_list]
t1 = time.perf_counter()
t2 = time.perf_counter()
squared_arr = np_arr * np_arr
t3 = time.perf_counter()
print(f"python list : {t1 - t0:.3f} s")
print(f"numpy array : {t3 - t2:.3f} s")
print(f"same answer : {squared_list[-1] == squared_arr[-1]}")
print("list memory :", (n * 8 + n * 28) // 1_000_000, "MB (approx)")
print("array memory:", np_arr.nbytes // 1_000_000, "MB")python list : 0.161 s numpy array : 0.008 s same answer : True list memory : 180 MB (approx) array memory: 40 MB
Those two timings come from one ordinary laptop. Your numbers will differ. The ratio depends on your CPU, your NumPy build, and what else is running. The shape of the result holds everywhere. The array version is roughly ten to fifty times faster, on a fraction of the memory.
The memory line is worth staring at. A Python list of five million floats stores five million pointers plus five million separate float objects. The array stores five million raw 8-byte values and nothing else.
Common mistakes
1. Thinking a slice is a copy.
import numpy as np
a = np.array([1, 2, 3, 4])
b = a[:2]
b[0] = 99
print(a)[99 2 3 4]
Writing to b changed a. A basic slice is a view, a window onto the same memory rather than a new array. Use a[:2].copy() when you want an independent array. Boolean and fancy indexing do return copies, which makes the rule feel inconsistent. It is inconsistent; memorise the two cases.
2. Integer arrays silently truncating.
import numpy as np
a = np.array([1, 2, 3]) # no decimal point anywhere, so this is an integer array
a[0] = 2.7
print(a)[2 2 3]
The .7 was thrown away without a warning. Fix: create with np.array([1, 2, 3], dtype=float) when decimals are coming.
3. Broadcasting the wrong way round.
A (3,) array added to a (3, 4) array fails. Alignment happens from the right, and 3 does not match 4.
ValueError: operands could not be broadcast together with shapes (3,4) (3,)
Fix: reshape to a column with bonus.reshape(3, 1) or bonus[:, None].
4. Comparing floats with ==.
0.1 + 0.2 == 0.3 is False in every language that uses binary floating point. Use np.isclose(a, b) or np.allclose(a, b) for arrays.
Try it yourself
Take the marks array and find each student's best subject without writing a loop. np.argmax(marks, axis=1) gives you the column position; use it to index into a list of subject names.
Then normalise every column so it has mean zero. Subtract marks.mean(axis=0) and print the new column means. They will be tiny values near zero rather than exact zeros, which is floating-point rounding doing its honest work.
What to learn next
- Pandas — labelled tables on top of these arrays.
- Vectors and matrices — the meaning behind shapes and dot products.
- Matplotlib — plotting arrays.
Researcher — Mathematics and papers.
The object
An ndarray is a thin Python wrapper over four things:
- a pointer to a single contiguous buffer of raw bytes,
- a
dtypegiving the size and interpretation of one element, - a
shapetuple(d0, d1, ..., d[k-1]), - a
stridestuple(s0, s1, ..., s[k-1])measured in bytes.
The element at index (i0, i1, ..., i[k-1]) lives at byte offset:
offset = i0*s0 + i1*s1 + ... + i[k-1]*s[k-1]kis the number of dimensions.i[j]is the index along dimensionj.s[j]is the stride: how many bytes you step to advance one position along dimensionj.d[j]is the extent: how many positions exist along dimensionj.
This single formula explains a great deal:
- Transpose is free.
a.Tpermutesshapeandstridesand reuses the buffer. Cost is O(1), not O(n). - Basic slicing is free. A slice adjusts the base offset, shape, and strides. It is a view.
reshapeis usually free. It must copy when the requested layout cannot be expressed with any stride pattern over the existing buffer. That happens after certain transposes.np.ascontiguousarrayforces the copy explicitly.- Fancy and boolean indexing must copy, because an arbitrary index list has no constant stride.
C order (row-major) means the last axis has the smallest stride. Fortran order means the first does. Iterating along the small-stride axis is cache-friendly. Iterating along the large-stride axis can be an order of magnitude slower, for identical arithmetic.
Broadcasting, stated precisely
Given operands with shapes A and B, align them from the trailing axis and prepend 1s to the shorter shape. The operation is defined when this holds on every axis j:
a[j] == b[j] or a[j] == 1 or b[j] == 1
output extent[j] = max(a[j], b[j])Here a[j] and b[j] are the extents of the two operands on axis j after alignment.
An axis of extent 1 is realised by setting its stride to 0. The same memory is then re-read on every step, and no data is duplicated. So broadcasting is cheap in memory but not free in bandwidth. A stride-0 axis still costs cache traffic.
Numerical behaviour
Machine epsilon, written eps, is the gap between 1.0 and the next representable number above it. It sets the relative accuracy of every operation in that format.
| dtype | significand bits | eps (approx) | typical use |
|---|---|---|---|
float64 | 53 | 2.2e-16 | scientific computing, NumPy default |
float32 | 24 | 1.2e-07 | neural network training and inference |
float16 | 11 | 9.8e-04 | mixed precision, needs loss scaling |
bfloat16 | 8 | 7.8e-03 | training; same exponent range as float32 |
Naive left-to-right summation of n values has a worst-case relative error bound that grows like n * eps. NumPy's reduction uses pairwise summation, which recursively splits the range in half, giving a bound of order eps * log2(n). Kahan compensated summation reduces the bound to order eps at roughly four times the arithmetic cost. For plain Python floats, math.fsum returns the exactly rounded sum.
Practical consequence: for large n, arr.sum() and sum(arr) can differ in the last few digits. The NumPy one is the more accurate.
Where the FLOPs go
Element-wise operations such as a * b are memory-bandwidth bound. The limit is how fast data reaches the CPU, not how fast it multiplies. For float64, each multiply reads 16 bytes and writes 8 bytes, for one floating-point operation. That is an arithmetic intensity near 0.04 FLOP per byte. A modern CPU can sustain tens of FLOP per byte, so the arithmetic units idle while memory catches up.
Consequence: chaining a * b + c as two statements materialises an intermediate array and pays that bandwidth twice. np.multiply(a, b, out=tmp) and the numexpr package exist to avoid the intermediate.
Matrix multiplication is the opposite case. Two n x n operands cost 2*n**3 - n**2 floating-point operations while touching only 3*n**2 elements. Intensity therefore grows like O(n), and np.dot reaches a large fraction of peak throughput.
NumPy dispatches it to a tuned BLAS (Basic Linear Algebra Subprograms) library, such as OpenBLAS or Intel MKL. Those libraries block the problem for cache and use SIMD registers. They also release the GIL and thread internally, so a single A @ B already uses every core.
Asymptotically faster algorithms exist, such as Strassen at O(n**2.807). Their large constants and weaker numerical stability keep them out of default BLAS kernels at practical sizes.
Universal functions
A ufunc is a compiled element-wise kernel with a broadcasting front end and a defined type-resolution table. np.add is a ufunc; so is np.exp. Useful properties:
out=writes in place and avoids an allocation.where=applies a boolean mask at the kernel level..reduce,.accumulate, and.reduceatgive folds, scans, and segmented folds without Python loops.np.errstatecontrols how overflow, underflow, division by zero, and invalid results are reported.
Custom kernels that cannot be expressed as ufunc compositions belong in Numba (@njit) or Cython, not in a Python loop.
References
- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. "Array programming with NumPy." Nature 585, 357–362, 2020.
- van der Walt, S., Colbert, S. C., Varoquaux, G. "The NumPy Array: A Structure for Efficient Numerical Computation." Computing in Science & Engineering 13(2), 22–30, 2011.
- Higham, N. J. Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002. Chapter 4 covers summation error bounds.
- Goldberg, D. "What Every Computer Scientist Should Know About Floating-Point Arithmetic." ACM Computing Surveys 23(1), 5–48, 1991.
- Williams, S., Waterman, A., Patterson, D. "Roofline: An Insightful Visual Performance Model for Multicore Architectures." CACM 52(4), 65–76, 2009. The arithmetic-intensity framing used above.
What to learn next
- Vectors and matrices — the algebra these shapes implement.
- Linear algebra — decompositions, rank, and conditioning.
- Pandas — labelled, heterogeneous data on the same foundation.