Numerical Computing with NumPy

Introduction

Numerical computing is the use of algorithms and data structures to represent, transform, and analyse numerical information. It supports scientific simulation, statistics, machine learning, computer vision, finance, engineering, and business analytics. Python is productive for numerical work because it combines readable syntax with a large ecosystem of libraries. The central library for dense numerical computation is NumPy, whose principal object is the multidimensional array, or ndarray.

A Python list is a general-purpose container. A NumPy array is a structured block of values with one data type, a fixed shape, and operations designed to process many values efficiently. This distinction enables vectorization: expressing an operation on an entire array rather than writing an explicit Python loop.

This chapter develops array construction and manipulation, indexing, broadcasting, vectorized operations, aggregation, linear algebra, random sampling, performance comparison, and responsible use of numerical results.

Learning Objectives

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

• explain the structure, shape, dimensions, data type, and memory properties of an ndarray;
• create, reshape, index, slice, combine, and copy multidimensional arrays;
• apply elementwise operations, universal functions, broadcasting, and vectorization;
• calculate descriptive aggregates and use axis semantics correctly;
• solve common linear-algebra tasks with NumPy;
• generate reproducible random samples and distinguish random generation from statistical inference;
• compare vectorized NumPy processing with conventional Python loops;
• identify numerical, performance, and governance risks in data-science applications.

1. The NumPy Array Model

An ndarray has several important properties.

For an array with shape (r, c), r is the number of rows and c is the number of columns. Its total number of elements is

where N is the element count. More generally, for shape (d₁, d₂, …, d_k),

The formula assumes that every axis has a finite, positive length. In practice, NumPy also permits zero-length dimensions in some operations.

import numpy as np

x = np.array([[10, 20, 30], [40, 50, 60]], dtype=np.int64)
print(x.ndim)   # 2
print(x.shape)  # (2, 3)
print(x.size)   # 6
print(x.dtype)  # int64

A two-dimensional array can be viewed as a table, but NumPy arrays may have any number of dimensions. For example, an RGB image commonly has shape (height, width, 3), while a batch of colour images may have shape (batch, height, width, 3).

Constructing arrays

Common constructors include:

np.array([1, 2, 3])
np.zeros((2, 4))
np.ones((3, 3), dtype=np.int32)
np.full((2, 3), 7)
np.eye(4)
np.arange(0, 10, 2)
np.linspace(0, 1, 5)

arange(start, stop, step) is useful for a step-based sequence, but floating-point steps can accumulate rounding effects. linspace(start, stop, num) is often preferable when the required number of equally spaced points is known. By default, linspace includes both endpoints.

The default floating-point type is usually float64, while integer arrays commonly use a platform-dependent integer type. Data type selection affects memory consumption, range, precision, and compatibility with downstream libraries. Converting a value to a narrower type can cause overflow or loss of precision.

2. Indexing, Slicing, and Array Shape

NumPy uses zero-based indexing. For a two-dimensional array a, a[i, j] selects row i and column j. A slice has the form start:stop:step, where stop is exclusive.

a = np.array([[4, 8, 12],
              [5, 10, 15],
              [6, 12, 18]])

first_row = a[0, :]
last_column = a[:, -1]
centre = a[1, 1]
block = a[:2, 1:]

Basic slices generally produce views rather than independent copies. A view shares memory with the original array; modifying one may modify the other.

b = a[:, :2]
b[0, 0] = 999       # may also change a[0, 0]
c = a[:, :2].copy() # independent array

This behaviour is efficient but requires care. Use .copy() when an independent working array is required.

Boolean and fancy indexing

A Boolean mask selects elements for which a condition is true.

scores = np.array([42, 67, 81, 55, 93])
passed = scores[scores >= 60]
scores[scores < 50] = 0

The mask must be compatible with the indexed array. For compound conditions, use elementwise operators &, |, and ~, with parentheses around each comparison.

mask = (scores >= 60) & (scores <= 90)
selected = scores[mask]

Fancy indexing uses integer arrays or lists of indices. Unlike basic slicing, it generally returns a copy.

rows = np.array([0, 2])
chosen = a[rows, :]

Reshaping and transposition

reshape changes the shape without changing the element order, provided the total number of elements remains constant. Thus, a reshape from (2, 6) to (3, 4) is valid because both contain 12 elements.

