12 - NumPy: Arrays, Vectorisation and Broadcasting

Machine Learning Engineer course · Module 01 · Python for data · Lesson 5

You have a loop that sums the squares of a million numbers and takes a tenth of a second, and a one-line expression that does the same job in less than a thousandth. This lesson, which comes after floating-point numbers, explains where that difference comes from and gives you the three ideas that hold up almost all machine learning code: views, which do not copy, broadcasting, which aligns different shapes, and axes, which say along which dimension you compute. By the end you can predict the shape of an operation's result before running it and you recognise the shape error that produces a wrong number without any message.

The commonplace says NumPy is fast because it is written in C. That is only half the truth: a C function called a million times from a Python loop stays slow. NumPy is fast because the data are homogeneous and contiguous in memory and because a single call works on all the elements. It pays to learn it well because almost every data library (Pandas, scikit-learn, Matplotlib) exchanges numbers in the form of NumPy arrays, or of objects that resemble them.

Why a Python list is not enough

The analogy: the mailbox versus the egg carton

A Python list is a row of tickets. Each ticket does not contain the number: it contains the address of a box in a warehouse, where a complete float object sits (with its reference counter, its type, its value). To sum a thousand numbers the interpreter has to go a thousand times to look for a box around the warehouse, open it, check what type it is.

A NumPy array is an egg carton: all of the same type, one next to the other, in a single block. To sum them you just walk along the block. There is nothing to check, because the carton has a single label on its lid: "8-byte eggs".

Let us measure it

import sys
import numpy as np

lista = [float(i) for i in range(1_000_000)]
array = np.arange(1_000_000, dtype=np.float64)

byte_lista = sys.getsizeof(lista) + sum(sys.getsizeof(x) for x in lista)
print("list: bytes of the pointers    :", sys.getsizeof(lista))
print("list: bytes of a single float  :", sys.getsizeof(lista[1]))
print("list: total with the objects   :", byte_lista)
print("array: bytes                   :", array.nbytes)
print("memory ratio                   :", round(byte_lista / array.nbytes, 1))

Output obtained:

list: bytes of the pointers    : 8448728
list: bytes of a single float  : 24
list: total with the objects   : 32448728
array: bytes                   : 8000000
memory ratio                   : 4.1

The list spends 8 bytes for every pointer and another 24 for every float object, so about 32 bytes per number against 8. Beyond the space, locality matters: a contiguous block is read in sequence, and a computer's memory is much faster at reading in sequence than at jumping from one address to another.

Trade-off. An array has a single type and a fixed size. If you have to mix strings, numbers and None, or keep adding elements, a list (or a DataFrame) is the right tool. NumPy pays for its speed with rigidity.

Anatomy of an array: shape, dtype, strides

An ndarray is, concretely, a block of memory plus three numbers that say how to read it:

• dtype: what kind of number each cell is (float64, int32, bool...). It also fixes how many bytes it takes (itemsize).
• shape: how many cells there are along each dimension, for example (3, 4) for 3 rows and 4 columns.
• strides: by how many bytes you move in memory to go to the next element along each dimension.

The first two numbers say what is in the block, the third says how to walk through it. Separating the two is the design choice that makes transposes, slices and reshapes cheap: none of these gestures touches the data, only the way of reading it changes. The link between the three is a single formula.

The formula that ties everything together is the address of the element a[i, j]:

address(i, j) = base + i · strides[0] + j · strides[1]

For a (3, 4) matrix of float64 stored by rows, moving to the next column shifts by 8 bytes (one float64), moving to the next row shifts by 4 · 8 = 32 bytes. So strides = (32, 8).

import numpy as np

x = np.arange(12, dtype=np.float64).reshape(3, 4)
print(x)
print("shape:", x.shape, "| dtype:", x.dtype, "| itemsize:", x.itemsize, "| strides:", x.strides)

t = x.T                      # the transpose
print("transpose strides:", t.strides, "| shares memory with x:", np.shares_memory(x, t))

Output obtained:

[[ 0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]]
shape: (3, 4) | dtype: float64 | itemsize: 8 | strides: (32, 8)
transpose strides: (8, 32) | shares memory with x: True

The transpose does not move a single number: it swaps the two strides from (32, 8) to (8, 32). The way of reading changes, not the content. This is why x.T costs practically nothing, and also why NumPy can afford so many "free" shape-changing operations.

