NumPy Array Arithmetic

Module 3: NumPy

NumPy Array Arithmetic

In the last session, you stored the marks of the midterm in a NumPy array: one row per student, one column per exercise, every mark between 0 and 1. You also stored the value of each exercise, in points. But the most important number is still missing: the final grade of each student.

Today we will compute it, and then use it to analyse the exam. Along the way, we will learn how NumPy does arithmetic with whole arrays at once. This array-based style is the foundation of vector programming.

import numpy as np

Here is the data, after the revisions of the last session.

students = ["Ana", "Bruno", "Chloe", "Diego", "Emma", "Farid"]

# Points of each exercise: 2 open problems, 3 debugging, 6 multiple choice
weights = np.array([20, 20, 10, 10, 10, 5, 5, 5, 5, 5, 5])

marks = np.array([
    [0.75, 1, 0.5, 1, 0.5, 1, 1, 0, 1, 1, 1],         # Ana
    [0.5, 0.25, 1, 0.5, 0, 1, 0, 1, 1, 0, 1],         # Bruno
    [1, 1, 1, 1, 0.5, 1, 1, 1, 0, 1, 1],              # Chloe
    [0.25, 0.5, 0.5, 0.5, 0.5, 0, 1, 0, 1, 0, 1],     # Diego
    [1, 1, 0.5, 1, 1, 1, 0, 1, 1, 1, 1],              # Emma
    [0.5, 0.5, 0.5, 0, 1, 1, 1, 0, 0, 1, 0],          # Farid
])

The Problem: Grading with Loops

The points of a student in one exercise are her mark times the value of the exercise. Her final score is the sum over all exercises, and dividing by 10 gives a grade out of 10. With what we know from Modules 1 and 2, that means two nested loops:

final_grades = []

for row in marks:
    total = 0
    for j in range(len(weights)):
        total = total + row[j] * weights[j]
    final_grades.append(float(total / 10))

print(final_grades)

[8.0, 5.0, 9.0, 4.5, 9.0, 5.0]

It works, but the code hides the idea behind indices and accumulators. Let’s see how NumPy expresses the same calculation.

Element-Wise Operations

When two arrays have the same shape, the arithmetic operators act element by element: the first element with the first element, the second with the second, and so on. Let’s take one student.

ana = marks[0]

# Multiplication element-wise: points earned in each exercise
ana_points = ana * weights
print(ana_points)

[15. 20.  5. 10.  5.  5.  5.  0.  5.  5.  5.]

The same holds for the other operators. Suppose that two professors graded Ana’s open problems and debugging exercises independently.

grader_1 = np.array([0.75, 1, 0.5, 1, 0.5])
grader_2 = np.array([0.5, 1, 0.5, 0.5, 0.5])

# Addition and division by a number: the average of both graders
print((grader_1 + grader_2) / 2)

[0.625 1.    0.5   0.75  0.5  ]

# Subtraction: where do the graders disagree?
print(grader_1 - grader_2)

[0.25 0.   0.   0.5  0.  ]

# Division: back from points to marks
print(ana_points / weights)

[0.75 1.   0.5  1.   0.5  1.   1.   0.   1.   1.   1.  ]

# Power element-wise
print(grader_1 ** 2)

[0.5625 1.     0.25   1.     0.25  ]

Notice that / 2 divides every element by 2. An operation between an array and a single number (a scalar) is applied to every element of the array.

To add up all the elements of an array, use np.sum(). Ana’s final grade is now one line:

ana_grade = np.sum(ana * weights) / 10
print(ana_grade)

8.0

Exercise: Points Lost

For each exercise, how many points did Ana lose? Compute it without a loop. Then check that the points lost plus the points earned add up to 100.

Solution

ana_lost = (1 - ana) * weights
print(ana_lost)

print(np.sum(ana_lost) + np.sum(ana * weights))

[5. 0. 5. 0. 5. 0. 0. 5. 0. 0. 0.]
100.0

1 - ana subtracts every mark from 1, which gives the fraction of each exercise that was missed.

Exercise: A Perfect Exam

Create an array with the marks of a perfect exam, without typing them, and compute its grade out of 10.

Solution

perfect = np.ones(11)
print(np.sum(perfect * weights) / 10)

10.0

Broadcasting in NumPy

We have the grade of one student. What about the whole class? marks has shape (6, 11), while weights has shape (11,). They do not have the same shape… and yet:

points = marks * weights
print(points)

[[15. 20.  5. 10.  5.  5.  5.  0.  5.  5.  5.]
 [10.  5. 10.  5.  0.  5.  0.  5.  5.  0.  5.]
 [20. 20. 10. 10.  5.  5.  5.  5.  0.  5.  5.]
 [ 5. 10.  5.  5.  5.  0.  5.  0.  5.  0.  5.]
 [20. 20.  5. 10. 10.  5.  0.  5.  5.  5.  5.]
 [10. 10.  5.  0. 10.  5.  5.  0.  0.  5.  0.]]

NumPy multiplied every row of marks by weights. This mechanism is called broadcasting: it allows NumPy to work with arrays of different shapes when performing arithmetic operations. We have already used its simplest form, a scalar with an array:

# Every mark, in percentage
print(marks * 100)

[[ 75. 100.  50. 100.  50. 100. 100.   0. 100. 100. 100.]
 [ 50.  25. 100.  50.   0. 100.   0. 100. 100.   0. 100.]
 [100. 100. 100. 100.  50. 100. 100. 100.   0. 100. 100.]
 [ 25.  50.  50.  50.  50.   0. 100.   0. 100.   0. 100.]
 [100. 100.  50. 100. 100. 100.   0. 100. 100. 100. 100.]
 [ 50.  50.  50.   0. 100. 100. 100.   0.   0. 100.   0.]]

How does broadcasting work?

NumPy compares the shapes of the arrays element-wise, starting from the trailing (rightmost) dimensions.
If the dimensions are equal, or one of them is 1, they are considered compatible.
NumPy then “stretches” the smaller array along the dimension with size 1 (or the missing dimension) to match the larger array.

