NumPy

What it is and why it exists

NumPy is Python’s numerical array library. Its core object is the N-dimensional array (ndarray) : homogeneous elements plus metadata that describes dimensions, data type, and memory layout. A two-dimensional array might mean “samples × features,” while a three-dimensional array might mean “image height × width × channels,” but NumPy doesn’t record those domain meanings for you.

Python lists are good at holding heterogeneous objects and changing collections, but they don’t have one set of array arithmetic rules. [1, 2] * 2 repeats a list, whereas multiplying a NumPy array by 2 computes element by element. NumPy gives regular numeric data a shared array protocol, so Pandas, SciPy, Matplotlib, and many machine learning tools can exchange data without a Python loop for every operation.

NumPy fits regular, dense, homogeneous numeric grids. Tables with column names and mixed types usually fit Pandas or Polars better; highly irregular object collections may still belong in Python containers. That boundary matters: don’t force irregular data into a dtype=object array merely to call the result “vectorized.”

The examples here were run with Python 3.14.3 and NumPy 2.5.2. Their displayed values use explicit rounding or list conversion, so they don’t depend on array truncation rules or terminal width.

How it works

Three contracts in one array

Every array has an array shape , a data type (dtype) , and a data buffer. shape is a tuple containing the length of each axis; ndim is the number of axes; size is the total element count. The dtype tells NumPy how to interpret each element, such as int64, float32, or bool.

An array variable also has an implicit domain contract: what each axis means. Two arrays with shape (32, 10) could represent 32 samples with 10 features or 32 time points from 10 sensors. NumPy checks whether lengths can participate in an operation, not whether the business meanings of the axes have been reversed.

An array stride is the number of bytes crossed when moving one position along an axis. Slicing and transposition can change only shape and strides while continuing to reference the same buffer. Two arrays that look different can therefore share data, which is why modifying a slice sometimes changes its source.

Creation fixes the initial type

np.array() builds from existing Python data; np.zeros(), np.ones(), and np.full() create initialized regular arrays; np.arange() and np.linspace() create numeric sequences. Specify dtype at business boundaries, especially when integer width and floating precision matter or when input may be empty. This prevents one batch of sample values from deciding the type by accident.

np.empty() allocates space without initializing its values. It suits internal code that will overwrite every element, not an “empty values” array. If any position could be read before assignment, use a creation function with a defined initial value.

reshape() requires the element count to stay unchanged. -1 asks NumPy to infer one axis length, but a shape-valid statement can still express the wrong data layout. Name the meaning of each axis before and after a reshape instead of checking only whether the element count divides evenly.

Common shape operations differ in useful ways:

Moving and joining axes

Transposition doesn’t mean “swap rows and columns for every array.” .T swaps the two axes of a 2D array, but it reverses the complete axis order of a 3D array. For batch, channel, and spatial dimensions, transpose() or moveaxis() states the intended order more clearly.

stack() and concatenate() solve different problems. stack() creates a new axis, so ten images shaped (64, 64) become (10, 64, 64); concatenate() extends an existing axis. Generated code often adds or loses one batch axis at this boundary.

Calling squeeze() without an axis removes every length-one axis. That seems convenient in sample data, but it can make a return value lose a dimension whenever the batch size happens to be one. Public APIs should name the axis to remove or keep a stable batch dimension.

np.newaxis is an alias for None that inserts a length-one axis in an index. values[:, np.newaxis] turns (n,) into (n, 1), while values[np.newaxis, :] produces (1, n). The elements are identical, but the broadcasting direction is completely different.

Before joining arrays, validate lengths outside the target axis. An error can report that dimensions conflict, but it can’t know that one input is channel-first while another is channel-last. Axis names in the interface matter more than repairing one ValueError.

Vectorization, universal functions, and reductions

Vectorization expresses a computation as operations on complete arrays. For example, prices * quantities means element-wise multiplication, while np.sqrt(values) applies one universal function to every element. Vectorization first changes the expression of the problem: NumPy handles the loop, and shape and type rules live in the array operation.

Reductions such as sum(), mean(), and max() use axis to name the axis they remove. For a “samples × features” matrix, axis=0 removes the sample axis and leaves one result per feature; axis=1 removes the feature axis and leaves one result per sample. keepdims=True retains a length-one axis and often makes later broadcasting easier to read.

Ordinary arithmetic is element-wise. * between two matrices is not matrix multiplication; use @ or np.matmul() for linear algebra multiplication. Use np.linalg.solve() for a square linear system rather than explicitly forming an inverse and multiplying by it.

Broadcasting starts at the trailing axes