v = np.arange(12)
m = v.reshape(3, 4)
column = v.reshape(-1, 1)
transposed = m.T

The value -1 asks NumPy to infer one dimension. Only one dimension may be inferred, and the known dimensions must divide the element count exactly. ravel() attempts to return a flattened view, whereas flatten() returns a copy.

concatenate, stack, vstack, and hstack combine arrays, but they have different shape requirements. concatenate joins along an existing axis; stack creates a new axis.

p = np.array([1, 2])
q = np.array([3, 4])
np.concatenate([p, q]) # shape (4,)
np.stack([p, q])       # shape (2, 2)

3. Elementwise Operations, Universal Functions, and Broadcasting

Arithmetic on arrays is normally elementwise.

x = np.array([1, 2, 3])
y = np.array([10, 20, 30])
x + y       # [11, 22, 33]
x * y       # [10, 40, 90]
x ** 2      # [1, 4, 9]

The symbol * denotes elementwise multiplication, not matrix multiplication. Matrix multiplication uses @ or np.matmul.

Universal functions, or ufuncs, apply compiled numerical operations elementwise:

np.sqrt(x)
np.exp(x)
np.log(x)
np.sin(x)
np.maximum(x, 2)

For logarithms and division, mathematical domain restrictions apply. For example, real-valued log(z) requires z > 0; division by zero may generate warnings, infinities, or nan, depending on the data type and operation.

Broadcasting

Broadcasting allows arrays with compatible shapes to participate in one operation without manually copying data. NumPy compares dimensions from right to left. Two dimensions are compatible when they are equal or one of them is 1; missing leading dimensions are treated as 1.

For example, an array of shape (3, 1) and one of shape (1, 4) broadcast to (3, 4). Conceptually, each value in the first array is used across a row and each value in the second across a column, although NumPy generally avoids physically replicating the data.

temperatures = np.array([[18.0], [21.0], [25.0]])  # (3, 1)
offset = np.array([0.5, 1.0, 1.5, 2.0])           # (4,)
result = temperatures + offset                      # (3, 4)

Broadcasting is powerful but can create unintended results if shapes are misunderstood. Inspect shape explicitly and use np.broadcast_shapes when validating an interface.

A common statistical transformation is standardization. For observations xᵢ, mean μ, and standard deviation σ, the standardized value is

This formula assumes σ > 0 and is applied to a defined population or sample convention. In NumPy, axis and ddof determine how μ and σ are calculated.

4. Aggregation and Axis Semantics

Aggregation reduces several values to one or fewer values. Important functions include sum, mean, median, min, max, std, var, argmin, and argmax.

For values x₁, x₂, …, xₙ, the arithmetic mean is

where n is the number of observations. The population variance is

where μ is the population mean. For an unbiased estimate of a population variance from a sample, NumPy uses ddof=1, giving a denominator of n − 1, provided that n > 1.

data = np.array([[70, 75, 80],
                 [60, 65, 70]])

row_means = data.mean(axis=1)
column_means = data.mean(axis=0)
population_sd = data.std()
sample_sd = data.std(ddof=1)

For a two-dimensional array, axis=0 reduces down rows and produces one result per column; axis=1 reduces across columns and produces one result per row. This wording is often confusing, so verify the resulting shape. keepdims=True retains reduced axes of length 1, which can make later broadcasting explicit.

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

NaN-aware variants such as np.nanmean ignore missing values, but they should not be used automatically: a result based on very few observations may be misleading. Missingness should be measured and governed as part of data preparation.

5. Vectorization and Performance

Consider computing squares for one million values. A conventional loop repeatedly executes Python bytecode:

out = []
for value in values:
    out.append(value * value)

A vectorized operation delegates the inner loop to compiled NumPy code:

out = values ** 2

If n values are processed, both approaches perform roughly O(n) arithmetic operations. The major practical difference is that vectorization reduces interpreter overhead and often accesses contiguous memory efficiently. It does not change every algorithm's asymptotic complexity, and it may allocate temporary arrays.

Benchmark representative workloads rather than relying on intuition.

import timeit

time_loop = timeit.timeit(
    'sum(v * v for v in values)',
    setup='values = range(100000)',
    number=10
)