reshape works in the same way: if the total number of elements adds up, it recomputes shape and strides without touching the data (when the layout in memory allows it).

Creating arrays

import numpy as np

print(np.zeros((2, 3)))
print(np.ones(3, dtype=np.int64))
print(np.eye(3))
print(np.arange(0, 1, 0.25))          # fixed step, right end excluded
print(np.linspace(0, 1, 5))           # fixed number of points, ends included
print(np.full((2, 2), 7.5))

Output obtained:

[[0. 0. 0.]
 [0. 0. 0.]]
[1 1 1]
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
[0.   0.25 0.5  0.75]
[0.   0.25 0.5  0.75 1.  ]
[[7.5 7.5]
 [7.5 7.5]]

A real trap: arange with a decimal step. np.arange(start, stop, step) produces ceil((stop − start) / step) elements (this is the formula in the documentation of numpy.arange, which also warns that the last element may turn out greater than stop). If the quotient is an integer "by a whisker" in floating point, the outcome is surprising:

import numpy as np

a = np.arange(0, 2.1, 0.3)
print("elements:", len(a), "| last:", a[-1], "| quotient 2.1/0.3 =", 2.1 / 0.3)
print("linspace, 8 exact points:", np.linspace(0, 2.1, 8))

Output obtained:

elements: 8 | last: 2.1 | quotient 2.1/0.3 = 7.000000000000001
linspace, 8 exact points: [0.  0.3 0.6 0.9 1.2 1.5 1.8 2.1]

On paper 2.1 / 0.3 = 7 exactly, so 7 elements (from 0 to 1.8) and the right end excluded. But 0.3 and 2.1 are not representable (previous lesson): the computed quotient is 7.000000000000001, the ceiling is 8, and arange returns eight elements, with the last one coinciding with the end that was supposed to be excluded. When you know how many points you want, use np.linspace(start, stop, n): the number of points is a parameter, not a result of floating-point arithmetic.

Numeric types and silent overflow

Every integer dtype has a finite range. With int8 (8 bits with sign) the values go from −128 to 127. What happens beyond? The counter wraps around, like the odometer of an old car.

import numpy as np

piccolo = np.array([127], dtype=np.int8)
print("127 + 1 in int8 :", piccolo + 1)
print("int8 range      :", np.iinfo(np.int8).min, np.iinfo(np.int8).max)

img_a = np.array([200], dtype=np.uint8)      # like the value of a pixel
img_b = np.array([100], dtype=np.uint8)
print("uint8: 200 + 100          :", img_a + img_b, "(wrapped: 300 - 256)")
print("uint8 converted first     :", img_a.astype(np.int64) + img_b.astype(np.int64))

print("float32 of 0.1  :", float(np.float32(0.1)))
print("float64 of 0.1  :", float(np.float64(0.1)))

Output obtained:

127 + 1 in int8 : [-128]
int8 range      : -128 127
uint8: 200 + 100          : [44] (wrapped: 300 - 256)
uint8 converted first     : [300]
float32 of 0.1  : 0.10000000149011612
float64 of 0.1  : 0.1

On arrays, no error and no warning: just a wrong number. (On single scalars NumPy may emit a RuntimeWarning, but computations on arrays, the ones that matter in data work, wrap around silently.) The derivation is simple: with n bits there are 2ⁿ values; in arithmetic modulo 2ⁿ, 300 mod 256 = 44. 8-bit images are the classic case: adding two uint8 images without converting them first produces dark pixels where they should be light.

float32 takes half the memory of float64 but has about 7 significant decimal digits against 15-16. For neural networks it is often enough (and halves memory and time); for summing millions of terms or computing variances often not, as seen in the lesson on floating-point numbers.

Practical rule: let NumPy use the default type (float64, int64) until you measure a memory or speed problem. Optimising the type before you have a problem adds risk with no demonstrated benefit.

Indexing: views versus copies

This is the section that avoids the most bugs. There are three ways of selecting, with different behaviours:

The complete rule is in the guide Copies and views: basic indexing (slices) always returns a view, advanced indexing (lists of indices, masks) always a copy.

Analogy. A view is a window onto a room: if someone moves a piece of furniture, you see it change from the window, and if you move it yourself through the window, it changes in the room. A copy is a photograph: from that moment the two lives are independent.

import numpy as np

