3  NumPy and Numerical Computing

3.1 Why arrays matter

A learning algorithm repeats numerical calculations over many observations. A table might contain one row per walking trial and one column per measurement. We need a representation in which the row and column meanings remain clear while arithmetic applies consistently to the whole collection.

NumPy provides the ndarray, a numerical array with an explicit shape and element type. An array is more structured than a Python list: ordinary numeric arrays store their elements in a common numerical type and support elementwise arithmetic. We will first expose calculations with loops, then express the same operations using arrays. Shorter code is useful only when the reader can explain which values it combines.

3.1.1 Learning objectives and prerequisites

After this chapter, you should be able to:

  • Distinguish a scalar, a one-dimensional vector, and a two-dimensional matrix.
  • Inspect an array’s shape, dimensions, numerical type, and finite values.
  • Select observations and features without confusing rows and columns.
  • Translate sums and weighted sums into loops and NumPy operations.
  • Distinguish elementwise multiplication from matrix multiplication.
  • Recognize broadcasting mistakes, shared-memory views, and floating-point limits.

The prerequisite is Chapter 2. Mathematical ideas are introduced locally; Chapter 6 develops their geometry later. Each executable block is independent. import numpy as np gives the library the conventional short name np, as used in the NumPy beginner guide.

3.2 Step 1: Describe the dimensions before calculating

A scalar is a single number. A vector is an ordered collection of numbers. A matrix is a rectangular arrangement of numbers with rows and columns. In NumPy, the number of array axes is ndim; shape lists the length of every axis; size counts all elements. A matrix with three observations and two features has shape (3, 2) and contains six numbers.

In this book, a feature matrix usually has shape \((n,p)\), where \(n\) is the number of observations and \(p\) is the number of features per observation. We write it mathematically as \(X\in\mathbb{R}^{n\times p}\). The entry \(x_{ij}\) is feature \(j\) of observation \(i\). Mathematical indices often start at one, while Python indices start at zero. We will state this distinction whenever it matters.

# Step 1: Import the numerical array library.
import numpy as np

# Step 2: Create data structures with increasing numbers of axes.
# A zero-dimensional array holds a single trial count.
trial_count_scalar = np.array(3)
# A one-dimensional array holds three durations in seconds.
durations_seconds = np.array([10.0, 12.0, 8.0], dtype=np.float64)
# A two-dimensional array has rows = trials and columns = distance, duration.
# Column units are meters and seconds, respectively.
trial_measurements = np.array(
    [[12.0, 10.0], [12.0, 12.0], [12.0, 8.0]],
    dtype=np.float64,
)

# Step 3: Inspect the structural information before using the values.
print("Scalar shape:", trial_count_scalar.shape)
print("Duration vector shape:", durations_seconds.shape)
print("Measurement matrix shape:", trial_measurements.shape)
print("Matrix axes:", trial_measurements.ndim)
print("Matrix elements:", trial_measurements.size)
print("Numerical type:", trial_measurements.dtype)
Scalar shape: ()
Duration vector shape: (3,)
Measurement matrix shape: (3, 2)
Matrix axes: 2
Matrix elements: 6
Numerical type: float64

The shapes are (), (3,), and (3, 2). The comma in (3,) marks a one-item Python tuple. A one-dimensional vector has no separate row or column orientation. Shape (3,) therefore differs from both (3, 1) and (1, 3). This difference can change the result of arithmetic even when all three arrays contain three values. The NumPy quickstart describes axes, reductions, and shape manipulation.

Arrays do not store physical units automatically. Both columns above are floating-point numbers, although they represent different units. Maintain a column description alongside the data. Adding the distance column to the duration column would be syntactically valid and physically meaningless.

3.3 Step 2: Select the intended rows and columns

For a matrix, measurements[row_index, column_index] selects one entry. : means all positions along an axis. measurements[:, 1] selects the duration column from every row. A slice 0:2 includes positions zero and one and excludes position two.

  1. Extract one entry and one row to inspect the data convention.
  2. Extract the duration feature across all trials.
  3. Form a Boolean condition and select only rows that satisfy it.
# Step 1: Define three synthetic rows with documented column meanings.
import numpy as np

# Columns contain distance in meters and duration in seconds.
trial_measurements = np.array(
    [[12.0, 10.0], [12.0, 12.0], [12.0, 8.0]],
    dtype=np.float64,
)