For NumPy comparisons, construct arrays outside the timed expression and report the Python and NumPy versions, hardware, data size, and number of repetitions. Avoid comparing array construction time with an already constructed list unless that is the actual application scenario.

Vectorization is not always best. A small operation, complex branch-heavy logic, object data, or a problem requiring sequential dependence may be clearer or faster with a loop, compiled extension, or specialised tool. Memory traffic can dominate arithmetic, and chained expressions can create temporary arrays. In high-performance workloads, tools such as numexpr, Numba, sparse arrays, or distributed systems may be appropriate, but they introduce additional complexity.

6. Linear Algebra with NumPy

Linear algebra represents vectors, matrices, and transformations. Let A be an m × n matrix and x an n-element vector. The matrix-vector product y = Ax has components

and produces an m-element vector. Matrix multiplication is defined only when the inner dimensions agree.

A = np.array([[2.0, 1.0],
              [1.0, 3.0]])
x = np.array([4.0, 5.0])
y = A @ x

Here, y = [13, 19], because 2×4 + 1×5 = 13 and 1×4 + 3×5 = 19.

Useful operations include:

A.T                  # transpose
np.linalg.det(A)     # determinant
np.linalg.inv(A)     # inverse, when nonsingular
np.linalg.norm(x)    # Euclidean norm by default
np.linalg.eig(A)     # eigenvalues and eigenvectors

For solving Ax = b, use np.linalg.solve(A, b) rather than explicitly calculating A⁻¹b. An inverse exists only for a square, nonsingular matrix, whereas solve is more numerically appropriate and generally more efficient.

b = np.array([13.0, 19.0])
solution = np.linalg.solve(A, b)

The solution is [4, 5]. Numerical results should be checked using a residual:

A small residual norm ‖r‖ indicates that the computed solution satisfies the equations approximately, although a small residual alone does not guarantee good sensitivity when the matrix is ill-conditioned. Use np.linalg.cond(A) to investigate conditioning. Singular-value decomposition is valuable for least squares, dimensionality reduction, and diagnosing rank deficiency.

For an overdetermined system, np.linalg.lstsq(X, y, rcond=None) computes a least-squares solution. The ordinary least-squares objective is

where X is the design matrix, y is the response vector, and β is the coefficient vector. This formulation assumes a linear model in the parameters; statistical interpretation additionally requires assumptions about errors, sampling, and model specification.

7. Random Sampling and Reproducibility

NumPy's modern random interface uses a Generator.

rng = np.random.default_rng(42)
coin_flips = rng.integers(0, 2, size=20)
noise = rng.normal(loc=0.0, scale=1.0, size=1000)
uniform_values = rng.uniform(0.0, 1.0, size=100)

A seed initializes a deterministic pseudorandom sequence. It supports reproducibility, debugging, and experimental comparison; it does not create cryptographic security. For passwords, tokens, or security-sensitive sampling, use an appropriate cryptographic random source.

The normal distribution with mean μ and standard deviation σ > 0 has density

This equation describes a model, not a guarantee that generated values will have exactly the specified sample mean and standard deviation. Random samples fluctuate. Sampling methods should match the problem: choice supports sampling from categories, permutation shuffles without replacement, and multivariate_normal generates correlated vectors when a valid covariance matrix is supplied.

Use separate generators or carefully controlled streams in parallel experiments. Record seeds, versions, data sources, and sampling assumptions in research and production pipelines.

8. Worked Application: Sensor-Based Energy Monitoring

The following is an illustrative case. A building-management team records hourly electricity use from three buildings over four hours. The goal is to calculate total consumption, compare buildings, identify unusual readings, and estimate the effect of a 10% efficiency reduction.

import numpy as np

energy = np.array([
    [120, 135, 128, 140],
    [90,  95,  92,  110],
    [160, 155, 170, 165]
], dtype=float)

building_totals = energy.sum(axis=1)
hourly_totals = energy.sum(axis=0)
building_means = energy.mean(axis=1)
threshold = 150
unusual = energy > threshold
reduced = energy * 0.90

Analysis

