What it actually does
A ratings matrix has millions of rows and columns and almost nothing in it — most users have touched almost nothing. The bet is that behaviour is driven by a small number of underlying factors, so the matrix is approximately low rank: describe each user by a short vector of tastes, each item by a short vector of attributes, and predict an interaction as the dot product between them. Nothing says what the factors mean, and they are not interpretable; what matters is that fitting them on the entries you observed lets you fill in the ones you did not. The crucial detail is that the sum runs only over observed entries — treating the blanks as zeros would be a completely different and much worse model.
- The interaction matrix is approximately low rank — a small number of latent factors explains most behaviour, which is the entire premise
- Entries are missing in a way the model can ignore, which is false in detail: users choose what to rate, so missingness is informative and the standard formulation assumes it away
- Observed interactions reflect preference, which explicit ratings support and implicit signals such as clicks support much more weakly
- Users and items have enough observations for their factors to be estimated; below a handful the regularizer is doing the predicting
- Preferences are static over the training window, since the model has no notion of time unless one is added
Loss function
Explain each observed interaction as a baseline plus an interaction between what the user likes and what the item is, and keep the factor vectors small. Two things deserve emphasis. The restriction to observed entries is the entire modelling decision: an SVD would require the missing values to be filled in, and filling them with zeros asserts that unrated means disliked, which is false. And the bias terms are not a technicality — a generous user rating a popular film is largely explained before any factor is consulted, and models that skip biases spend their factor capacity relearning them.
- the set of OBSERVED entries — the sum runs over these only, which is the difference between this and an SVD
- latent factor vectors for user u and item i, typically 10 to 200 dimensions
- global mean and per-user, per-item biases — these explain more of the variance than the factors do
- regularization, and not optional: a user with three ratings can otherwise fit them exactly
How it is minimized
Alternating least squares, or stochastic gradient descent over observed entries
- the factor vectors of the items user u actually rated, stacked — only those rows participate
- ridge term added to the diagonal, which is also what keeps the tiny k-by-k system invertible for a user with few ratings
- alternating: solve every user with items fixed, then every item with users fixed, and repeat
- the observed ratings of user u, with the biases already subtracted
The objective is bilinear: not convex in P and Q jointly, but exactly a ridge regression in either one with the other held fixed. So alternating minimization applies, and each half-step has a closed form — a k-by-k solve per user, then per item — with k typically under two hundred, which makes those solves trivially cheap. This is the same alternating-minimization structure as k-means and EM, and it inherits the same property: monotone descent, no guarantee of a global optimum, results that depend on initialization. The reason to prefer ALS over SGD is parallelism, and it is decisive: within a half-step every user is independent of every other, so the work distributes across cores or machines with no coordination. SGD converges in fewer passes on a single machine and does not partition that way.
k (rank)
Number of latent factors, and the capacity dial. Diminishing returns arrive quickly and overfitting arrives with them (10 to 200; the Netflix-era consensus was that gains flatten well before 200)
lambda
Regularization. Not optional — a user with three observed entries will otherwise reproduce them exactly and generalize to nothing (0.01 to 1.0, cross-validated; often scaled by the number of ratings per user)
iterations
ALS sweeps. Convergence is fast, usually 10 to 20, because each half-step is an exact solve rather than a gradient step
alpha (implicit feedback)
Confidence scaling on observed interactions when there are no explicit ratings. Turns a click into a weighted positive rather than a rating
biases
Whether to fit per-user and per-item offsets. Almost always yes; they explain more variance than the factors and cost k=0 parameters each
ALS descends monotonically and converges quickly — ten to twenty sweeps is typical, because each half-step solves its subproblem exactly rather than stepping toward it. What it converges to depends on the initialization, since the joint objective is non-convex. The failure modes are about data rather than optimization. Cold start is structural: a user or item with no observed entries has no equation to solve and no factor vector, and no amount of training fixes it — content features or a popularity fallback are the only answers. Popularity bias is self-reinforcing: popular items appear in more equations, get better-estimated factors, get recommended more, and accumulate more interactions. And the factors are not identified — any invertible transformation of P with its inverse applied to Q gives the same predictions — so interpreting individual factors is unfounded no matter how suggestive they look.
ALS: O(|O|·k^2 + (m + n)·k^3) per sweep, where |O| is the number of observed entries — linear in the data and cubic only in the small rank. For implicit feedback the naive form is O(m·n·k^2) because every unobserved pair is a weak negative, and the Gramian trick reduces it back to O(|O|·k^2 + (m + n)·k^3), which is what makes implicit ALS feasible at all. Memory is O((m + n)·k) for the factors, which is small — a million users at k=64 is a few hundred megabytes. Inference is a single dot product, or a nearest-neighbour search over item factors when the task is top-N retrieval.
How, where, and why it is used
Time-Series Forecasting
AdaptedMatrix completion over a series-by-time matrix. With many related series observed irregularly — sensors dropping out, stores reporting late — the matrix of series against periods is low rank because the series share drivers, and factorizing it fills the gaps. Temporal regularized variants add a penalty tying adjacent time factors together, which is what turns completion into something that can extrapolate one step rather than only interpolate.
- Imputing missing observations across a panel of related series before a forecaster consumes them
- Traffic and sensor networks with systematic dropout, where the spatial correlation is exactly the low-rank structure
- Extracting shared latent factors from a wide panel, as a variant of the factor-model approach PCA takes
It handles the situation PCA cannot: a panel with holes in it. PCA needs a complete matrix, and imputing before decomposing biases the result toward whatever was imputed; this fits only what was observed. The limitation is that plain matrix factorization has no notion of time — the columns could be shuffled and the fit would be identical — so it interpolates well and extrapolates not at all unless a temporal regularizer is added, at which point it is a different and more complex model.
- Centre by series and by period before factorizing, so the biases carry the level and the factors carry the structure
- Add a temporal smoothness penalty on the time factors if forecasting rather than imputing is the goal
- Hold out observed entries at random for validation, and separately hold out a time block, since the two measure different things
- Keep the rank low — a wide panel with few observations per series overfits quickly
For imputation, error on held-out observed entries. For forecasting, a rolling-origin backtest of whatever consumes the completed matrix. Reconstruction error on the training entries says nothing about either.
- Column order is irrelevant to the model, so nothing stops it from "predicting" a period using information from later ones unless the split is temporal
- Series with almost no observations get factors driven entirely by the regularizer, which is a shrunk guess rather than a prediction
- Treating an imputed matrix as observed data downstream, so the uncertainty of the completion disappears from every subsequent estimate
Anomaly Detection
Not applicableThe model reconstructs the entries it was fitted on and has no notion of normality for a whole observation; a large residual on a single cell usually means that user or item had too few observations, not that anything anomalous occurred.
Optimization
AdaptedThe objective is bilinear — non-convex jointly, exactly a ridge regression in each factor alone — which makes it the cleanest example of why alternating minimization is reached for. Each half-step has a closed form, so there is no step size and no line search, and the descent is monotone by construction. It is the same pattern as EM and as k-means, with the advantage that here the subproblem is a textbook regularized least squares rather than something bespoke.
- The reference example of a bilinear objective and why alternating least squares beats gradient descent on it
- Low-rank matrix completion as a relaxation of a rank-constrained problem, which is where nuclear-norm methods come from
- Recommendation as a constrained allocation problem, where the predicted preference surface feeds a slate or budget optimizer
The interesting content is the contrast with gradient descent, and it comes out in favour of exploiting structure: because each subproblem is solved exactly rather than stepped toward, ALS converges in tens of sweeps where SGD needs many passes, and it parallelizes because the subproblems within a half-step are independent. The counterweight is the honest one — monotone descent on a non-convex surface says nothing about the global optimum, and different initializations reach genuinely different factorizations with similar loss.
- Add lambda to the diagonal of every small solve; it is simultaneously the regularizer and what keeps a sparse user’s system invertible
- Scale the regularizer by the number of observations per user or item, so a user with three ratings is shrunk harder than one with three hundred
Verify the loss decreases monotonically across sweeps — an increase is a bug, almost always a bias term updated inconsistently between the two half-steps. Compare across random initializations to see how non-convex the surface is on this data.
- Reading monotone descent as convergence to the global optimum
- Inverting the k-by-k matrix explicitly instead of solving with a Cholesky factorization, which is both slower and worse conditioned
- Regularizing biases at the same strength as factors, which shrinks the baseline the model most depends on
Recommendation & Ranking
Primary fitThe canonical application and the one it was developed for. Factorize the user-item interaction matrix into latent factors, then score a candidate by the dot product of the two vectors plus the biases. For top-N recommendation the item factors form a vector index and retrieval becomes an approximate nearest-neighbour search over them, which is what makes serving a catalogue of millions tractable.
- The core of collaborative filtering systems since the Netflix Prize, and still a production baseline
- Implicit-feedback recommendation from clicks, plays, and purchases, using the confidence-weighted variant
- Candidate generation ahead of a learned ranker, where item factors are indexed and queried by user vector
- Learning item embeddings for downstream use, where the factorization is a means to the vectors rather than to the predictions
It captures the fact that similar users like similar things without anyone specifying what similar means, and it does so from the interaction matrix alone. Two properties keep it deployed: the model is small — factors, not data — and scoring is a dot product, so it fits inside a request path where a neighbourhood method would not. Its structural limits are equally clear and are what two-tower models exist to fix: no way to use content or context, and nothing to say about a user or item with no history.
- Always fit biases; user generosity and item popularity explain more than the factors do and cost two parameters each
- For implicit feedback use confidence weighting rather than treating unobserved as negative — a click is evidence of interest, and a non-click is mostly evidence of not having looked
- Scale the regularizer by observation count per user and item, so sparse rows are shrunk harder
- Plan the cold-start fallback explicitly: a content model or popularity, because the factorization simply has no vector to offer
Recall@k and NDCG on a temporal split, not RMSE. The Netflix Prize optimized RMSE and the field has since largely concluded that was the wrong target — a model can improve rating prediction and rank worse. Measure catalogue coverage alongside, or a popularity baseline in disguise will score well.
- Optimizing rating prediction when the product is a ranked list, which are different objectives that disagree
- Popularity bias compounding through the feedback loop, so recommendations narrow over time without any metric moving
- Random rather than temporal splits, which let the model see a user’s future interactions and inflate every offline number
- Cold start treated as an edge case when it is most of the catalogue on any real platform
Natural Language
AdaptedWord embeddings are matrix factorization, whether or not they are presented that way. Build a word-context co-occurrence matrix, reweight it — pointwise mutual information is the usual choice — and factorize it, and the resulting word vectors carry the same distributional semantics that a neural embedding model learns. Skip-gram with negative sampling was shown to be implicitly factorizing a shifted PMI matrix, which makes the connection exact rather than analogical.
- Explicit count-based embeddings, of which GloVe is a weighted least-squares factorization of log co-occurrence counts
- Topic-adjacent decomposition of term-document matrices, where non-negative variants give parts-based components
- Small-corpus embedding where a neural model has too little data to train and a factorization does not
It is worth knowing because it demystifies the neural version: the semantics come from the co-occurrence statistics and the low-rank assumption, not from the network. On a small corpus the explicit factorization is often the better engineering choice, since it is deterministic, has no training loop, and needs no hyperparameter search. What it cannot do is context — one vector per word type, so every sense of a polysemous word is averaged into a single point, which is exactly the limitation contextual embeddings removed.
- Reweight raw counts with PMI or a log transform before factorizing; raw co-occurrence counts are dominated by frequent words
- Apply a context window and subsample frequent tokens, which is what the neural variants do implicitly
- Truncate rare words rather than factorizing a matrix mostly composed of noise rows
Downstream task performance, or word similarity and analogy benchmarks with the caveat that those correlate loosely with real use. Reconstruction error of the co-occurrence matrix is not the objective anyone cares about.
- Factorizing raw counts, where frequency swamps every semantic signal
- One vector per word type, so polysemy is silently averaged away
- Treating dimensions as interpretable when the factorization is identified only up to an invertible transformation
What it costs to run
Linear in the number of observed entries and cubic only in the small rank, which makes it one of the better-scaling models here. ALS parallelizes within each half-step with no coordination, so a billion interactions is a cluster job rather than an impossible one.
A dot product per candidate — microseconds — and for top-N the item factors go into a vector index so retrieval is an approximate nearest-neighbour query rather than a scan. The model is factors only, so a million users at rank 64 is a few hundred megabytes and the item side is far smaller.
Typically nightly or weekly in batch. The awkward part is between refits: a new user accumulates interactions the model cannot use, so production systems usually fold new users in by a cheap approximate solve against fixed item factors rather than waiting for the next full fit.
- Track catalogue coverage of the recommendations, since popularity bias narrows the served set gradually and no accuracy metric registers it
- Monitor the share of requests falling back to the cold-start path — a rise means the refit cadence is no longer keeping up with catalogue turnover
- Watch the distribution of factor norms; a user or item whose norm is dominated by the regularizer has effectively no personalization
- Compare held-out ranking metrics on a temporal split at each refit, never a random one, which flatters every model here
- Factors are identified only up to an invertible transformation, so they are not comparable between refits and nothing downstream may cache or interpret an individual dimension
- Cold start is structural rather than an edge case: a new user or item has no vector at all, and the fallback path is part of the system design, not an afterthought
- Treating unobserved entries as zero rather than as unknown is the most common and most damaging implementation error, and it produces a model that trains fine and recommends badly
- Item factors and the serving index must be rebuilt together; a refreshed factorization against a stale index scores users against vectors that no longer exist
- Biases must be applied consistently at scoring time — omitting them shifts every prediction by the global mean and quietly inverts the ranking near the top
Pros and cons, in context
Learns latent structure from interactions alone, with no feature engineering
Discovers that two items are similar without anyone describing either of them, which is what made collaborative filtering work at all. Worth nothing when there is no interaction history, which is exactly the cold-start case.
The model is small and scoring is a dot product
Factors rather than data, so serving a catalogue of millions is a vector index and a dot product. This cost profile is why it survived long after more accurate models existed.
ALS parallelizes with no coordination
Within a half-step every user is independent, so training distributes across cores or machines cleanly. SGD converges in fewer passes on one machine and does not partition this way — the choice between them is about deployment, not accuracy.
Fits only the observed entries
The distinction from an SVD, and the reason it works on a 99% empty matrix. Imputing the blanks first — which a plain SVD requires — biases the whole decomposition toward whatever was imputed.
Cold start is structural and unfixable within the model
No history means no vector, and no amount of data about other users helps. Every production system needs a separate content or popularity path, which is a design constraint rather than a tuning issue.
Cannot use content, context, or side information
The matrix is all it sees, so time of day, device, item text and user demographics are unavailable. This is precisely the gap two-tower and feature-based models fill, at the cost of a much larger model.
Factors are not interpretable and not identified
Any invertible transformation of one side with its inverse on the other gives identical predictions, so a factor has no meaning independent of its partner. Interpretations of individual dimensions get published anyway and are not defensible.
Popularity bias compounds through the feedback loop
Popular items appear in more equations, get better factors, get recommended more, and gather more interactions. Nothing in the objective resists this, and offline metrics computed on logged data will not reveal it.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Matrix factorization by SGD - the objective, transcribed.
For every OBSERVED entry, compute the error, then nudge both factor vectors
and both biases along their gradients. The loop runs over observations, not
over the matrix - a blank is never visited, which is the entire difference
between this and an SVD.
"""
import random
def fit(observations, n_users, n_items, k=10, lr=0.01, reg=0.05, epochs=20, seed=0):
"""observations is a list of (user, item, rating) triples."""
rng = random.Random(seed)
# Small random init: zeros would make every gradient zero, since the
# gradient with respect to p_u is proportional to q_i and vice versa.
P = [[rng.gauss(0.0, 0.1) for _ in range(k)] for _ in range(n_users)]
Q = [[rng.gauss(0.0, 0.1) for _ in range(k)] for _ in range(n_items)]
user_bias = [0.0] * n_users
item_bias = [0.0] * n_items
mu = sum(rating for _, _, rating in observations) / len(observations)
for _ in range(epochs):
rng.shuffle(observations)
for user, item, rating in observations:
# prediction = mu + b_u + b_i + p_u . q_i
dot = 0.0
for f in range(k):
dot += P[user][f] * Q[item][f]
prediction = mu + user_bias[user] + item_bias[item] + dot
error = rating - prediction
# Biases first: they explain more of the variance than the
# factors do, and they are two parameters rather than k.
user_bias[user] += lr * (error - reg * user_bias[user])
item_bias[item] += lr * (error - reg * item_bias[item])
# Both factor vectors move together, each toward the other
# scaled by the error - which is why they must be updated from
# saved copies rather than in sequence.
for f in range(k):
p_uf = P[user][f]
q_if = Q[item][f]
P[user][f] += lr * (error * q_if - reg * p_uf)
Q[item][f] += lr * (error * p_uf - reg * q_if)
return P, Q, user_bias, item_bias, mu
def predict(P, Q, user_bias, item_bias, mu, user, item):
dot = 0.0
for f in range(len(P[user])):
dot += P[user][f] * Q[item][f]
return mu + user_bias[user] + item_bias[item] + dotO(|O| * k) per epoch in interpreter loops, and inherently sequential — two observations sharing a user or item race on the same vector.
How this language actually wants it written.
"""Matrix factorization by ALS - typed, CSR-indexed, closed-form half-steps."""
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
@dataclass(frozen=True)
class SparseMatrix:
"""CSR: for row r, the entries live at indices[indptr[r]:indptr[r + 1]].
Both orientations are needed - one half-step iterates users, the other
iterates items - so the transpose is built once rather than searched for
repeatedly.
"""
indptr: NDArray
indices: NDArray
values: Vector
shape: tuple[int, int]
def row(self, index: int) -> tuple[NDArray, Vector]:
start, stop = self.indptr[index], self.indptr[index + 1]
return self.indices[start:stop], self.values[start:stop]
@dataclass(frozen=True)
class FactorModel:
user_factors: Matrix # (m, k)
item_factors: Matrix # (n, k)
user_bias: Vector
item_bias: Vector
global_mean: float
def predict(self, user: int, item: int) -> float:
return float(
self.global_mean
+ self.user_bias[user]
+ self.item_bias[item]
+ self.user_factors[user] @ self.item_factors[item]
)
def top_n(self, user: int, n: int) -> NDArray:
"""Scores every item for one user - one matvec, not a loop."""
scores = self.item_factors @ self.user_factors[user] + self.item_bias
return np.argpartition(scores, -n)[-n:]
def _solve_side(
matrix: SparseMatrix,
fixed: Matrix,
bias_self: Vector,
bias_other: Vector,
global_mean: float,
regularization: float,
) -> Matrix:
"""One ALS half-step: solve every row independently.
Each row is a ridge regression of its observed residuals on the fixed
side's factor vectors - a k-by-k system, so the solve is trivial. The
rows are independent of one another, which is the property that makes ALS
parallelizable and SGD not.
"""
rows, k = matrix.shape[0], fixed.shape[1]
updated = np.zeros((rows, k), dtype=np.float64)
identity = np.eye(k)
for row in range(rows):
columns, values = matrix.row(row)
if columns.size == 0:
continue # cold start: no equation to solve
partner = fixed[columns] # (nnz, k)
residual = values - global_mean - bias_self[row] - bias_other[columns]
# lambda scaled by observation count: a row with three entries should
# be shrunk harder than one with three hundred.
gram = partner.T @ partner + regularization * columns.size * identity
updated[row] = np.linalg.solve(gram, partner.T @ residual)
return updated
def fit(
matrix: SparseMatrix,
transposed: SparseMatrix,
k: int = 32,
regularization: float = 0.05,
iterations: int = 15,
seed: int = 0,
) -> FactorModel:
"""Fit by ALS. Raises ValueError on malformed input."""
if matrix.shape[::-1] != transposed.shape:
raise ValueError(f"transpose shape {transposed.shape} does not match {matrix.shape}")
if k < 1:
raise ValueError(f"k must be at least 1, got {k}")
if regularization <= 0.0:
raise ValueError("regularization must be positive; without it sparse rows overfit exactly")
m, n = matrix.shape
rng = np.random.default_rng(seed)
global_mean = float(matrix.values.mean())
user_factors = rng.normal(0.0, 0.1, size=(m, k))
item_factors = rng.normal(0.0, 0.1, size=(n, k))
user_bias = np.zeros(m, dtype=np.float64)
item_bias = np.zeros(n, dtype=np.float64)
for _ in range(iterations):
# Biases first, and closed form: the optimal offset for a row is the
# mean of its residuals, shrunk by the regularizer.
for user in range(m):
columns, values = matrix.row(user)
if columns.size == 0:
continue
residual = values - global_mean - item_bias[columns] - (
user_factors[user] @ item_factors[columns].T
)
user_bias[user] = residual.sum() / (regularization + columns.size)
for item in range(n):
rows, values = transposed.row(item)
if rows.size == 0:
continue
residual = values - global_mean - user_bias[rows] - (
item_factors[item] @ user_factors[rows].T
)
item_bias[item] = residual.sum() / (regularization + rows.size)
user_factors = _solve_side(
matrix, item_factors, user_bias, item_bias, global_mean, regularization
)
item_factors = _solve_side(
transposed, user_factors, item_bias, user_bias, global_mean, regularization
)
return FactorModel(
user_factors=user_factors,
item_factors=item_factors,
user_bias=user_bias,
item_bias=item_bias,
global_mean=global_mean,
)What changed: The optimizer changes from SGD to alternating least squares, which is a change of kind rather than of style. The objective is a ridge regression in either factor with the other fixed, so each half-step has a closed form — a k-by-k solve per row — and converges in tens of sweeps rather than many noisy passes. It also makes every row independent within a half-step, which is the property the next stage exploits and which SGD structurally lacks. Around that: CSR indexing in both orientations so neither half-step searches for its entries, biases fitted in closed form rather than by gradient, regularization scaled by observation count so sparse rows shrink harder, and the cold-start case handled as a branch rather than as a division by zero.
O(|O| * k^2 + (m + n) * k^3) per sweep, with the per-row solves still driven from a Python loop.
Faster, and the justification for each change.
"""Implicit-feedback ALS - the Gramian trick, Cholesky solves, preallocated."""
import numpy as np
from numpy.typing import NDArray
from scipy.linalg import cho_factor, cho_solve
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
def solve_side_implicit(
indptr: NDArray,
indices: NDArray,
values: Vector,
fixed: Matrix,
regularization: float,
alpha: float,
) -> Matrix:
"""One implicit-feedback ALS half-step.
Implicit feedback has no negatives: every unobserved pair is a weak
negative with confidence 1, and every observed one a strong positive with
confidence 1 + alpha*r. Written naively that means each row's system sums
over ALL n items - O(m * n * k^2), which is impossible at catalogue scale.
The Gramian trick: split the confidence into a baseline of 1 everywhere
plus a correction on the observed entries only. The baseline part is
Y^T Y, which is the SAME k-by-k matrix for every row and is computed ONCE
per half-step. Only the correction depends on the row, and it touches
just that row's observed entries. Cost drops to O(|O| * k^2 + m * k^3).
"""
rows = indptr.size - 1
k = fixed.shape[1]
# Computed once, reused by every row: this single line is the whole trick.
gramian = fixed.T @ fixed + regularization * np.eye(k)
updated = np.empty((rows, k), dtype=np.float64) # allocated once
accumulator = np.empty((k, k), dtype=np.float64) # reused per row
target = np.empty(k, dtype=np.float64)
for row in range(rows):
start, stop = indptr[row], indptr[row + 1]
if start == stop:
updated[row] = 0.0 # cold start
continue
columns = indices[start:stop]
confidence = alpha * values[start:stop] # the correction, c - 1
partner = fixed[columns] # (nnz, k)
# accumulator = Y^T Y + Y_o^T diag(c - 1) Y_o + lambda I
np.copyto(accumulator, gramian)
accumulator += (partner * confidence[:, None]).T @ partner
# Target is (c - 1 + 1) on observed entries and 0 elsewhere, so the
# right-hand side touches only the observed rows.
np.matmul(partner.T, confidence + 1.0, out=target)
# Cholesky, not an explicit inverse: the system is symmetric positive
# definite by construction, and inverting it would be both slower and
# worse conditioned for a row with few observations.
updated[row] = cho_solve(cho_factor(accumulator, lower=True), target)
return updated
def fit_implicit(
indptr: NDArray,
indices: NDArray,
values: Vector,
t_indptr: NDArray,
t_indices: NDArray,
t_values: Vector,
n_users: int,
n_items: int,
k: int = 64,
regularization: float = 0.01,
alpha: float = 40.0,
iterations: int = 15,
seed: int = 0,
) -> tuple[Matrix, Matrix]:
rng = np.random.default_rng(seed)
# C-contiguous float64: every half-step forms Y^T Y over these, and a
# non-contiguous or float32 array would force a copy on each one.
user_factors = np.ascontiguousarray(rng.normal(0.0, 0.01, size=(n_users, k)))
item_factors = np.ascontiguousarray(rng.normal(0.0, 0.01, size=(n_items, k)))
for _ in range(iterations):
user_factors = solve_side_implicit(
indptr, indices, values, item_factors, regularization, alpha
)
item_factors = solve_side_implicit(
t_indptr, t_indices, t_values, user_factors, regularization, alpha
)
return user_factors, item_factorsWhat changed: The model changes to implicit feedback, which is what most real systems actually have, and that change introduces a cost problem the trick then solves. With no negatives, every unobserved pair is a weak negative, so each row’s normal equations naively sum over the entire catalogue — O(m·n·k²), which is impossible at scale. Splitting the confidence into a baseline of one everywhere plus a correction on observed entries only makes the baseline term the same k-by-k Gramian for every row: computed once per half-step, reused by all of them, with only the correction depending on the row. That single line takes the cost back to linear in the observed entries. The solves become Cholesky rather than explicit inverses, and the per-row buffers are allocated once.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
The shared Gramian is one GEMM per half-step rather than a per-row sum over the catalogue, and each row’s correction is a small weighted product over its observed entries only.
Trade-off: The identity depends on the confidence baseline being exactly uniform; any per-pair weighting scheme that breaks that assumption invalidates the shared term and returns the cost to O(m·n·k^2).
Replace a closed-form solve with a numerically stabler factorization
The k-by-k system is symmetric positive definite by construction, so a Cholesky factorization solves it at half the cost of general elimination and far better conditioned than forming an inverse.
Trade-off: Cholesky fails outright if the accumulated matrix loses positive definiteness through round-off on a row with very few observations, so the regularizer is load-bearing rather than merely helpful.
Pre-allocate output arrays and use in-place operations
The output factors, the per-row accumulator and the right-hand side are allocated once and rewritten per row, rather than allocating three arrays per row across millions of rows.
Trade-off: The buffers are shared mutable state, so the function is not reentrant and the row loop cannot be parallelized as written without giving each worker its own copies.
O(|O| * k^2 + (m + n) * k^3) per sweep, versus O(m * n * k^2) naively. Illustrative, not a measured benchmark.