Vectors & Matrices from Scratch

Why vectors and matrices?

A spreadsheet row describing a house (area, bedrooms, age), the red/green/blue values of a pixel, the position and velocity of a satellite, the word counts of a document: each is an ordered list of numbers. Mathematicians call such a list a vector. A table of such rows, or a rule that turns one vector into another, is a matrix. Once data is in this shape, a small set of operations (adding, scaling, dot products, multiplication) does a surprising amount of work: it measures similarity between documents, rotates 3-D models, fits lines to data, and powers every layer of a neural network.

This module uses only high-school math: coordinates on a grid, Pythagoras' theorem, and a little trigonometry (cosine). Everything else is built here, from scratch, in pure Python. In M48 you will swap the hand-written loops for NumPy, and you will understand exactly what NumPy is doing for you.

Remember: a vector is an ordered list of numbers; a matrix is a rectangular grid of them. The rest of this module gives those lists a geometric meaning.

Vectors as arrows and as lists

Draw a grid. Start at the origin (0, 0), walk 3 units right and 2 units up. The arrow from where you started to where you stopped is the vector v = (3, 2). Its numbers are called components (or coordinates). The number of components is the vector's dimension: (3, 2) is 2-dimensional, a pixel (255, 128, 0) is 3-dimensional, and a document's word counts over a 10,000-word vocabulary is a 10,000-dimensional vector. We can't draw 10,000 dimensions, but every rule we discover in 2-D works the same way there.

Three basic operations

• Addition v + w: add matching components. Only defined when both vectors have the same dimension. Geometrically: tip to tail.
• Scalar multiplication c · v: multiply every component by a plain number c (a scalar). 2·v points the same way but is twice as long; -1·v points the opposite way; 0.5·v is half as long.
• Magnitude (also called length or norm) |v|: by Pythagoras, the arrow (3, 2) is the hypotenuse of a right triangle with sides 3 and 2, so its length is √(3² + 2²). In any dimension: square every component, add, take the square root.

A vector with length 1 is a unit vector. Dividing any non-zero vector by its own length ("normalizing" it) keeps its direction and makes its length exactly 1. Unit vectors are how we talk about pure direction.

In Python, the most direct representation is a list (or tuple) of floats. Python has no built-in "add these lists element-wise" operator: [3, 2] + [1, 3] concatenates to [3, 2, 1, 3]. So we write small functions with zip.

Remember: add and scale component by component; length is Pythagoras in any number of dimensions; a unit vector is pure direction.

The dot product, angles, and cosine similarity

The dot product of two vectors of the same dimension is: multiply matching components, then add the results. For v = (3, 2) and w = (1, 3): v · w = 3·1 + 2·3 = 9. It takes two vectors and returns one number. That number has a beautiful geometric meaning:

v · w = |v| × |w| × cos(θ)        # θ is the angle between the arrows

You don't need to prove it to use it, but here is the intuition. cos(θ) is 1 when the arrows point the same way (θ = 0°), 0 when they are perpendicular (θ = 90°), and -1 when they point in opposite directions (θ = 180°). So the sign of the dot product tells you whether two vectors broadly agree (positive), are unrelated (zero), or disagree (negative). Note also that v · v = |v|²: the dot product of a vector with itself is its squared length, which is exactly Pythagoras again.

Rearranging the formula gives the angle: cos θ = (v · w) / (|v| |w|). That ratio is called cosine similarity. It ignores length entirely and keeps only direction, which is exactly what you want when comparing documents: a long article and a short summary about the same topic have very different word counts (lengths) but point in a similar direction. Cosine similarity ranges from -1 (opposite) through 0 (unrelated) to 1 (same direction). It is the default way search engines, recommender systems, and embedding models (see M57) compare items.

Remember: v · w = Σ vᵢwᵢ = |v||w|cos θ. Positive means "same general direction", zero means perpendicular, negative means opposed. Cosine similarity is the dot product of the two unit vectors.

Matrices as lists of lists

