Numerical Computing with NumPy

Introduction

Data science depends on repeated numerical operations: transforming observations, calculating descriptive statistics, solving equations, simulating uncertainty, and preparing data for machine-learning algorithms. Python provides convenient syntax, but ordinary Python lists and loops are often inefficient for large numerical workloads. NumPy—short for Numerical Python—addresses this problem through compact, homogeneous multidimensional arrays and highly optimized operations implemented largely in compiled code.

NumPy is not merely a faster list. It provides a computational model based on array shape, data type, mathematical operations, memory layout, and broadcasting rules. Understanding these ideas is essential because many scientific Python libraries, including pandas, SciPy, scikit-learn, and parts of PyTorch, use NumPy arrays or concepts derived from them.

This chapter treats NumPy as both a programming library and a numerical-computing framework. Examples assume that NumPy has been installed and imported as follows:

import numpy as np

Learning Objectives

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

1. Explain the structure, data types, dimensions, shape, and memory characteristics of NumPy arrays.
2. Create, inspect, reshape, index, slice, and filter arrays safely.
3. Apply vectorized operations and broadcasting instead of unnecessary Python loops.
4. Perform aggregation, statistical calculation, and missing-value-aware analysis.
5. Use NumPy for matrix operations, linear systems, norms, and least-squares estimation.
6. Generate reproducible random samples using modern random-number generators.
7. Diagnose numerical, memory, and performance problems in array-based programs.
8. Select appropriate NumPy techniques for a realistic data-science problem.

1. The NumPy Array Model

A NumPy array, commonly called an ndarray, is an indexed collection of values of a declared data type. Unlike a standard Python list, an ordinary NumPy array normally stores values of one primary type in a compact memory block. This regularity permits low-level optimized operations.

Important array attributes include:

For an array with shape (r, c), r is the number of rows and c is the number of columns, so its total number of elements is N = r × c. More generally, if the shape is (d₁, d₂, …, dₖ), then N = ∏_i=1^k d_i.

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

The array has two axes. Axis 0 moves between rows; axis 1 moves between columns. This distinction is central to aggregation and reshaping.

Creating arrays

Common construction methods include:

np.array([2, 4, 6])
np.zeros((2, 3))
np.ones(4, dtype=np.int32)
np.full((2, 2), 7.5)
np.arange(0, 10, 2)       # 0, 2, 4, 6, 8
np.linspace(0, 1, 5)      # five evenly spaced values
np.eye(3)                 # 3 × 3 identity matrix

arange(start, stop, step) is useful for integer-like steps, but floating-point steps can accumulate rounding differences. When the number of desired points is known, linspace(start, stop, num) is often preferable. In linspace, the endpoint is included by default.

Data types and conversion

NumPy data types include signed and unsigned integers, floating-point numbers, Boolean values, complex numbers, and fixed-width strings. The choice of type affects memory use, range, precision, and operation results.

values = np.array([1, 2, 3])
real_values = values.astype(np.float64)

Converting a floating-point array to an integer array discards the fractional part rather than rounding it:

np.array([1.9, -1.9]).astype(int)   # [1, -1]

This conversion may therefore lose information. Also, integer arithmetic can overflow when values exceed the representable range of the chosen integer type. For scientific and financial calculations, dtype selection should be deliberate rather than left to chance.

2. Indexing, Slicing, and Selection

For a one-dimensional array, indexing begins at zero:

scores = np.array([61, 74, 88, 93])
scores[0]       # 61
scores[-1]      # 93
scores[1:3]     # [74, 88]

The slice a[start:stop:step] includes start but excludes stop. In two dimensions, use one index per axis:

matrix = np.array([[1, 2, 3], [4, 5, 6]])
matrix[1, 2]     # 6
matrix[:, 1]     # second column: [2, 5]
matrix[0, :]     # first row: [1, 2, 3]

A basic slice usually produces a view, not an independent copy. Changes to the view can change the original array:

row = matrix[0, :]
row[0] = 99
# matrix is now [[99, 2, 3], [4, 5, 6]]

