Six questions, three languages

Everything here is used in every lecture that follows. Most of it is used far beyond this course.

Every Idea, Written Twice

This session is a translation exercise. Every symbol you meet has a line of code, and every line of code has a symbol.

And a third language, the first one: the sentence you would say out loud. If you cannot say it in English, the symbols will not save you. The next slide says one idea all three ways.

Python in six linesEnough to read every code slide in this lecture

import numpy as np              # a library: arrays and linear algebra
x = np.array([2.0, 3.0, 1.0])   # a vector
print(x.shape)                  # (3,): one axis holding three numbers
A = np.array([[1, 2], [3, 4]])  # a matrix: a list of rows
print(A.shape)                  # (2, 2)
print(A @ A, A * A)             # matrix product, then elementwise
    

Press Run, then change a number and predict the output before you run it again. (3,), with its trailing comma, is how Python writes a shape with one axis.

The loop that predicts the future

Try it: the loopWorld population, 1950 to 2023

import pandas as pd

P = pd.read_csv('world_population.csv').set_index('year').population
r = (P[1970] / P[1950]) ** (1 / 20) - 1   # growth per year, 1950 to 1970

p = P[1970]
for year in range(1970, 2023):
    p = p * (1 + r)                        # next year: this year, grown by r
print(f'r = {r:.2%}: the loop says {p / 1e9:.2f} billion in 2023')
print(f'the world had {P[2023] / 1e9:.2f} billion')
    

About 1.99%, 10.48 billion predicted, 8.09 billion real. The loop was fine; the rate changed.

Vectors

Vector Operations

import numpy as np

x = np.array([2, 3, 1])
y = np.array([1, -2, 4])
print(x.shape, len(x))

# addition and scaling
print(x + y)
print(3 * x)

# the inner product
print(np.dot(x, y))
print(np.linalg.norm(x))

The inner product \( x^\top y = \sum_i x_i y_i \) is the single most common operation in this course.

How Long, How Far, How Similar

Try it: vectors and matricesChange a number, predict the output, then press Run

import numpy as np

x = np.array([2, 3, 1])
y = np.array([1, -2, 4])
print(x + y, 3 * x, np.dot(x, y), np.linalg.norm(x))

A = np.array([[1, 2, 3],
              [4, 5, 6]])
print(A * A)      # elementwise
print(A @ A.T)    # matrix product: (2, 3) times (3, 2)
    

Try A @ A. The error message is the shape rule speaking: the inner dimensions must match.

Matrices

Matrices in NumPy

A = np.array([[1, 2, 3],
              [4, 5, 6]])
print(A.shape)   # (2, 3)
print(A.size)    # 6

B = np.zeros((3, 3))
C = np.ones((2, 4))
D = np.eye(3)
E = np.random.rand(2, 3)

# the two products people confuse
print(A * A)        # elementwise
print(A @ A.T)      # matrix product

A matrix is a machine that moves space

Two machines in a row are one machine

Two Ways to Read \( Ax \)

Span, Independence, Rank

• Span. Everything you can build from a set of vectors by scaling and adding. Two independent vectors in the plane span the plane; two parallel ones span only a line.
• Linear independence. No vector in the set is a combination of the others, which is the same as saying only the zero recipe reaches zero. \( Xv = 0 \;\Longrightarrow\; v = 0 \)
• Rank. The number of independent columns, so the dimension of the column space. It can never exceed \( \min(m, n) \).

Inverse, Determinant, and When Things Break

• Inverse. The matrix that, when multiplied by \(A\), is the identity \(I\), the matrix version of 1. It solves a linear system. \( A A^{-1} = A^{-1} A = I, \qquad Ax = b \;\Rightarrow\; x = A^{-1} b \)
• Determinant. The factor by which \(A\) scales volume. Zero means the map flattens space, and nothing can unflatten it. \( \det(A) = 0 \;\Longleftrightarrow\; A \text{ has no inverse} \)
• When an inverse exists. Square, and full rank: no column is a combination of the others.
• Why you will care. The normal equations need this inverse. Duplicate a feature and it stops existing. \( \theta = (X^\top X)^{-1} X^\top y \)