A matrix is a rectangular grid of numbers with m rows and n columns; we say its shape is m × n (rows first, always). In pure Python, the natural representation is a list of rows, each row a list:

A = [[1, 2, 3],
     [4, 5, 6]]          # shape 2 × 3: 2 rows, 3 columns
rows, cols = len(A), len(A[0])
A[1][2]                  # row 1, column 2 (zero-based) -> 6

Mathematicians write the entry in row i, column j as Aᵢⱼ and count from 1; Python counts from 0. A matrix with one row is a row vector; one with one column is a column vector. You can read a matrix two ways: as a stack of rows (each row a data record, like a spreadsheet), or as a side-by-side list of columns (each column a feature). Both views are useful, and switching between them is what the transpose does.

Transpose

The transpose Aᵀ flips a matrix across its main diagonal: row i becomes column i, so Aᵀ[j][i] == A[i][j] and a 2 × 3 matrix becomes 3 × 2. Python's zip(*A) idiom does exactly this: *A passes each row as a separate argument, and zip groups the first items, then the second items, and so on.

Remember: shape is (rows, columns). M[i][j] is row i, column j. [list(c) for c in zip(*M)] transposes. Never build rows with [row] * n.

Matrix multiplication

Matrix times vector

Multiplying a matrix A (shape m × n) by a vector x (dimension n) gives a new vector of dimension m. Each output entry is the dot product of one row of A with x:

[[1, 2],      [5]     [1·5 + 2·6]     [17]
 [3, 4]]  ×   [6]  =  [3·5 + 4·6]  =  [39]

There is a second, equally important reading: A x is a weighted sum of the columns of A, with weights taken from x. Above, 5·(1, 3) + 6·(2, 4) = (17, 39). The row view is how you compute it; the column view is how you think about it, especially for transformations.

Matrix times matrix

For C = A B, the entry C[i][j] is the dot product of row i of A with column j of B. That requires each row of A to be as long as each column of B: the inner dimensions must match. Shapes combine like dominoes: (m × n) · (n × p) → (m × p). Step through it below and watch which row and column meet for each output cell.

The identity matrix

The identity matrix I has 1s on the main diagonal and 0s elsewhere. It is the matrix version of the number 1: I A = A I = A and I x = x. As a transformation, it leaves every point where it is.

Why AB usually differs from BA

With numbers, 3 × 4 = 4 × 3. With matrices that's almost never true. The run above shows it: B swaps things, and "swap the rows of A" (BA) is a different result from "swap the columns of A" (AB). Shapes make it even clearer: (1×3)(3×1) is 1×1, while (3×1)(1×3) is 3×3. Geometrically, "rotate then stretch" lands points in different places from "stretch then rotate", which you will see in the next section. Matrix multiplication is associative, though: (AB)C = A(BC), and the transpose reverses order: (AB)ᵀ = BᵀAᵀ.

The cost: O(n³)

