Numerical Computing with NumPy

Introduction

Numerical computing is the use of algorithms and data structures to represent, transform, and analyse numbers efficiently. In data science, numerical computing supports data cleaning, statistical analysis, machine learning, simulation, optimisation, and scientific modelling. Python’s built-in lists are flexible, but they are not designed for large-scale numerical operations. NumPy (Numerical Python) provides multidimensional arrays and highly optimised operations implemented largely in compiled code.

The central idea is to express an operation on an entire array rather than repeatedly processing individual elements in Python. This practice, called vectorization, usually produces code that is shorter, clearer, and faster.

Learning Objectives

By the end of this chapter, you should be able to:

• create and inspect one- and multidimensional NumPy arrays;
• select and modify data using indexing and slicing;
• apply vectorized operations and understand broadcasting;
• calculate descriptive aggregates along selected axes;
• perform basic matrix and linear algebra operations;
• generate reproducible random samples and simulations;
• compare NumPy computation with native Python approaches; and
• identify numerical, performance, and data-governance limitations.

1. The NumPy Array Model

NumPy’s main structure is the ndarray. An array contains elements of one data type, arranged over zero or more dimensions. Important properties include:

A scalar has zero dimensions, a vector commonly has one dimension, a table has two dimensions, and a data cube may have three or more.

import numpy as np

scores = np.array([72, 81, 69, 94])
print(scores.shape)   # (4,)
print(scores.ndim)    # 1
print(scores.dtype)

marks = np.array([[72, 81, 69],
                  [94, 77, 88]])
print(marks.shape)    # (2, 3)

The two-dimensional array marks has two rows and three columns. NumPy uses zero-based indexing, so marks[0, 1] is 81.

Creating Arrays

Common constructors include:

np.zeros((2, 3))              # two rows, three columns of zeros
np.ones(4)                    # four ones
np.full((2, 2), 7)            # every element is 7
np.arange(0, 10, 2)           # 0, 2, 4, 6, 8
np.linspace(0, 1, 5)          # five equally spaced values
np.eye(3)                     # 3 × 3 identity matrix

arange is convenient for integer-like steps but can produce unexpected endpoint behaviour with floating-point steps. linspace is preferable when the number of values, rather than the step size, is the primary requirement.

Data Types and Conversion

NumPy selects a common data type when constructing an array. Explicit conversion is often useful:

values = np.array([1, 2, 3], dtype=np.float64)
integers = values.astype(np.int32)

Converting floating-point values to integers discards the fractional part; it does not round to the nearest integer. Numeric overflow is also possible with fixed-width integer types, so the chosen type must be suitable for the expected range.

2. Indexing, Slicing, and Reshaping

For a two-dimensional array, use array[row, column].

A = np.array([[10, 20, 30],
              [40, 50, 60],
              [70, 80, 90]])

A[1, 2]       # 60
A[0, :]       # first row
A[:, 1]       # second column
A[:2, :2]     # upper-left 2 × 2 region
A[::2, ::2]   # every second row and column

A slice normally creates a view rather than an independent copy. Therefore, changing a slice may change the original array:

part = A[:2, :2]
part[0, 0] = 999

Use A[:2, :2].copy() when independent storage is required. This distinction is important in data pipelines because an apparently local modification can otherwise alter a source dataset.

Boolean indexing selects elements satisfying a condition:

temperatures = np.array([18.5, 21.0, 29.2, 16.8])
hot = temperatures[temperatures > 20]
temperatures[temperatures < 18] = 18

Integer or fancy indexing selects specified positions:

data = np.array([11, 22, 33, 44, 55])
selected = data[[0, 3, 4]]

Reshaping changes the organisation of elements without changing their total number. If an array has m × n elements, a valid reshape must also contain m × n positions.

x = np.arange(12)
B = x.reshape(3, 4)
C = B.T                 # transpose: shape (4, 3)
flat = B.ravel()        # flattened view where possible

The -1 placeholder allows NumPy to infer one dimension: x.reshape(3, -1) produces shape (3, 4).

3. Vectorization and Performance

Consider adding two sequences. A native Python approach is:

result = [a + b for a, b in zip(list_a, list_b)]

The NumPy version is:

result = array_a + array_b

The second expression applies compiled element-wise addition to the entire array. This is vectorization. It reduces Python interpreter overhead and often uses contiguous memory more efficiently.

For an array x with n elements, element-wise addition with another compatible array performs n additions. The mathematical operation is:

Here xᵢ, yᵢ, and zᵢ are corresponding elements, and both arrays must have compatible shapes and numeric types.

A simple comparison is:

import time