a = np.arange(10)
vista = a[2:5]
vista[0] = 99                       # I modify the "window"
print("original after changing the view:", a)

copia = a[[2, 3, 4]]
copia[0] = -1                       # I modify the "photograph"
print("original after changing the copy:", a)

print("view shares memory:", np.shares_memory(a, vista))
print("copy shares memory:", np.shares_memory(a, copia))

indipendente = a[2:5].copy()        # to get an explicit copy
print("explicit copy shares:", np.shares_memory(a, indipendente))

Output obtained:

original after changing the view: [ 0  1 99  3  4  5  6  7  8  9]
original after changing the copy: [ 0  1 99  3  4  5  6  7  8  9]
view shares memory: True
copy shares memory: False
explicit copy shares: False

When you need which. A view is very fast and takes no memory, perfect for reading. If you want to modify a subset without touching the original, call .copy() explicitly. When in doubt, np.shares_memory(a, b) answers with certainty.

Vectorisation: one instruction, a million elements

Vectorising means replacing a Python loop with a single expression that NumPy executes internally in compiled code, on all the elements.

Let us compare two ways of computing the sum of the squares of a million numbers.

import numpy as np

import time

v = np.random.default_rng(0).random(1_000_000)

t0 = time.perf_counter()
s_ciclo = 0.0
for e in v.tolist():
    s_ciclo += e * e
t_ciclo = time.perf_counter() - t0

t0 = time.perf_counter()
s_vett = np.dot(v, v)
t_vett = time.perf_counter() - t0

print(f"Python loop : {t_ciclo:.4f} s")
print(f"vectorised  : {t_vett:.5f} s")
print(f"ratio       : about {t_ciclo / t_vett:.0f} times")
print("same result (within tolerance):", np.isclose(s_ciclo, s_vett))

Output obtained:

Python loop : 1.2632 s
vectorised  : 0.02296 s
ratio       : about 55 times
same result (within tolerance): True

Times change from one run to the next and from one machine to another, so the output above is a single run. It was taken while the machine was heavily loaded (load average above 40), which is why the loop is slower than usual and the ratio lower. In five successive runs on a Ryzen 5 5600H with Python 3.14.4 and NumPy 2.3.5 (the measurement of the Italian edition of this article, taken on a quiet machine), the loop took from 0.10 to 0.15 s, the vectorised version about 0.0009 s, and the ratio was between 105 and 168 times: the order of magnitude (a factor around a hundred) is stable, the exact figure is not. Note that the loop already starts from v.tolist(), that is with the conversion into Python objects left out of the measurement: it is a comparison favourable to the loop.

The reason is the one in the opening analogy. For each element, the Python loop has to read the float object, check its type, find the multiplication function, create a new object for the result and repeat for the sum. NumPy, with the dtype fixed in advance, runs a compiled loop that does only the multiplications and the sums.

The functions that operate element by element and that NumPy calls ufuncs (universal functions) are everywhere: np.exp, np.log, np.sqrt, +, *, the comparisons. They compose freely:

import numpy as np

x = np.linspace(-3, 3, 7)
sigmoide = 1 / (1 + np.exp(-x))          # definition: 1 / (1 + e^{-x})
print(np.round(x, 2))
print(np.round(sigmoide, 4))

Output obtained:

[-3. -2. -1.  0.  1.  2.  3.]
[0.0474 0.1192 0.2689 0.5    0.7311 0.8808 0.9526]

A single line, no loop, works for 7 elements as for 7 million.

When NOT to vectorise. If the computation of each element depends on the previous one (a recurrence), it cannot be parallelised as a block and a loop may be unavoidable (tools such as Numba exist to compile it, but they are outside this lesson). And a vectorisation that creates huge intermediate arrays can exhaust memory: a loop over blocks is sometimes the right answer.

Axes and aggregations

An aggregation function (sum, mean, max, std...) reduces an array. Along which dimension? axis decides. There is only one rule to remember:

The axis you indicate is the one that disappears.

For a matrix with 3 rows and 4 columns, axis=0 collapses the rows (you get one value per column: shape (4,)), axis=1 collapses the columns (one value per row: shape (3,)).

import numpy as np

m = np.array([[1, 2, 3, 4],
              [5, 6, 7, 8],
              [9, 10, 11, 12]])