# Step 2: Use zero-based indices to inspect individual values and rows.
# The first row and second column contain the first duration.
first_duration_seconds = trial_measurements[0, 1]
# Selecting one row returns its two feature values.
second_trial = trial_measurements[1, :]
# Selecting one column returns a length-three vector.
durations_seconds = trial_measurements[:, 1]

# Step 3: Build a row-selection mask using the duration feature.
# True indicates a duration strictly below 11 seconds.
short_trial_mask = durations_seconds < 11.0
# Apply the same mask to complete rows to preserve feature alignment.
short_trials = trial_measurements[short_trial_mask, :]

print("First duration (s):", first_duration_seconds)
print("Second trial [m, s]:", second_trial)
print("Row mask:", short_trial_mask)
print("Selected rows [m, s]:\n", short_trials)
First duration (s): 10.0
Second trial [m, s]: [12. 12.]
Row mask: [ True False  True]
Selected rows [m, s]:
 [[12. 10.]
 [12.  8.]]

The mask is [True, False, True]. Two rows remain. If targets or participant identifiers accompany the feature matrix, apply the same row selection to them. Filtering inputs and targets independently can leave equally long arrays whose rows describe different observations.

For multiple array conditions, use parentheses with & for elementwise “and” and | for elementwise “or”; for example, (durations_seconds > 0) & (durations_seconds < 11). Python’s scalar and asks for one truth value and cannot represent this row-by-row decision. See NumPy indexing.

3.4 Step 3: Understand whether a selection shares data

A view is another array object referring to the same underlying data. A copy has independent data. Basic slicing commonly creates a view; Boolean indexing creates a copy. This matters when a preprocessing step changes a selected subset: modifying a view can also modify the original data. Use .copy() when you intentionally need independence, and do not infer ownership merely because the selection has a new variable name. The behavior is documented in NumPy copies and views.

# Step 1: Create a small original signal in arbitrary sensor units.
import numpy as np

# The first value is intentionally distinct so mutation is easy to detect.
original_signal = np.array([2.0, 4.0, 6.0], dtype=np.float64)
# A basic slice shares the first two values with the original array.
signal_view = original_signal[0:2]
# An explicit copy preserves the first two values independently.
signal_copy = original_signal[0:2].copy()

# Step 2: Change the view and inspect both underlying relationships.
signal_view[0] = 99.0
print("Original after editing view:", original_signal)
print("Independent copy:", signal_copy)
print("View shares memory:", np.shares_memory(original_signal, signal_view))
Original after editing view: [99.  4.  6.]
Independent copy: [2. 4.]
View shares memory: True

The original now begins with 99.0; the copy still begins with 2.0. This is a memory relationship, not a random failure. In model development, preserve raw measurements so that transformed values can always be traced to their source.

3.5 Step 4: Reduce values along the correct axis

A reduction combines several values into fewer values, for example a sum, mean, or maximum. For \(n\) duration values \(t_1,\ldots,t_n\), their mean is

\[ \bar{t}=\frac{1}{n}\sum_{i=1}^{n}t_i. \]

Here \(t_i\) is measured in seconds, the sum also has units of seconds, and \(\bar{t}\) is a typical duration in seconds. In a two-dimensional feature matrix, axis=0 combines observations separately for each feature. The output has one value per column. axis=1 combines features separately for each observation. Whether that second calculation makes sense depends on the feature units.

# Step 1: Store two same-unit sensor channels for three time samples.
import numpy as np

# Rows are time samples; both columns are acceleration in m/s^2.
acceleration_samples = np.array(
    [[1.0, 3.0], [2.0, 4.0], [3.0, 5.0]],
    dtype=np.float64,
)

# Step 2: Average over rows to obtain one mean per sensor channel.
# axis=0 removes the time-sample axis from the result.
channel_mean_acceleration = np.mean(acceleration_samples, axis=0)

# Step 3: Average over columns to obtain one mean per time sample.
# This average is meaningful here because both columns have the same units.
sample_mean_acceleration = np.mean(acceleration_samples, axis=1)

print("Channel means (m/s^2):", channel_mean_acceleration)
print("Time-sample means (m/s^2):", sample_mean_acceleration)
Channel means (m/s^2): [2. 4.]
Time-sample means (m/s^2): [2. 3. 4.]

The channel means are [2, 4]; the sample means are [2, 3, 4]. Neither result is “the mean” without specifying the intended axis. Applying np.mean without an axis combines all entries; that is rarely the right choice for a feature table with different physical quantities.