To see a zero determinant, press "flatten" on "A matrix is a machine that moves space": the plane collapses to a line, and nothing can unflatten it.

When the Inverse Quietly Lies

Real measurements, 342 penguins, four features plus an intercept column. Now add one more column: bill length again, in centimetres.

X = np.column_stack([np.ones(342), bill_len, bill_dep,
                     flipper, mass])

np.linalg.matrix_rank(X)      # 5 of 5 columns
np.linalg.cond(X.T @ X)       # 1.8e+10

# the same feature, in different units
X2 = np.column_stack([X, bill_len / 10])

np.linalg.matrix_rank(X2)     # 5 of 6 columns
np.linalg.cond(X2.T @ X2)     # 7.1e+21

# and this does not raise
theta = np.linalg.inv(X2.T @ X2) @ X2.T @ y
np.linalg.det(X2.T @ X2)      # -7.9e+09, not 0

Singular: no inverse. Condition number: how close to singular (1 is perfect; past 1e+16, numerically singular). lstsq: least squares without an inverse.

Try it: the inverse that lies342 penguins, and the same column twice

X = np.column_stack([np.ones(342), bill_len, bill_dep, flipper, mass])
X2 = np.column_stack([X, bill_len / 10])       # bill length again, in cm

print(np.linalg.matrix_rank(X), np.linalg.matrix_rank(X2))
print(f'{np.linalg.cond(X.T @ X):.1e}  {np.linalg.cond(X2.T @ X2):.1e}')
print(f'det = {np.linalg.det(X2.T @ X2):.1e}')  # not zero, and not trustworthy

theta = np.linalg.inv(X2.T @ X2) @ X2.T @ y     # no error is raised
print(theta.round(2))
    

Rank 5, then 5 of 6. The second condition number is astronomical and the determinant, which should be zero, is a large number whose exact value depends on your machine: on a laptop 7.1e+21 and −7.9e+09, in this browser something else. That instability is the lesson. The inverse runs anyway and returns nonsense; replace the last line with np.linalg.lstsq(X2, y, rcond=None) and compare.

A derivative is what a curve looks like up close

A rate and its totalBeirut's rain, one year at a time

Try it: a derivative and an integral, in two wordsBeirut, September 2019 to August 2020

import json, numpy as np

d = json.load(open('beirut_weather.json'))['daily']
dates = np.array(d['time'], dtype='datetime64[D]')
rain = np.array(d['precipitation_sum'])              # mm per day: a rate

year = (dates >= np.datetime64('2019-09-01')) & (dates < np.datetime64('2020-09-01'))
total = np.cumsum(rain[year])                        # the integral: rain so far
print(round(total[-1]))                              # the whole year, in mm
print(np.allclose(np.diff(total), rain[year][1:]))   # its derivative is the rate
    

You should see 932 and True. Change the two dates on the year line to see another rain year; then load d["temperature_2m_max"] the same way and find, with np.diff, the day with the biggest jump.

The Calculus You Actually Need

• Derivative. For \( J: \mathbb{R} \to \mathbb{R} \), the slope: how much the output moves when the input nudges.
• Gradient. For \( J: \mathbb{R}^d \to \mathbb{R} \), one partial derivative per parameter, collected into a vector that points uphill. \( \nabla_\theta J = \left[ \frac{\partial J}{\partial \theta_1}, \dots, \frac{\partial J}{\partial \theta_d} \right]^\top \in \mathbb{R}^d \)

The One Gradient You Will Use All Semester

Least squares in matrix form, differentiated once, slowly. Here \( X \in \mathbb{R}^{m \times d} \), \( y \in \mathbb{R}^{m} \), \( \theta \in \mathbb{R}^{d} \), and \( J(\theta) \) is a single number.

Read it as a sentence: \( X\theta - y \) is the residual, and \( X^\top \) asks each feature how much of that residual it can still explain.

Least Squares Is a Projection

Written with an inverse; computed with np.linalg.lstsq.

The chain rule, as a pictureHow every neural network computes its gradient