print("sum of everything   :", m.sum())
print("axis=0 (per column) :", m.sum(axis=0), "shape", m.sum(axis=0).shape)
print("axis=1 (per row)    :", m.sum(axis=1), "shape", m.sum(axis=1).shape)
print("keepdims            :", m.sum(axis=1, keepdims=True).shape)

Output obtained:

sum of everything   : 78
axis=0 (per column) : [15 18 21 24] shape (4,)
axis=1 (per row)    : [10 26 42] shape (3,)
keepdims            : (3, 1)

keepdims=True keeps the axis with length 1: it is often needed for the broadcasting we will see shortly. In tabular data, the convention is rows = examples, columns = features: X.mean(axis=0) is the mean of each feature.

Broadcasting: the rules, derived

The analogy: the printer that repeats the label

You have 100 folders and a single label with today's date. You do not print 100 labels: you apply it to all the folders. Broadcasting does the same: when two arrays have different shapes, NumPy virtually "stretches" the smaller one, repeating it, without copying the data.

The rules

To combine two arrays, NumPy compares the shapes from the last dimension towards the first. Two dimensions are compatible if:

1. they are equal, or
2. one of the two is 1 (which is stretched until it matches the other).

If a dimension is missing (an array has fewer dimensions), it is considered of length 1 on the left. If some pair respects neither rule: error (ValueError, as the Broadcasting guide describes). The shape of the result is, in each position, the larger of the two values, with a single exception: a dimension of length 0 against one of length 1 gives 0 (for example (3, 1) and (0,) give (3, 0), verified).

Let us apply them to three cases, reading the shapes from the right:

import numpy as np

X = np.array([[1., 2., 3., 4.],
              [5., 6., 7., 8.],
              [9., 10., 11., 12.]])
media = X.mean(axis=0)                  # shape (4,)
print("X - mean:\n", X - media)         # (3,4) - (4,) -> (3,4)

colonna = np.array([[10], [20], [30]])  # (3,1)
riga = np.array([[1, 2, 3, 4]])         # (1,4)
print("table of sums:\n", colonna + riga)

v = np.array([1., 2., 3.])              # (3,)
try:
    X - v
except ValueError as e:
    print("expected error:", e)
print("fixed with v[:, None]:\n", X - v[:, None])   # (3,4) - (3,1)

Output obtained:

X - mean:
 [[-4. -4. -4. -4.]
 [ 0.  0.  0.  0.]
 [ 4.  4.  4.  4.]]
table of sums:
 [[11 12 13 14]
 [21 22 23 24]
 [31 32 33 34]]
expected error: operands could not be broadcast together with shapes (3,4) (3,) 
fixed with v[:, None]:
 [[0. 1. 2. 3.]
 [3. 4. 5. 6.]
 [6. 7. 8. 9.]]

The v[:, None] trick (equivalent to v.reshape(3, 1)) inserts an axis of length 1 and turns a row vector into a column. It is the most frequent gesture in broadcasting.

Why it costs no memory

The secret is in the strides: along a "stretched" axis, the stride is zero. The step in memory is null, so the same element is always read again.

import numpy as np

stirato = np.broadcast_to(np.arange(4), (3, 4))
print("shape:", stirato.shape, "| strides:", stirato.strides, "(zero bytes to move to the next row)")

Output obtained:

shape: (3, 4) | strides: (0, 8) (zero bytes to move to the next row)

Standardising a matrix

An operation you will meet in almost every model is standardisation: subtract from each column its mean and divide by its standard deviation, z = (x − μ) / σ. With broadcasting it is one line:

import numpy as np

rng = np.random.default_rng(1)
dati = rng.normal(loc=[50, 1000, 0.5], scale=[10, 300, 0.1], size=(500, 3))

z = (dati - dati.mean(axis=0)) / dati.std(axis=0)
print("mean of each column ~ 0:", np.allclose(z.mean(axis=0), 0))
print("std of each column  = 1:", np.allclose(z.std(axis=0), 1))

Output obtained:

mean of each column ~ 0: True
std of each column  = 1: True

The three columns started from very different scales (tens, thousands, tenths). After standardisation they all have mean 0 and standard deviation 1.

Pairwise distances without loops, and the hidden cost

Given n points in d dimensions, the Euclidean distance between every pair is obtained by extending the axes:

import numpy as np