Use .copy() when independent storage is required:

safe_row = matrix[0, :].copy()

Boolean and fancy indexing

Boolean indexing selects values satisfying a condition:

scores = np.array([61, 74, 88, 93])
passed = scores[scores >= 70]

For multiple conditions, use elementwise operators and parentheses. Python’s and and or are not suitable for array-wise comparisons:

selected = scores[(scores >= 70) & (scores < 90)]

Fancy indexing uses an array of integer positions:

scores[[0, 3]]       # values at positions 0 and 3

Boolean indexing creates a new array. It is therefore convenient and safe for selection, but repeated large selections may increase memory use.

3. Vectorization and Broadcasting

Vectorization means expressing an operation on an entire array rather than writing an explicit Python loop. For example:

temperatures_c = np.array([0.0, 10.0, 20.0])
temperatures_f = temperatures_c * 9 / 5 + 32

The mathematical transformation is applied element by element. If C is temperature in degrees Celsius and F is temperature in degrees Fahrenheit, then:

The formula assumes the Celsius and Fahrenheit scales are being converted as temperature readings, not temperature differences. A difference of 1 °C corresponds to 1.8 °F, without the additive 32 offset.

Vectorization generally improves clarity and often improves speed because NumPy delegates loops to optimized compiled implementations. It does not guarantee that every vectorized expression is optimal: several temporary arrays may be created, increasing memory traffic.

Broadcasting

Broadcasting allows arrays with compatible shapes to participate in one operation. NumPy compares dimensions from right to left. Two dimensions are compatible when they are equal or when 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 shape (3, 4):

row = np.array([[10, 20, 30, 40]])       # (1, 4)
column = np.array([[1], [2], [3]])       # (3, 1)
result = column + row                    # (3, 4)

Conceptually, each value in column is paired with every value in row. Broadcasting does not necessarily create literal copies of the repeated data, which can save memory.

A frequent application is standardizing columns. Given observations X_ij, column means μ_j, and column standard deviations σ_j, the standardized value is:

Here σ_j > 0 is required. A column with zero variance cannot be standardized using this formula because division by zero is undefined.

X = np.array([[10., 100.], [12., 110.], [14., 120.]])
mean = X.mean(axis=0)
std = X.std(axis=0)
Z = (X - mean) / std

The vectors mean and std have shape (2,) and broadcast across the rows.

4. Reshaping and Combining Arrays

Reshaping changes the organization of an array without changing its values, provided the number of elements remains constant. Thus a six-element array can be reshaped to (2, 3) or (3, 2), but not (4, 2).

x = np.arange(6)
x.reshape(2, 3)
x.reshape(-1, 2)     # NumPy infers the first dimension as 3

The inferred dimension must be uniquely determined. reshape may return a view when the memory layout permits; otherwise it may return a copy. Code should not rely on whether a view or copy is returned unless this behavior has been explicitly checked.

Related operations include:

x.ravel()            # flattened view when possible
x.flatten()          # flattened copy
x.T                  # transpose for a 2-D array
np.concatenate([a, b], axis=0)
np.vstack([a, b])
np.hstack([a, b])

concatenate requires compatible dimensions except along the chosen axis. Explicit shape checks are valuable in production data pipelines because a one-dimensional shape (n,) is not identical to a column-vector shape (n, 1).

5. Aggregation and Descriptive Statistics

Aggregations reduce multiple values to summaries. Common methods include sum, mean, median, min, max, std, var, and argmax.

For values x₁, x₂, …, x_n, the arithmetic mean is:

The population variance is:

The sample variance, used when estimating a population variance from a sample, is:

NumPy’s np.var and np.std use ddof=0 by default, corresponding to the population-style denominator n. For the sample estimate, use ddof=1, provided the sample contains at least two observations.

x = np.array([2., 4., 6., 8.])
x.mean()                 # 5.0
x.std()                  # population standard deviation
x.std(ddof=1)            # sample standard deviation

The axis argument