Try it: backpropagation in twenty linesA number that remembers how it was made

class Value:
    def __init__(self, v, parents=(), local=()):
        self.v, self.grad = v, 0.0
        self.parents, self.local = parents, local   # d(self)/d(parent)

    def __add__(self, o): return Value(self.v + o.v, (self, o), (1, 1))
    def __sub__(self, o): return Value(self.v - o.v, (self, o), (1, -1))
    def __mul__(self, o): return Value(self.v * o.v, (self, o), (o.v, self.v))

    def backward(self, g=1.0):
        # chain rule: multiply along each path back, add over paths
        self.grad += g
        for p, d in zip(self.parents, self.local):
            p.backward(g * d)
        

x, y = Value(2.0), Value(5.0)
w, b = Value(1.0), Value(0.5)
e = w * x + b - y
L = e * e
L.backward()
print('L =', L.v, ' dL/dw =', w.grad, ' dL/db =', b.grad)
print('by hand: 2e*x =', 2 * e.v * x.v, ' 2e =', 2 * e.v)
h = 1e-6                                     # and by nudging w a little
print('nudge:', ((((1 + h) * 2 + 0.5 - 5) ** 2) - L.v) / h)
        

A variable is a label, not a box

Broadcasting, and the Bug It Causes

mu = X.mean(axis=0)     # (4,)
sd = X.std(axis=0)      # (4,)
Z  = (X - mu) / sd      # (342, 4)

Broadcasting, drawn

Drawn with 3 rows; in the penguin data the 3 would be 342 and nothing else changes. Solid cells are the data you have; faint ones are the copies NumPy pretends to make.

Avoid Loops if You Can! Vectorize.

Loops are not the enemy: a loop over time, like the population loop at the start, is how a model steps forward. The loop to avoid is the one over the rows of your data.

Try it: the raceTimed on your own machine, right now

import numpy as np, time

n = 1_000_000
a, b = np.random.rand(n), np.random.rand(n)

t = time.perf_counter()
s = 0.0
for i in range(n):
    s += a[i] * b[i]
loop = time.perf_counter() - t

t = time.perf_counter()
s2 = a @ b
vec = time.perf_counter() - t
print(f'loop {loop * 1000:.0f} ms, a @ b {vec * 1000:.2f} ms: {loop / vec:.0f} times faster')
    

Python here runs inside your browser, so the loop is slower than on a laptop; the gap is the lesson. Run it twice: the second time is faster, because the first run warms things up.

Reading and Exploring Data

import pandas as pd

df = pd.read_csv("penguins.csv")
print(df.shape)

print(df.head(3))
print(df.describe())
print(df.isnull().sum())

# the four numeric columns, as one matrix
cols = ["bill_length_mm", "bill_depth_mm",
        "flipper_length_mm", "body_mass_g"]
X = df[cols].dropna().values
print(X.shape, X.dtype)

Always describe() and check for missing values before modelling anything. Two of these 344 penguins were never measured, and dropna() is a decision, not a formality.

Try it: read the penguinsThe first four lines of every analysis

import pandas as pd

df = pd.read_csv('penguins.csv')
print(df.shape)
print(df.isnull().sum())
print(df.describe().round(1))
    

Two birds have no measurements and eleven have no sex recorded. dropna() is a decision, not a formality: which rows would it remove here?

The Design Matrix

Every dataset in this course arrives in this one shape: \( m \) examples down, \( d \) features across, and the thing you want to predict pulled out into its own column.

Statistics Is Linear Algebra

Palmer Station penguins (Gorman et al.), CC0, vendored with the deck.

Look at It Before You Model It

Anscombe's quartet, 1973

Two Measurements, 342 Birds

The Same Data, Two Answers

The Lab: Five Notebooks

• 00 Python intro: syntax, collections, functions
• 01 Vectors: creation, arithmetic, inner products, norms
• 02 Matrices: indexing, slicing, the two products
• 03 Matrix properties and calculus: inverse, determinant, rank, gradients
• 04 Data reading and plotting: pandas, summary statistics, figures

Summary

Next: supervised learning and linear regression, where all of this becomes one line of code.