3.6 Step 5: Move from a loop to elementwise arithmetic

Suppose each row contains distance \(d_i\) in meters and duration \(t_i\) in seconds. The average speed of trial \(i\) is \(v_i=d_i/t_i\). We can implement the same formula once per row or apply it to two aligned vectors.

Vectorization expresses a repeated numerical operation as an array operation. NumPy’s underlying routines perform the elementwise work without a Python-level loop for each value. It does not remove the arithmetic, and it does not automatically make every possible expression memory efficient. Understand the loop before replacing it.

# Step 1: Store aligned distances and durations for three synthetic trials.
import numpy as np

# Distance traveled in each trial, measured in meters; shape (3,).
distances_meters = np.array([12.0, 12.0, 12.0], dtype=np.float64)
# Matching trial durations in seconds; shape (3,).
durations_seconds = np.array([10.0, 12.0, 8.0], dtype=np.float64)
# Count the trials so the loop output can be allocated explicitly.
number_of_trials = len(distances_meters)

# Step 2: Check alignment and measurement validity before dividing.
assert distances_meters.shape == durations_seconds.shape
assert np.all(np.isfinite(distances_meters))
assert np.all(np.isfinite(durations_seconds))
assert np.all(distances_meters >= 0)
assert np.all(durations_seconds > 0)

# Step 3: Calculate speed with an explicit educational loop.
# One output position is needed for each trial.
loop_speeds_meters_per_second = np.zeros(number_of_trials, dtype=np.float64)
for trial_index in range(number_of_trials):
    # Select both measurements from the same trial before dividing.
    distance_meters = distances_meters[trial_index]
    # Each duration is positive in this controlled example.
    duration_seconds = durations_seconds[trial_index]
    loop_speeds_meters_per_second[trial_index] = distance_meters / duration_seconds

# Step 4: Perform the identical elementwise division on aligned vectors.
# Each output is distances_meters[i] divided by durations_seconds[i].
array_speeds_meters_per_second = distances_meters / durations_seconds

# Step 5: Verify agreement rather than assuming shorter code is correct.
assert np.allclose(loop_speeds_meters_per_second, array_speeds_meters_per_second)
print("Loop speeds (m/s):", loop_speeds_meters_per_second)
print("Array speeds (m/s):", array_speeds_meters_per_second)
Loop speeds (m/s): [1.2 1.  1.5]
Array speeds (m/s): [1.2 1.  1.5]

The speeds are [1.2, 1.0, 1.5] m/s. This is an agreement check, not a timing experiment. For \(n\) trials, both methods perform order-\(n\) arithmetic and store \(n\) outputs. Array operations can reduce Python overhead for large numeric collections, but performance also depends on allocation, data layout, and the operation. Do not claim a universal speedup from this tiny example.

3.7 Step 6: Distinguish products and check dimensions

Elementwise multiplication first_array * second_array multiplies corresponding entries after any allowed broadcasting. Matrix multiplication first_array @ second_array forms weighted sums according to matrix-product rules. They solve different tasks.

For a feature matrix \(X\) with shape \((n,p)\) and coefficient vector \(\boldsymbol{w}\) with shape \((p,)\), the product \(X\boldsymbol{w}\) has shape \((n,)\). Its \(i\)th entry is

\[ \widehat{y}_i=\sum_{j=1}^{p}x_{ij}w_j. \]

Each observation combines \(p\) feature contributions into one output. This is the basic calculation behind a linear prediction. The following coefficients are chosen for arithmetic practice; they have not been fitted to observations.

# Step 1: Define a dimensionless feature matrix and chosen coefficients.
import numpy as np

# Three observations have two dimensionless features each; shape (3, 2).
feature_matrix = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])
# Each coefficient weights one feature; shape (2,).
feature_weights = np.array([0.5, 2.0])

# Step 2: Calculate the first weighted sum by hand-sized operations.
# Feature zero contributes 1.0 * 0.5 to observation zero.
first_contribution = feature_matrix[0, 0] * feature_weights[0]
# Feature one contributes 2.0 * 2.0 to the same observation.
second_contribution = feature_matrix[0, 1] * feature_weights[1]
# Summing contributions gives the first predicted output.
first_prediction = first_contribution + second_contribution