n = 1_000_000
python_values = list(range(n))
numpy_values = np.arange(n)

start = time.perf_counter()
python_result = [v * 2 for v in python_values]
python_time = time.perf_counter() - start

start = time.perf_counter()
numpy_result = numpy_values * 2
numpy_time = time.perf_counter() - start

The exact timings depend on hardware, data type, memory, and program environment. NumPy is not automatically faster for every small operation; array creation and conversion costs can dominate tiny workloads. Vectorization also does not eliminate algorithmic complexity: an inefficient vectorized algorithm may still be inefficient.

4. Element-Wise Operations and Broadcasting

Arithmetic operators on arrays are generally element-wise:

x = np.array([1, 4, 9])
y = np.array([2, 3, 5])
x + y
x * y
x ** 2
np.sqrt(x)
np.log(x)

The expression x * y is not matrix multiplication. It multiplies corresponding elements. Matrix multiplication uses @ or np.matmul.

Broadcasting allows arrays with different but compatible shapes to participate in element-wise operations. NumPy compares shapes from the trailing dimension backwards. Two dimensions are compatible when they are equal or one of them is 1; missing leading dimensions are treated as 1.

For example, a matrix with shape (3, 2) can be added to a vector with shape (2,). The vector is conceptually used for every row:

M = np.array([[1, 2], [3, 4], [5, 6]])
offset = np.array([10, 100])
M + offset

The result is [[11, 102], [13, 104], [15, 106]]. No three-row copy of offset is necessarily created, which saves memory.

A column vector with shape (3, 1) can combine with a row vector with shape (2,) to produce a (3, 2) result:

rows = np.array([[1], [2], [3]])
columns = np.array([10, 20])
rows + columns

Broadcasting can also cause accidental results. A shape (3,) array and a shape (3, 1) array produce a (3, 3) array, not a pairwise one-dimensional operation. Inspect shapes before calculations.

5. Aggregation and Axes

Aggregations reduce multiple values to summaries:

marks = np.array([[72, 81, 69],
                  [94, 77, 88]])

marks.sum()
marks.mean()
marks.min()
marks.max()
marks.std()

For observations arranged in rows and variables in columns, axis=0 reduces down the rows and returns one value per column. axis=1 reduces across columns and returns one value per row.

column_means = marks.mean(axis=0)
student_means = marks.mean(axis=1)

The arithmetic mean of n values x₁, …, xₙ is:

The mean is appropriate for quantitative measurements when each observation has equal importance and extreme values are not misleading. NumPy’s mean uses floating-point arithmetic and can be sensitive to precision for some data types.

The standard deviation describes dispersion. NumPy’s default is the population form:

Here μ is the population mean and n is the number of observations. For a sample estimate, use np.std(x, ddof=1), which divides by n − 1 when n > 1. Missing values require special treatment: np.mean propagates NaN, whereas np.nanmean ignores NaN values, provided that at least one valid value exists.

The keepdims=True option preserves reduced dimensions and can simplify later broadcasting:

row_means = marks.mean(axis=1, keepdims=True)
centred = marks - row_means

6. Linear Algebra

