9. Numerical Computing with NumPy and Scientific Python

Introduction and Context

Artificial intelligence and machine learning depend on numerical computation. Data is represented as vectors, matrices, tensors, probability distributions, and sequences of measurements. Model training repeatedly performs operations such as matrix multiplication, gradient calculation, normalization, sampling, and optimization. Although Python provides an expressive programming interface, ordinary Python loops are often inefficient for large numerical workloads. NumPy addresses this problem through compact, typed, multidimensional arrays and vectorized operations implemented largely in optimized compiled code.

NumPy is not merely a faster list library. It provides a computational model based on homogeneous arrays, predictable memory layouts, broadcasting rules, linear algebra routines, random-number generation, and numerical functions. The wider scientific Python ecosystem extends this model through libraries such as SciPy for numerical methods and scientific algorithms, pandas for tabular data, Matplotlib for visualization, and scikit-learn for classical machine learning.

For postgraduate AI and ML students, numerical computing requires more than learning syntax. It requires understanding shapes, data types, memory, numerical stability, algorithmic complexity, and the difference between mathematically equivalent expressions and computationally equivalent implementations.

Learning Objectives

After completing this chapter, the learner should be able to:

1. Explain the structure and computational advantages of NumPy arrays.

2. Create, inspect, index, slice, reshape, and combine multidimensional arrays.

3. Apply vectorization and broadcasting correctly.

4. Compare vectorized expressions with explicit Python loops.

5. Perform common linear algebra operations and interpret their shapes.

6. Generate reproducible random samples using modern NumPy random generators.

7. Use SciPy for numerical integration, root finding, optimization, and related tasks.

8. Identify numerical precision, memory, performance, and reproducibility risks.

9. Select appropriate numerical methods for AI, ML, and scientific applications.

10. Design efficient and reliable numerical workflows.

The NumPy Array as a Computational Object

The central object in NumPy is the `ndarray`, an n-dimensional array whose elements generally share one data type. A Python list may contain integers, strings, and objects simultaneously. An ndarray normally stores values of a uniform type in a structured memory region. This uniformity enables low-level optimized operations and reduces the overhead associated with processing individual Python objects.

```python

import numpy as np

x = np.array([2, 4, 6, 8])

print(x.shape) # (4,)

print(x.ndim) # 1

print(x.dtype) # commonly int64 or int32

print(x.size) # 4

```

The `shape` describes the length of each axis, `ndim` gives the number of axes, `size` gives the total number of elements, and `dtype` describes the element representation. These properties are fundamental when debugging machine-learning pipelines.

Arrays can be constructed from sequences or generated directly:

```python

A = np.zeros((3, 4))

B = np.ones((2, 2), dtype=np.float32)

C = np.eye(4)

D = np.arange(0, 10, 2)

E = np.linspace(0, 1, 5)

```

`arange` uses a step size and is convenient for integer-like sequences. `linspace` specifies the number of evenly spaced values and is usually preferable when endpoints and sampling density matter. In floating-point computations, `linspace` avoids some of the endpoint uncertainty associated with repeatedly adding a non-integer step.

Shape, Axis, and Dimension

An array of shape `(m, n)` contains `m` rows and `n` columns. In NumPy, axis 0 refers to the first dimension and axis 1 to the second. For a matrix, `A.mean(axis=0)` computes a mean for each column, while `A.mean(axis=1)` computes a mean for each row.

```python

A = np.array([[1, 2, 3],

[4, 5, 6]])

column_means = A.mean(axis=0) # [2.5, 3.5, 4.5]

row_means = A.mean(axis=1) # [2., 5.]

```

Many errors in AI code are shape errors rather than mathematical errors. A programmer may intend to normalize each feature but accidentally normalize each observation. Explicitly documenting the meaning of each axis is therefore a professional practice.

Indexing, Slicing, and Views

NumPy uses zero-based indexing. For a two-dimensional array, `A[i, j]` selects one element, `A[i, :]` selects a row, and `A[:, j]` selects a column.

```python

A = np.arange(12).reshape(3, 4)

first_row = A[0, :]

last_column = A[:, -1]

block = A[1:, 1:3]

```

Basic slicing commonly returns a view into the original array rather than an independent copy. Consequently, modifying a slice may modify the source array:

```python

row = A[0, :]

row[:] = 99

```

If independent data is required, use `A[0, :].copy()`. Confusing views and copies can cause subtle data corruption in preprocessing pipelines.

Boolean indexing selects elements according to a condition:

```python

x = np.array([3, -1, 7, 0])

positive = x[x > 0]

x[x < 0] = 0

```

