NumPy Interview Questions (2026): Most Asked, With Answers

An ndarray is a fixed-size, homogeneous, multidimensional array backed by a single contiguous block of memory, with metadata like shape, dtype, and strides describing how to read that memory (NumPy, ndarray reference). A Python list stores pointers to separately allocated Python objects scattered across memory, each carrying its own type info and reference count.

That difference is why NumPy is fast. Iterating a Python list means dereferencing a pointer and unboxing a PyObject on every single access. Operating on an ndarray means walking one contiguous buffer of raw bytes with a known stride, which is exactly the access pattern CPUs and vectorized instructions are built to handle well. The speedup isn't about NumPy being "written in C" so much as it's about the data layout letting the loop skip Python's per-object overhead entirely.

dtype describes the exact binary format of every element in the array, byte order, size, and whether it's an integer, float, boolean, or something else, and every element in an ndarray shares the same one. That's what "homogeneous" means in practice.

import numpy as np

a = np.array([1, 2, 3.5])
print(a.dtype) # float64, every element upcast to float

NumPy silently upcasts the whole array to the type that can represent every value without losing information, here float64. It won't split the difference or keep a mixed-type array unless you explicitly ask for dtype=object, which throws away most of the performance benefit of using NumPy at all.

arange takes a step size and generates values up to, but not including, a stop value, the same mental model as Python's range but supporting floats. linspace takes a count of points and includes both endpoints by default, dividing the interval evenly.

np.arange(0, 1, 0.1)    # step-based, may not land exactly on 1.0 due to float error
np.linspace(0, 1, 10)    # count-based, always exactly 10 evenly spaced points

I reach for linspace whenever the count of points matters more than the step size, plotting axes, sampling a fixed number of test points. arange is right when the step itself is the meaningful quantity, and floating point accumulation error is the reason I avoid it for anything where landing exactly on the endpoint matters.

np.zeros allocates memory and writes zero into every element. np.empty allocates the same memory but skips initialization entirely, so the array starts with whatever garbage values happened to already be sitting in that memory.

a = np.empty((1000, 1000)) # fast allocation, uninitialized values
a[:] = 0 # you still have to fill it yourself if you need known values

The only reason to use empty is speed when you're about to overwrite every element anyway, inside a loop that fills the array from a computation, for instance. If there's any chance you read from the array before writing to every index, empty will hand you silent garbage instead of an error, which is a genuinely nasty bug to track down.

.shape is a tuple giving the length of each axis..ndim is the number of axes, just len(shape)..size is the total number of elements across every axis, the product of all the values in shape.

a = np.zeros((3, 4, 5))
a.shape  # (3, 4, 5)
a.ndim  # 3
a.size  # 60, which is 3 * 4 * 5

Mixing these up is a small mistake but a telling one, confusing size with shape usually means someone hasn't had to reason carefully about multidimensional arrays yet, only ever worked with 1D ones where the distinction barely shows up.

Basic slicing uses a start:stop:step syntax and always returns a view into the original array's memory, no data is copied. Fancy indexing means indexing with a list or array of integers, or with a boolean array, and it always returns a new array with copied data.

a = np.arange(10)
view = a[2:5]     # slicing, view
copy = a[[2, 3, 4]]  # fancy indexing, copy

view[0] = 99
print(a[2]) # 99, original changed

copy[0] = -1
print(a[2]) # still 99, original untouched

Negative indices count from the end of the array, -1 is the last element, and a negative step walks the array backward instead of forward, which is exactly the mechanism NumPy uses to reverse an array without allocating a new one.

a = np.arange(10)
a[-1]   # 9, last element
a[::-1]  # reversed array, still a view, no copy
a[::-2]  # every other element, reversed

Broadcasting is the set of rules NumPy uses to let operations work on arrays of different shapes without you manually replicating the smaller one to match the larger one's size. Instead of allocating a full-size copy of the smaller array, NumPy conceptually stretches it across the missing dimensions during the operation itself.

a = np.array([[1, 2, 3], [4, 5, 6]])
b = np.array([10, 20, 30])
a + b
# [[11, 22, 33],
# [14, 25, 36]]

Yes, a scalar is treated as a shape-() array with zero dimensions, and it satisfies the broadcasting rules trivially against any shape since there's nothing to compare it against except itself being stretched everywhere. It's the simplest possible case of the same general mechanism that handles the more complex shape combinations.

