NumPy and Vectorization

Let's take a break from the trilogy. This post has no new model to train, no cost to compute, no hill to walk down. It's a tool post, like pulling into a gas station to swap a bald tire before it blows out on the highway. It exists because starting in Week 2 of the course the model gains several features at once (size, bedroom count, house age, all together), and at that point computing with a for loop stops being viable. Let's understand the tool that solves this before we actually need it.

Python lists work, so why switch

Python already ships with lists, which store numbers just fine. So why does the entire ML community use a separate library for this? Short answer: memory.

A Python list stores each number as a separate object, scattered around memory, with a bunch of extra bookkeeping stuck to it (reference count, and so on). Doing math on top of that means the processor keeps jumping from address to address, unwrapping each object to get the value inside. A NumPy array, on the other hand, stores raw numbers of the same type, packed right next to each other in one contiguous block of memory. The processor can crunch a whole bunch of them at once, using special instructions that exist exactly for this.

That difference in memory layout is the real reason behind NumPy's speed. It's not magic, it's low-level engineering.

Vectors: the basic building block

A vector is a bunch of numbers arranged in order, all of the same type. In the course's notation, a vector is a bold lowercase letter, like $\mathbf{x}$. A vector's size is called its shape, written $(n,)$ for a vector with $n$ elements. Notice the lone comma inside the parentheses: that's Python's way of saying "this is a one-element tuple", and it's easy to confuse $(n,)$ with $(n,1)$ later on, so keep an eye on that.

import numpy as np

a = np.zeros(4)
print(f"np.zeros(4): a = {a}, shape = {a.shape}, dtype = {a.dtype}")
# a = [0. 0. 0. 0.], shape = (4,), dtype = float64

a = np.array([5, 4, 3, 2])
print(f"np.array: a = {a}, shape = {a.shape}, dtype = {a.dtype}")
# a = [5 4 3 2], shape = (4,), dtype = int64

A detail that tripped me up as a newcomer: a single decimal-point value in the list is enough for the whole array to become float64. Makes sense, since every element of an array has to share the same type.

Indexing and slicing

Indexing (a[2]) grabs one element, while slicing (a[2:7]) grabs a chunk. The rules are the same as Python lists: counting starts at zero, negative indices count from the end (a[-1] is the last element), and the end of a slice is not included (a[2:7] grabs indices 2, 3, 4, 5, and 6, five elements, not six).

a = np.arange(10)          # [0 1 2 3 4 5 6 7 8 9]

print(a[2])                 # 2, grabbing one element = a scalar
print(a[-1])                 # 9, the last element
print(a[2:7:1])              # [2 3 4 5 6], from index 2 through 6
print(a[3:])                  # [3 4 5 6 7 8 9], from index 3 to the end
print(a[:3])                   # [0 1 2], from the start through index 2

Operations with no loop at all

This is where the good part lives. Vector operations run over the entire array at once, with no for written by you:

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

print(-a)          # [-1 -2 -3 -4], flips the sign of everyone at once
print(np.sum(a))   # 10, sums it all up
print(a ** 2)       # [1 4 9 16], squares every element
print(5 * a)          # [5 10 15 20], multiplies every element by 5

Notice that last example: the 5 was just a loose number, but NumPy "stretched" it on its own to match the array's 4 elements. That has a name, broadcasting, and it'll show up a lot from here on.

That's the whole point of this lab: you describe the operation over the entire array, and let NumPy figure out how to apply it to each element under the hood. If you catch yourself writing a for to walk through an array, there's probably a vectorized equivalent waiting to be used instead.

The dot product: the operation linear regression runs on

The dot product is this lab's most important operation, because it's literally what's going to replace that w * x[i] + b we wrote with a loop in the previous posts, once the model gains several features.

$\mathbf{a} \cdot \mathbf{b} = \sum_{i=0}^{n-1} a_i b_i$

Multiply pairwise, sum it all up. That simple. Before reaching for NumPy's built-in version, I wrote my own, just to make what's happening under the hood explicit:

def my_dot(a, b):
    x = 0
    for i in range(a.shape[0]):
        x = x + a[i] * b[i]
    return x

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

print(my_dot(a, b))       # 24
print(np.dot(a, b))       # 24, same result, NumPy's built-in function

The speed gap, live in your browser

The original notebook tests this with 10-million-element arrays, comparing np.dot (vectorized) against a hand-written for loop, in Python. The gain there runs dozens to hundreds of times over.

Below, you can run a similar test, except in your own browser, right now, in JavaScript. I need to be honest with you about one thing first, though: the browser is not going to show a 100x gain. And it's not because the test is broken.

