Linear Algebra for Machine Learning from Scratch

01A vector is a list that means something

A language model can't work with the word "cat" directly. It looks the word up in a table and gets back a list of numbers, 768 of them in GPT-2 small, called an embedding. Training nudges those numbers until words used in similar ways end up with similar lists. In math, a list of numbers like this is a vector, and the first thing a model needs is a way to measure how similar two vectors are.

The tool for that is the dot product: multiply the matching entries of two vectors and add up the results. That's all it is.

Draw each vector as an arrow from the origin and the result has a geometric meaning. The dot product is large and positive when the arrows point the same way, zero when they're at right angles, and negative when they point in opposite directions. It also grows with the arrows' lengths. A vector's length, its norm, comes from the same operation: √(a · a), Pythagoras' theorem for any number of entries. For a = (2, 1) that's √5 ≈ 2.24.

In words: multiply matching entries and add them up. The result always equals the length of a, times the length of b, times the cosine of the angle between them.^[1] That second form is why the dot product measures direction.

For the example, ‖a‖‖b‖ = √5 × √10 ≈ 7.07, so cos θ = 5 / 7.07 ≈ 0.71 and the angle is 45°. You can't picture an angle between two 768-number vectors, but the arithmetic is identical, so "angle" still means "how closely do these point the same way".

Often you want direction only, without length. Divide the dot product by both lengths and you get cosine similarity, a score from −1 to 1. It's the standard way to compare embeddings, because a long vector and a short one pointing the same way should count as alike.^[2] The dot product itself shows up everywhere else in a model:

• A neuron computes w · x + bias: a dot product of its weights with its input. Each weight says how much one input matters.
• Attention, the core of a transformer, scores how relevant each earlier word is by taking a dot product between a "query" vector and a "key" vector for every pair of words.
• Vector search ranks documents by the cosine similarity of their embeddings to the question's embedding.

02A matrix turns one list into another

A layer in a neural network takes a vector in and produces a different vector out: in GPT-2 small, one layer turns each 768-number vector into a 3072-number one. The layer's weights are stored as a grid with one row per output and one column per input. That grid is a matrix, and running the layer means multiplying the matrix by the input vector.

Each output number is the dot product of one row with the input. So a layer is a stack of dot products, and each output asks the same question as a neuron: how much does this input line up with my row of weights?

The column reading is the one that makes a matrix easy to picture. The first column, (1, 3), is where the input (1, 0) ends up; the second, (2, 4), is where (0, 1) ends up. Every other input is a mix of those two, so its output is the same mix of the columns.^[1] That's why mathematicians call a matrix a linear map: it moves every point of space in a way that's fully decided by where a few basic directions land.

In words: the output is a blend of the matrix's columns, and the input's entries say how much of each column to use. Reading row by row gives exactly the same numbers.^[2]

The widget below lets you edit a 2×2 matrix and see what it does to every point in the plane at once. Two numbers describe the result. The determinant says how much the matrix scales area; if it's zero, the matrix squashes the whole plane onto a line and the original input can't be recovered. The singular values say how much the matrix stretches along its strongest and weakest directions: a circle of inputs becomes an ellipse, and the ellipse's half-widths are the singular values.

Rotations are a special case worth remembering: they turn vectors without changing any length or any dot product, so a rotation matrix never loses or amplifies information.^[1] Matrices with that property are called orthogonal. Transformers use exactly this property to tell the model where each word sits in the sentence: rotary position embeddings rotate each word's vectors by an angle that depends on its position, which changes their direction but never their size.

03Matrix multiplication is the cost of a model

Models process many inputs at once. Stack 32 input vectors as the rows of one matrix and a layer becomes matrix times matrix. For large language models, these matrix multiplications are where most of the arithmetic goes, so it pays to know their rules and their cost.

The rule: each entry of the result is the dot product of a row of the left matrix with a column of the right one. For that to work, a row of the left matrix must be as long as a column of the right one. The shared size gets summed away, and the outer sizes give the result's shape: (m×k) times (k×n) gives (m×n).^[2]

In words: to fill in row i, column j of the result, walk along row i of A and down column j of B, multiplying pairs and adding them up. Swapping the order, BA instead of AB, gives a different result and often a different shape.

The cost follows directly. The result has m·n entries and each one is a dot product of length k, so the whole thing takes m·n·k multiply-adds. Hardware specs usually count a multiply and an add as two FLOPs, giving 2mnk.

The arithmetic is fixed by the shapes. What the programmer controls is the order the loops visit the numbers in, and that decides how fast the memory can keep up. Most languages and libraries, including C, C++, Rust, NumPy and PyTorch, store a matrix row by row by default.^[3] A CPU fetches memory in chunks (64-byte cache lines on typical x86 chips), so reading the next entry along a row is nearly free, while jumping down a column means a fresh fetch every step.