a = np.array([[1, 2], [3, 4]])
a * 2
# [[2, 4],
# [6, 8]]

People rarely think of scalar multiplication as "broadcasting" because it's so intuitive, but it's worth naming explicitly in an interview since it shows you understand broadcasting as one general rule rather than a handful of separate special cases you memorized.

A Python for loop over an array pulls one element out as a Python object on every iteration, applies the operation through Python's interpreter, and pays that overhead every single time. A vectorized operation dispatches once into a compiled C loop that walks the underlying memory directly, with no per-element Python interpreter overhead at all.

a = np.arange(1_000_000)

# slow: Python-level loop, one interpreter step per element
total = 0
for x in a:
  total += x

# fast: single call into a compiled C loop
total = a.sum()

axis names the dimension that gets collapsed, not the dimension that survives. axis=0 collapses rows, summing down each column, and axis=1 collapses columns, summing across each row. It's a common source of confusion because the intuitive reading, "sum along axis 0 means go across row 0," is backward from what actually happens.

a = np.array([[1, 2, 3], [4, 5, 6]])
a.sum(axis=0)  # [5, 7, 9], one value per column, rows collapsed
a.sum(axis=1)  # [6, 15], one value per row, columns collapsed

max and min return the value itself, the largest or smallest element. argmax and argmin return the index of that element instead of the value. On a multidimensional array without an axis argument, argmax returns the index into the flattened array, which is a common trip-up if you forget to unravel it back into a multidimensional index.

a = np.array([[1, 9, 3], [7, 2, 8]])
a.max()     # 9, the value
a.argmax()    # 1, index into the flattened array
np.unravel_index(a.argmax(), a.shape) # (0, 1), the real row/col position

sum and prod collapse an array down to one final total or product. cumsum and cumprod return an array the same size as the input, where each position holds the running total or running product up through that point rather than only the final answer.

a = np.array([1, 2, 3, 4])
a.sum()   # 10, single value
a.cumsum()  # [1, 3, 6, 10], running total at each step
a.cumprod() # [1, 2, 6, 24], running product at each step

These come up constantly for anything involving running totals over time, cumulative revenue, a running balance, since the alternative is writing a manual loop that keeps a running variable, which is exactly the pattern NumPy is built to replace.

np.random.seed sets the state of NumPy's older, global legacy random module, RandomState, which every call to np.random.rand or np.random.randint reads from implicitly. default_rng() constructs a new, local Generator object with its own independent state that you pass around explicitly instead of relying on hidden global state.

np.random.seed(42)
np.random.rand(3)       # legacy global state

rng = np.random.default_rng(42)
rng.random(3)          # modern, explicit local state

default_rng is the recommended approach in current NumPy documentation (NumPy, Random Generator). It avoids the global-state footgun where two unrelated pieces of code both reseed or draw from the same global generator and quietly affect each other's output.

rng.random() draws floats uniformly from [0, 1). rng.integers(low, high) draws integers, and by default high is exclusive, matching Python's range convention, so integers(0, 10) never returns 10 itself.

rng = np.random.default_rng(0)
rng.random(3)       # e.g. [0.63, 0.82, 0.02], floats in [0, 1)
rng.integers(0, 10, 3)   # integers from 0 to 9 inclusive, 10 excluded
rng.integers(0, 10, 3, endpoint=True) # 0 to 10 inclusive, opt in explicitly

Forgetting that high is exclusive by default is a genuinely common off-by-one bug, especially for anyone coming from a language or library where the upper bound is inclusive by convention.

np.array always makes a copy of the input by default. np.asarray only copies if the input isn't already an ndarray of the requested dtype, otherwise it returns the original array unchanged, no copy at all.

a = np.array([1, 2, 3])
b = np.asarray(a)
print(b is a) # True, no copy since a is already an ndarray with a compatible dtype

c = np.array(a)
print(c is a) # False, np.array always copies

Use asarray inside a function that accepts array-like input but shouldn't force an unnecessary copy when the caller already passed a proper ndarray. It's a small detail that matters a lot once you're writing library code other people call in a hot path.

ndarray requires every dimension to have a consistent size, so a list of lists where the inner lists have different lengths can't form a proper multidimensional array. Recent NumPy versions raise a ValueError for this instead of silently falling back to an object array like older versions used to.