rng = np.random.default_rng(1)
punti = rng.normal(size=(4, 2))
diff = punti[:, None, :] - punti[None, :, :]       # (4,1,2) - (1,4,2) -> (4,4,2)
distanze = np.sqrt((diff ** 2).sum(axis=-1))        # sum over the last axis -> (4,4)
print(np.round(distanze, 3))
print("diagonal (distance of a point from itself):", distanze.diagonal())

Output obtained:

[[0.    2.125 0.674 0.915]
 [2.125 0.    1.842 2.074]
 [0.674 1.842 0.    1.449]
 [0.915 2.074 1.449 0.   ]]
diagonal (distance of a point from itself): [0. 0. 0. 0.]

Trade-off. The intermediate array diff has shape (n, n, d): for n = 100,000 points and d = 10 you would need 10¹⁰ · 10 · 8 bytes, that is 800 GB. Broadcasting is elegant, but memory is not infinite. For large sets dedicated algorithms are used (for example scipy.spatial.distance.cdist, or computing in blocks).

Matrix product versus element-wise product

The * operator multiplies element by element (with broadcasting). The matrix product is written @ (or np.matmul, np.dot). For A of shape (m, n) and B of shape (n, p):

(A @ B)ᵢⱼ = Σₖ Aᵢₖ · Bₖⱼ, with a result of shape (m, p).

The number of columns of A must coincide with the number of rows of B: it is the condition for being able to do the sum over k.

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 1], [1, 0]])
print("element-wise A * B:\n", A * B)
print("matrix product A @ B:\n", A @ B)
print("B @ A differs from A @ B:", not np.array_equal(A @ B, B @ A))

Output obtained:

element-wise A * B:
 [[0 2]
 [3 0]]
matrix product A @ B:
 [[2 1]
 [4 3]]
B @ A differs from A @ B: True

Confusing * with @ is a very common mistake: no exception if the shapes are compatible with both, completely different results. The matrix product is not commutative: A @ B and B @ A are in general different matrices. It is the central operation of linear algebra; here we use it without developing its theory.

Masks, where and missing values

A boolean mask is an array of True/False of the same shape, and it selects the true elements.

import numpy as np

temperature = np.array([18.5, 22.0, np.nan, 30.1, 15.2, np.nan, 27.8])
print("nan == nan:", np.nan == np.nan)               # NaN is not equal to anything, not even itself
validi = ~np.isnan(temperature)
print("naive mean       :", temperature.mean())      # the NaN contaminates everything
print("mean of the valid:", temperature[validi].mean())
print("equivalent       :", np.nanmean(temperature))

caldo = np.where(temperature > 25, "hot", "mild")  # (NaN > 25 is False)
print(caldo)
print("clip between 20 and 28:", np.clip(temperature, 20, 28))

Output obtained:

nan == nan: False
naive mean       : nan
mean of the valid: 22.72
equivalent       : 22.72
['mild' 'mild' 'mild' 'hot' 'mild' 'mild' 'hot']
clip between 20 and 28: [20.  22.   nan 28.  20.   nan 27.8]

NaN ("not a number") is the value IEEE 754 reserves for undefined results. Its most insidious property: any comparison with NaN is false (including NaN == NaN), so you test only with np.isnan. And it propagates: the mean of an array that contains a single NaN is NaN. Deciding what to do with missing values (discard them, fill them, flag them) is an analysis choice, not a technical detail: we come back to it with Pandas.

Reproducible random numbers

For simulations and experiments you need a reproducible source of pseudo-random numbers: with the same seed you get the same sequence, hence the same result.

import numpy as np

rng = np.random.default_rng(seed=42)       # the modern way
print(rng.integers(0, 10, size=5))
print(np.round(rng.normal(size=3), 4))
print(rng.choice(["a", "b", "c"], size=4, replace=True))

rng2 = np.random.default_rng(seed=42)
print("same sequence with the same seed:", rng2.integers(0, 10, size=5))

Output obtained:

[0 7 6 4 4]
[ 0.9406 -1.951  -1.3022]
['c' 'c' 'c' 'c']
same sequence with the same seed: [0 7 6 4 4]

You use np.random.default_rng (a Generator object) and pass it to the functions that require it, instead of the global state np.random.seed, which is the older interface. The local generator prevents two parts of the program from disturbing each other; the NumPy documentation recommends it for new code. With different NumPy versions the sequence for the same seed is not guaranteed to be identical: the generator's compatibility policy reserves to the designers the possibility of changing the algorithm, and even calling rng.random() five times may not give the same numbers as rng.random(5). This is why in this course every experiment declares the versions used.