For a matrix, axis=0 reduces down rows and produces one result per column; axis=1 reduces across columns and produces one result per row.

sales = np.array([[10, 12, 8], [9, 15, 11]])
sales.sum(axis=0)     # totals for each product/column
sales.sum(axis=1)     # totals for each period/row

keepdims=True retains reduced axes as length-one dimensions, which is useful for broadcasting:

column_mean = sales.mean(axis=0, keepdims=True)
centered = sales - column_mean

Missing values

A regular float array may contain np.nan, representing a missing or undefined floating-point value. Most ordinary reductions propagate it:

x = np.array([2.0, np.nan, 6.0])
x.mean()             # nan
np.nanmean(x)       # 4.0

nanmean, nanmedian, nansum, and related functions ignore NaN values, but they do not determine whether missingness is random, systematic, or meaningful. Missing-data treatment is a modelling and governance decision, not merely a programming detail. An all-NaN slice may produce a warning and a NaN result.

6. Universal Functions and Numerical Stability

NumPy universal functions, or ufuncs, operate elementwise and include np.sqrt, np.exp, np.log, np.sin, np.abs, and np.clip.

x = np.array([1., 4., 9.])
roots = np.sqrt(x)

The natural logarithm np.log(x) is defined for positive real x. For zero it returns negative infinity, and for negative real inputs it generally produces NaN in a real-valued array. Domain assumptions should therefore be checked before applying mathematical functions.

Numerical stability concerns how rounding and finite precision affect results. For example, computing np.exp(x) for very large positive x can overflow. Subtracting nearly equal large numbers can cause loss of significant digits. Practical safeguards include:

• selecting an appropriate dtype;
• checking inputs and outputs with np.isfinite;
• avoiding unnecessary subtraction of nearly equal quantities;
• using stable library routines such as np.logaddexp where applicable;
• tolerating small floating-point discrepancies with np.isclose rather than exact equality.

Floating-point representation is approximate. Consequently, a calculation such as 0.1 + 0.2 == 0.3 may be false. Use a stated tolerance when comparing numerical results:

np.isclose(a, b, rtol=1e-5, atol=1e-8)

The absolute tolerance atol matters when values are near zero; the relative tolerance rtol scales with the magnitude of the reference value.

7. Linear Algebra with numpy.linalg

Linear algebra is fundamental to regression, dimensionality reduction, optimization, and scientific modelling. NumPy provides a dedicated linalg module.

For vectors a and b of equal length, the dot product is:

In NumPy, use np.dot(a, b) or, for suitable arrays, a @ b. The @ operator expresses matrix multiplication, whereas * performs elementwise multiplication.

A = np.array([[2., 1.], [1., 3.]])
b = np.array([5., 7.])
product = A @ b

A linear system is written as:

where A is an m × n coefficient matrix, x is an unknown vector, and b is a known vector. If A is square and nonsingular, a unique solution exists. In computation, use np.linalg.solve(A, b) rather than explicitly calculating A^-1, because solving is usually more efficient and numerically preferable.

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

A matrix inverse exists only for a square, nonsingular matrix. np.linalg.inv(A) should be reserved for situations that genuinely require the inverse as an object.

The Euclidean norm of a vector is:

np.linalg.norm(x)

Least squares

When a system is overdetermined, an exact solution may not exist. Least squares chooses coefficients β that minimize the sum of squared residuals:

Here X is the design matrix and y is the response vector. Use np.linalg.lstsq:

X = np.array([[1., 1.], [1., 2.], [1., 3.]])
y = np.array([2., 3., 5.])
beta, residuals, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)

The first column of ones represents an intercept. Least squares assumes that the chosen model form is appropriate for the analysis; NumPy computes the solution but does not establish causality, validate assumptions, or prevent leakage between training and test data.

8. Random Sampling and Reproducibility

Randomness supports simulation, resampling, randomized algorithms, and experimental design. NumPy’s modern interface uses a Generator object:

rng = np.random.default_rng(42)
values = rng.normal(loc=0.0, scale=1.0, size=5)