Broadcasting lets arrays with different shapes participate in element-wise operations. NumPy compares axes from the right: a pair is compatible when the lengths are equal or one length is 1; a missing leading axis behaves as length 1. Otherwise, the operation raises a shape-compatibility error.

For example, a sample matrix shaped (4, 3) can subtract feature means shaped (3,) because both trailing axes have length 3. Per-sample means shaped (4,) can’t be subtracted directly because trailing lengths 3 and 4 conflict. Retaining those means as (4, 1) gives them a length-one trailing axis that broadcasts across features.

Broadcasting describes a conceptual expansion; it doesn’t mean NumPy first creates repeated input data. The arithmetic result will usually still allocate its own array. Predict the result shape before writing the expression. “It didn’t raise” proves only that the shapes are compatible, not that the axes mean the right things.

Indexing decides sharing or copying

Basic slicing, such as array[1:4] or matrix[:, ::2], usually returns an array view that shares data. Element changes through that view are visible in the source. view.base helps exploration, but the base of a chain of views isn’t guaranteed to be the array in your hand, so use np.shares_memory() to test whether two arrays overlap.

Integer-array and Boolean indexing are advanced indexing and return copies. You can modify array[[0, 2]] and array[array > 0] independently afterward, but the copy allocates a new buffer. When basic and advanced indexing are mixed, the result shape can also depart from a simple axis-by-axis intuition, so confirm it with a small input.

Assignment has a separate rule: array[mask] = 0 writes into the source because the index is the assignment target. In contrast, assigning selected = array[mask] and then changing selected changes only that copy. A review must follow the complete data flow; one index expression isn’t enough to tell where the final write lands.

Masks, conditions, and order

Comparing arrays produces Boolean arrays. Combine conditions with the element-wise operators &, |, and ~, usually with parentheses around each comparison. Python’s and and or require one truth value for the complete array and therefore raise an ambiguity error.

Boolean indexing retains only matching elements and often flattens the original structure. np.where(condition, left, right) instead chooses element by element between two equally shaped or broadcastable values and keeps the broadcasted shape. The right choice depends on whether downstream code still needs the original axes.

Missing data has more than one representation. Floating arrays can store NaN directly, while integer arrays cannot; creating an array from integers and NaN will generally promote it to a floating type. If missingness is mixed with zero, an empty string, or a sentinel, normalize that meaning before a reduction.

np.argsort() returns indexes that can reorder an array, not the sorted values themselves. To reorder labels and features together, compute one order and apply it to every associated array. Sorting only one column breaks the correspondence within each sample.

Mask and order-index lengths are part of the array contract too. A stale mask produced before one filter can’t safely be reused after that filter. A pipeline should carry a mask alongside its matching data or recompute it at each stage.

Errors and numerical warnings

Incompatible shapes usually raise ValueError. Read the input shapes in the exception and compare them from the trailing axes before inserting a blind reshape(). A shape patch that makes code run doesn’t necessarily preserve the data’s meaning.

Unsafe in-place casts fail. For example, a floating division result can’t be written back into an integer array under the same_kind casting rule. The usual fix is to create a result with the required floating type, not force a conversion that truncates information.

Floating divide-by-zero, invalid operations, and overflow may emit RuntimeWarning by default while producing inf or NaN. Those values can continue into later calculations if the caller ignores warnings, so “the function didn’t raise” isn’t a sufficient check.

np.errstate() changes floating-error handling inside one explicit block. Tests can temporarily promote relevant warnings to exceptions on critical numerical paths. An algorithm that expects special values can handle them locally and validate the result immediately.

np.seterr() changes NumPy’s error settings for the current thread and shouldn’t be changed casually inside an ordinary library function. A local context is easier to reason about than hidden environmental state and doesn’t make the caller restore settings afterward.

np.isfinite() finds NaN and both infinities together. Input validation and important intermediate results can use it, but the domain contract still decides whether those values are rejected, replaced, or preserved.

Warning policy doesn’t replace boundary tests. Cover division by zero, very large and small values, and all-missing slices separately because their promotion, warnings, and result shapes can differ.

Examples

Read shape, type, and axes

Rows in this matrix are days and columns are sensors. axis=1 removes the sensor axis to produce daily means; axis=0 removes the day axis to produce a peak for each sensor.

The shape (2, 3) doesn’t say which axis is the day. Variable names, nearby comments, and interface documentation must preserve that meaning. After a reduction, check that the result shape matches what its consumer expects.

Standardize features with broadcasting

This example treats each column as one feature. keepdims=True keeps the mean and standard deviation shaped (1, 3), so subtraction and division clearly broadcast across samples. The third column is constant, and the code replaces its zero standard deviation with 1.0 to avoid division by zero.