np.array([[1, 2, 3], [4, 5]])
# ValueError: setting an array element with a sequence.
# The requested array has an inhomogeneous shape

If you actually need a container of differently-sized sequences, the honest answer is that NumPy is the wrong tool, a Python list of arrays, or an object-dtype array created with dtype=object explicitly, is what you want, and you give up vectorized math the moment you do that.

A boolean mask is an array of True/False values the same shape as the array you're indexing, or at least a compatible shape. Indexing with it returns only the elements where the mask is True, flattened into a new one-dimensional array.

a = np.array([1, -2, 3, -4, 5])
mask = a > 0
print(a[mask])  # [1, 3, 5]

a[a < 0] = 0   # in-place conditional replacement, no explicit loop
print(a)     # [1, 0, 3, 0, 5]

The reason this comes up constantly in interviews is that it's the direct replacement for a Python if-check inside a for loop, and combined with assignment it lets you do conditional replacement in a single vectorized line instead of iterating.

Called with three arguments, np.where(condition, x, y) is a vectorized ternary, it returns x where the condition is True and y where it's False, elementwise. Called with just a condition, it returns the indices where that condition is True, as a tuple of arrays, one per dimension.

a = np.array([1, -2, 3, -4])
np.where(a > 0, a, 0)     # [1, 0, 3, 0], replace negatives with 0
np.where(a > 0)        # (array([0, 2]),), indices where True

Candidates who only know the three-argument form get tripped up the first time they hit code that calls np.where with one argument and expects it to return values instead of indices.

Ellipsis expands to as many full slices (:) as needed to fill out the remaining dimensions, so you don't have to write one : per axis when you only care about indexing the first or last dimension of a high-dimensional array.

a = np.zeros((3, 4, 5, 6))
a[..., 0].shape  # (3, 4, 5), same as a[:, :, :, 0]
a[0,...].shape  # (4, 5, 6), same as a[0, :, :, :]

It's a small convenience but it shows up constantly in real image and tensor code where arrays commonly carry four or five dimensions and writing out every colon by hand gets unreadable fast.

Functionally, np.take(a, indices, axis=...) does the same gather operation as a[indices] along a given axis, it's a copy either way. The practical differences are that take lets you specify the axis explicitly as an argument rather than needing a full tuple of slices, and it accepts an out parameter and a mode argument controlling what happens with out-of-bounds indices, wrap or clip instead of raising.

a = np.arange(12).reshape(3, 4)
np.take(a, [0, 2], axis=1)     # same values as a[:, [0, 2]]
np.take(a, [5], mode='wrap')    # wraps out-of-bounds index instead of erroring

Check the base attribute. If an array is a view, its.base points back to the original array it shares memory with. If it's an independent array, whether created fresh or via.copy(),.base is None.

a = np.arange(10)
view = a[2:5]
copy = a[[2, 3, 4]]

print(view.base is a)  # True, view shares memory with a
print(copy.base is None) # True, copy owns its own memory

flags.owndata tells you roughly the same thing from the other direction, but base is usually the more direct answer when someone asks you to prove it in an interview.

reshape returns a view whenever the new shape is compatible with the array's existing strides, which is the common case for a contiguous array. It falls back to returning a copy when the requested shape can't be expressed with any valid stride pattern given the current memory layout, most often after a non-contiguous slice or a transpose.

a = np.arange(12).reshape(3, 4)
b = a.reshape(4, 3)     # view, a is contiguous
c = a.T.reshape(12)     # copy, a.T isn't contiguous in the needed order
print(c.base is None)    # True

If you specifically need a guaranteed view and want an error rather than a silent copy, use.reshape() cautiously or check contiguity with.flags first. If you need a guaranteed flat copy regardless,.flatten() always copies, while.ravel() behaves like reshape and returns a view when it can.

NumPy compares shapes dimension by dimension starting from the trailing (rightmost) axis and working backward. Two dimensions are compatible if they're equal, or if one of them is 1. If a shape has fewer dimensions than the other, it's treated as if padded with size-1 dimensions on the left (NumPy, Broadcasting rules).

# shape (256, 256, 3) and shape (3,)
# compare trailing dims: 3 == 3, ok
# missing dims on the smaller array padded with 1: (1, 1, 3)
# result shape: (256, 256, 3)

