6  Arrays, Shapes, and Vectorization

6.1 Learning objectives

After completing this unit, you will be able to:

• create NumPy arrays with a stated shape and dtype;
• explain the relationship between indices, axes, and array shape;
• replace simple componentwise loops with vectorized operations;
• use broadcasting only after checking shape compatibility;
• distinguish a memory-sharing view from an independent copy;
• formulate invariants for shapes, values, and the state of inputs; and
• limit experimental conclusions so they are not mistaken for general proofs.

Local prerequisites: Units 1-2, simple linear functions, lists and tuples, and coordinates in the plane. This unit assumes neither calculus nor formal linear algebra.

6.2 Arrays as data with a contract

An array is more than a collection of numbers. Its minimal contract includes:

1. shape (shape): the number of positions along each axis;
2. number of dimensions (ndim);
3. representation type (dtype); and
4. the meaning of each axis, which we must document ourselves.

import numpy as np

temperature = np.array(
    [[28.0, 29.5, 30.0], [27.5, 29.0, 31.0]],
    dtype=np.float64,
)
(temperature.shape, temperature.ndim, temperature.dtype)

Listing 6.1

((2, 3), 2, dtype('float64'))

The shape (2, 3) says only that there are two rows and three columns. It does not say whether the rows represent cities, days, or repeated trials. That meaning is part of the data specification.

A one-dimensional array with shape (3,) differs from a single-row table with shape (1, 3) and a single-column table with shape (3, 1). All three contain three numbers, but operations along their axes can behave differently.

6.3 dtype is a representation, not a mathematical set

NumPy stores the elements of an ordinary array in a single dtype. An explicit choice makes the contract easier to check:

integers = np.array([1, 2, 3], dtype=np.int64)
approximations = np.array([0.1, 0.2, 0.3], dtype=np.float64)
(integers.dtype, approximations.dtype, approximations.sum())

Listing 6.2

(dtype('int64'), dtype('float64'), np.float64(0.6000000000000001))

int64 stores integers within a finite range. Operations that go outside that range can overflow; it is not the same as the set of all mathematical integers. float64 stores binary approximations with finite precision; it is not the same as the real numbers.

NumPy can use dtype=object for Python objects such as integers of arbitrary precision, subject to practical limits, or Fraction values, but many benefits of numerical array computation are lost. If the question requires exact arithmetic, deliberately choose an exact tool (such as Fraction, SymPy, or Sage) instead of calling a float64 result exact.

6.4 Vectorization and a scalar reference

Suppose each data value $x_i$ is transformed by the function

$$y_i=3x_i-2.$$

A clear scalar loop can serve as a reference:

x = np.array([1, 2, 3, 4], dtype=np.int64)
y_reference = np.array([3 * int(value) - 2 for value in x])
y_reference

Listing 6.3

array([ 1,  4,  7, 10])

NumPy lets us write the same rule for the entire array:

y = 3 * x - 2
y

Listing 6.4

array([ 1,  4,  7, 10])

Vectorization expresses an array operation without a Python loop that we write ourselves. It is often more concise and can be faster, but brevity is not a proof of correctness. Compare the result with a simple reference and check invariants derived from the formula:

• y.shape == x.shape;
• x is unchanged; and
• $\sum_i y_i=3\sum_i x_i-2n$, where $n$ is the number of elements.

The scalar reference must also be checked. Two implementations that copy the same conceptual error can agree and still be wrong.

6.5 Broadcasting with shape checks

Broadcasting extends operations between arrays without requiring us to manually create tiled copies. Suppose each column of a table receives a different offset:

table = np.array([[10, 20, 30], [40, 50, 60]], dtype=np.int64)
offset = np.array([1, -2, 3], dtype=np.int64)
table + offset

Listing 6.5

array([[11, 18, 33],
       [41, 48, 63]])

The shapes (2, 3) and (3,) are compatible because their last axes both have size three. The offsets are applied to every row. Do not rely on successful execution as the only check: an unintended shape may still broadcast and produce a plausible-looking answer.

The add_column_offsets function in the Unit 3 code specifies a stricter contract: the table must be two-dimensional, the offsets must be one-dimensional, and the number of offsets must equal the number of columns.

6.6 Copies and views

A NumPy array slice is usually a view: a new object that accesses the original array’s memory. Changing the view can change the original data.

original = np.arange(6)
view = original[1:4]
view[0] = -10
original

Listing 6.6

array([  0, -10,   2,   3,   4,   5])

Use .copy() when an experiment needs an independent value:

original = np.arange(6)
independent_copy = original[1:4].copy()
independent_copy[0] = -10
original

Listing 6.7

array([0, 1, 2, 3, 4, 5])

Make this choice deliberately. Views are useful for avoiding large copies; copies are useful for preventing unexpected changes. np.shares_memory(a, b) can help check the relationship in tests.