import numpy as np

samples = np.array(
    [[1.0, 10.0, 7.0], [2.0, 20.0, 7.0], [3.0, 30.0, 7.0]]
)

means = samples.mean(axis=0, keepdims=True)
scales = samples.std(axis=0, keepdims=True)
safe_scales = np.where(scales == 0, 1.0, scales)
standardized = (samples - means) / safe_scales

np.set_printoptions(precision=2, suppress=True)
print("mean shape:", means.shape)
print(standardized)
print("constant feature:", standardized[:, 2].tolist())

mean shape: (1, 3)
[[-1.22 -1.22  0.  ]
 [ 0.    0.    0.  ]
 [ 1.22  1.22  0.  ]]
constant feature: [0.0, 0.0, 0.0]

Replacing a zero scale with 1.0 is an explicit policy in this example, not a universal rule. Another system might drop constant features or reject the input. Numerical code needs to put that domain choice in its interface and tests.

Distinguish a view from a copy

The basic slice scores[:, 1:] shares data with the source, while Boolean indexing creates a copy. The assignments look similar, but only the first one writes back into scores.

import numpy as np

scores = np.array(
    [[70, 80, 90], [60, 75, 85], [88, 92, 95]],
    dtype=np.int64,
)

last_two = scores[:, 1:]
selected = scores[scores[:, 0] >= 70]

print("slice shares memory:", np.shares_memory(scores, last_two))
print("mask shares memory:", np.shares_memory(scores, selected))

last_two[0, 0] = 999
selected[0, 0] = -1
print("source first row:", scores[0].tolist())
print("selected first row:", selected[0].tolist())

slice shares memory: True
mask shares memory: False
source first row: [70, 999, 90]
selected first row: [-1, 80, 90]

Call .copy() when a modification must be isolated in a new buffer. When a function updates in place, say so in its name or documentation. Making callers guess about ownership lets an apparently local cleaning step corrupt upstream data.

Isolate random state and solve least squares

default_rng() returns an independent random generator rather than resetting module-level random state. This example uses it to produce training and test indexes, then fits a one-feature linear model with an intercept through np.linalg.lstsq().

import numpy as np

rng = np.random.default_rng(42)
feature = np.array([0, 1, 2, 3, 4, 5], dtype=np.float64)
target = 1.0 + 2.0 * feature

order = rng.permutation(feature.size)
train_index = order[:4]
test_index = order[4:]

design = np.column_stack(
    [np.ones(train_index.size), feature[train_index]]
)
intercept, slope = np.linalg.lstsq(
    design, target[train_index], rcond=None
)[0]
predicted = intercept + slope * feature[test_index]

print("train index:", train_index.tolist())
print(f"model: y = {intercept:.1f} + {slope:.1f}x")
print("test index:", test_index.tolist())
print("prediction:", np.round(predicted, 6).tolist())

train index: [3, 2, 5, 4]
model: y = 1.0 + 2.0x
test index: [1, 0]
prediction: [3.0, 1.0]

A fixed seed makes this run repeatable, but a long-lived experiment should also save the actual indexes, NumPy version, and selected bit generator. Don’t treat “same seed” as complete data provenance across every version and algorithm. Real model splits must also define constraints such as stratification, grouping, and time order before shuffling.

Pitfalls

Fix: Assert shapes at boundaries and document what every axis means. Use keepdims=True, np.newaxis, or an explicit reshape() to name the expanded axis, then check the result with a tiny hand-calculated input.

Fix: Call .copy() when isolation is required and np.shares_memory() when overlap must be confirmed. Tests for in-place functions should check both the returned value and whether the input changed according to contract.

Fix: Choose a wide enough dtype at the input boundary and use np.result_type() to inspect the result type of a mixed operation. Check types before division, normalization, and cumulative sums, and test maxima, minima, and negative values.

Fix: Use empty() only when you can prove every element is assigned before it is read. Otherwise, use zeros(), full(), or an explicit initial value that communicates the missing-data policy.

Fix: Accept a Generator parameter or create a local Generator with an explicit seed inside the function. When an experiment must replay across processes, save the samples and bit-generator information too.

Fix: Use np.isclose() or np.allclose() with tolerances from the domain. Decide whether NaN should be rejected, ignored, or propagated before choosing mean(), nanmean(), or explicit validation; don’t ignore missing values merely to obtain a number.

Further reading

• NumPy: the absolute basics for beginners
• NumPy: broadcasting
• NumPy: copies and views
• NumPy: random generator
• NumPy: linear algebra