If neither condition holds for some pair of dimensions, NumPy raises ValueError instead of guessing what you meant. Knowing that it's a strict, mechanical rule rather than a heuristic is the part that actually helps you predict failures before running the code.

Shapes (4, 3) and (4,) look like they should line up since both have a 4 somewhere, but broadcasting compares from the trailing axis, so it compares 3 against 4 first, not 4 against 4.

a = np.zeros((4, 3))
b = np.array([1, 2, 3, 4]) # shape (4,)
a + b
# ValueError: operands could not be broadcast together with shapes (4,3) (4,)

The fix, if the intent really was to add a per-row value, is to reshape b to (4, 1) explicitly so the trailing comparison becomes 3 vs 1, which is compatible. That reshape step is exactly the part candidates skip when they just guess at shapes instead of working through the rule.

np.newaxis inserts a new axis of size 1 at a specific position, which is exactly the tool you reach for when broadcasting's automatic left-padding doesn't line dimensions up the way you actually need, like the row-vector-versus-column-vector case from a shape-mismatch error.

a = np.array([1, 2, 3])     # shape (3,)
col = a[:, np.newaxis]      # shape (3, 1)
row = a[np.newaxis, :]      # shape (1, 3)

col + row # shape (3, 3), outer sum

A ufunc, universal function, operates elementwise on ndarrays and supports broadcasting automatically, and most of them are written in C and compiled, which is where the speed comes from. np.add, np.multiply, np.sqrt, and the operators they back (+, *) are all ufuncs. Not every NumPy function is a ufunc, np.sum and np.concatenate, for instance, are regular Python functions that call into compiled code but don't follow the ufunc protocol themselves.

np.add(np.array([1, 2, 3]), np.array([10, 20, 30])) # ufunc
np.add.reduce(np.array([1, 2, 3]))          # ufuncs expose reduce, accumulate, outer

out tells the ufunc to write its result directly into an existing array's memory instead of allocating a brand new array for the output. It avoids the allocation and the garbage collection of the temporary result that you'd otherwise create and immediately discard.

a = np.arange(1_000_000, dtype=np.float64)
result = np.empty_like(a)
np.multiply(a, 2, out=result) # no new allocation for the output

This matters most inside a loop that runs the same operation repeatedly on same-shaped arrays, streaming data processing or an optimization loop, where the repeated allocation and deallocation adds up to real, measurable overhead over thousands of iterations.

Calling a ufunc directly applies it elementwise between two arrays. reduce applies it repeatedly along an axis to collapse an array down to a single aggregated value per remaining axis, the same idea as sum but generalized to any ufunc. accumulate does the same repeated application but keeps every intermediate result instead of only the final one.

a = np.array([1, 2, 3, 4])
np.add.reduce(a)    # 10, same as a.sum()
np.add.accumulate(a)  # [1, 3, 6, 10], running total, same as np.cumsum(a)
np.multiply.reduce(a) # 24, same as a.prod()

Knowing that sum and cumsum are really just add.reduce and add.accumulate under a friendlier name is a good signal that you understand ufuncs as a general mechanism rather than a list of unrelated built-in functions.

By default, reducing along an axis removes that axis from the result's shape entirely. keepdims=True keeps it as a size-1 axis instead, so the result still has the same number of dimensions as the original.

a = np.ones((4, 5))
a.sum(axis=1).shape         # (4,), axis removed
a.sum(axis=1, keepdims=True).shape # (4, 1), axis kept as size 1

You need keepdims whenever you're about to broadcast the reduced result back against the original array, normalizing each row by its own sum, say. Without keepdims, the shapes won't line up for broadcasting without an extra reshape or newaxis step.

For plain 2D matrix multiplication they behave the same, and @ is really just operator syntax for np.matmul. They diverge on higher-dimensional arrays and on 1D arrays. np.dot on two 1D arrays computes a scalar inner product, same as matmul, but np.dot treats extra dimensions differently from matmul when arrays have more than 2 dimensions, matmul specifically treats the trailing two dimensions as stacked matrices and broadcasts over the leading ones, while dot's behavior with N-D arrays follows a different, less intuitive rule involving sums over the last axis of one and the second-to-last of the other.

a = np.random.default_rng(0).random((10, 3, 4))
b = np.random.default_rng(1).random((10, 4, 5))
a @ b      # shape (10, 3, 5), 10 independent matmuls, one per stack
np.matmul(a, b) # identical result to @