In our case:

marks     (6, 11)
weights       (11,)   -> treated as (1, 11), then stretched to (6, 11)
result    (6, 11)

Now consider a different situation. Chloe and Farid delivered their exams 10 minutes late, and the rules say their marks are multiplied by 0.9. We store one factor per student:

late_factor = np.array([1, 1, 0.9, 1, 1, 0.9])
print(late_factor.shape)

(6,)

penalized = marks * late_factor

ValueError: operands could not be broadcast together with shapes (6,11) (6,)

The trailing dimensions are 11 and 6: they are different and neither is 1, so the shapes are incompatible. We need late_factor to be a column, with shape (6, 1), so that it is stretched along the exercises. A reshape is enough:

late_column = late_factor.reshape(6, 1)
print(late_column)

[[1. ]
 [1. ]
 [0.9]
 [1. ]
 [1. ]
 [0.9]]

penalized = marks * late_column
print(penalized)

[[0.75 1.   0.5  1.   0.5  1.   1.   0.   1.   1.   1.  ]
 [0.5  0.25 1.   0.5  0.   1.   0.   1.   1.   0.   1.  ]
 [0.9  0.9  0.9  0.9  0.45 0.9  0.9  0.9  0.   0.9  0.9 ]
 [0.25 0.5  0.5  0.5  0.5  0.   1.   0.   1.   0.   1.  ]
 [1.   1.   0.5  1.   1.   1.   0.   1.   1.   1.   1.  ]
 [0.45 0.45 0.45 0.   0.9  0.9  0.9  0.   0.   0.9  0.  ]]

marks        (6, 11)
late_column  (6, 1)   -> stretched to (6, 11)
result       (6, 11)

Exercise: Which Shapes Broadcast?

Without running any code, decide which of these arrays can be added to marks (shape (6, 11)). For those that work, what is the shape of the result?

A scalar, 0.1
An array with shape (11,)
An array with shape (6,)
An array with shape (1, 11)
An array with shape (6, 1)
An array with shape (2, 11)

Solution

Works, (6, 11): a scalar is applied to every element.
Works, (6, 11): one value per exercise, repeated for every student.
Fails: 6 and 11 are compared and are different.
Works, (6, 11): the 1 is stretched to 6 rows.
Works, (6, 11): the 1 is stretched to 11 columns.
Fails: 2 and 6 are different and neither is 1.

Practice

Deactivate AI assistant tools, and try the following exercises.

Exercise: Scale and Broadcast

Compute the following operation, where the products are element-wise:

\[5 \cdot A \cdot v\]

Where \[A=\begin{pmatrix} 1 & 2 & 3 & 4\\ 5 & 6 & 7 & 8\\ 9 & 10 & 11 & 12 \end{pmatrix}\] and $v = [-2, -1, 0, 1]$

Solution

A = np.array([
  [1, 2, 3, 4],
  [5, 6, 7, 8],
  [9, 10, 11, 12]
])
v = np.array([-2, -1, 0, 1])

result = 5 * A * v
print(result)

[[-10 -10   0  20]
 [-50 -30   0  40]
 [-90 -50   0  60]]

Exercise: Column Weights

Given any matrix $A$, multiply the elements of its first column times 1, the second column times 2, the third times 3, and so on.

For instance:

\[A = \begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ \end{pmatrix}\]