6.7 The Unit 3 experiment

Run from the project root:

python source/code/unit03_arrays.py --output output/unit03-results.json

The script records the NumPy version, shapes, dtype, vectorized results, broadcasting results, view and copy traces, floating-point limitations, and the limits of what the conclusions establish. JSON is written as UTF-8 with sorted keys, fixed indentation, and LF line endings. In the same environment, the same invocation produces identical bytes.

6.8 Computation and proof

For the array [1, 2, 3, 4], the program checks that the vectorized operation agrees with the scalar reference. This result supports a bounded claim about the case and implementation that were run. It does not prove equality for every array, every dtype, or every size.

For arithmetic without overflow, the general argument proceeds component by component. At index $i$, both implementations compute the same expression, $3x_i-2$. Since the index $i$ is arbitrary, equality holds at every component; the shape is also preserved because exactly one output is created for each input.

The argument must be qualified by the representation domain. With int64, overflow can change the relationship to mathematical integers. With float64, rounding can cause algebraically equivalent arrangements of operations to produce different bits. Tests help expose these effects; the specification and argument determine what is actually being claimed.

6.9 Exercises

6.9.1 Exercise 1 - shape and axis meanings

An array stores measurements from three sensors at four times and has shape (4, 3). Explain the meaning of both axes. What shapes do data[2, :] and data[:, 1] produce?

NoteHint

A single index removes the selected axis, whereas : retains all positions along that axis.

TipAnswer and solution

Axis 0 represents the four times, and axis 1 represents the three sensors. data[2, :] selects all sensors at the third time and therefore has shape (3,). data[:, 1] selects the second sensor at all times and therefore has shape (4,). This meaning comes from the data contract, not the shape alone.

6.9.2 Exercise 2 - vectorization and invariants

Vectorize a loop that computes $z_i=5x_i+1$. State two testable invariants and derive the sum invariant.

NoteHint

Apply the operation directly to x, then sum the component equations.

TipAnswer and solution

Write z = 5*x + 1. Two useful invariants are z.shape == x.shape and that x is unchanged. If there are $n$ elements, then

$$\sum_i z_i=\sum_i(5x_i+1)=5\sum_i x_i+n.$$

Tests can compare the vectorized result with a scalar reference and check all three properties on selected cases.

6.9.3 Exercise 3 - a broadcasting contract

Determine whether the following shape pairs are compatible for addition intended to mean “one offset per column”: (5, 3) with (3,), (5, 3) with (5,), and (5, 3) with (1, 3). Explain each answer.

NoteHint

Compare dimensions from the right. Their sizes must be equal, or one of them must be one.

TipAnswer and solution

(5, 3) with (3,) is compatible and applies three offsets to every row. (5, 3) with (5,) is incompatible because the last sizes, 3 and 5, differ. (5, 3) with (1, 3) is compatible; the first axis, of size one, is expanded to five. Although two pairs are compatible under NumPy’s rules, the function contract must still specify that the second input really represents column offsets.

6.9.4 Exercise 4 - tracing a view

For the following code, determine the final values of a, b, and c without running it.

a = np.array([0, 1, 2, 3])
b = a[1:3]
c = a[1:3].copy()
b[0] = 8
c[1] = 9

NoteHint

b shares memory with a, whereas c does not.

TipAnswer and solution

The final arrays are a: [0, 8, 2, 3], b: [8, 2], and c: [1, 9]. The change through b affects the same element in a. The copy c was created before the change, and changes to c do not propagate back to a. To test an entire array for equality, use a check such as np.array_equal(a, [0, 8, 2, 3]). In contrast, a == [0, 8, 2, 3] returns an array of elementwise Boolean comparisons, not a single Boolean.

6.9.5 Exercise 5 - experiments are not general proofs

You generate a million arrays of small numbers and find that 3*x-2 always agrees with the loop reference. State the strongest justified conclusion, then outline a general proof and give one limitation of the representation.

NoteHint

Separate the sample’s coverage, the componentwise argument, and the behavior of finite dtype representations.

TipAnswer and solution

The strongest empirical conclusion is that no difference was found in the million cases and environment tested. For a general proof in the specified arithmetic domain, choose an arbitrary index $i$: both methods compute $3x_i-2$, so all their components are equal and the shape is preserved. A limitation is that int64 can overflow and float64 rounds; therefore, any relationship to integer or real arithmetic must state the range and representation model.

6.10 Summary

• shape, ndim, dtype, and axis meanings together form the array contract.
• Vectorization expresses componentwise operations but still needs a reference and invariants.
• Broadcasting is safe to teach and use when both shape and axis intent are checked.
• Slices are usually views; .copy() breaks the memory-sharing relationship.
• Canonical output records experiments that can be compared byte for byte.
• Checking many cases provides computational evidence, not a replacement for a general argument; dtype limitations must be part of the claim.