My rule of thumb: use @ for anything that's genuinely matrix multiplication, including batched, and reach for np.dot only for the simple scalar or vector inner-product case where the extra behavior doesn't matter.

Computing the explicit inverse and multiplying, np.linalg.inv(A) @ b, does more arithmetic than necessary and is numerically less stable, small floating point errors in computing the full inverse get amplified before you even use it. np.linalg.solve is built specifically to solve Ax = b directly using a decomposition method, without ever forming the inverse at all.

x = np.linalg.solve(A, b)     # preferred, faster and more numerically stable
x_bad = np.linalg.inv(A) @ b    # works, but slower and less accurate

The only real reason to compute the inverse explicitly is if you actually need the inverse matrix itself for something beyond solving one system, otherwise solve is strictly the better call.

A single NaN in an array poisons any ordinary reduction that touches it, np.mean and np.sum both propagate NaN through the whole result rather than skipping it, because NaN compares unequal to everything, including itself, and any arithmetic involving it produces NaN. NumPy provides nan-aware versions specifically to opt out of that behavior.

a = np.array([1.0, 2.0, np.nan, 4.0])
np.mean(a)   # nan, the whole result is poisoned
np.nanmean(a) # 2.333..., NaN entries skipped

The trade-off worth knowing: nanmean silently changes the denominator, since it divides by the count of non-NaN values, not the original array length. That's usually what you want, but it can quietly mask a data quality problem if you never check how many NaNs got dropped.

The new design splits responsibilities. A BitGenerator, PCG64 by default, is responsible purely for producing a stream of random bits, and a Generator sits on top of it and transforms those bits into useful distributions, normal, binomial, uniform, and so on. That separation lets NumPy swap in better underlying bit-generation algorithms over time without touching the distribution-sampling code built on top.

RandomState's older Mersenne Twister algorithm is locked in place for backward compatibility, since changing it would silently change the output of every seeded RandomState call ever written. Generator makes no such backward-compatibility promise about exact output across NumPy versions, which is precisely what frees it up to keep improving.

rng.choice with replace=False samples without putting drawn elements back, which is what you want any time duplicates would be wrong, picking a random validation subset, shuffling a deck without repeats. For parallel work across multiple processes, each worker needs its own independent stream rather than all of them sharing or accidentally reusing the same seed.

rng = np.random.default_rng(0)
rng.choice(100, size=10, replace=False)  # 10 distinct values from 0-99

# independent, non-overlapping streams for parallel workers
streams = np.random.default_rng(0).spawn(4)

spawn() is the modern answer to the parallel-seeding problem, it hands each worker a child generator that's statistically independent from the others, instead of the older and riskier pattern of manually adding a worker index to a shared seed.

Floating point arithmetic accumulates small rounding error, so two values that are mathematically identical, like 0.1 + 0.2 and 0.3, can differ in their last bits and fail an exact == comparison even though they represent the same intended quantity.

0.1 + 0.2 == 0.3      # False
np.isclose(0.1 + 0.2, 0.3) # True, compares within a tolerance

np.isclose and np.allclose compare within a configurable relative and absolute tolerance instead of demanding bit-exact equality, and that's the right default for basically any comparison involving computed floats rather than values you set directly and never touched.

concatenate joins arrays along an existing axis, so the input arrays need to already agree on every other dimension, and the result has the same number of dimensions as the inputs. stack introduces a brand new axis and requires all input arrays to have exactly the same shape, increasing the dimensionality by one. vstack and hstack are convenience wrappers around concatenate for the common 1D and 2D row or column cases.

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.concatenate([a, b])  # shape (6,), same dimensionality, joined along axis 0
np.stack([a, b])     # shape (2, 3), new axis introduced
np.vstack([a, b])     # shape (2, 3), stacks as rows

The mistake I see most often is reaching for stack when concatenate was needed, or the reverse, then being confused why the output has one dimension too many or too few. Checking the resulting.shape immediately after either call catches it fast.

np.unique returns the sorted distinct values in an array. With return_counts=True it also returns a second array, the same length, giving how many times each of those distinct values appeared. It's the direct NumPy equivalent of a value-counts operation.

labels = np.array([0, 1, 1, 2, 2, 2, 0])
values, counts = np.unique(labels, return_counts=True)
# values: [0, 1, 2], counts: [2, 2, 3]