The textbook loop order (row of A, column of B, then the inner sum) walks B down a column: bad. Swapping the two inner loops walks along rows of both B and the result: good. Same multiplications, same answer, different memory pattern. The memory hierarchy section of the containers page explains why it matters so much; the widget counts it.

Production libraries go much further: they cut the matrices into tiles small enough to stay in cache, reuse each tile many times, and use SIMD or GPU tensor-core instructions that do many multiply-adds at once. Loop order is the first step on that ladder, and on matrices bigger than the cache it alone commonly changes run time several-fold.

04Shapes and the bugs that don't crash

Real model code juggles arrays with more than two axes. A batch of 32 sentences, each 128 tokens long, each token a 768-number vector, is one array of shape (32, 128, 768). Machine learning code calls an array like that a tensor, and its shape is the first thing to check when the numbers look wrong.

A very common operation is adding one vector to every row: a layer adds the same bias vector of 768 numbers to every token's output. Writing a loop for that is tedious, so array libraries let you write outputs + bias with shapes (32, 768) and (768,) and stretch the smaller array to fit. That stretching is broadcasting, and NumPy, PyTorch and JAX all use the same rule:^[4]

1. Line the two shapes up from the right.
2. Compare sizes pair by pair. They must be equal, or one of them must be 1. A missing axis counts as 1.
3. Any size-1 axis is stretched to match the other. Anything else is an error.

The dangerous cases are the ones that don't raise an error. Say you have 4 examples with 4 features each, shape (4, 4), and you want to scale each example by its own weight, a vector of shape (4,). Broadcasting lines the vector up with the features axis, so it silently scales each column instead of each row. Or add a (32, 1) column to a (1, 32) row and you get a 32×32 grid nobody asked for. The fix is cheap: reshape with explicit axes, such as (4, 1) for per-row weights, and assert shapes at function boundaries while developing.

Weight files have a shape trap too. PyTorch's linear layer stores its weight matrix as (outputs, inputs) and computes y = xWᵀ + b, the transpose of the Ax convention used above.^[5] Other libraries store the other orientation. Load one into the other without transposing and the model runs, produces numbers, and is wrong.

05What a matrix stretches, and how to shrink it

Two practical questions lead to the same idea. Why do signals in a deep network sometimes blow up or fade to nothing as they pass through layer after layer? And how can you fine-tune a model with billions of weights by training only a tiny fraction of them? Both answers come from asking which directions a matrix stretches, and by how much.

Start with a square matrix. Most input directions get both stretched and turned. A few special directions only get stretched: the output points the same way as the input. Those directions are the matrix's eigenvectors, and the stretch factor for each is its eigenvalue.^[1]

Now multiply some other vector by A again and again. Every pass triples the part of it that lies along (1, 1) and leaves the part along (1, −1) alone, so after a few passes the vector points almost exactly along (1, 1). That repeated multiplication is called power iteration, and it's a small model of a deep network. A signal passing through many similar layers gets pulled toward the direction those layers stretch most, and grows or shrinks by roughly that stretch factor at every layer. Stretch factors above 1 make signals explode; below 1 they fade. The training page shows how careful initialization keeps them near 1.

Eigenvectors only exist for square matrices, and most weight matrices aren't square. The tool that works for every matrix is the singular value decomposition (SVD). It says that any matrix does three simple things in a row: turn the input, stretch it along the axes by a set of factors called the singular values, then turn the result.^[1]

The useful reading for machine learning is the second form in the equation below. It splits the matrix into a sum of simple layers, each one a single column times a single row, weighted by its singular value and sorted biggest first. The big singular values carry most of the matrix; the small ones carry fine detail.

In words: the matrix equals turn, stretch, turn. Equivalently, it equals a stack of single-column-times-single-row pieces, each scaled by how much it matters. The number of pieces with a non-zero weight is the matrix's rank, the number of independent directions it really uses.

Keep only the k biggest pieces and you get a rank-k matrix that is the closest possible rank-k approximation of the original. That's the Eckart–Young theorem.^[6] Storing it takes k(m + n + 1) numbers instead of m × n. For a 4096×4096 matrix at k = 16, that's about 131 thousand numbers instead of 16.8 million.

This is the idea behind LoRA, the standard way to fine-tune a large model on modest hardware. Instead of updating a full 4096×4096 weight matrix, LoRA freezes it and trains a correction made of two thin matrices, 4096×16 and 16×4096. Their product has the full shape but only rank 16, and the two factors hold under 1% of the numbers.^[7] The fine-tuning page measures what that small rank costs. The same decomposition also gives PCA, the classic technique for finding the few directions along which a dataset varies most.^[1]

06A checked matrix kernel