Your browser's JavaScript engine (V8, if you're on Chrome or Edge) already optimizes a plain for loop in a way pretty similar to what it does for built-in methods like .reduce(). When I compared "loop" against ".reduce()" on the same kind of array, it came out close to a tie, and that was expected, not a bug in the test.

The comparison that actually shows a real, honest gap is a different one: a plain JavaScript Array (which stores numbers somewhat "boxed up", similar to the Python list we talked about above) versus a Float64Array (contiguous memory, single type, no boxing). That's the same underlying reason NumPy is fast, just reproduced here in the browser with a smaller, honest gain (usually 2 to 5x, sometimes more), instead of me pretending I could recreate Python's giant number in an environment that wasn't built for that.

Run it a few times with different sizes. Notice that "Array comum + laço" (plain Array + loop) tends to be visibly slower than the other two, and that "Float64Array + laço" (typed array + loop) and "Float64Array + .reduce()" land close to each other. That's exactly the pattern the explanation above predicted.

Matrices: when one vector isn't enough

A matrix is a two-dimensional array, written with a bold uppercase letter ($\mathbf{X}$), with shape $(m, n)$. In the course's context the convention is always the same: row is a training example, column is a feature. One house per row, one attribute (size, bedrooms, age) per column.

X = np.array([[1, 5], [2, 3], [3, 1]])   # shape (3, 2): 3 examples, 2 features each

print(X.shape)      # (3, 2)
print(X[1])          # [2 3], the entire row 1, becomes a 1-D vector
print(X[1, 0])         # 2, one specific element, becomes a scalar
print(X[:, 0])           # [1 2 3], the entire column 0

The detail that confuses newcomers the most here: indexing a matrix with just one index (X[1]) returns an array with one fewer dimension, not a single-row matrix. X[1] has shape (2,), not (1, 2). This specific gotcha is behind a good chunk of the dimension errors you'll run into later in the course, and I've learned to flag it mentally every time I touch shapes.

And this is exactly why the dot product comes back into play: once the model has $n$ features, each row of the matrix $\mathbf{X}$ becomes a vector of shape $(n,)$, ready to take a dot product directly with the weight vector $\mathbf{w}$, also of shape $(n,)$. Example $i$'s prediction becomes np.dot(w, X[i]) + b, a single line of code, no loop at all, whether $n$ is 1 or 100.

Wrapping up

How this connects to the rest of the playlist: with these tools in hand, the model you built in Lab 02 with a hand-written w * x + b becomes np.dot(w, x) + b, working for any number of features without rewriting anything. This exact foundation is what the course uses to extend everything you've already seen (model, cost, gradient descent) to linear regression with multiple variables.

Practical application

Same real housing dataset from the previous posts (Housing Prices Regression, Kaggle). This post doesn't train anything, so here the idea is just applying the matrix/vector mechanics on top of real column names, instead of the toy a and b examples.

import pandas as pd
import numpy as np

df = pd.read_csv("real_estate_dataset.csv")

X = df[["Square_Feet", "Num_Bedrooms"]].to_numpy()[:3]   # first 3 houses, 2 features

print("X.shape:", X.shape)        # (3, 2): 3 examples, 2 features
print("X[1]:", X[1])              # the entire second house, becomes a 1-D vector
print("X[:, 0]:", X[:, 0])        # the entire Square_Feet column, all houses

Output: X.shape: (3, 2) / X[1]: [55.15 5.] / X[:, 0]: [143.64 55.15 202.96]

# A real dot product, with arbitrary weights just to illustrate the mechanics
# (we're not training anything in this post, this isn't the "right" fit):
w = np.array([2000, 10000])

for i in range(X.shape[0]):
    print(f"np.dot(w, X[{i}]) = {np.dot(w, X[i]):,.0f}")

Output: np.dot(w, X[0]) = 297,280 / np.dot(w, X[1]) = 160,300 / np.dot(w, X[2]) = 455,920

Same math, same mechanics from the rest of the post, just running on top of real Square_Feet and Num_Bedrooms instead of [1, 2, 3, 4]. Actually training these weights for real (finding the $\mathbf{w}$ that makes sense for the data) is content further ahead in the course, with multiple features at once.

Picking the matrix columns, live

Remember $\mathbf{X}$ is just a matrix, one row per example, one column per feature? Pick which two columns become the $x$ and $z$ axes below (height always stays price) and rotate to feel how each column pairing relates to the price of the 50 houses:

Loading real data...