A common real use is checking class balance in a labeled dataset before training a model, one line tells you exactly how skewed the classes are without writing a manual counting loop.

A slice like 0:2 describes a regular, evenly-spaced pattern through memory, a fixed start, stop, and stride, which NumPy can express entirely through the existing strides metadata without touching the underlying bytes. Fancy indexing with an arbitrary index list has no such regular pattern, the indices could be [0, 1] or they could be [7, 2, 9, 0], so there's no fixed stride that can describe jumping between them.

Because there's no stride that generalizes to arbitrary index lists, NumPy has no way to build a view for fancy indexing, it has to gather the selected elements into a freshly allocated, contiguous array instead. This is also why fancy indexing is slower than a comparable slice for large arrays, it always pays a full memory copy.

strides is a tuple giving the number of bytes to step in memory to move one position along each axis. For a C-contiguous (row-major) 2D array of float64, moving one column over costs 8 bytes, moving one row down costs 8 times the number of columns.

a = np.zeros((3, 4), dtype=np.float64)
print(a.strides)     # (32, 8), 4 cols * 8 bytes per row, 8 bytes per element

b = a.T
print(b.shape, b.strides) # (4, 3), (8, 32), swapped, no data moved

Transpose doesn't touch the underlying data at all, it just returns a new view with the shape and strides tuples swapped. That's why transposing a huge array is instant regardless of size, and it's also exactly why it's a view, writing to a.T writes to a itself.