These complete programs build a small matrix type stored row by row, multiply in the row-friendly loop order, and check three claims from this page: the fast order gives the same answer as the textbook order, a rotation keeps a vector's length, and a column times a row gives a matrix where every row is a multiple of the same vector (rank one). Both listings use the same names so you can read them side by side.

#include <cassert>
#include <cmath>
#include <cstddef>
#include <vector>
struct Matrix {
    std::size_t rows, columns;
    std::vector<double> data; // Row-major: element (row, column) sits at row*columns + column.
    Matrix(std::size_t rowCount, std::size_t columnCount)
        : rows(rowCount), columns(columnCount), data(rowCount * columnCount, 0.0) {}
    double& at(std::size_t row, std::size_t column) { return data[row * columns + column]; }
    double at(std::size_t row, std::size_t column) const { return data[row * columns + column]; }
};
// Row-friendly order: the inner loop walks along a row of right and a row of product,
// so consecutive reads and writes are adjacent in memory.
Matrix multiply(const Matrix& left, const Matrix& right) {
    assert(left.columns == right.rows); // Inner sizes must agree.
    Matrix product(left.rows, right.columns);
    for (std::size_t row = 0; row < left.rows; ++row)
        for (std::size_t inner = 0; inner < left.columns; ++inner) {
            const double leftValue = left.at(row, inner); // Constant across the inner loop.
            for (std::size_t column = 0; column < right.columns; ++column)
                product.at(row, column) += leftValue * right.at(inner, column);
        }
    return product;
}
Matrix multiplyTextbook(const Matrix& left, const Matrix& right) {
    Matrix product(left.rows, right.columns);
    for (std::size_t row = 0; row < left.rows; ++row)
        for (std::size_t column = 0; column < right.columns; ++column)
            for (std::size_t inner = 0; inner < left.columns; ++inner)
                product.at(row, column) += left.at(row, inner) * right.at(inner, column);
    return product;
}
double norm(const Matrix& vector) {
    double sum = 0;
    for (double value : vector.data) sum += value * value;
    return std::sqrt(sum);
}
int main() {
    Matrix left(3, 4), right(4, 5);
    for (std::size_t index = 0; index < left.data.size(); ++index) left.data[index] = double(index % 7) - 3;
    for (std::size_t index = 0; index < right.data.size(); ++index) right.data[index] = double(index % 5) * 0.5;
    const Matrix fast = multiply(left, right), slow = multiplyTextbook(left, right);
    for (std::size_t index = 0; index < fast.data.size(); ++index)
        assert(std::abs(fast.data[index] - slow.data[index]) < 1e-12); // Same arithmetic, different order.
    Matrix rotation(2, 2), point(2, 1);
    const double angle = 0.6;
    rotation.at(0, 0) = std::cos(angle); rotation.at(0, 1) = -std::sin(angle);
    rotation.at(1, 0) = std::sin(angle); rotation.at(1, 1) = std::cos(angle);
    point.at(0, 0) = 3; point.at(1, 0) = -4;
    assert(std::abs(norm(multiply(rotation, point)) - 5.0) < 1e-12); // Orthogonal maps preserve length.
    Matrix columnU(3, 1), rowV(1, 4);
    columnU.data = {1, 2, -1}; rowV.data = {2, 0, 1, 3};
    const Matrix outer = multiply(columnU, rowV); // Rank one: every row is a multiple of rowV.
    for (std::size_t row = 0; row < 3; ++row)
        for (std::size_t column = 0; column < 4; ++column)
            assert(outer.at(row, column) == columnU.at(row, 0) * rowV.at(0, column));
}