For two n × n matrices, there are n² output cells, and each one is a dot product of length n: n multiplications and about n additions. Total: n · n · n = n³ multiply-adds. Doubling n makes the work 8 times bigger. A 1000 × 1000 product is a billion multiply-adds. (Clever algorithms such as Strassen's reach about O(n^2.81), but real libraries mostly win with the hardware tricks you will meet below, not by a better exponent.)

Remember: C[i][j] = row i of A · column j of B; shapes (m×n)(n×p) → (m×p); AB ≠ BA in general; I is the do-nothing matrix; naive matmul is O(n³).

Matrices as transformations of the plane

Here's the idea that makes linear algebra click. A 2 × 2 matrix is a machine that moves every point of the plane. Feed it a point (x, y) as a vector, and A (x, y) is where that point goes. Because of the column view, you can read the whole transformation straight off the matrix: the first column is where (1, 0) lands, and the second column is where (0, 1) lands. Every other point follows along, because (x, y) = x·(1, 0) + y·(0, 1).

Such machines are called linear transformations: they keep the origin fixed, keep straight lines straight, and keep evenly spaced grid lines evenly spaced. The common ones:

Now draw it for real. The cell below transforms the corners of a unit square (and an "F" shape, which shows reflections, because an F looks different in a mirror) with several matrices and plots each result. The first import matplotlib in the browser downloads the library, so it takes a few seconds; later runs are fast.

Composition is multiplication

Applying R first and then S to a point is S (R p). Associativity says that equals (S R) p, so the single matrix SR does both steps at once. Read products right to left: in SR, R happens first. This is why games and graphics engines multiply all the transformation matrices together once and then apply one matrix to millions of vertices. And because SR ≠ RS, "rotate then stretch" is a different shape from "stretch then rotate". Change the last transform in the cell to RS and compare.

Remember: the columns of a 2 × 2 matrix are the images of (1, 0) and (0, 1). Multiplying matrices composes transformations, applied right to left.

Memory layout: row-major vs column-major

Computer memory is one long line of numbered bytes, not a grid. To store a 2-D matrix you must choose an order. Row-major order (used by C and by NumPy by default) stores row 0 completely, then row 1, and so on. Column-major order (used by Fortran, MATLAB, R, Julia, and much of the LAPACK library that NumPy and SciPy call) stores column 0 completely, then column 1.

Why it matters: the cache

A CPU can do arithmetic far faster than it can fetch data from main memory. To hide that gap, it keeps small, very fast memories called caches, and it fetches memory in chunks (a cache line, typically 64 bytes, or eight 64-bit floats). If your loop reads numbers that sit next to each other, one fetch serves eight reads. If it jumps around, almost every read misses the cache and waits for main memory, which can be tens of times slower. So: loop in the order the data is stored. For row-major data, the innermost loop should walk along a row.

Pure Python lists blur this picture: a list of lists stores pointers to separate float objects scattered around the heap, so even "consecutive" numbers are not adjacent in memory. You can still see a small difference below (partly cache, partly the extra indexing work), but the dramatic effect appears with contiguous arrays in M48, where the numbers really do sit side by side.

Remember: memory is 1-D; row-major puts rows together, column-major puts columns together. Loop in storage order so the cache can help you.

Why pure Python is slow for this, and what comes next

Your matmul is correct, but time it. Every a * b in a Python loop does far more than one multiplication: the interpreter fetches the next bytecode, looks up the types of both objects, finds the right __mul__, allocates a brand-new float object for the result, updates reference counts, and later frees the old object. Each number is a full Python object (24 bytes for a float, plus an 8-byte pointer in the list), scattered across memory. The actual multiplication is perhaps 1% of the work.

Optimized numerical libraries remove all of that overhead. They store raw 8-byte floats contiguously, run the loops in compiled C or Fortran, use SIMD instructions that multiply 4 to 16 numbers at once, split work across CPU cores, and arrange the computation in cache-sized blocks. That's the BLAS (Basic Linear Algebra Subprograms) library that NumPy calls for @. The cell below compares the two. (It imports NumPy, which takes a few seconds to download the first time in the browser.)

Notice two things. First, the Python time grows roughly 8× when n doubles: that's O(n³) in action. Second, NumPy is typically hundreds of times faster, and the gap grows with n. In the browser (WebAssembly, single thread) the gap is smaller than on your own machine; run it locally with pixi run python timing.py to see the full effect.

pixi add numpy matplotlib     # conda-forge packages; NumPy links an optimized BLAS
pixi run python timing.py

Remember: interpreted loops pay object and dispatch overhead on every number. Use pure Python to understand; use NumPy to compute.

Practice

These exercises are graded in the browser. Use only plain Python (no NumPy), and the math module.

import math

def dot(v, w):
    if len(v) != len(w):
        raise ValueError(f"length mismatch: {len(v)} vs {len(w)}")
    return sum(a * b for a, b in zip(v, w))

def norm(v):
    return math.sqrt(dot(v, v))

import math

def cosine_similarity(v, w):
    nv = math.sqrt(sum(a * a for a in v))
    nw = math.sqrt(sum(b * b for b in w))
    if nv == 0 or nw == 0:
        raise ValueError("cosine similarity is undefined for the zero vector")
    return sum(a * b for a, b in zip(v, w, strict=True)) / (nv * nw)

def most_similar(query, items):
    return max(items, key=lambda name: cosine_similarity(query, items[name]))

def shape(M):
    if not M:
        return (0, 0)
    return (len(M), len(M[0]))

def transpose(M):
    rows, cols = shape(M)
    return [[M[i][j] for i in range(rows)] for j in range(cols)]

def identity(n):
    return [[1 if i == j else 0 for j in range(n)] for i in range(n)]

def matmul(A, B):
    if len(A[0]) != len(B):
        raise ValueError(f"inner dimensions differ: {len(A[0])} vs {len(B)}")
    n, p = len(B), len(B[0])
    return [[sum(A[i][k] * B[k][j] for k in range(n)) for j in range(p)] for i in range(len(A))]

def matpow(A, k):
    result = identity(len(A))
    for _ in range(k):
        result = matmul(result, A)
    return result

import math
from numbers import Real

class Vector:
    __slots__ = ("_c",)

    def __init__(self, *components):
        self._c = tuple(components)

    def __repr__(self):
        return f"Vector({', '.join(repr(c) for c in self._c)})"

    def __len__(self):
        return len(self._c)

    def __iter__(self):
        return iter(self._c)

    def __getitem__(self, i):
        return self._c[i]

    def __eq__(self, other):
        if not isinstance(other, Vector):
            return NotImplemented
        return self._c == other._c

    def __hash__(self):
        return hash(self._c)

    def _check(self, other):
        if len(self) != len(other):
            raise ValueError(f"dimension mismatch: {len(self)} vs {len(other)}")

    def __add__(self, other):
        if not isinstance(other, Vector):
            return NotImplemented
        self._check(other)
        return Vector(*(a + b for a, b in zip(self, other)))

    def __sub__(self, other):
        if not isinstance(other, Vector):
            return NotImplemented
        self._check(other)
        return Vector(*(a - b for a, b in zip(self, other)))

    def __mul__(self, scalar):
        if not isinstance(scalar, Real):
            return NotImplemented
        return Vector(*(scalar * a for a in self))

    __rmul__ = __mul__

    def __matmul__(self, other):
        if not isinstance(other, Vector):
            return NotImplemented
        self._check(other)
        return sum(a * b for a, b in zip(self, other))

    def __abs__(self):
        return math.sqrt(self @ self)

    def __neg__(self):
        return self * -1

Mini-project: a tiny 2-D transformation toolkit

pixi init transform2d && cd transform2d
pixi add python=3.14 matplotlib pytest
pixi run python demo.py          # writes shapes.png
pixi run pytest -q

Quiz

Recap

• A vector is an ordered list of numbers and an arrow. Add and scale component-wise; length is √(Σ vᵢ²).
• The dot product Σ vᵢwᵢ = |v||w|cos θ measures agreement; cosine similarity is its length-free version.
• A matrix has shape (rows, cols). The transpose swaps rows and columns: Aᵀ[j][i] = A[i][j].
• (AB)[i][j] is row i of A dotted with column j of B. Shapes chain as (m×n)(n×p) → (m×p). Multiplication is associative but not commutative, and costs O(n³).
• A 2 × 2 matrix transforms the plane; its columns are where the basis vectors land; products compose transformations right to left.
• Memory is linear: row-major vs column-major order decides which loop order is cache friendly.
• Pure Python pays object and interpreter overhead per number, which is why NumPy comes next.

Going further

• 3Blue1Brown: Essence of Linear Algebra, the best visual introduction to everything in this module
• Python docs: emulating numeric types (__add__, __rmul__, __matmul__)
• PEP 465: why Python added the @ operator
• math.sumprod, math.hypot, and math.dist
• Next: M48 · NumPy Fundamentals