Would end up being: \[A' = \begin{pmatrix} 1 & 4 & 9\\ 4 & 10 & 18\\ \end{pmatrix}\]

Your code should work for a matrix with any number of columns.

Solution

A = np.array([
  [1, 2, 3],
  [4, 5, 6]
])
column_weights = np.arange(1, A.shape[1] + 1)

result = A * column_weights
print(result)

[[ 1  4  9]
 [ 4 10 18]]

This is exactly what we did with marks * weights.

Aggregation Functions

We have the table of points, with shape (6, 11). np.sum() adds all its elements, which is not what we want: that would be the total points of the whole class.

print(np.sum(points))

405.0

The axis argument chooses the dimension that is “collapsed”. To add the points of each student, we sum along the exercises, which is dimension 1 (the columns). One line, and we have the final grades of the class:

final_grades = np.sum(points, axis=1) / 10
print(final_grades)

[8.  5.  9.  4.5 9.  5. ]

Summing along dimension 0 (the rows) adds the points of every student in each exercise:

points_per_exercise = np.sum(points, axis=0)
print(points_per_exercise)

[80. 85. 40. 40. 35. 25. 20. 15. 20. 20. 25.]

A good way to remember it: the axis you pass is the one that disappears from the shape. points has shape (6, 11); axis=1 leaves (6,), one value per student, and axis=0 leaves (11,), one value per exercise.

The cumulative sum keeps every partial total. How did Ana’s score grow as she went through the exam?

print(np.cumsum(points[0]))

[15. 35. 40. 50. 55. 60. 65. 65. 70. 75. 80.]

Other aggregations work the same way. For instance, np.prod() multiplies all the elements:

print(np.prod(np.array([[1, 2], [3, 4]])))

Exercise: The Three Parts

Compute, for every student, the points obtained in each part of the exam: open problems (exercises 1 and 2), debugging (3 to 5), and multiple choice (6 to 11). Each result should be an array with one value per student.

Solution

open_points = np.sum(points[:, :2], axis=1)
debugging_points = np.sum(points[:, 2:5], axis=1)
mc_points = np.sum(points[:, 5:], axis=1)

print(open_points)
print(debugging_points)
print(mc_points)

[35. 15. 40. 15. 40. 20.]
[20. 15. 25. 15. 25. 15.]
[25. 20. 25. 15. 25. 15.]

We slice the columns of each part first, and then collapse the exercises with axis=1.

Exercise: Passing Point

Using the cumulative sum, find after which exercise Ana had already reached 50 points (a pass). You can read it from the output.

Solution

print(np.cumsum(points[0]))

[15. 35. 40. 50. 55. 60. 65. 65. 70. 75. 80.]

The fourth value is 50, so Ana had reached a pass after exercise 4. In a later session, we will learn how to let NumPy find that position for us.

Sorting and Searching

Which student got the highest grade? np.argmax() returns the index of the largest value, which we can use with our list of names.

best = np.argmax(final_grades)
print(best, students[best])

2 Chloe

Chloe and Emma both got a 9. When there is a tie, np.argmax() returns the first one. np.argmin() works the same way for the smallest value.

worst = np.argmin(final_grades)
print(worst, students[worst])

3 Diego

np.argsort() returns the indices that would sort the array, from the lowest to the highest value:

order = np.argsort(final_grades)
print(order)

[3 1 5 0 2 4]

Exercise: The Hardest Exercise

Which exercise was the hardest for the class? Compute the average mark of each exercise, then find the exercise number (not the index) with the lowest average.

Solution

average_per_exercise = np.mean(marks, axis=0)
print(average_per_exercise)

hardest = np.argmin(average_per_exercise) + 1
print(f"The hardest exercise was exercise {hardest}")

[0.66666667 0.70833333 0.66666667 0.66666667 0.58333333 0.83333333
 0.66666667 0.5        0.66666667 0.66666667 0.83333333]
The hardest exercise was exercise 8

We add 1 because exercise 1 is at index 0.

Exercise: Ranking

Print the names of the students from the highest to the lowest grade.

Solution

ranking = np.argsort(final_grades)[::-1]

for i in ranking:
    print(students[i], final_grades[i])

Emma 9.0
Chloe 9.0
Ana 8.0
Farid 5.0
Bruno 5.0
Diego 4.5

np.argsort() sorts from low to high, so we reverse it with the slice [::-1].

Statistical Functions

How did the class do as a whole? NumPy computes the usual descriptive statistics.

The mean is the average of a data set.

print("Mean:", np.mean(final_grades))

Mean: 6.75

The median is the middle value when the data set is ordered.

print("Median:", np.median(final_grades))

Median: 6.5

The standard deviation measures the dispersion of data from its mean.

print("Standard Deviation:", np.std(final_grades))

Standard Deviation: 1.9525624189766635

Variance is the square of the standard deviation.

print("Variance:", np.var(final_grades))

Variance: 3.8125

They accept the axis argument too. This is the average mark of each exercise we used before:

print(np.mean(marks, axis=0))

[0.66666667 0.70833333 0.66666667 0.66666667 0.58333333 0.83333333
 0.66666667 0.5        0.66666667 0.66666667 0.83333333]

Exercise: Grade Report

Compute the mean, median, mode, standard deviation, and variance of the final grades. NumPy has no mode function: look at np.unique() with the argument return_counts=True.

Solution

print("Mean:", np.mean(final_grades))
print("Median:", np.median(final_grades))
# Count each unique value and pick the most frequent
values, counts = np.unique(final_grades, return_counts=True)
print("Mode:", values[np.argmax(counts)])
print("Standard Deviation:", np.std(final_grades))
print("Variance:", np.var(final_grades))

Mean: 6.75
Median: 6.5
Mode: 5.0
Standard Deviation: 1.9525624189766635
Variance: 3.8125

Two values appear twice (5 and 9). np.argmax() returns the first of them, so the mode reported is 5.0.

Mathematical Functions

NumPy also provides the usual mathematical functions, applied element by element. Some professors curve the grades with a square root, which raises low grades more than high ones:

\[\text{curved} = 10 \cdot \sqrt{\frac{\text{grade}}{10}}\]

curved = 10 * np.sqrt(final_grades / 10)
print(curved)

[8.94427191 7.07106781 9.48683298 6.70820393 9.48683298 7.07106781]

# For higher roots, we can just operate as we do with arrays
print(final_grades ** (1 / 3))

[2.         1.70997595 2.08008382 1.65096362 2.08008382 1.70997595]

So why use np.sqrt if we can do **(1/2)?

np.sqrt() clearly communicates the intended operation and uses NumPy’s optimized array implementation. The expression arr ** (1 / 2) is also valid. The following optional timing experiment should be run manually rather than whenever the notes are rendered.

arr = np.arange(1, 10000)  # We use a very big array

for _ in range(10000):
  np.sqrt(arr)

# We are not printing anything, because we do not want to fill the output
# We only want to see how much time this takes!

for _ in range(10000):
  arr**(1/2)

# We are not printing anything, because we do not want to fill the output
# We only want to see how much time this takes!

The same logic applies to exponentiation!

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

# Exponentiation
print(np.exp(arr))

[  2.71828183   7.3890561   20.08553692  54.59815003 148.4131591 ]

Geometry Functions

NumPy can do more than just element-wise operations. It also supports matrix multiplication, transposition, and other matrix math.

Dot products: For one-dimensional vectors, np.dot() computes the usual dot product.

\[\mathbf{u} \cdot \mathbf{v} = \sum_i u_i v_i\]

Look closely at the formula: multiply element by element, then add. That is exactly how we computed Ana’s score!

print(np.dot(ana, weights))
print(np.sum(ana * weights))

80.0
80.0

# Vector dot product
u = np.array([1, 2, 3])
v = np.array([4, 5, 6])
print(np.dot(u, v))

Matrix multiplication:

\[A\ B = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{pmatrix} = \begin{pmatrix} a_{11}b_{11} + a_{12}b_{21} & a_{11}b_{12} + a_{12}b_{22} \\ a_{21}b_{11} + a_{22}b_{21} & a_{21}b_{12} + a_{22}b_{22} \end{pmatrix}\]

# Matrix multiplication (option A)
A = np.array([[1, 2], [3, 4]])
B = np.array([[4, 3], [2, 1]])

C = np.matmul(A, B)
print(C)

[[ 8  5]
 [20 13]]

# Matrix multiplication (option B)
C = A @ B
print(C)

[[ 8  5]
 [20 13]]

Every element of a matrix product is the dot product of a row of the first matrix with a column of the second. So the product of the table of marks, (6, 11), with the vector of weights, (11,), gives one dot product per student: the scores of the whole class.

print(marks @ weights / 10)

[8.  5.  9.  4.5 9.  5. ]

For two-dimensional arrays, np.dot(A, B) performs matrix multiplication. Prefer the @ operator when that is what you mean because it makes the operation explicit.

If you want the sum of element-wise products of two matrices—the Frobenius inner product—write it directly:

# Sum of element-wise products
frobenius_product = np.sum(A * B)
print(frobenius_product)

Both dot and matmul are used for matrix multiplication in NumPy, but they behave differently, especially when dealing with higher-dimensional arrays (i.e., tensors). The primary difference emerges in multi-dimensional array multiplication.

The official documentation says: matmul differs from dot in two important ways.

Multiplication by scalars is not allowed.
Stacks of matrices are broadcast together as if the matrices were elements.

Cross product: A cross product is defined for vectors with two or three components. With a stack of vectors, NumPy applies the operation to corresponding rows.

For two-dimensional vectors, the result is the scalar $u_1v_2-u_2v_1$. For three-dimensional vectors, the result is another three-dimensional vector perpendicular to both inputs.

# Create two stacks of 2D vectors
u = np.array([[1, 2], [3, 4]])
v = np.array([[4, 3], [2, 1]])

# One scalar cross product per row
cross_products = np.cross(u, v)
print(cross_products)

ValueError: Both input arrays must be (arrays of) 3-dimensional vectors, but they are 2 and 2 dimensional instead.

Transposition:

\[A^T = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}^T = \begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \end{pmatrix}\]