The seed initializes a pseudorandom sequence. A fixed seed makes a computation reproducible on the same software environment and algorithmic configuration; it does not make the generator cryptographically secure.

Examples:

rng.integers(1, 7, size=10)                 # integer outcomes 1 through 6
rng.choice(np.array(['A', 'B', 'C']), size=5)
rng.uniform(0, 1, size=1000)
rng.normal(100, 15, size=1000)

For a normal distribution with mean μ and standard deviation σ > 0, the density is:

A simulation is meaningful only when the distribution, parameters, independence assumptions, and random seed policy are justified. Do not use predictable pseudorandom numbers for passwords, tokens, or other security-sensitive purposes; use a cryptographic source such as Python’s secrets module.

9. Performance-Aware Computation

Performance has at least two dimensions: execution time and memory consumption. A vectorized operation may be faster than a Python loop because it reduces interpreter overhead, but a large expression can create several temporary arrays.

Measuring performance

Use representative data and repeat measurements:

import timeit

a = np.random.default_rng(1).random(1_000_000)
time = timeit.timeit('a * a + 2 * a', globals=globals(), number=10)

Big-O notation describes growth as input size increases, but constant factors, memory bandwidth, cache behavior, and algorithm implementation also matter. A single benchmark on a small array is not a reliable performance conclusion.

Practical principles

1. Prefer vectorized NumPy operations to Python loops for regular numerical work.
2. Choose an appropriate dtype; reducing precision can reduce memory but may harm accuracy.
3. Avoid unnecessary copies and inspect whether .copy() is required.
4. Reuse arrays or use the out= argument where safe and supported.
5. Process data in chunks when the full dataset cannot fit in memory.
6. Use contiguous, well-shaped arrays where downstream libraries benefit from them.
7. Profile before optimizing; bottlenecks are often different from assumptions.
8. Use compiled or parallel tools only after a clear baseline has been established.

A loop may remain preferable when operations are highly conditional, sequential, irregular, or difficult to express with array operations. Vectorization is a tool, not an absolute rule.

10. Applied Case: Detecting Anomalous Sensor Readings

Background

The following is an illustrative case. A manufacturing facility collects temperature readings from four sensors every hour. The data-science team must calculate sensor summaries, standardize readings, identify unusually high observations, and estimate the average temperature profile.

Problem

The team receives a matrix whose rows represent hours and columns represent sensors. One reading is missing. The team needs a transparent and reproducible NumPy workflow that does not silently mix sensor scales.

readings = np.array([
    [20.1, 20.4, 19.9, 20.3],
    [20.3, 20.2, 20.1, 20.4],
    [20.2, np.nan, 20.0, 20.2],
    [25.8, 20.5, 20.2, 20.3]
])

Analysis and approach

First inspect structure and missingness:

print(readings.shape)                 # (4, 4)
print(np.isnan(readings).sum(axis=0))  # missing count per sensor

Compute sensor means while ignoring the missing value:

sensor_means = np.nanmean(readings, axis=0)

For a simple illustrative analysis, impute the missing second-sensor value with that sensor’s mean. This is not automatically the best scientific method; it is used here to demonstrate broadcasting and must be documented:

filled = readings.copy()
missing = np.isnan(filled)
filled[missing] = np.take(sensor_means, np.where(missing)[1])

Standardize each sensor using its own mean and standard deviation:

means = filled.mean(axis=0)
stds = filled.std(axis=0)
z = (filled - means) / stds

Flag observations with absolute standardized value above 2:

flags = np.abs(z) > 2

The threshold 2 is a screening rule, not a universal law. Its usefulness depends on distributional assumptions, sample size, sensor reliability, and the consequences of false alarms. In this small dataset, the fourth hour’s first sensor reading will be conspicuous relative to the other readings, but a production system should use validated limits and temporal context.

The average across sensors for each hour is:

hourly_mean = filled.mean(axis=1)

Outcome