Typical mistakes (symptom, cause, remedy)

Silent broadcasting

It is the costliest trap, because it produces no error. A vector of predictions of shape (n, 1) and one of true values of shape (n,) do not give an error: they give an (n, n) matrix, and the mean of that matrix is an arbitrary number.

import numpy as np

y_vero = np.array([1.0, 2.0, 3.0])
y_prev = np.array([[1.0], [2.0], [4.0]])          # shape (3,1) instead of (3,)

mse_sbagliato = ((y_prev - y_vero) ** 2).mean()
mse_giusto = ((y_prev[:, 0] - y_vero) ** 2).mean()
print("(3,1) - (3,) has shape:", (y_prev - y_vero).shape)
print("WRONG MSE:", round(float(mse_sbagliato), 4), "| RIGHT MSE:", round(float(mse_giusto), 4))

Output obtained:

(3,1) - (3,) has shape: (3, 3)
WRONG MSE: 2.3333 | RIGHT MSE: 0.3333

Remedy: check the shapes with assert y_prev.shape == y_vero.shape at the start of error functions. A line that costs nothing.

The other traps

• I modify a view without meaning to. Symptom: the original data change. Remedy: .copy() when you need independence.
• Wrong axis in an aggregation. Symptom: unexpected result shape. Remedy: print the shape, remember "the axis indicated disappears".
• Overflow of small integers (uint8, int16). Remedy: convert before computing (astype).
• == between floats. Remedy: np.isclose, np.allclose (see the lesson on floating-point numbers).
• A Python loop over an array. Symptom: slowness. Remedy: look for the vectorised expression; if none exists, assess whether the loop is really unavoidable.

The first two traps produce a result of unexpected shape or content and are found by printing shape; overflow and float comparison instead produce plausible, wrong numbers, and must be prevented by construction. The loop over an array is the only one that gives no wrong answers, only slow ones: that is why it is the last to fix.

Best practices

• Print shape and dtype when something does not add up: it is the first debugging step, and it solves the majority of cases.
• Check shapes with assert in functions that take them for granted.
• Prefer array operations to loops, but do not sacrifice readability: an incomprehensible line is not a gain.
• Use keepdims=True or [:, None] explicitly to make broadcasting visible.
• Fix the seed and pass the Generator as an argument.
• Keep float64 until you measure a problem.

The six rules serve the same purpose: making visible what NumPy does in silence. Shapes and types printed or checked with assert turn a wrong number into an error with a name, and the fixed seed turns an unrepeatable experiment into one a colleague can redo. They cost a few lines at the start of a function and save hours of searching.

Mini-project

Write a function standardizza(X) that:

1. checks that X has two dimensions (otherwise it raises ValueError with a clear message);
2. returns (Z, media, dev_std), where Z is the standardised matrix;
3. handles columns with zero standard deviation (constant columns) without producing NaN: choose and document what to return.

The three requirements cover the input check, the output contract and the edge case, which is where numerical functions usually break. Test with an X of shape (100, 5) in which one column is constant, and check that Z.mean(axis=0) is close to zero with np.allclose. Try writing the function first, then compare it with the annotated solution below.

The annotated solution of the mini-project: standardizza

Let us put together views, axes and broadcasting in the function the mini-project asks for, paying attention to the edge cases.

import numpy as np

def standardizza(X):
    """Return (Z, mean, std). Constant columns stay at zero instead of becoming NaN."""
    X = np.asarray(X, dtype=float)
    if X.ndim != 2:
        raise ValueError(f"a matrix with 2 dimensions is required, got {X.ndim}")
    media = X.mean(axis=0)
    dev = X.std(axis=0)
    dev_sicura = np.where(dev == 0, 1.0, dev)           # avoids division by zero in constant columns
    return (X - media) / dev_sicura, media, dev

rng = np.random.default_rng(3)
X = rng.normal(loc=[10.0, -5.0, 3.0, 100.0, 0.0], scale=[2.0, 0.5, 1.0, 30.0, 1.0], size=(100, 5))
X[:, 2] = 7.0                                            # a constant column