Integer-array indexing supports selection by specified positions. Unlike basic slicing, advanced indexing generally creates a copy. Boolean masks are especially useful for filtering invalid observations, applying thresholds, and implementing piecewise rules.

Vectorization and Broadcasting

Vectorization means expressing an operation on an entire array rather than writing a Python loop that processes one element at a time.

```python

x = np.arange(1_000_000)

y = 3 * x + 2

```

This expression delegates iteration to optimized compiled code. It is usually shorter, clearer, and faster than an equivalent loop. The benefit arises not because the mathematical operation changes, but because repeated interpreter-level overhead is reduced and contiguous numerical data can be processed efficiently.

Universal functions, or ufuncs, apply elementwise operations:

```python

z = np.sqrt(x)

probabilities = 1 / (1 + np.exp(-z))

```

Common ufuncs include arithmetic, trigonometric, exponential, logarithmic, comparison, and rounding operations. Reductions such as `sum`, `mean`, `max`, and `argmin` summarize arrays.

Broadcasting is NumPy’s rule for combining arrays with different shapes. NumPy compares dimensions from the rightmost side. Two dimensions are compatible when they are equal or one of them is 1. Missing leading dimensions are treated as 1.

```python

X = np.array([[1., 2.],

[3., 4.],

[5., 6.]])

mu = X.mean(axis=0) # shape (2,)

sigma = X.std(axis=0) # shape (2,)

X_standardized = (X - mu) / sigma

```

Here, `mu` and `sigma` are broadcast across the rows. Broadcasting does not normally replicate data physically; it provides a compatible interpretation of the smaller array. Nevertheless, careless broadcasting can produce a mathematically valid but unintended result. Shape assertions and explicit reshaping improve reliability.

```python

weights = np.array([0.2, 0.5, 0.3])

weights = weights.reshape(1, -1) # row vector

```

The distinction between `(n,)`, `(n, 1)`, and `(1, n)` is important. A one-dimensional array has no explicit row or column orientation, while the latter two are two-dimensional arrays.

Reshaping and Memory Layout

`reshape` changes the apparent dimensions without changing the data order when possible. The number of elements must remain constant.

```python

x = np.arange(24)

T = x.reshape(2, 3, 4)

flat = T.reshape(-1)

```

The value `-1` allows NumPy to infer one dimension. `transpose` or `.T` exchanges axes for two-dimensional arrays, while `np.moveaxis` and `np.swapaxes` provide more general control for tensors.

Array layout affects performance. C-order arrays store the last axis contiguously, whereas Fortran-order arrays prioritize the first axis. Most users need not manipulate layout directly, but performance-sensitive applications should understand that non-contiguous views may require copying before passing data to external compiled routines.

Data types also affect memory and accuracy. A `float32` array uses approximately half the storage of a `float64` array, but has lower precision and a smaller numerical range. Integer overflow is possible when computations remain in a fixed-width integer type. Silent or implicit type conversion should not be assumed safe.

Linear Algebra for Machine Learning

Vectors and matrices are the language of many ML algorithms. NumPy uses `@` for matrix multiplication:

```python

X = np.array([[1., 2.], [3., 4.]])

w = np.array([0.5, -1.0])

b = 0.2

scores = X @ w + b

```

The `*` operator performs elementwise multiplication, not matrix multiplication. Confusing these operations is one of the most common numerical programming errors.

Important operations include:

```python

A.T # transpose

np.linalg.norm(w) # vector norm

np.linalg.det(A) # determinant

np.linalg.eig(A) # eigenvalues and eigenvectors

np.linalg.svd(A) # singular value decomposition

```

Although `np.linalg.inv(A) @ b` represents a solution to a linear system, explicitly computing an inverse is often less stable and less efficient than solving the system directly:

```python

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

```

For least-squares problems, use `np.linalg.lstsq`. For symmetric positive-definite systems, Cholesky-based methods may provide efficient and stable solutions. In practical ML systems, specialized libraries may offer better hardware utilization, automatic differentiation, or sparse-matrix support, but the underlying numerical principles remain the same.

Conditioning is a key concept. A well-conditioned problem changes its answer modestly when inputs are perturbed. An ill-conditioned problem can amplify small measurement or rounding errors. A large condition number often signals that numerical results require caution.

Random Sampling and Reproducibility

Randomness supports simulation, train-test splitting, initialization, bootstrapping, and uncertainty estimation. Modern NumPy recommends creating a random generator rather than relying on a global random state.

```python

rng = np.random.default_rng(42)

normal_samples = rng.normal(loc=0, scale=1, size=(1000, 3))

indices = rng.integers(0, 100, size=10)

choices = rng.choice([0, 1], size=20, p=[0.7, 0.3])

```

