Matrix Factorization

The Story That Explains Everything

That compression — from one enormous sparse table into two small dense ones — is the entire idea. The rest of this tutorial unpacks how it works, why it's mathematically sound, and how to implement it from scratch and with libraries.

The Rating Matrix R — Our Starting Point

Before factorizing anything, you need to understand what R looks like. It is a matrix of shape m × n where m = number of users and n = number of items. Each observed entry R[u, i] is the rating user u gave item i. Most entries are missing — we use a ? to denote unobserved cells.

Only 14 of 30 possible cells are filled — that's a 53% sparsity rate. In real systems this is 99.99%+. Our job: predict every ? accurately.

Latent Factors — The Hidden Language of Taste

A latent factor is a hidden dimension that explains observed behaviour without being directly measured. We choose how many latent factors to use (call it k). Typical values: 10, 20, 50, 100 for large systems.

The Core Formula — R ≈ PQᵀ

The fundamental idea: decompose the rating matrix R (shape m×n) into two smaller matrices P (shape m×k) and Q (shape n×k), such that their product approximates R.

Low-Rank Approximation — The Mathematical Backbone

The phrase "low-rank approximation" is the formal mathematical description of what we are doing. It comes from linear algebra — let's build intuition carefully.

SVD — Singular Value Decomposition

The theoretically optimal low-rank approximation for a complete matrix comes from Singular Value Decomposition (SVD). SVD decomposes any matrix R into three matrices:

Hidden Embeddings — The Secret Geometry of Taste

Learning P and Q — Gradient Descent on the Loss

We initialise P and Q with small random values, then optimise them to minimise prediction error on the observed ratings only. The loss function is:

Matrix Factorization from Scratch in Python

The cleanest way to build intuition is to implement the SGD update loop yourself. The entire algorithm fits in about 30 lines of NumPy.

import numpy as np

# ── Toy rating matrix (0 = not rated) ────────────────────────────────────
R = np.array([
    [5, 4, 0, 0, 0, 5],
    [0, 0, 4, 5, 4, 0],
    [4, 5, 0, 0, 0, 0],
    [0, 0, 3, 0, 5, 0],
    [2, 0, 0, 3, 0, 1],
], dtype=float)

m, n = R.shape   # 5 users, 6 movies
k     = 2         # latent factors
eta   = 0.01      # learning rate
lam   = 0.01      # regularisation

# ── Initialise P and Q randomly ────────────────────────────────────────
np.random.seed(42)
P = np.random.normal(0, 0.1, (m, k))
Q = np.random.normal(0, 0.1, (n, k))

# ── Index of observed ratings ──────────────────────────────────────────
observed = [(u, i) for u in range(m) for i in range(n) if R[u, i] > 0]

# ── SGD training loop ──────────────────────────────────────────────────
for epoch in range(500):
    for u, i in observed:
        pred  = P[u] @ Q[i]                   # dot product
        err   = R[u, i] - pred                 # residual error
        P[u] += eta * (err * Q[i] - lam * P[u])
        Q[i] += eta * (err * P[u] - lam * Q[i])

    if epoch % 100 == 0:
        total_loss = sum((R[u,i] - P[u] @ Q[i])**2 for u,i in observed)
        print(f"Epoch {epoch:4d} | Loss: {total_loss:.4f}")

# ── Reconstruct full matrix ────────────────────────────────────────────
R_hat = P @ Q.T
print("\nPredicted Rating Matrix R̂:")
print(np.round(R_hat, 2))

Truncated SVD in Practice

When your matrix is reasonably complete (e.g. after mean-imputing missing values for a quick baseline), scipy.sparse.linalg.svds gives you the rank-k truncated SVD directly without computing all singular values.

import numpy as np
from scipy.sparse.linalg import svds

# ── Impute missing values with row means (simple baseline) ─────────────
R_full = R.copy()
row_means = np.true_divide(R_full.sum(1), (R_full != 0).sum(1))
for u in range(m):
    R_full[u, R_full[u] == 0] = row_means[u]

# ── Truncated SVD with k=2 factors ─────────────────────────────────────
k = 2
U, sigma, Vt = svds(R_full, k=k)

# sigma is a 1D array; make it a diagonal matrix
Sigma = np.diag(sigma)

# ── Reconstruct ─────────────────────────────────────────────────────────
R_svd = U @ Sigma @ Vt
print("SVD Low-rank Reconstruction:")
print(np.round(R_svd, 2))

# ── Singular values show variance captured ──────────────────────────────
print(f"\nSingular values: {np.round(sigma, 3)}")
total_var = np.sum(sigma**2)
print(f"Variance explained by k={k}: {100*total_var/np.sum(np.linalg.svd(R_full,compute_uv=False)**2):.1f}%")

Production Implementation — Surprise Library

For real datasets, the Surprise library provides optimised, cross-validated matrix factorization via its SVD class (which is actually the SGD-optimised MF, not classical SVD — a common naming quirk).

from surprise import SVD, Dataset, Reader
from surprise.model_selection import cross_validate, GridSearchCV
import pandas as pd

# ── Load MovieLens 100k ───────────────────────────────────────────────
data = Dataset.load_builtin('ml-100k')

# ── Baseline SVD — default hyperparams ───────────────────────────────
algo = SVD(n_factors=50, n_epochs=20, lr_all=0.005, reg_all=0.02, random_state=42)
results = cross_validate(algo, data, measures=['RMSE', 'MAE'], cv=5, verbose=True)

print(f"\nMean RMSE: {results['test_rmse'].mean():.4f}")
print(f"Mean MAE:  {results['test_mae'].mean():.4f}")

# ── Grid search over k and regularisation ────────────────────────────
param_grid = {
    'n_factors': [20, 50, 100],
    'n_epochs':  [20, 40],
    'lr_all':    [0.003, 0.005],
    'reg_all':   [0.02, 0.1]
}

gs = GridSearchCV(SVD, param_grid, measures=['rmse'], cv=3)
gs.fit(data)
print(f"\nBest RMSE:   {gs.best_score['rmse']:.4f}")
print(f"Best params: {gs.best_params['rmse']}")

# ── Predict a specific user-item rating ──────────────────────────────
trainset = data.build_full_trainset()
algo.fit(trainset)
pred = algo.predict(uid='196', iid='302', r_ui=4, verbose=True)