C order (row-major, NumPy's default) stores a row's elements consecutively in memory, so the last axis varies fastest. Fortran order (column-major) stores a column's elements consecutively, so the first axis varies fastest. Both describe the same logical array, just a different physical byte layout.

a = np.zeros((1000, 1000), order='C')
b = np.zeros((1000, 1000), order='F')

It matters for performance whenever you're looping over an array's memory in a specific order, iterating column-by-column over a C-contiguous array walks memory with a large stride and thrashes the CPU cache, while row-by-row iteration walks it sequentially. It also matters at library boundaries, some Fortran-backed linear algebra routines expect Fortran order and will silently copy your array to convert it if you hand them the wrong layout, which is a real, measurable slowdown on large matrices.

.view(dtype) reinterprets the exact same underlying bytes as a different dtype, without converting any values, so the byte count has to be compatible or the array's last dimension changes size..astype(dtype) actually converts each value to the new type, allocating a new array with newly computed data.

a = np.array([1, 2, 3], dtype=np.int32)
a.view(np.float32)  # reinterprets the same 4 bytes per element as float32, garbage values
a.astype(np.float32) # [1.0, 2.0, 3.0], actual conversion

This distinction trips people up constantly because both return an array of a different dtype..view is a memory-layout trick, useful for things like inspecting the raw bits of a float or packing/unpacking binary data..astype is a value conversion, useful for basically everything else.

No, broadcasting itself doesn't materialize a full-size copy of the smaller array in memory, it uses a stride of 0 along the stretched dimension so the same underlying values get reread repeatedly without duplication. That part genuinely is close to free.

What isn't free is the output. The result of a broadcast operation is always allocated at full size, so adding a (10000, 1) array to a (1, 10000) array produces a (10000, 10000) output, which is 100 million elements from two arrays that together only held 20000 values. The broadcasting mechanism saved memory on the inputs, but the output size can still blow up in a way that's easy to miss until the process runs out of RAM.

No, and this is a common trap. np.vectorize is a convenience wrapper that takes a plain Python function and applies it elementwise using a Python-level loop internally, it doesn't compile anything or push the loop into C. It exists purely for broadcasting convenience and readability, not performance.

def custom_logic(x):
  return x ** 2 if x > 0 else 0

slow = np.vectorize(custom_logic)(a)  # still a Python loop underneath
fast = np.where(a > 0, a ** 2, 0)   # actually vectorized, no per-element Python call

If you genuinely need custom conditional logic per element and want real speed, np.where, boolean masking, or reformulating the logic as arithmetic almost always beats np.vectorize. I've seen candidates confidently reach for np.vectorize as "the fast way to apply a function" and it's the opposite of that.

matmul treats the last two axes of each array as the matrix dimensions and treats every axis before that as a batch dimension that gets broadcast normally. So (10, 3, 4) @ (10, 4, 5) runs 10 independent 3x4 by 4x5 matrix multiplications and stacks the 10 results, giving shape (10, 3, 5).

# leading dims (10,) match on both sides, so they broadcast together
# trailing dims (3,4) and (4,5) are matmul-compatible: inner dims 4 == 4
# result: (10, 3, 5)

This is exactly the mechanism batched neural network layers rely on, one weight matrix applied across a whole batch of inputs at once, without writing an explicit Python loop over the batch dimension.

A structured array is an ndarray whose dtype defines multiple named fields, each with its own type, packed into a single record per element, closer to a C struct than a plain numeric array.

dt = np.dtype([('name', 'U10'), ('age', 'i4'), ('weight', 'f8')])
people = np.array([('Ann', 32, 68.5), ('Bo', 41, 81.2)], dtype=dt)
people['age']  # array([32, 41])
people[0]     # ('Ann', 32, 68.5)

Reach for a structured array when you need a fixed-layout, low-level binary record, reading a fixed-format binary file directly with np.fromfile, or interoperating with a C struct without a copy. For anything involving heterogeneous, labeled tabular data with real analysis work, joins, groupby, missing-value handling, pandas is almost always the better tool. Structured arrays are a lower-level building block, not a pandas replacement.

Operations between two arrays of different but compatible dtypes promote to a common type that can represent both, int plus float promotes to float, for instance. That's usually safe. The dangerous edge case is fixed-width integer overflow within the same dtype, since NumPy integers, unlike Python's arbitrary-precision int, wrap around silently instead of raising an error.

a = np.array([120, 100], dtype=np.int8)
a + a
# array([-16, -56], dtype=int8) # wrapped around 127, no warning, no error

This is a real production bug pattern, someone picks a small dtype to save memory on a large array, values grow past the type's range during a later computation, and the result silently wraps instead of failing loudly. Checking the dtype's actual range with np.iinfo before choosing a narrow integer type is the cheap insurance against it.

np.ma.MaskedArray pairs the data array with a separate boolean mask array marking which entries are invalid, and masked-aware functions automatically exclude those entries from computations. Using NaN as a sentinel is simpler, it needs no separate structure, but NaN only exists for floating point dtypes, you can't represent a missing integer or missing boolean with NaN without upcasting the whole array to float first.

data = np.ma.masked_array([1, 2, -999, 4], mask=[False, False, True, False])
data.mean()  # 2.333..., -999 excluded because masked, no need to filter first

ints = np.array([1, 2, 3])
# can't represent a missing int value as NaN without converting to float

In practice most teams reach for NaN because it's simpler and plays well with pandas, and reserve masked arrays for cases involving integer or boolean data where NaN genuinely isn't representable, or where you need the mask itself to survive as separate metadata rather than being baked into the values.

einsum expresses arbitrary sums-of-products over array axes using Einstein summation notation, a compact string describing which axes to multiply together and which to sum over. It generalizes matmul, transpose, trace, and outer product into one function, and for operations that don't map cleanly onto an existing NumPy function, it's often both clearer and faster than writing the equivalent nested loop or chain of reshapes and transposes yourself.

A = np.random.default_rng(0).random((3, 4))
B = np.random.default_rng(1).random((4, 5))

np.einsum('ij,jk->ik', A, B)  # equivalent to A @ B
np.einsum('ii', A[:3, :3])   # trace, sum of the diagonal

The honest caveat is that einsum's notation has a real learning curve, and for a plain matrix multiply, A @ B is more readable and just as fast. Where einsum earns its place is custom contractions across three or more tensors that would otherwise need an ugly chain of transpose and reshape calls to express with ordinary functions.

Comparing trailing dimensions, the second axis is 1 versus 4, compatible because one side is 1, and it stretches to 4. Comparing the first axis, 3 versus 1, also compatible for the same reason, and it stretches to 3. Both arrays end up broadcast to shape (3, 4) before the operation applies.

a = np.array([[1], [2], [3]])    # shape (3, 1)
b = np.array([[10, 20, 30, 40]])  # shape (1, 4)
a + b
# [[11, 21, 31, 41],
# [12, 22, 32, 42],
# [13, 23, 33, 43]]

This is the pattern behind computing an outer sum or outer product without an explicit loop, and it's the same trick that underlies a lot of distance-matrix and pairwise-comparison code in real numerical work.