# Step 3: Use matrix multiplication to calculate all three weighted sums.
# The inner dimensions match: (3, 2) @ (2,) produces (3,).
predicted_outputs = feature_matrix @ feature_weights
print("First weighted sum:", first_prediction)
print("All weighted sums:", predicted_outputs)
print("Prediction shape:", predicted_outputs.shape)
First weighted sum: 4.5
All weighted sums: [ 4.5  9.5 14.5]
Prediction shape: (3,)

The first result is \(1(0.5)+2(2)=4.5\). The three results are [4.5, 9.5, 14.5]. feature_matrix * feature_weights would produce a (3, 2) array of separate contributions, not their row sums. See the numpy.matmul documentation for the @ operation. Later, Chapter 14 explains how coefficients are learned and how an intercept is included.

3.8 Step 7: Use broadcasting deliberately

Broadcasting allows operations between certain different shapes. Compare dimensions starting from the right: they must be equal or one of them must be one; absent leading dimensions behave as one. A scalar can therefore be added to every array entry. A (p,) vector can be subtracted from each row of an (n,p) matrix. These are useful rules, but an unintended compatible shape can create an incorrect result without raising an error. See NumPy broadcasting.

A particularly dangerous case is subtracting predictions shaped (n, 1) from targets shaped (n,). NumPy creates an (n,n) array of pairwise differences, even though the intention was \(n\) matched errors.

# Step 1: Store two matched targets and predictions in different shapes.
import numpy as np

# Actual outputs for two observations; shape (2,).
observed_targets = np.array([1.0, 3.0])
# Predictions stored as a column matrix; shape (2, 1).
column_predictions = np.array([[0.5], [2.5]])

# Step 2: Inspect the unintended pairwise result.
# Broadcasting compares every target with each prediction row.
pairwise_differences = observed_targets - column_predictions
print("Unintended difference shape:", pairwise_differences.shape)
print(pairwise_differences)

# Step 3: Make the intended one-prediction-per-observation shape explicit.
# reshape changes the organization of two entries, not their meaning or order.
predicted_targets = column_predictions.reshape(2)
# Check matching structure before computing matched errors.
assert observed_targets.shape == predicted_targets.shape
# Each difference now belongs to the same observation in both arrays.
matched_errors = observed_targets - predicted_targets
print("Matched errors:", matched_errors)
Unintended difference shape: (2, 2)
[[ 0.5  2.5]
 [-1.5  0.5]]
Matched errors: [0.5 0.5]

The intended errors are [0.5, 0.5]. reshape is appropriate here because each column entry already corresponds to one target in the correct order. It cannot repair mismatched observation identities. Transposing also does not make a one-dimensional (n,) vector into a column; use an explicitly understood reshape such as .reshape(n, 1) when a column is required.

3.9 Step 8: Account for precision and invalid numbers

Floating-point numbers have finite precision. Values such as one tenth generally do not have exact finite binary representations. Thus a comparison of computed quantities using == can fail even when their mathematical values agree. np.allclose checks whether differences are small relative to chosen absolute and relative tolerances. Choose tolerances according to the scale of the calculation, not merely because a check passes. See numpy.allclose.

# Step 1: Compare a short decimal calculation with its intended value.
import numpy as np

# Binary floating-point calculation to inspect.
computed_value = np.float64(0.1) + np.float64(0.2)
# Intended mathematical result, also stored approximately.
reference_value = np.float64(0.3)
# Allow only a tiny absolute rounding difference for this unit-scale example.
absolute_tolerance = 1e-12

print("Stored sum:", repr(computed_value))
print("Exact stored equality:", computed_value == reference_value)
print(
    "Within chosen tolerance:",
    np.allclose(computed_value, reference_value, rtol=0.0, atol=absolute_tolerance),
)

# Step 2: Identify nonfinite values before using measurements in a model.
# NaN represents a missing or undefined numeric entry; infinity is also invalid here.
candidate_measurements = np.array([1.0, np.nan, np.inf])
# One Boolean value marks the validity of each numeric entry.
finite_value_mask = np.isfinite(candidate_measurements)
print("Finite entries:", finite_value_mask)
Stored sum: np.float64(0.30000000000000004)
Exact stored equality: False
Within chosen tolerance: True
Finite entries: [ True False False]

float64 uses more storage and offers greater precision than float32, but no fixed-width type gives unlimited accuracy. Fixed-width integer arrays can overflow; assigning a decimal into an integer array can discard its fractional part. Inspect dtype and select an appropriate representation before calculating. The NumPy data types guide explains these limits.