A matrix is a rectangular arrangement of numbers. If A has shape (m, n) and B has shape (n, p), their matrix product C = AB has shape `(m, p)), with entries:

The inner dimensions must match. In NumPy:

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = A @ B

Here C is [[19, 22], [43, 50]]. In contrast, A * B is element-wise multiplication.

The dot product of vectors x and y of equal length n is:

x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
x @ y       # 32

A linear system asks for x in:

If A is square and nonsingular, a unique solution exists. Numerically, use np.linalg.solve(A, b) rather than explicitly calculating A⁻¹ and multiplying, because solving is generally more stable and efficient.

A = np.array([[3., 1.], [1., 2.]])
b = np.array([9., 8.])
solution = np.linalg.solve(A, b)

The solution is x = 2 and y = 3. The determinant, inverse, eigenvalues, and norm are available through np.linalg, but these tools must be interpreted carefully. A matrix can be invertible in theory yet poorly conditioned in finite-precision arithmetic. Residual checking is good practice:

np.allclose(A @ solution, b)

7. Random Numbers and Simulation

Random simulation is useful for risk analysis, sampling, testing, and Monte Carlo estimation. NumPy’s modern interface uses a Generator:

rng = np.random.default_rng(42)
rolls = rng.integers(1, 7, size=10)
normal_values = rng.normal(loc=0, scale=1, size=1000)
uniform_values = rng.uniform(0, 1, size=1000)

The seed makes a pseudorandom sequence reproducible on the same software implementation. A seed does not make a process genuinely random or secure. Do not use ordinary simulation generators for passwords, tokens, or security keys.

For independent Bernoulli trials with success probability p, the expected number of successes in n trials is:

and the variance is:

These formulas assume independent trials with the same success probability p. Simulation results fluctuate around theoretical values.

Worked Simulation: Estimating π

Draw points uniformly in the unit square [0, 1] × [0, 1]. Points inside the quarter-circle satisfy:

The quarter-circle area is π/4, while the square area is 1. Therefore, if q is the proportion of points inside the quarter-circle, π can be estimated by:

rng = np.random.default_rng(7)
N = 200_000
points = rng.random((N, 2))
inside = (points[:, 0]**2 + points[:, 1]**2) <= 1
pi_estimate = 4 * inside.mean()

The estimate will not exactly equal π. Increasing N generally reduces random error, but it increases computation. A simulation should report its seed, sample size, assumptions, and uncertainty where decisions depend on the result.

8. Applied Case: Monitoring a Manufacturing Process

Background

This illustrative case concerns a factory measuring the temperature of three production lines over four hourly readings. The quality team wants line averages, deviations from each line’s mean, and a warning for readings above 80 °C.

temps = np.array([
    [76.0, 78.0, 79.0, 77.0],
    [81.0, 83.0, 82.0, 84.0],
    [74.0, 75.0, 73.0, 76.0]
])

Rows represent lines and columns represent hours. The row means are:

• Line 1: (76 + 78 + 79 + 77) ÷ 4 = 77.5 °C
• Line 2: (81 + 83 + 82 + 84) ÷ 4 = 82.5 °C
• Line 3: (74 + 75 + 73 + 76) ÷ 4 = 74.5 °C

line_means = temps.mean(axis=1)
line_deviation = temps - line_means[:, None]
alerts = temps > 80

line_means[:, None] changes the shape from (3,) to (3, 1), allowing broadcasting across the four hourly columns. The alert matrix identifies all readings above 80 °C; Line 2 triggers four alerts, while Lines 1 and 3 trigger none.

Analysis and Outcome

The mean indicates typical temperature, while the deviations show whether individual readings are unusually high or low relative to their line. However, an operational decision should not rely on the mean alone. Engineers may also require control limits, sensor calibration records, time ordering, and knowledge of safe operating ranges.

Learning Points

The case demonstrates a complete NumPy workflow: represent structured measurements as an array, aggregate along the correct axis, use broadcasting for centering, and apply Boolean indexing for threshold detection. It also shows that technical results require domain interpretation and validated measurement procedures.

9. Best Practices, Limitations, and Responsible Use

Use descriptive variable names, document array shapes, and validate inputs before computation. Prefer np.random.default_rng() for new work. Use np.isfinite to detect invalid numeric values and select np.nan* functions only when ignoring missing values is justified. Use np.allclose for floating-point comparisons rather than exact equality when small rounding differences are expected.

Floating-point numbers approximate real numbers. For example, decimal values such as 0.1 may not have an exact binary representation. Results can therefore contain small errors. Unit consistency is also essential: do not combine temperatures in Celsius with quantities requiring absolute temperature in kelvins without applying the correct transformation.

NumPy arrays are homogeneous, unlike Python lists, so mixed textual and numeric data should generally be handled with a tabular tool such as pandas. NumPy also does not provide complete data validation, provenance, access control, or audit trails. In professional settings, retain source information, transformation steps, version details, and random seeds where reproducibility matters.

Numerical software can support AI and machine learning, but a fast calculation is not automatically a valid decision. Biased samples, inappropriate aggregation, leakage of sensitive data, and unexamined thresholds can produce harmful outcomes. Personal or operational data should be minimised, protected, and processed according to applicable privacy and organisational requirements.

Chapter Summary

NumPy provides homogeneous, multidimensional arrays for efficient numerical computing. Array shape, dimensionality, data type, and memory behaviour determine how operations work. Indexing and slicing support precise selection, while views and copies must be distinguished to avoid unintended mutation. Vectorization replaces many explicit Python loops with array-level operations. Broadcasting permits compatible arrays of different shapes to interact, but shape errors can produce invalid or surprising results. Aggregations such as sums, means, and standard deviations become especially powerful when the axis is selected correctly. Matrix multiplication, dot products, and linear-system solvers support foundational linear algebra. NumPy’s random generators enable reproducible simulation, subject to clearly stated assumptions and responsible use. Finally, performance gains must be balanced with numerical precision, validation, interpretation, privacy, and domain requirements.

Key Terms

Array: A structured collection of elements, usually having a common data type.

Axis: A dimension along which an array operation is performed.

Broadcasting: NumPy’s rule-based alignment of compatible array shapes.

Boolean indexing: Selecting or modifying elements using True/False conditions.

Copy: An independent array with separate storage.

Data type (dtype): The numerical or logical representation used for array elements.

Dimension: The number of axes in an array.

Fancy indexing: Selection using arrays or lists of integer positions.

Generator: NumPy object that produces pseudorandom values.

Matrix: A two-dimensional rectangular array of values.

Monte Carlo simulation: Numerical estimation using repeated random sampling.

Shape: The size of an array along each axis.

Slicing: Selecting a range or pattern of positions.

Standard deviation: A measure of dispersion around a mean.

Transpose: An operation that exchanges rows and columns.

Vectorization: Applying an operation to whole arrays rather than explicit Python loops.

View: An array representation sharing storage with another array.

ndarray: NumPy’s multidimensional array data structure.

Review Questions

Short Answer Questions

1. What information is provided by an array’s shape and dtype?
2. Distinguish between element-wise multiplication and matrix multiplication in NumPy.
3. Why can modifying a slice modify the original array?
4. What does axis=0 generally mean when aggregating a two-dimensional array?
5. Why is a random seed useful in simulation?

Long Answer Questions

1. Explain vectorization and compare it with a Python list comprehension, including situations where NumPy may not provide a substantial advantage.
2. Describe NumPy broadcasting rules and illustrate them with arrays of shapes (3, 2) and (2,).
3. Explain the difference between population and sample standard deviation and how ddof=1 changes NumPy’s calculation.
4. Discuss why np.linalg.solve(A, b) is usually preferable to explicitly computing A⁻¹.
5. Explain the main numerical and governance risks that should be considered when using NumPy in a professional data-science workflow.

Analytical / Application Questions

1. A dataset has shape (500, 12), representing 500 observations and 12 variables. Write NumPy expressions to calculate one mean per variable, centre every variable, and identify observations containing at least one value above a chosen threshold.
2. Design a simulation using NumPy to estimate the probability that at least 6 of 10 independent trials succeed when the success probability is 0.4. State how you would assess reproducibility and simulation error.
3. A developer uses A * B to solve a matrix model and obtains plausible-looking values. Explain how you would diagnose the error, replace the operation correctly, and validate the result.

Multiple Choice Questions with Answers

Question 1

Which NumPy attribute gives the number of dimensions of an array?

A. size

B. shape

C. ndim

D. dtype

Correct Answer: C. ndim

Question 2

What is the result of np.array([1, 2, 3]) * 2?

A. The matrix product with 2

B. An array containing [2, 4, 6]

C. An error because 2 is a scalar

D. The array [1, 2, 3, 2]

Correct Answer: B. An array containing [2, 4, 6]

Question 3

For an array with shape (3, 2), which shape can be broadcast for element-wise addition?

A. (3,)

B. (2,)

C. (4,)

D. (3, 3)

Correct Answer: B. (2,)

Question 4

If M has shape (4, 3), what does M.mean(axis=0) return?

A. One mean for each of four rows

B. One mean for all elements

C. One mean for each of three columns

D. A matrix with shape (4, 3)

Correct Answer: C. One mean for each of three columns

Question 5

What does A @ B represent for compatible two-dimensional arrays?

A. Element-wise multiplication

B. Matrix multiplication

C. Boolean indexing

D. Array concatenation

Correct Answer: B. Matrix multiplication

Question 6

Why is np.linalg.solve(A, b) generally preferred to calculating np.linalg.inv(A) @ b?

A. It changes all values to integers

B. It usually provides a more direct and numerically preferable solution

C. It works only when A is nonsquare

D. It avoids the need for compatible dimensions

Correct Answer: B. It usually provides a more direct and numerically preferable solution

Question 7

Which approach creates a modern NumPy pseudorandom generator with a fixed seed?

A. np.random.default_rng(42)

B. np.seed(42)

C. np.array.random(42)

D. np.random.fixed(42)

Correct Answer: A. np.random.default_rng(42)

Question 8

A slice of a NumPy array may be dangerous to modify because it often is a:

A. Scalar

B. View sharing storage

C. Random copy

D. String object

Correct Answer: B. View sharing storage

Question 9

What condition identifies points inside the unit quarter-circle in the simulation described?

A. x + y <= 1

B. x2 + y2 <= 1

C. x2 - y2 <= 1

D. x * y == 1

Correct Answer: B. x2 + y2 <= 1

Question 10

If two vectors are [1, 2, 3] and [4, 5, 6], what is their dot product?

A. 21

B. 30

C. 32

D. 36

Correct Answer: C. 32