A seed does not make a process intrinsically random or non-random; it initializes a deterministic pseudorandom sequence. Reproducibility requires controlling all relevant sources of randomness, including data shuffling, model libraries, parallel execution, and hardware-specific behavior.

Sampling distributions must match the intended model. Uniform, normal, binomial, Poisson, exponential, and multivariate normal distributions each represent different assumptions. Poorly chosen sampling assumptions can produce misleading experiments even when code executes correctly.

Scientific Python and Numerical Methods

NumPy supplies the array foundation, while SciPy provides algorithms for scientific computation. Typical modules include `scipy.optimize`, `scipy.integrate`, `scipy.stats`, `scipy.linalg`, `scipy.signal`, and `scipy.sparse`.

Root finding seeks a value satisfying `f(x) = 0`:

```python

from scipy.optimize import brentq

root = brentq(lambda x: x**2 - 2, 0, 2)

```

Optimization finds a minimum or maximum of an objective function. The choice of method depends on differentiability, dimensionality, constraints, convexity, and noise. Numerical integration approximates an integral from function evaluations:

```python

from scipy.integrate import quad

area, error = quad(lambda x: np.exp(-x**2), 0, 1)

```

Numerical methods involve approximation. Tolerance settings express acceptable error, but tighter tolerances increase computational cost and may not improve results when input data itself is inaccurate. Algorithms should therefore be selected based on the problem’s mathematical structure rather than by default preference.

Sparse matrices are appropriate when most entries are zero, as in graph and text representations. Storing such data densely wastes memory and computation. SciPy sparse formats, including CSR and CSC, support specialized operations while requiring awareness that not every dense-array operation has a sparse equivalent.

Efficiency: From Loops to Reliable Vectorized Programs

A vectorized solution is not automatically optimal. It may allocate several large temporary arrays, increasing memory traffic. For example, an expression containing multiple intermediate results may be slower than a fused or in-place implementation.

```python

x = np.random.default_rng(1).random(10_000_000)

y = x.copy()

y *= 2

y += 1

```

The accidental leading space in the final line would cause an indentation error; more importantly, in-place operations should be used only when overwriting the original values is safe. In-place computation reduces allocations but can harm clarity or create aliasing problems.

Performance should be measured, not guessed. Use `timeit` for small operations and profiling tools for complete workflows. Compare equivalent implementations under realistic array sizes. Consider:

• algorithmic complexity;
• number and size of temporary arrays;
• memory bandwidth and cache behavior;
• data transfer between CPU and accelerator;
• use of compiled libraries;
• parallelism and thread configuration;
• numerical precision requirements.

Loops are not always wrong. A loop may be appropriate when each step depends on the previous step, when control flow is irregular, when arrays are small, or when vectorization would obscure a critical algorithm. Numba, Cython, JAX, PyTorch, and other tools can compile or transform suitable numerical code, but their semantics and performance characteristics differ.

Applied Case Study: Efficient Feature Standardization for a Predictive Model

Background

Consider an illustrative case involving a hospital analytics team developing a model to predict readmission risk. The data contains 2 million patient records and 80 numerical features. The team must standardize the features before model training while preserving missing-value handling and ensuring reproducibility.

Situation and Problem

An initial implementation loops through records and features in Python. It is slow, uses inconsistent means between training and validation data, and occasionally generates invalid values when a feature has zero variance.

Analysis

The desired transformation for feature j is

`z_ij = (x_ij - μ_j) / σ_j`,

where `μ_j` and `σ_j` are computed from the training set only. Computing statistics over the complete dataset would leak information from validation or test data. A zero standard deviation indicates a constant feature and requires a defined policy, such as dropping the feature or using a denominator of one after centering.

Approach and Solution

The team stores the feature matrix as a two-dimensional floating-point NumPy array. Missing values are represented consistently and handled before standardization. Training statistics are computed along axis 0:

```python

mu = np.nanmean(X_train, axis=0)

sigma = np.nanstd(X_train, axis=0)

sigma_safe = np.where(sigma == 0, 1.0, sigma)

X_train_scaled = (X_train - mu) / sigma_safe

X_valid_scaled = (X_valid - mu) / sigma_safe

```

The same `mu` and `sigma_safe` are applied to validation data. The operation is vectorized and broadcasts feature statistics across rows. The team verifies shapes, checks for non-finite outputs, records the data type, and fixes a random seed for any sampling used during evaluation.

Outcome