Do not silently replace every nonfinite value by zero. Determine why it exists. A failed sensor channel, an absent record, and division by zero require different investigation. More decimal places in a printout also do not create more measurement accuracy.

3.10 Step 9: Generate repeatable synthetic values

Synthetic data let us practice an operation without downloading a dataset or implying a real scientific result. A random generator produces a reproducible sequence when initialized and used consistently. A seed is a label for that initial state, not a quality score. Fixing it helps compare code changes; it does not show that a conclusion is robust to alternative samples.

# Step 1: Create a local random generator for a small synthetic example.
import numpy as np

# A fixed seed makes this example repeatable in the same software setting.
random_seed = 42
# Keep generator state local rather than changing a global random sequence.
random_generator = np.random.default_rng(random_seed)

# Step 2: Simulate four dimensionless noise values around zero.
# loc is the distribution mean; scale is its standard deviation.
synthetic_noise = random_generator.normal(loc=0.0, scale=0.1, size=4)
print("Synthetic noise:", synthetic_noise)
print("All finite:", np.all(np.isfinite(synthetic_noise)))
Synthetic noise: [ 0.03047171 -0.10399841  0.07504512  0.09405647]
All finite: True

Calling the same generator again advances its state and produces new draws. Recreating it with the same seed repeats the initial sequence in the same compatible setting. Record software versions when exact reproduction matters; NumPy documents the generator interface and reproducibility considerations in its random generator reference.

3.11 Exercises and worked answers

3.11.1 1. Explain the shape

There are 20 trials and four measurements per trial. What shape should the feature matrix have, and what shape results from taking its mean along axis=0?

Answer: The matrix has shape (20, 4). Reducing over observations produces shape (4,), one mean for each feature. Check units separately for each feature.

3.11.2 2. Identify the operation

For a matrix [[1, 2], [3, 4]] and weights [2, 1], calculate the elementwise contributions and matrix-product predictions.

Answer: Elementwise multiplication gives [[2, 2], [6, 4]]. Summing each row, or using @, gives [4, 10]. An array of contributions and a vector of predictions have different meanings and shapes.

3.11.3 3. Why does subtraction produce nine values?

Targets have shape (3,) and predictions have shape (3, 1). What shape does their difference have?

Answer: (3, 3). Broadcasting compares an implicit (1, 3) target row with a (3, 1) prediction column. If there is one prediction per target, deliberately use matching shapes and confirm their observation order.

3.11.4 4. Does a view protect the original?

You take selected_values = original_values[1:3] and change selected_values[0]. What should you expect?

Answer: With a basic slice, the selected array is a view and changing it changes the corresponding original entry. Use .copy() if the selection must be independent.

3.11.5 5. Implement mean squared error

Observed outputs are [1, 2, 3] and predictions are [1, 3, 2]. Calculate the mean squared error using explicit steps.

Answer: Matched observed-minus-predicted errors are [0, -1, 1]. Squared errors are [0, 1, 1]. Their sum is \(2\) and the mean is \(2/3\). In NumPy, first verify identical shapes, then store errors, squared_errors = errors ** 2, and mean_squared_error = np.mean(squared_errors). If outputs have units of meters, squared error has units of square meters.

3.11.6 6. Interview question: does vectorization change complexity?

Does replacing an explicit loop over \(n\) values by one elementwise array division change linear time into constant time?

Answer: No. It still performs one division per entry, so arithmetic work scales with \(n\). It changes where the loop executes and may reduce interpreter overhead. Broadcasting can even create larger temporary arrays if shapes are chosen poorly.

3.12 One-page revision summary

Question Habit
What does each axis mean? State observations, features, channels, or time explicitly
Are dimensions correct? Inspect shape, ndim, and operation outputs
Are rows still matched? Apply one selection to every aligned variable
Can a selection change the original? Understand views; copy intentionally
What is being reduced? Name the intended axis and output units
Is this a weighted sum? Use @; distinguish it from *
Could broadcasting hide an error? Check matching target/prediction shapes
Are numeric values trustworthy? Inspect dtype, finiteness, units, and tolerances
Is a random example repeatable? Use a local seeded generator and record the setup

Arrays provide numerical structure, but they do not preserve a rich table schema. Continue to Chapter 4 to keep column names, identifiers, missingness, and plots connected to the measurements.