struct Matrix { rows: usize, columns: usize, data: Vec<f64> } // Row-major storage.
impl Matrix {
    fn new(rows: usize, columns: usize) -> Matrix { Matrix { rows, columns, data: vec![0.0; rows * columns] } }
    fn at(&self, row: usize, column: usize) -> f64 { self.data[row * self.columns + column] }
    fn set(&mut self, row: usize, column: usize, value: f64) { self.data[row * self.columns + column] = value; }
}
// Row-friendly order: the inner loop walks along a row of right and a row of product,
// so consecutive reads and writes are adjacent in memory.
fn multiply(left: &Matrix, right: &Matrix) -> Matrix {
    assert!(left.columns == right.rows); // Inner sizes must agree.
    let mut product = Matrix::new(left.rows, right.columns);
    for row in 0..left.rows {
        for inner in 0..left.columns {
            let left_value = left.at(row, inner); // Constant across the inner loop.
            for column in 0..right.columns {
                product.data[row * product.columns + column] += left_value * right.at(inner, column);
            }
        }
    }
    product
}
fn multiply_textbook(left: &Matrix, right: &Matrix) -> Matrix {
    let mut product = Matrix::new(left.rows, right.columns);
    for row in 0..left.rows { for column in 0..right.columns { for inner in 0..left.columns {
        product.data[row * product.columns + column] += left.at(row, inner) * right.at(inner, column);
    }}}
    product
}
fn norm(vector: &Matrix) -> f64 { vector.data.iter().map(|value| value * value).sum::<f64>().sqrt() }
fn main() {
    let mut left = Matrix::new(3, 4);
    let mut right = Matrix::new(4, 5);
    for (index, value) in left.data.iter_mut().enumerate() { *value = (index % 7) as f64 - 3.0; }
    for (index, value) in right.data.iter_mut().enumerate() { *value = (index % 5) as f64 * 0.5; }
    let fast = multiply(&left, &right);
    let slow = multiply_textbook(&left, &right);
    for (fast_value, slow_value) in fast.data.iter().zip(&slow.data) {
        assert!((fast_value - slow_value).abs() < 1e-12); // Same arithmetic, different order.
    }
    let angle: f64 = 0.6;
    let mut rotation = Matrix::new(2, 2);
    rotation.set(0, 0, angle.cos()); rotation.set(0, 1, -angle.sin());
    rotation.set(1, 0, angle.sin()); rotation.set(1, 1, angle.cos());
    let mut point = Matrix::new(2, 1);
    point.set(0, 0, 3.0); point.set(1, 0, -4.0);
    assert!((norm(&multiply(&rotation, &point)) - 5.0).abs() < 1e-12); // Orthogonal maps preserve length.
    let column_u = Matrix { rows: 3, columns: 1, data: vec![1.0, 2.0, -1.0] };
    let row_v = Matrix { rows: 1, columns: 4, data: vec![2.0, 0.0, 1.0, 3.0] };
    let outer = multiply(&column_u, &row_v); // Rank one: every row is a multiple of row_v.
    for row in 0..3 { for column in 0..4 {
        assert!(outer.at(row, column) == column_u.at(row, 0) * row_v.at(0, column));
    }}
}

07Where the shapes lie

Cosine similarity breaks near zero. It divides by both lengths, so it's undefined for an all-zero vector, and for a nearly-zero one a tiny change can swing the angle anywhere. Check that embedding lengths stay well away from zero before normalizing them, and keep the raw dot product when length carries meaning.

Floating-point dot products aren't exact. Adding up a long vector piles up rounding error, and the error depends on the order of the additions, so two correct matrix kernels can disagree in the last few digits. The floating point page explains why. In practice, compare results with a tolerance scaled to the size of the numbers, as the programs above do.

Determinants and eigenvalues need square matrices. For a rectangular weight matrix, the singular values tell you how much a layer can stretch its input; the largest one is called the spectral norm.^[1] Weight matrices with large spectral norms are one common contributor to activations that explode in deep stacks, which the training page returns to.

08What's next

A model's last layer produces a vector of scores, one per possible next word. Probability and Information Theory explains how those scores become probabilities, how to grade a prediction against what actually happened, and why that grade, the cross-entropy loss, is what nearly every model in this series is trained to minimize.

09Sources

The cited texts define the objects used here. Widget readouts report values computed in the page from the displayed matrices.

1. Marc Peter Deisenroth, A. Aldo Faisal, Cheng Soon Ong, 2020. Mathematics for Machine Learning, Cambridge University Press. Norms (Definition 3.1), the Cauchy–Schwarz inequality and angles (Sections 3.3–3.4), orthogonal matrices (Definition 3.8), eigenvalues (Definition 4.6), the SVD theorem (4.22), spectral norm (Theorem 4.24), the Eckart–Young theorem (4.25), and PCA as truncated SVD (Section 10.4).
2. Aston Zhang, Zachary C. Lipton, Mu Li, Alexander J. Smola, 2023. Dive into Deep Learning, Section 2.3: Linear Algebra, Cambridge University Press. Dot products as weighted sums and cosines, matrix–vector and matrix–matrix products as dot products.
3. NumPy developers. The N-dimensional array: memory layout. Row-major (C order) storage and strides.
4. NumPy developers. Broadcasting. The trailing-axis alignment rule, size-1 stretching, and the broadcastable examples.
5. PyTorch contributors. torch.nn.Linear. The y = xAᵀ + b convention and the (out_features, in_features) weight shape.
6. Carl Eckart, Gale Young, 1936. The approximation of one matrix by another of lower rank, Psychometrika 1(3). The optimality of truncated SVD as a least-squares (Frobenius-norm) approximation; source 1 states the spectral-norm form.
7. Edward J. Hu et al., 2021. LoRA: Low-Rank Adaptation of Large Language Models; reference implementation at github.com/microsoft/LoRA. Frozen weights with learned rank-decomposition updates scaled by α/r.