The workflow produces: a documented missing-value count, sensor-specific summary statistics, standardized values comparable across sensors, anomaly flags, and hourly aggregate values. Because the fourth row contains an extreme first-sensor reading, the team should investigate calibration, environmental conditions, and whether the reading reflects a real event before deleting or correcting it.

Learning points

• Shape and axis semantics determine whether summaries are calculated by sensor or by hour.
• Broadcasting enables column-wise standardization without explicit loops.
• Imputation changes the data and must be recorded.
• An anomaly flag is an investigative signal, not proof of an error.
• Reproducibility requires preserving code, versions, input data, and documented decisions.

11. Limitations, Governance, and Contemporary Practice

NumPy performs numerical operations; it does not automatically provide data provenance, access control, statistical validity, or domain interpretation. Analysts must validate units, ranges, timestamps, identifiers, and assumptions before calculation.

Important governance practices include:

• retain the original data and distinguish it from transformed arrays;
• document imputation, filtering, dtype conversion, and random seeds;
• avoid exposing personally identifiable information in demonstrations or logs;
• test edge cases such as empty arrays, all-missing columns, zero variance, overflow, and incompatible shapes;
• record library versions when reproducibility is important;
• validate numerical outputs with independent calculations or known test cases;
• distinguish exploratory thresholds from operational or regulated decision rules.

Contemporary numerical workflows increasingly combine NumPy with pandas for labelled tables, SciPy for specialized scientific algorithms, scikit-learn for machine learning, and GPU-oriented frameworks for large parallel workloads. NumPy remains important because it supplies a clear array abstraction and interoperability conventions. Array APIs and alternative implementations may execute similar code on CPUs, GPUs, or other accelerators, but portability is not guaranteed for every function or dtype. Professional practice therefore requires testing on the actual deployment environment.

Chapter Summary

NumPy provides efficient homogeneous multidimensional arrays and a powerful foundation for numerical data science. Array attributes such as shape, ndim, size, and dtype describe the computational object. Indexing and slicing provide structured access, while Boolean and fancy indexing support selection. Vectorization replaces many explicit loops with optimized array operations, and broadcasting aligns compatible shapes for elementwise computation. Aggregations depend critically on the selected axis, and ddof distinguishes population-style from sample-style variance and standard deviation. Missing values require explicit, documented treatment.

NumPy’s universal functions support mathematical transformation, but domain restrictions and floating-point limitations must be considered. The linalg module supports matrix multiplication, norms, linear-system solving, and least-squares estimation. The modern random Generator interface supports reproducible simulation, although pseudorandom generation is not a substitute for cryptographic security. Finally, performance-aware programming balances vectorization, memory use, numerical accuracy, profiling, and maintainability. NumPy can execute calculations efficiently, but professional data science still requires valid assumptions, transparent governance, and critical interpretation.

Key Terms

Array: An indexed NumPy data structure containing values arranged along one or more axes.

Axis: A dimension of an array along which indexing or reduction can occur.

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

Data type (dtype): The type and precision used to store array elements.

Dimension (ndim): The number of axes in an array.

Vectorization: Applying an operation to an entire array rather than iterating in Python.

View: An array object that shares data storage with another array.

Copy: An independent array allocation whose values can change without changing the source.

Boolean indexing: Selecting array elements using a Boolean condition array.

Fancy indexing: Selecting elements using integer index arrays.

Aggregation: Reducing multiple values to a summary such as a sum or mean.

keepdims: An option that retains reduced axes as dimensions of length one.

NaN: A floating-point marker representing a missing or undefined numerical value.

Universal function (ufunc): An elementwise NumPy function such as sqrt or exp.

Matrix multiplication: Algebraic multiplication represented in NumPy by @.

Least squares: Estimation that minimizes the squared Euclidean residual norm.

Random generator: A NumPy object used to produce pseudorandom samples.

Broadcast-compatible shape: A pair of shapes whose aligned dimensions are equal or include 1.

Numerical stability: Resistance of a computation to unacceptable effects of finite precision and rounding.

Profiling: Measuring execution time or resource use to locate computational bottlenecks.

Review Questions