The building totals are [523, 387, 650] energy units, and the hourly totals are [370, 385, 390, 415]. The total across all observations is 1,552 units. Building 3 has the greatest total and average consumption; its mean is 650 ÷ 4 = 162.5 units per hour. The Boolean mask identifies readings above 150: all four readings for Building 3, and none for Buildings 1 or 2.

The reduction model multiplies every reading by 0.90. This assumes the intervention produces a uniform 10% reduction, an assumption that should be validated rather than treated as fact. The estimated new total is 1,552 × 0.90 = 1,396.8 units, a reduction of 155.2 units.

Professional interpretation

The array representation makes the dimensions explicit: rows represent buildings and columns represent hours. Aggregation along different axes answers different operational questions. However, the analysis does not establish causation, account for weather or occupancy, or identify whether sensors were calibrated. Production use should include timestamps, units, missing-value rules, calibration records, access controls, and an audit trail. Energy data can also reveal occupancy patterns, so privacy and appropriate access are relevant governance concerns.

9. Numerical Reliability and Best Practices

Floating-point numbers approximate real numbers using finite binary representations. Consequently, some decimal values cannot be represented exactly.

0.1 + 0.2 == 0.3        # commonly False
np.isclose(0.1 + 0.2, 0.3) # True under default tolerances

Use tolerances deliberately. A relative-and-absolute comparison can be expressed conceptually as

where atol is absolute tolerance and rtol is relative tolerance. The appropriate tolerances depend on scale, units, and error requirements. For financial decimal arithmetic, binary floating point may be inappropriate; use decimal or integer minor units where required.

Other best practices include:

1. Validate shapes, units, ranges, and data types at system boundaries.
2. Prefer explicit axis arguments over relying on defaults.
3. Avoid silent casting and inspect warnings for overflow, invalid operations, and division by zero.
4. Preserve raw data and document transformations.
5. Use copies when mutation could create hidden side effects.
6. Profile before optimising, and measure memory as well as execution time.
7. Test edge cases: empty arrays, NaNs, repeated values, extreme magnitudes, and singular matrices.
8. Keep random seeds and software versions for reproducible experiments.
9. Never infer more precision than the measurements support.
10. Separate computational correctness from statistical validity and business interpretation.

In machine learning, NumPy arrays frequently form the interface between data preparation, model training, evaluation, and deployment. Shape errors can produce immediate exceptions, but semantic errors—such as mixing rows and columns or scaling training and test data together—may remain hidden. Numerical pipelines therefore require both automated tests and domain review.

Contemporary Developments

NumPy remains a foundational array interface, but modern numerical computing increasingly includes array-compatible libraries for GPUs, distributed systems, automatic differentiation, and sparse computation. The Array API standard seeks to improve portability across array libraries. These developments do not eliminate the need to understand shape, dtype, numerical stability, and memory layout. They make those concepts more important because the same expression may execute on different hardware with different performance and precision characteristics.

Chapter Summary

NumPy provides dense, typed, multidimensional arrays for efficient numerical computing. Array shape, axis, dtype, and memory-sharing behaviour determine how operations work. Indexing and Boolean masks support precise selection; reshaping changes organisation without changing element count; and broadcasting enables compatible arrays to interact without explicit replication. Vectorized operations reduce Python interpreter overhead, although performance must be measured and memory use considered. Aggregation uses axis semantics to produce summaries, while np.linalg supports matrix products, equation solving, least squares, norms, and decompositions. NumPy's Generator provides reproducible pseudorandom sampling, but seeds are not security mechanisms. Reliable numerical work requires attention to floating-point error, conditioning, missing data, units, reproducibility, privacy, and interpretation.

Key Terms

ndarray: NumPy's multidimensional homogeneous array object.

Shape: The length of an array along each axis.

Axis: A dimension of an array used to specify operation direction.

Dtype: The data type used to store array elements.

View: An array representation that may share memory with another array.

Copy: An independent array with separate storage.

Broadcasting: NumPy's rules for applying operations to compatible shapes.

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

Ufunc: A compiled universal function that operates elementwise on arrays.

Boolean mask: An array of Boolean values used to select elements.

Fancy indexing: Indexing with integer sequences or arrays.

Aggregation: Reduction of multiple values into summaries such as sums or means.

Matrix multiplication: A product defined when the inner dimensions of two arrays agree.

Residual: The difference between an observed right-hand side and a computed model result.

