Background

This blog post was spurred by an interesting discussion with a coworker of mine. The question was relatively simple, “why do we use graphs for recommender systems, why not just use tabular ML?”. This question caused a sort of existential crisis in me, somewhat more than the average Murakami character (in the first 10 chapters).
We will use a strategy that is commonly used across mathematics, we will show that something is a simplified version of something else. In order to do this we need to have a principled way of thinking about recommender systems.

Basic Math of Recommender Systems

The most general way to think about a recommender system is to actually think of it in terms of products x users. This is either the adjacency matrix for the graph or the matrix to be factorized for the matrix factorization approach. So let us create this,

A 1 means the user watched the movie; a 0 means they were shown it and did not watch; NA means we have no information.

Matrix factorization approach

Here $U$ is a $5 \times 2$ matrix of movie embeddings and $V^\top$ is a $2 \times 4$ matrix of user embeddings, where 2 is the number of latent factors. The entries of $\hat{A}$ will not necessarily be 1s and 0s but rather as close as possible, minimising:

where $u_i$ and $v_j$ are the rows of $U$ and $V$ corresponding to movie $m_i$ and user $u_j$. Defining $e_{ij} = A_{ij} - u_i \cdot v_j$, gradient descent updates each embedding by nudging it in the direction that reduces the error:

where $\alpha$ is the learning rate.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

import numpy as np
from sklearn.decomposition import NMF


A_observed = np.array([
    [1,   np.nan, 1,   0  ],
    [np.nan, 1,   0,   1  ],
    [1,   0,   np.nan, 1  ],
    [0,   np.nan, 1,   1  ],
    [1,   1,   1,   np.nan],
], dtype=float)

# NAs are unobserved; fill with 0 for factorisation
A = np.nan_to_num(A_observed, nan=0.0)

rng = np.random.default_rng(42)
noise = rng.uniform(-0.15, 0.15, A.shape)
A = np.clip(A + noise, 0, 1)

model = NMF(n_components=2, init='nndsvda', random_state=0, max_iter=1000)
U  = model.fit_transform(A)
VT = model.components_

# Normalize columns of U to [0,1] and absorb the scale into VT
col_max = U.max(axis=0)
U  = U  / col_max
VT = VT * col_max[:, None]

A_hat = U @ VT

print("U =\n", np.round(U, 2))
print("\nV^T =\n", np.round(VT, 2))
print("\nU @ V^T =\n", np.round(A_hat, 2))
print("\nMax reconstruction error:", np.round(np.max(np.abs(A - A_hat)), 4))

Graph Convolutional Network Approach

In the GCN approach we treat movies and users as nodes in a bipartite graph. An edge connects movie $m_i$ to user $u_j$ whenever $A_{ij} = 1$. We stack these into a single $(5+4) \times (5+4)$ adjacency matrix:

The top-right block is the movie-to-user interactions from $A$; the bottom-left is its transpose. There are no movie-movie or user-user edges.

Before propagating information we symmetrically normalise $R$ by the degree matrix $D$, where $D_{ii} = \sum_j R_{ij}$:

A single GCN aggregation layer then updates every node’s embedding by averaging its neighbours’ and applying a learnable weight matrix $W^{(0)}$:

where $E^{(0)} \in \mathbb{R}^{9 \times d}$ is a matrix of initial node embeddings (one row per node, $d$ dimensions). After one pass, a movie node’s new embedding is a weighted average of the embeddings of the users who watched it, and vice versa. This is the key structural difference from matrix factorisation: instead of fitting two flat factor matrices independently, the GCN lets each node gather information from its neighbours before a prediction is made.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

import numpy as np
import torch
import torch.nn as nn

A_observed = np.array([
    [1,   np.nan, 1,   0  ],
    [np.nan, 1,   0,   1  ],
    [1,   0,   np.nan, 1  ],
    [0,   np.nan, 1,   1  ],
    [1,   1,   1,   np.nan],
], dtype=np.float32)

# Mask of observed entries (both 0s and 1s)
mask = ~np.isnan(A_observed)
A_clean = np.nan_to_num(A_observed, nan=0.0)

n_movies, n_users = A_clean.shape

# Build the (9 x 9) bipartite adjacency matrix
R = np.block([
    [np.zeros((n_movies, n_movies)), A_clean],
    [A_clean.T, np.zeros((n_users, n_users))],
]).astype(np.float32)

# Symmetric normalisation: D^{-1/2} R D^{-1/2}
deg = np.where(R.sum(axis=1) == 0, 1, R.sum(axis=1))
R_hat = np.diag(1 / np.sqrt(deg)) @ R @ np.diag(1 / np.sqrt(deg))
R_hat_t = torch.tensor(R_hat)