5 Short Answer Questions

1. What information is provided by an array’s shape and dtype?
2. Why can modifying a slice modify the original NumPy array?
3. Distinguish elementwise multiplication * from matrix multiplication @.
4. What is the purpose of ddof=1 when calculating standard deviation?
5. Why is a fixed random seed useful in a data-science experiment?

5 Long Answer Questions

1. Explain NumPy broadcasting, including the rules for determining whether two shapes are compatible. Illustrate your answer with a column-wise standardization example.
2. Compare basic slicing, Boolean indexing, and fancy indexing in terms of syntax, use, and whether they commonly return views or copies.
3. Discuss the assumptions and computational considerations involved in solving a linear system with np.linalg.solve and fitting a least-squares model with np.linalg.lstsq.
4. Explain how axis selection affects aggregation for a two-dimensional dataset. Include an example of a row-wise and a column-wise calculation.
5. Discuss at least five performance-aware NumPy practices and explain why profiling is necessary before optimization.

3 Analytical / Application Questions

1. A matrix of sales has shape (12, 5), representing 12 months and 5 products. Design NumPy operations to calculate monthly totals, product totals, product-wise standardized sales, and the month with the highest total. State the axis used for each operation.
2. A sensor array contains NaN values, one constant-variance-zero column, and a very large outlier. Propose a robust analysis workflow. Explain how you would handle missingness, avoid division by zero during standardization, investigate the outlier, and document decisions.
3. An algorithm computes a large expression involving several temporary arrays and is slower than expected. Develop a performance investigation plan involving correctness tests, timing, memory assessment, dtype review, possible chunking, and comparison with a loop or alternative implementation.

Multiple Choice Questions with Answers

Question 1

Which attribute gives the tuple of lengths along an array’s axes?

A. dtype

B. shape

C. itemsize

D. ndim

Correct Answer: B. shape

Question 2

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

A. The scalar 12

B. [1, 2, 3, 1, 2, 3]

C. [2, 4, 6]

D. An error because arrays cannot use scalars

Correct Answer: C. [2, 4, 6]

Question 3

Which expression correctly selects values of x that are at least 5 and less than 10?

A. x[x >= 5 and x < 10]

B. x[(x >= 5) & (x < 10)]

C. x[x >= 5 or x < 10]

D. x[5 <= x <= 10]

Correct Answer: B. x[(x >= 5) & (x < 10)]

Question 4

For an array with shape (4, 3), what is the shape of array.mean(axis=0)?

A. (4,)

B. (3,)

C. (1, 3)

D. `()

Correct Answer: B. (3,)

Question 5

Which NumPy function is generally preferred for solving Ax = b without explicitly forming the inverse of A?

A. np.linalg.solve

B. np.linalg.det

C. np.linalg.diagonal

D. np.sort

Correct Answer: A. np.linalg.solve

Question 6

What does np.nanmean([2.0, np.nan, 6.0]) return?

A. nan

B. 2.0

C. 4.0

D. 8.0

Correct Answer: C. 4.0

Question 7

Which condition is required for two aligned dimensions to be broadcast-compatible?

A. They must always be equal

B. Their sum must equal the number of dimensions

C. They must both be greater than 1

D. They must be equal or one must be 1

Correct Answer: D. They must be equal or one must be 1

Question 8

In NumPy, what does a * b usually perform for two same-shaped numeric arrays?

A. Elementwise multiplication

B. Matrix inversion

C. A dot product only

D. Concatenation

Correct Answer: A. Elementwise multiplication

Question 9

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

A. It creates a cryptographically secure generator

B. It creates a reproducible pseudorandom generator when used consistently

C. It guarantees identical results across every possible NumPy implementation

D. It returns only normally distributed values

Correct Answer: B. It creates a reproducible pseudorandom generator when used consistently

Question 10

For x = np.array([1.9, -1.9]).astype(int), what is the usual result?

A. [2, -2]

B. [1, -1]

C. [1.9, -1.9]

D. An exception due to negative values

Correct Answer: B. [1, -1]