The revised pipeline reduces interpreter overhead, applies a consistent transformation, avoids division by zero, and prevents validation leakage. It also becomes easier to integrate with a formal preprocessing component and to test independently.

Learning Points

The case demonstrates that performance and correctness are connected. Vectorization improves speed, but axis selection, data leakage prevention, missing-value policy, numerical safeguards, and reproducibility determine whether the result is trustworthy.

Common Errors, Risks, and Limitations

Several mistakes recur in numerical Python programs. Using `*` instead of `@` changes matrix multiplication into elementwise multiplication. Omitting `axis` may reduce an entire array when per-feature statistics were intended. Unintended broadcasting may silently produce a result with a plausible shape. Integer division or integer overflow may damage calculations. Comparing floating-point numbers with exact equality can fail because of rounding; use tolerances such as `np.isclose` where appropriate.

NaN and infinity values require explicit policies. Functions such as `np.nanmean` ignore NaNs, but ignoring missing values is not always statistically appropriate. Overflow in exponential functions can be addressed with stable formulations or specialized functions. Subtracting nearly equal large numbers can cause loss of significance.

Memory is another limiting resource. Loading a dataset into a dense array may be impossible, and creating multiple temporary arrays may exhaust available memory. Chunking, memory mapping, sparse representations, and out-of-core processing can help. However, chunked processing requires careful treatment of global statistics and order-dependent algorithms.

Professional and Ethical Considerations

Numerical code used in finance, medicine, public services, or safety-critical systems requires traceability and validation. A reproducible seed does not validate a model. Developers should document data provenance, transformations, assumptions, precision, software versions, and tolerance settings. Numerical errors can create unfair or unsafe decisions when they disproportionately affect particular groups or when missing data is handled inconsistently.

Sensitive datasets must be protected during experimentation. Randomized test data should not be mistaken for anonymization, and computational efficiency should never justify bypassing access controls or privacy requirements. Numerical results should be communicated with appropriate uncertainty rather than presented as exact facts.

Current Developments and Best Practices

Modern scientific computing increasingly combines NumPy-compatible interfaces with accelerators, automatic differentiation, just-in-time compilation, distributed arrays, and specialized tensor libraries. Array APIs are improving portability across CPU, GPU, and other devices. Nonetheless, portability does not eliminate the need to understand shape semantics, data types, memory transfer, and numerical stability.

Recommended practices include:

1. Define the mathematical meaning of every axis.

2. Inspect shapes and dtypes at system boundaries.

3. Prefer vectorized operations when they improve clarity and performance.

4. Use `solve` or least-squares routines instead of explicit matrix inversion.

5. Use local random generators and record seeds for experiments.

6. Measure performance with realistic data.

7. Check finite values, tolerances, and conditioning.

8. Avoid unnecessary copies and temporary arrays.

9. Separate training statistics from evaluation statistics.

10. Test both typical and pathological inputs.

Chapter Summary

NumPy provides the array-based foundation of scientific Python and much of AI and ML computation. Its multidimensional arrays, ufuncs, reductions, indexing mechanisms, broadcasting rules, and linear algebra routines allow mathematical operations to be expressed compactly and executed efficiently. Correct numerical programming depends on understanding shape, axis, dtype, views, copies, memory layout, and precision.

Vectorization generally reduces Python interpreter overhead, but it must be used with attention to memory and semantics. Broadcasting is powerful but can conceal shape errors. Linear algebra should use numerically appropriate routines such as `solve` and least squares rather than unnecessary matrix inversion. NumPy’s random generators support controlled sampling and reproducible experiments, while SciPy extends the ecosystem with optimization, integration, statistics, sparse computation, and other numerical methods. In professional AI and ML work, efficiency must be balanced with numerical stability, data governance, reproducibility, and ethical responsibility.

Key Terms

Array: A structured multidimensional collection of elements, usually sharing one data type.

Axis: A dimension of an array along which an operation may be performed.

Broadcasting: Rules that allow compatible arrays of different shapes to participate in one operation.

Vectorization: Applying an operation to an entire array without an explicit Python-level element loop.

Universal function (ufunc): A NumPy function that performs elementwise operations on array values.

Dtype: The data type used to represent array elements.

Shape: The size of an array along each of its dimensions.

Stride: The number of bytes used to move along an array dimension in memory.

View: An array representation sharing data with another array.

Copy: An independent array containing separately stored data.

Reduction: An operation, such as sum or mean, that aggregates values along one or more axes.

Matrix multiplication: Algebraic multiplication of compatible matrices, represented by `@`.

Condition number: A measure of how sensitively a problem’s solution responds to input perturbations.

Floating-point precision: The degree of numerical detail representable by a floating-point type.