Z, media, dev = standardizza(X)
print("shape of Z:", Z.shape, "| no NaN:", not np.isnan(Z).any())
print("column means ~ 0:", np.allclose(Z.mean(axis=0), 0))
print("std (non-constant columns) = 1:", np.allclose(Z.std(axis=0)[[0, 1, 3, 4]], 1))
print("constant column -> all 0:", np.allclose(Z[:, 2], 0), "| original std:", dev[2])

try:
    standardizza(np.array([1.0, 2.0, 3.0]))
except ValueError as e:
    print("ValueError:", e)

# the transformation can be inverted: Z * std + mean recovers X (where the std was not zero)
X_ricostruita = Z * np.where(dev == 0, 1.0, dev) + media
print("X rebuilt from Z:", np.allclose(X_ricostruita, X))

Output obtained:

shape of Z: (100, 5) | no NaN: True
column means ~ 0: True
std (non-constant columns) = 1: True
constant column -> all 0: True | original std: 0.0
ValueError: a matrix with 2 dimensions is required, got 1
X rebuilt from Z: True

How to read the solution. np.asarray(X, dtype=float) also accepts lists and converts them into an array of decimals without copying if it is already a suitable array. The X.ndim != 2 check rejects a vector at once: better a clear error than a result that is wrong because of broadcasting. X.mean(axis=0) and X.std(axis=0) have shape (5,), and (X - media) / dev_sicura uses the broadcasting of (5,) vectors over a (100, 5) matrix.

The delicate point is the constant column: its standard deviation is zero and 0/0 would give NaN over the whole column (with a warning). np.where(dev == 0, 1.0, dev) replaces the zeros with 1 only for the division: the constant column becomes all zeros, which is the honest choice ("it carries no information"), and the original dev is returned intact so the caller knows that column was constant.

The last check shows that the transformation is reversible: saving media and dev lets you bring new data to the same scale or go back to the original units (this is exactly what a scikit-learn scaler does, where media and dev are computed on the training set only and then applied, unchanged, to the test data). A note of style: the function does not modify X (it builds a new array) and returns everything needed to invert the operation; these are the two habits that make code on data predictable.

Challenge

• Easy. Compute the moving average with window 3 of an array of 10 numbers, without loops, using only slicing and sums.
• Medium. Given two matrices A (n, d) and B (m, d), compute with broadcasting the (n, m) matrix of the squared distance. Then prove the identity ‖a − b‖² = ‖a‖² − 2 a·b + ‖b‖² and rewrite it using a matrix product to avoid the intermediate (n, m, d) array.
• Hard. Measure with time.perf_counter the time of the sum of squares for n = 10³, 10⁴, 10⁵, 10⁶, 10⁷ in the loop version and in the vectorised one. From which n is vectorising worth it? Why does the advantage shrink for small n?
• Expert. With the formula ceil((stop − start) / step), search by trial for a decimal pair (stop, step) for which len(np.arange(0, stop, step)) differs from the exact count (use fractions.Fraction for the exact value) and explain why by looking at the binary representation of the quotient.
• Enterprise. Take a real table from your organisation with at least five numeric columns. Standardise it with NumPy and document, for each column, what a value z = 2 means.

The levels rise from practice with slices (Easy) to algebra (Medium), to measurement (Hard) and to the search for a numerical counterexample (Expert); the last has no code solution but an explanation for whoever will read your data. For Hard remember to repeat each measurement several times and to report an interval, as above.

What you can do now

You can read an array as a block of memory with dtype, shape and strides, and you can say when a selection is a view and when a copy. You can choose the axis of an aggregation remembering that the one indicated disappears, predict the shape of an operation with broadcasting by aligning shapes from the right, and recognise the (n, 1) against (n,) case that produces a wrong number without warning. The next step is Pandas, which puts labels and mixed types on top of arrays to work with real tables.

The same series also includes: 15 - scikit-learn: Estimators, Pipelines and Validation, 23 - Linear Algebra 1: Vectors, Matrices, Systems, 56 - Recurrent Neural Networks (RNN): Giving the Network a Memory (which refer back to this lesson).

Sources

Bibliographic reference not consulted in this revision: C. R. Harris et al., Array programming with NumPy, Nature 585, 2020. All the outputs were recomputed with Python 3.14.4 and NumPy 2.3.5 (each block run twice, identical output, except the timing block, reported as an interval).