a = np.array([[1, 2, 3], [4, 5, 6]])
print(a)

# Matrix transposition
b = a.T
print(b)

[[1 2 3]
 [4 5 6]]
[[1 4]
 [2 5]
 [3 6]]

Determinants provide information about the matrix’s invertibility, while the inverse of a matrix is crucial for solving systems of linear equations.

Exercise: Matrix Product with a Transpose

Let $A$ be the $3 \times 4$ matrix containing the integers from 1 to 12, ordered by rows. Compute the following matrix operation:

\[A B^T / 6\]

Where

\[B=\begin{pmatrix} -1 & 0 & 0 & -1\\ 0 & 1 & 1 & 0\\ 1 & -1 & 1 & -1 \end{pmatrix}\]

Solution

A = np.arange(1, 13).reshape(3, 4)
B = np.array([
    [-1, 0, 0, -1],
    [0, 1, 1, 0],
    [1, -1, 1, -1]
])

result = (A @ B.T) / 6
print(result)

[[-0.83333333  0.83333333 -0.33333333]
 [-2.16666667  2.16666667 -0.33333333]
 [-3.5         3.5        -0.33333333]]

Linear Algebra Functions

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

det_a = np.linalg.det(a)

print(det_a)

-2.0000000000000004

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

inverse_a = np.linalg.inv(a)

print(inverse_a)

[[-2.   1. ]
 [ 1.5 -0.5]]

Trigonometric Functions

NumPy operates in radians. If you are using degrees, make sure to convert them first!

angles = np.array([0, 30, 45, 60, 90])  # in degrees

# NumPy operates in radians!
angles_rad = np.radians(angles)

print(angles_rad)

[0.         0.52359878 0.78539816 1.04719755 1.57079633]

sines = np.sin(angles_rad)

print(sines)

[0.         0.5        0.70710678 0.8660254  1.        ]

cosines = np.cos(angles_rad)

print(cosines)

[1.00000000e+00 8.66025404e-01 7.07106781e-01 5.00000000e-01
 6.12323400e-17]

tangents = np.tan(angles_rad)

print(tangents)

[0.00000000e+00 5.77350269e-01 1.00000000e+00 1.73205081e+00
 1.63312394e+16]

Practice

Deactivate AI assistant tools, and try the following exercises.

Exercise: Normalize a Matrix

Given a “magic square” matrix $A$:

\[A = \begin{pmatrix} 23 & 28 & 21\\ 22 & 24 & 26\\ 27 & 20 & 25 \end{pmatrix}\]

Normalize its values (subtract mean and divide by standard deviation). We call the new matrix $B$.

\[B[i,j] = \frac{A[i,j] - \text{mean}(A)}{\text{std}(A)}\]

Solution

arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]])
b = (arr - np.mean(arr)) / np.std(arr)
print(b)

[[-0.38729833  1.54919334 -1.161895  ]
 [-0.77459667  0.          0.77459667]
 [ 1.161895   -1.54919334  0.38729833]]

Exercise: Transpose the Result

Compute $B^T$, using the matrix $B$ from the previous exercise.

Solution

arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]])
b = (arr - np.mean(arr)) / np.std(arr)
print(b.T)

[[-0.38729833 -0.77459667  1.161895  ]
 [ 1.54919334  0.         -1.54919334]
 [-1.161895    0.77459667  0.38729833]]

Exercise: Matrix Product

Compute the matrix product $A B^T$.

Solution

arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]])
b = (arr - np.mean(arr)) / np.std(arr)
print(arr @ b.T)

[[10.0697567  -1.54919334 -8.52056336]
 [-1.54919334  3.09838668 -1.54919334]
 [-8.52056336 -1.54919334 10.0697567 ]]

Exercise: Magic Square Properties

What makes $A$ a “magic square”? Compute the sum of each row, the sum of each column, and the sum of the main diagonal.

Solution

Every row, every column, and the diagonal add up to the same number: 72.

arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]])

print(np.sum(arr, axis=0))
print(np.sum(arr, axis=1))
print(np.trace(arr))

[72 72 72]
[72 72 72]
72

Challenge: The Ambiguous Question

After publishing the grades, several students complain that exercise 8, the hardest one, had an ambiguous statement. You agree, and decide that everyone gets full marks in that exercise. Starting from marks and weights, and without writing any loop:

Create a copy of the marks with fixed = marks.copy() and give everyone a 1 in exercise 8.
Compute the new final grades with a matrix product.
Compute how much each student’s grade went up, and the name of the student who gained the most.
Compare the class mean before and after the correction.
Apply the square-root curve to the new grades and print each student’s name with their curved grade, rounded to two decimals with np.round().

Possible solution

fixed = marks.copy()
fixed[:, 7] = 1  # exercise 8 is at index 7

new_grades = fixed @ weights / 10
print(new_grades)

gain = new_grades - final_grades
print(gain)
print(students[np.argmax(gain)])

print(np.mean(final_grades), np.mean(new_grades))

curved = np.round(10 * np.sqrt(new_grades / 10), 2)
for i in range(len(students)):
    print(students[i], curved[i])

[8.5 5.  9.  5.  9.  5.5]
[0.5 0.  0.  0.5 0.  0.5]
Ana
6.75 7.0
Ana 9.22
Bruno 7.07
Chloe 9.49
Diego 7.07
Emma 9.49
Farid 7.42

Several students gain 0.5 points; np.argmax() reports the first one. The class mean goes up by the average gain.

Further Practice

You will find more exercises here!

Summary and Reflection

Arithmetic operators act element by element on compatible arrays; an operation with a scalar applies to every element.
Broadcasting compares shapes from right to left and allows dimensions that match or equal 1.
The axis argument controls the direction of an aggregation: the axis you pass disappears from the shape.
* is element-wise multiplication; @ is matrix multiplication. A weighted sum is a dot product.

Given arrays with shapes (4, 3) and (3,), explain why addition works. Would addition still work if the second shape were (4,)?

We now know the grade of each student, but not yet who passed. In a later session, we will learn to select the elements of an array that satisfy a condition.

Homework and Review

Work through the homework exercises available here.

Homework