Conditioning: Sensitivity of a numerical problem to perturbations in input.

Pseudorandom generator: A deterministic algorithm producing sequences with random-like properties.

Floating-point error: Approximation error arising from finite numerical representation.

Reproducibility: The ability to obtain consistent results under documented conditions.

Review Questions

Short Answer Questions

1. What information is represented by an array's shape and dtype?
2. Explain the difference between a basic slice and a copied array.
3. State the two compatibility rules used by broadcasting.
4. Why is np.linalg.solve(A, b) generally preferred to explicitly computing an inverse?
5. What is the purpose of supplying a seed to a pseudorandom generator?

Long Answer Questions

1. Explain NumPy's array model and discuss how shape, axis, dtype, and memory layout affect computation.
2. Describe vectorization and broadcasting, including their benefits, risks, and appropriate use cases.
3. Explain mean, population standard deviation, and sample standard deviation, including the role of ddof.
4. Discuss numerical reliability in floating-point computation, matrix conditioning, and tolerance-based comparisons.
5. Design a NumPy workflow for cleaning, transforming, aggregating, and validating a multidimensional business dataset.

Analytical / Application Questions

1. An array has shape (24, 7, 3) representing hourly observations, days, and sensors. Determine its total number of elements and describe which axis should be reduced to obtain one value per sensor.
2. A dataset contains a matrix X of shape (500, 8) and a coefficient vector β of shape (8,). Explain why X @ β is valid, state the output shape, and contrast it with X * β.
3. A production model gives different results after a software upgrade. Propose a diagnostic plan involving dtype inspection, random seeds, tolerances, shape checks, version recording, and comparison of intermediate arrays.

Multiple Choice Questions with Answers

Question 1

Which NumPy attribute gives the number of axes in an array?

A. shape

B. ndim

C. size

D. dtype

Correct Answer: B. ndim

Question 2

What is the shape of the result of broadcasting an array with shape (3, 1) and an array with shape (4,)?

A. (3, 1)

B. (4, 1)

C. (3, 4)

D. (7,)

Correct Answer: C. (3, 4)

Question 3

For a two-dimensional array, what does axis=0 mean in a reduction such as mean(axis=0)?

A. Reduce across rows to produce one value per column

B. Reduce across columns to produce one value per row

C. Flatten the array before averaging

D. Add a new axis of length zero

Correct Answer: A. Reduce across rows to produce one value per column

Question 4

Which operation performs matrix multiplication in NumPy?

A. A * B

B. A + B

C. A / B

D. A @ B

Correct Answer: D. A @ B

Question 5

What is the preferred NumPy approach for solving Ax=b when A is an appropriate square nonsingular matrix?

A. np.linalg.solve(A, b)

B. np.mean(A, b)

C. np.linalg.det(A, b)

D. np.concatenate([A, b])

Correct Answer: A. np.linalg.solve(A, b)

Question 6

An array contains 12 elements. Which reshape is valid?

A. (2, 5)

B. (3, 4)

C. (4, 4)

D. (1, 13)

Correct Answer: B. (3, 4)

Question 7

Which statement about np.random.default_rng(42) is correct?

A. It provides cryptographically secure random numbers

B. It guarantees every sample has the exact theoretical mean

C. It creates a reproducible pseudorandom generator sequence

D. It removes all sampling error

Correct Answer: C. It creates a reproducible pseudorandom generator sequence

Question 8

What is a likely advantage of vectorized NumPy computation over a Python loop for large numeric arrays?

A. It changes every algorithm from O(n) to O(1)

B. It reduces Python interpreter overhead in the inner operation

C. It guarantees no temporary memory is allocated

D. It automatically corrects invalid input units

Correct Answer: B. It reduces Python interpreter overhead in the inner operation

Question 9

For a sample standard deviation in NumPy, which argument commonly uses the n−1 denominator?

A. axis=1

B. keepdims=True

C. ddof=1

D. dtype=float32

Correct Answer: C. ddof=1

Question 10

Why can np.isclose be preferable to == for comparing floating-point results?

A. It compares values using specified absolute and relative tolerances

B. It converts all values to strings

C. It guarantees exact mathematical equality

D. It ignores the magnitude of both values

Correct Answer: A. It compares values using specified absolute and relative tolerances