Pseudorandom generator: A deterministic algorithm producing sequences with statistical properties resembling randomness.

Numerical stability: The ability of an algorithm to limit the amplification of computational errors.

Sparse matrix: A matrix in which most entries are zero and stored using a specialized representation.

Broadcast-compatible: Describing shapes that satisfy NumPy’s dimension compatibility rules.

Vectorized reduction: Aggregation of array values along a specified axis or axes.

Numerical integration: Approximation of a definite integral using sampled function values.

Review Questions

Short Answer Questions

1. What information is provided by an ndarray’s `shape`, `ndim`, `size`, and `dtype`?

2. Distinguish between a view and a copy in NumPy.

3. State the basic compatibility rule used by broadcasting.

4. Why is `np.linalg.solve` generally preferable to explicitly computing a matrix inverse?

5. What role does a random seed play in computational experiments?

Long Answer Questions

1. Explain vectorization, including its performance advantages and situations in which an explicit loop may remain appropriate.

2. Describe NumPy indexing, slicing, Boolean indexing, and advanced indexing, including their implications for copying data.

3. Explain broadcasting with reference to feature standardization and the shapes `(n,)`, `(n, 1)`, and `(1, n)`.

4. Discuss the major numerical and computational factors that influence the reliability of linear algebra operations.

5. Describe how NumPy and SciPy complement one another in a scientific-computing workflow.

Analytical / Application Questions

1. A feature matrix has shape `(5000, 20)`. Its mean is computed with `axis=1`, but the intended operation is feature-wise standardization. Diagnose the error and propose the correct computation.

2. A vectorized expression is slower than a loop for a small array but faster for a large array. Explain why this can occur and identify the measurements needed before choosing an implementation.

3. Design a reproducible preprocessing strategy for a train-validation dataset containing missing values, constant features, and randomly sampled observations. Include the treatment of statistics, random generators, and numerical checks.

22. Multiple Choice Questions with Answers

Question 1

Which expression performs matrix multiplication in NumPy?

A. A * B

B. A @ B

C. A / B

D. A + B

Correct Answer: B. A @ B

Question 2

For an array X with shape (100, 5), which operation computes one mean for each feature?

A. X.mean(axis=0)

B. X.mean(axis=1)

C. X.mean(axis=2)

D. X.mean(axis=None, keepdims=True)

Correct Answer: A. X.mean(axis=0)

Question 3

What is the main performance advantage of vectorization?

A. It eliminates all memory use

B. It changes floating-point arithmetic into exact arithmetic

C. It moves repeated operations from Python-level loops to optimized array code

D. It guarantees that every algorithm becomes linear time

Correct Answer: C. It moves repeated operations from Python-level loops to optimized array code

Question 4

Which pair of shapes is broadcasting-compatible?

A. (3, 4) and (3, 2)

B. (2, 3) and (4,)

C. (3, 1) and (4,)

D. (5, 2) and (3, 2)

Correct Answer: C. (3, 1) and (4,)

Question 5

Why might modifying a basic slice unexpectedly modify the original NumPy array?

A. Basic slices commonly return views sharing the original data

B. NumPy automatically broadcasts modifications globally

C. All NumPy arrays are immutable

D. Slicing always converts arrays to lists

Correct Answer: A. Basic slices commonly return views sharing the original data

Question 6

Which method is generally preferable for solving Ax = b?

A. np.linalg.inv(A) @ b

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

C. A * b

D. np.linalg.det(A) + b

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

Question 7

In a train-validation standardization pipeline, which statistics should normally be used to transform validation data?

A. Statistics computed separately from validation data

B. Statistics computed from the combined train and validation data

C. Statistics computed from the training data

D. Random statistics generated for each batch

Correct Answer: C. Statistics computed from the training data

Question 8

What does setting a seed for a pseudorandom generator primarily provide?

A. A guarantee that samples are statistically perfect

B. A reproducible sequence when the same conditions are used

C. True physical randomness

D. Higher numerical precision

Correct Answer: B. A reproducible sequence when the same conditions are used

Question 9

Which representation is usually appropriate when most entries of a large matrix are zero?

A. A dense floating-point matrix only

B. A sparse matrix format

C. A Python string

D. A scalar array

Correct Answer: B. A sparse matrix format

Question 10

Why can explicit matrix inversion be problematic in numerical computing?

A. It is always impossible for square matrices

B. It can be less stable and less efficient than solving the system directly

C. It changes every matrix into a sparse matrix

D. It prevents all use of floating-point numbers

Correct Answer: B. It can be less stable and less efficient than solving the system directly