--- title: "NumPy Array Arithmetic" subtitle: "Module 3: NumPy" format: html execute: error: true --- # NumPy Array Arithmetic In the last session, you stored the marks of the midterm in a NumPy array: one row per student, one column per exercise, every mark between `0` and `1`. You also stored the value of each exercise, in points. But the most important number is still missing: the **final grade** of each student. Today we will compute it, and then use it to analyse the exam. Along the way, we will learn how NumPy does arithmetic with whole arrays at once. This array-based style is the foundation of vector programming. ```{python} import numpy as np ``` Here is the data, after the revisions of the last session. ```{python} students = ["Ana", "Bruno", "Chloe", "Diego", "Emma", "Farid"] # Points of each exercise: 2 open problems, 3 debugging, 6 multiple choice weights = np.array([20, 20, 10, 10, 10, 5, 5, 5, 5, 5, 5]) marks = np.array([ [0.75, 1, 0.5, 1, 0.5, 1, 1, 0, 1, 1, 1], # Ana [0.5, 0.25, 1, 0.5, 0, 1, 0, 1, 1, 0, 1], # Bruno [1, 1, 1, 1, 0.5, 1, 1, 1, 0, 1, 1], # Chloe [0.25, 0.5, 0.5, 0.5, 0.5, 0, 1, 0, 1, 0, 1], # Diego [1, 1, 0.5, 1, 1, 1, 0, 1, 1, 1, 1], # Emma [0.5, 0.5, 0.5, 0, 1, 1, 1, 0, 0, 1, 0], # Farid ]) ``` ## The Problem: Grading with Loops The points of a student in one exercise are her mark times the value of the exercise. Her final score is the sum over all exercises, and dividing by 10 gives a grade out of 10. With what we know from Modules 1 and 2, that means two nested loops: ```{python} final_grades = [] for row in marks: total = 0 for j in range(len(weights)): total = total + row[j] * weights[j] final_grades.append(float(total / 10)) print(final_grades) ``` It works, but the code hides the idea behind indices and accumulators. Let's see how NumPy expresses the same calculation. ## Element-Wise Operations When two arrays have the same shape, the arithmetic operators act **element by element**: the first element with the first element, the second with the second, and so on. Let's take one student. ```{python} ana = marks[0] # Multiplication element-wise: points earned in each exercise ana_points = ana * weights print(ana_points) ``` The same holds for the other operators. Suppose that two professors graded Ana's open problems and debugging exercises independently. ```{python} grader_1 = np.array([0.75, 1, 0.5, 1, 0.5]) grader_2 = np.array([0.5, 1, 0.5, 0.5, 0.5]) ``` ```{python} # Addition and division by a number: the average of both graders print((grader_1 + grader_2) / 2) ``` ```{python} # Subtraction: where do the graders disagree? print(grader_1 - grader_2) ``` ```{python} # Division: back from points to marks print(ana_points / weights) ``` ```{python} # Power element-wise print(grader_1 ** 2) ``` Notice that `/ 2` divides **every** element by 2. An operation between an array and a single number (a scalar) is applied to every element of the array. To add up all the elements of an array, use `np.sum()`. Ana's final grade is now one line: ```{python} ana_grade = np.sum(ana * weights) / 10 print(ana_grade) ``` ::: {.callout-note} ## Exercise: Points Lost For each exercise, how many points did Ana **lose**? Compute it without a loop. Then check that the points lost plus the points earned add up to 100. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} ana_lost = (1 - ana) * weights print(ana_lost) print(np.sum(ana_lost) + np.sum(ana * weights)) ``` `1 - ana` subtracts every mark from 1, which gives the fraction of each exercise that was missed. ::: ::: {.callout-note} ## Exercise: A Perfect Exam Create an array with the marks of a perfect exam, without typing them, and compute its grade out of 10. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} perfect = np.ones(11) print(np.sum(perfect * weights) / 10) ``` ::: --- ## Broadcasting in NumPy We have the grade of one student. What about the whole class? `marks` has shape `(6, 11)`, while `weights` has shape `(11,)`. They do not have the same shape... and yet: ```{python} points = marks * weights print(points) ``` NumPy multiplied **every row** of `marks` by `weights`. This mechanism is called **broadcasting**: it allows NumPy to work with arrays of different shapes when performing arithmetic operations. We have already used its simplest form, a scalar with an array: ```{python} # Every mark, in percentage print(marks * 100) ``` How does broadcasting work? - NumPy compares the shapes of the arrays element-wise, starting from the trailing (rightmost) dimensions. - If the dimensions are equal, or one of them is 1, they are considered compatible. - NumPy then "stretches" the smaller array along the dimension with size 1 (or the missing dimension) to match the larger array. In our case: ```text marks (6, 11) weights (11,) -> treated as (1, 11), then stretched to (6, 11) result (6, 11) ``` Now consider a different situation. Chloe and Farid delivered their exams 10 minutes late, and the rules say their marks are multiplied by `0.9`. We store one factor per student: ```{python} late_factor = np.array([1, 1, 0.9, 1, 1, 0.9]) print(late_factor.shape) ``` ```{python} penalized = marks * late_factor ``` The trailing dimensions are `11` and `6`: they are different and neither is `1`, so the shapes are incompatible. We need `late_factor` to be a **column**, with shape `(6, 1)`, so that it is stretched along the exercises. A reshape is enough: ```{python} late_column = late_factor.reshape(6, 1) print(late_column) ``` ```{python} penalized = marks * late_column print(penalized) ``` ```text marks (6, 11) late_column (6, 1) -> stretched to (6, 11) result (6, 11) ``` ::: {.callout-note} ## Exercise: Which Shapes Broadcast? Without running any code, decide which of these arrays can be added to `marks` (shape `(6, 11)`). For those that work, what is the shape of the result? 1. A scalar, `0.1` 2. An array with shape `(11,)` 3. An array with shape `(6,)` 4. An array with shape `(1, 11)` 5. An array with shape `(6, 1)` 6. An array with shape `(2, 11)` ::: ::: {.callout-tip collapse="true"} ## Solution 1. Works, `(6, 11)`: a scalar is applied to every element. 2. Works, `(6, 11)`: one value per exercise, repeated for every student. 3. Fails: `6` and `11` are compared and are different. 4. Works, `(6, 11)`: the `1` is stretched to 6 rows. 5. Works, `(6, 11)`: the `1` is stretched to 11 columns. 6. Fails: `2` and `6` are different and neither is `1`. ::: ### Practice Deactivate AI assistant tools, and try the following exercises. ::: {.callout-note} ## Exercise: Scale and Broadcast Compute the following operation, where the products are element-wise: $$5 \cdot A \cdot v$$ Where $$A=\begin{pmatrix} 1 & 2 & 3 & 4\\ 5 & 6 & 7 & 8\\ 9 & 10 & 11 & 12 \end{pmatrix}$$ and $v = [-2, -1, 0, 1]$ ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} A = np.array([ [1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12] ]) v = np.array([-2, -1, 0, 1]) result = 5 * A * v print(result) ``` ::: ::: {.callout-note} ## Exercise: Column Weights Given any matrix $A$, multiply the elements of its first column times 1, the second column times 2, the third times 3, and so on. For instance: $$A = \begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ \end{pmatrix}$$ Would end up being: $$A' = \begin{pmatrix} 1 & 4 & 9\\ 4 & 10 & 18\\ \end{pmatrix}$$ Your code should work for a matrix with any number of columns. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} A = np.array([ [1, 2, 3], [4, 5, 6] ]) column_weights = np.arange(1, A.shape[1] + 1) result = A * column_weights print(result) ``` This is exactly what we did with `marks * weights`. ::: --- ## Aggregation Functions We have the table of `points`, with shape `(6, 11)`. `np.sum()` adds all its elements, which is not what we want: that would be the total points of the whole class. ```{python} print(np.sum(points)) ``` The `axis` argument chooses the dimension that is "collapsed". To add the points of each student, we sum **along the exercises**, which is dimension 1 (the columns). One line, and we have the final grades of the class: ```{python} final_grades = np.sum(points, axis=1) / 10 print(final_grades) ``` Summing along dimension 0 (the rows) adds the points of every student in each exercise: ```{python} points_per_exercise = np.sum(points, axis=0) print(points_per_exercise) ``` A good way to remember it: the `axis` you pass is the one that **disappears** from the shape. `points` has shape `(6, 11)`; `axis=1` leaves `(6,)`, one value per student, and `axis=0` leaves `(11,)`, one value per exercise. The cumulative sum keeps every partial total. How did Ana's score grow as she went through the exam? ```{python} print(np.cumsum(points[0])) ``` Other aggregations work the same way. For instance, `np.prod()` multiplies all the elements: ```{python} print(np.prod(np.array([[1, 2], [3, 4]]))) ``` ::: {.callout-note} ## Exercise: The Three Parts Compute, for every student, the points obtained in each part of the exam: open problems (exercises 1 and 2), debugging (3 to 5), and multiple choice (6 to 11). Each result should be an array with one value per student. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} open_points = np.sum(points[:, :2], axis=1) debugging_points = np.sum(points[:, 2:5], axis=1) mc_points = np.sum(points[:, 5:], axis=1) print(open_points) print(debugging_points) print(mc_points) ``` We slice the columns of each part first, and then collapse the exercises with `axis=1`. ::: ::: {.callout-note} ## Exercise: Passing Point Using the cumulative sum, find after which exercise Ana had already reached 50 points (a pass). You can read it from the output. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} print(np.cumsum(points[0])) ``` The fourth value is `50`, so Ana had reached a pass after exercise 4. In a later session, we will learn how to let NumPy find that position for us. ::: --- ## Sorting and Searching Which student got the highest grade? `np.argmax()` returns the **index** of the largest value, which we can use with our list of names. ```{python} best = np.argmax(final_grades) print(best, students[best]) ``` Chloe and Emma both got a 9. When there is a tie, `np.argmax()` returns the first one. `np.argmin()` works the same way for the smallest value. ```{python} worst = np.argmin(final_grades) print(worst, students[worst]) ``` `np.argsort()` returns the indices that would sort the array, from the lowest to the highest value: ```{python} order = np.argsort(final_grades) print(order) ``` ::: {.callout-note} ## Exercise: The Hardest Exercise Which exercise was the hardest for the class? Compute the average mark of each exercise, then find the **exercise number** (not the index) with the lowest average. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} average_per_exercise = np.mean(marks, axis=0) print(average_per_exercise) hardest = np.argmin(average_per_exercise) + 1 print(f"The hardest exercise was exercise {hardest}") ``` We add `1` because exercise 1 is at index 0. ::: ::: {.callout-note} ## Exercise: Ranking Print the names of the students from the highest to the lowest grade. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} ranking = np.argsort(final_grades)[::-1] for i in ranking: print(students[i], final_grades[i]) ``` `np.argsort()` sorts from low to high, so we reverse it with the slice `[::-1]`. ::: --- ## Statistical Functions How did the class do as a whole? NumPy computes the usual descriptive statistics. The **mean** is the average of a data set. ```{python} print("Mean:", np.mean(final_grades)) ``` The **median** is the middle value when the data set is ordered. ```{python} print("Median:", np.median(final_grades)) ``` The **standard deviation** measures the dispersion of data from its mean. ```{python} print("Standard Deviation:", np.std(final_grades)) ``` **Variance** is the square of the standard deviation. ```{python} print("Variance:", np.var(final_grades)) ``` They accept the `axis` argument too. This is the average mark of each exercise we used before: ```{python} print(np.mean(marks, axis=0)) ``` ::: {.callout-note} ## Exercise: Grade Report Compute the mean, median, mode, standard deviation, and variance of the final grades. NumPy has no `mode` function: look at `np.unique()` with the argument `return_counts=True`. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} print("Mean:", np.mean(final_grades)) print("Median:", np.median(final_grades)) # Count each unique value and pick the most frequent values, counts = np.unique(final_grades, return_counts=True) print("Mode:", values[np.argmax(counts)]) print("Standard Deviation:", np.std(final_grades)) print("Variance:", np.var(final_grades)) ``` Two values appear twice (5 and 9). `np.argmax()` returns the first of them, so the mode reported is `5.0`. ::: --- ## Mathematical Functions NumPy also provides the usual mathematical functions, applied element by element. Some professors **curve** the grades with a square root, which raises low grades more than high ones: $$\text{curved} = 10 \cdot \sqrt{\frac{\text{grade}}{10}}$$ ```{python} curved = 10 * np.sqrt(final_grades / 10) print(curved) ``` ```{python} # For higher roots, we can just operate as we do with arrays print(final_grades ** (1 / 3)) ``` So why use `np.sqrt` if we can do `**(1/2)`? `np.sqrt()` clearly communicates the intended operation and uses NumPy's optimized array implementation. The expression `arr ** (1 / 2)` is also valid. The following optional timing experiment should be run manually rather than whenever the notes are rendered. ```{python} #| eval: false arr = np.arange(1, 10000) # We use a very big array for _ in range(10000): np.sqrt(arr) # We are not printing anything, because we do not want to fill the output # We only want to see how much time this takes! ``` ```{python} #| eval: false for _ in range(10000): arr**(1/2) # We are not printing anything, because we do not want to fill the output # We only want to see how much time this takes! ``` The same logic applies to exponentiation! ```{python} arr = np.array([1, 2, 3, 4, 5]) # Exponentiation print(np.exp(arr)) ``` --- ## Geometry Functions NumPy can do more than just element-wise operations. It also supports matrix multiplication, transposition, and other matrix math. **Dot products:** For one-dimensional vectors, `np.dot()` computes the usual dot product. $$\mathbf{u} \cdot \mathbf{v} = \sum_i u_i v_i$$ Look closely at the formula: multiply element by element, then add. That is exactly how we computed Ana's score! ```{python} print(np.dot(ana, weights)) print(np.sum(ana * weights)) ``` ```{python} # Vector dot product u = np.array([1, 2, 3]) v = np.array([4, 5, 6]) print(np.dot(u, v)) ``` **Matrix multiplication:** $$A\ B = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{pmatrix} = \begin{pmatrix} a_{11}b_{11} + a_{12}b_{21} & a_{11}b_{12} + a_{12}b_{22} \\ a_{21}b_{11} + a_{22}b_{21} & a_{21}b_{12} + a_{22}b_{22} \end{pmatrix}$$ ```{python} # Matrix multiplication (option A) A = np.array([[1, 2], [3, 4]]) B = np.array([[4, 3], [2, 1]]) C = np.matmul(A, B) print(C) ``` ```{python} # Matrix multiplication (option B) C = A @ B print(C) ``` Every element of a matrix product is the dot product of a row of the first matrix with a column of the second. So the product of the table of marks, `(6, 11)`, with the vector of weights, `(11,)`, gives one dot product per student: the scores of the whole class. ```{python} print(marks @ weights / 10) ``` For two-dimensional arrays, `np.dot(A, B)` performs matrix multiplication. Prefer the `@` operator when that is what you mean because it makes the operation explicit. If you want the sum of element-wise products of two matrices—the Frobenius inner product—write it directly: ```{python} # Sum of element-wise products frobenius_product = np.sum(A * B) print(frobenius_product) ``` Both `dot` and `matmul` are used for matrix multiplication in NumPy, but they behave differently, especially when dealing with higher-dimensional arrays (i.e., tensors). The primary difference emerges in multi-dimensional array multiplication. The [official documentation](https://numpy.org/doc/stable/reference/generated/numpy.matmul.html) says: `matmul` differs from `dot` in two important ways. - Multiplication by scalars is not allowed. - Stacks of matrices are broadcast together as if the matrices were elements. **Cross product:** A cross product is defined for vectors with two or three components. With a stack of vectors, NumPy applies the operation to corresponding rows. For two-dimensional vectors, the result is the scalar $u_1v_2-u_2v_1$. For three-dimensional vectors, the result is another three-dimensional vector perpendicular to both inputs. ```{python} # Create two stacks of 2D vectors u = np.array([[1, 2], [3, 4]]) v = np.array([[4, 3], [2, 1]]) # One scalar cross product per row cross_products = np.cross(u, v) print(cross_products) ``` **Transposition:** $$A^T = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}^T = \begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \end{pmatrix}$$ ```{python} a = np.array([[1, 2, 3], [4, 5, 6]]) print(a) # Matrix transposition b = a.T print(b) ``` Determinants provide information about the matrix's invertibility, while the inverse of a matrix is crucial for solving systems of linear equations. ::: {.callout-note} ## Exercise: Matrix Product with a Transpose Let $A$ be the $3 \times 4$ matrix containing the integers from 1 to 12, ordered by rows. Compute the following matrix operation: $$A B^T / 6$$ Where $$B=\begin{pmatrix} -1 & 0 & 0 & -1\\ 0 & 1 & 1 & 0\\ 1 & -1 & 1 & -1 \end{pmatrix}$$ ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} A = np.arange(1, 13).reshape(3, 4) B = np.array([ [-1, 0, 0, -1], [0, 1, 1, 0], [1, -1, 1, -1] ]) result = (A @ B.T) / 6 print(result) ``` ::: --- ## Linear Algebra Functions ```{python} a = np.array([[1, 2], [3, 4]]) det_a = np.linalg.det(a) print(det_a) ``` ```{python} a = np.array([[1, 2], [3, 4]]) inverse_a = np.linalg.inv(a) print(inverse_a) ``` --- ## Trigonometric Functions NumPy operates in radians. If you are using degrees, make sure to convert them first! ```{python} angles = np.array([0, 30, 45, 60, 90]) # in degrees # NumPy operates in radians! angles_rad = np.radians(angles) print(angles_rad) ``` ```{python} sines = np.sin(angles_rad) print(sines) ``` ```{python} cosines = np.cos(angles_rad) print(cosines) ``` ```{python} tangents = np.tan(angles_rad) print(tangents) ``` --- ## Practice Deactivate AI assistant tools, and try the following exercises. ::: {.callout-note} ## Exercise: Normalize a Matrix Given a "magic square" matrix $A$: $$A = \begin{pmatrix} 23 & 28 & 21\\ 22 & 24 & 26\\ 27 & 20 & 25 \end{pmatrix}$$ Normalize its values (subtract mean and divide by standard deviation). We call the new matrix $B$. $$B[i,j] = \frac{A[i,j] - \text{mean}(A)}{\text{std}(A)}$$ ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]]) b = (arr - np.mean(arr)) / np.std(arr) print(b) ``` ::: ::: {.callout-note} ## Exercise: Transpose the Result Compute $B^T$, using the matrix $B$ from the previous exercise. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]]) b = (arr - np.mean(arr)) / np.std(arr) print(b.T) ``` ::: ::: {.callout-note} ## Exercise: Matrix Product Compute the matrix product $A B^T$. ::: ::: {.callout-tip collapse="true"} ## Solution ```{python} arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]]) b = (arr - np.mean(arr)) / np.std(arr) print(arr @ b.T) ``` ::: ::: {.callout-note} ## Exercise: Magic Square Properties What makes $A$ a "magic square"? Compute the sum of each row, the sum of each column, and the sum of the main diagonal. ::: ::: {.callout-tip collapse="true"} ## Solution Every row, every column, and the diagonal add up to the same number: 72. ```{python} arr = np.array([[23, 28, 21], [22, 24, 26], [27, 20, 25]]) print(np.sum(arr, axis=0)) print(np.sum(arr, axis=1)) print(np.trace(arr)) ``` ::: --- ## Challenge: The Ambiguous Question After publishing the grades, several students complain that **exercise 8**, the hardest one, had an ambiguous statement. You agree, and decide that everyone gets full marks in that exercise. Starting from `marks` and `weights`, and without writing any loop: 1. Create a copy of the marks with `fixed = marks.copy()` and give everyone a `1` in exercise 8. 2. Compute the new final grades with a matrix product. 3. Compute how much each student's grade went up, and the name of the student who gained the most. 4. Compare the class mean before and after the correction. 5. Apply the square-root curve to the new grades and print each student's name with their curved grade, rounded to two decimals with `np.round()`. ::: {.callout-tip collapse="true"} ## Possible solution ```{python} fixed = marks.copy() fixed[:, 7] = 1 # exercise 8 is at index 7 new_grades = fixed @ weights / 10 print(new_grades) gain = new_grades - final_grades print(gain) print(students[np.argmax(gain)]) print(np.mean(final_grades), np.mean(new_grades)) curved = np.round(10 * np.sqrt(new_grades / 10), 2) for i in range(len(students)): print(students[i], curved[i]) ``` Several students gain 0.5 points; `np.argmax()` reports the first one. The class mean goes up by the average gain. ::: --- # Further Practice [You will find more exercises here!](https://drive.google.com/file/d/1280Tk_Nw6JXPJ_bD4qQhCtyFKKnjpsAK/view?usp=drive_link) --- ## Summary and Reflection - Arithmetic operators act element by element on compatible arrays; an operation with a scalar applies to every element. - Broadcasting compares shapes from right to left and allows dimensions that match or equal `1`. - The `axis` argument controls the direction of an aggregation: the axis you pass disappears from the shape. - `*` is element-wise multiplication; `@` is matrix multiplication. A weighted sum is a dot product. Given arrays with shapes `(4, 3)` and `(3,)`, explain why addition works. Would addition still work if the second shape were `(4,)`? We now know the grade of each student, but not yet **who passed**. In a later session, we will learn to select the elements of an array that satisfy a condition. --- # Homework and Review Work through the homework exercises available here. [Homework](homework_numpy.qmd#array-arithmetic-and-broadcasting){.btn .btn-primary}