# Observed (movie, user) index pairs and their labels
rows, cols = np.where(mask)
y = torch.tensor(A_observed[rows, cols])

# Learnable parameters
torch.manual_seed(0)
E0 = nn.Parameter(torch.randn(n_movies + n_users, 2))
W  = nn.Parameter(torch.randn(2, 2))

optimizer = torch.optim.Adam([E0, W], lr=0.01)

for step in range(500):
    E1 = R_hat_t @ E0 @ W
    e_movie = E1[rows]
    e_user  = E1[n_movies + cols]
    y_hat   = (e_movie * e_user).sum(dim=1)
    loss    = ((y - y_hat) ** 2).mean()
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    if step % 100 == 0:
        print(f"step {step:4d}  loss={loss.item():.4f}")

print("\nE0 (learned):\n", np.round(E0.detach().numpy(), 3))
print("\nW  (learned):\n", np.round(W.detach().numpy(), 3))

What does the GCN result in?

In the matrix factorisation case we get $U$ and $V$, separate movie and user embedding matrices. The GCN produces the same thing, but because of the structure of $R$ both sets of embeddings live in a single matrix $E^{(0)}$: movie rows stacked on top of user rows. The rank of this matrix, like the number of latent factors in matrix factorisation, is the number of columns $d$. Both $E^{(0)}$ and $W^{(0)}$ are learned via backpropagation by minimising the error between the predicted and observed interactions:

where $e_{u_j}$ and $e_{m_i}$ are the rows of $E^{(1)}$ corresponding to user $u_j$ and movie $m_i$.

Recap

In both cases you get a matrix of vectors of shape $(m+n) \times d$, where each row is a latent representation of either a movie or a user. The idea is that if we can represent movies and users in a shared latent space such that similar movies and users are closer together, then predicting an interaction between a user and a movie they have not seen before is trivial, simply take the dot product of their latent representations and you have a number that tells you how likely they are to interact. Notice that you are learning these representations solely by looking at who watched which movies.

What if I have a representation that I want to use?

After learning that all the recommender system tries to do is find a latent representation of a user and a movie, you might say, hang on a minute, I actually know of an overt (as in NOT latent) representation that already places similar users close together. For example, suppose a user filled out a form asking what genres of movies they are interested in, and answered honestly. Would that not also place users with the same preferences closer together? With the added advantage that this is something they have explicitly told you?

That is a fair point. It turns out you can incorporate those explicit features alongside the latent ones, you learn weights for both the latent embeddings and the explicit features, and you also learn the embeddings themselves. The only thing you do not update are the feature columns in $X$, since those are representations you already trust. Think of it as freezing the columns you know are good and letting the model learn everything else around them.

Matrix case

Suppose each movie $m_i$ has a feature vector $x_{m_i} \in \mathbb{R}^p$ (e.g. genre indicators) and each user $u_j$ has a feature vector $x_{u_j} \in \mathbb{R}^q$ (e.g. stated genre preferences). We can augment the latent factor prediction with explicit features:

where $\beta_m \in \mathbb{R}^p$ and $\beta_u \in \mathbb{R}^q$ are learned coefficient vectors. The loss sums over all observed pairs (both 1s and 0s), where $\Omega$ denotes the set of observed entries in $A$:

and the gradient descent updates are:

where $e_{ij} = y_{ij} - \hat{y}_{ij}$. The structure is identical to the latent factor case. The only difference is that $x_{m_i}$ and $x_{u_j}$ are fixed observed features rather than vectors being learned from scratch.

GCN case

The GCN equivalent keeps the latent and feature streams separate, each with its own weight matrix:

The graph propagation enriches both streams: after one aggregation step, each node’s embedding is a weighted average of its neighbours’ latent and feature representations independently transformed before being summed.

Conclusion

The two prediction equations summarise everything, and the pure latent models fall out as special cases:

Using only the interaction matrix is not a different class of model. It is the full model with the feature coefficients constrained to zero ($\beta_m, \beta_u = 0$ in the MF case and $W_{\text{features}} = 0$ in the GCN case).

Appendix

What is $E^{(0)}$?

How does the loss function use $R$?

1. Graph structure for propagation. $\hat{R}$ defines which nodes aggregate from which neighbours at every forward pass, exactly as shown above.

2. Source of training signal. The entries of $R$ are the labels. For a given user $u_j$ and movie $m_i$, the predicted score is the dot product of their final embeddings:

and the ground truth is simply $y_{ij} = R_{ij} = 1$ for every observed interaction. The loss penalises the difference between the two across all observed pairs:

Minimising this loss via backpropagation shapes $E^{(0)}$ so that users end up close in embedding space to the items they have interacted with. $R$ is therefore both the adjacency matrix that drives message passing and the labelled dataset that supervises learning.