Building a Movie Recommender System

Dataset

As a dataset we chose the MovieLens 100K Dataset. It consists of 100.000 movie ratings from 943 different users on 1682 movies from the late 1990s. In the following, we will often refer to movies as items. The data is stored in different files:

• u.data (contains user_id, item_id (movie), rating and timestamp),

• u.user (contains user_id, age, gender, occupation and zip code),

• u.item (contains item_id, movie name, release date, ImDb link, genre in a hot-one format).

# Load the MovieLens 100k dataset
header = ['user_id', 'item_id', 'rating', 'timestamp']
ratings_df = pd.read_csv('movie_data/ml-100k/u.data', sep='\t', names=header)

pivot_df = ratings_df.pivot(index='user_id', columns='item_id', values='rating')

Cosine Similarity

Collaborative Filtering is an umbrella term for algorithms that generate recommendations by identifying patterns in user ratings. One of the most common techniques is User-Based Collaborative Filtering, where we calculate the similarity between users based on their rating history. To do this, we treat each user as a vector in a high-dimensional space and measure the “distance” between them using Cosine Similarity.

To calculate the cosine similarity between all users (rows) in a matrix \(X\), we normalize this matrix and then multiply it with its transposose.

If \(\hat{X}\) is the row-normalized version of \(X\) such that each row \(i\) is defined as \(\hat{x}_i = \frac{x_i}{\|x_i\|}\), then the entire Cosine Similarity matrix \(S\) is calculated via the gramian matrix:

Using NumPy, we perform these operations using vectorization for efficiency:

# L2 norm of each user vector
norms = np.linalg.norm(X, axis=1, keepdims=True)

# Normalize user vectors
X_norm = X / norms

# Cosine similarity = dot product of normalized vectors
user_similarity = X_norm @ X_norm.T

In SystemDS, we follow a similar logic. First, we import and initialize the SystemDSContext.

from systemds.context import SystemDSContext

with SystemDSContext() as sds:

In this context window, we load the data into SystemDS and do our calculations, using simple matrix functions: Matrix.

# Load into SystemDS
X = sds.from_numpy(X_np)

# Compute L2 norms
row_sums = (X * X).sum(axis=1)
norms = row_sums.sqrt()

# Normalize user vectors
X_norm = X / norms

# Cosine similarity = dot product of normalized vectors
user_similarity = X_norm @ X_norm.t()

In SystemDS, the line user_similarity_op = X_norm @ X_norm.t() does not execute any math. Instead, it creates an execution plan. The actual computation only occurs when we call .compute(), allowing SystemDS to optimize the entire operation.

user_similarity = user_similarity.compute()

In both cases user_similarity gives us a diagonal matrix that shows the similarity for every user-user pair.

While both methods produce the same results, SystemDS takes slightly longer for this specific dataset.

Matrix Factorization

Another powerful method for generating movie recommendations is Matrix Factorization. Instead of looking at surface-level data, this technique uncovers latent factors, the hidden patterns that represent a user’s specific tastes (like a preference for 90s rom-coms) and a movie’s unique characteristics (like its level of whimsy).

In a real-world scenario, our user-item interaction matrix \(R\) is incredibly sparse because most users have only rated a tiny fraction of the available movies. Matrix factorization solves this by decomposing \(R\) into two much smaller, lower-dimensional matrices:

• \(P\): Representing user preferences.

• \(Q\): Representing item characteristics.

By multiplying these two matrices back together, we can estimate the missing values in our original matrix:

To find \(P\) and \(Q\), we use the optimization algorithm called Alternating Least Squares (ALS).

In NumPy, we manually iterate through users and items, solving a least-squares problem for each. This gives us full control but can be computationally expensive as the dataset grows. We can compute \(\hat{R} = P \cdot Q^T\) like this (cf. CodeSignal):

# Random initialization of user and item factors
P = np.random.rand(num_users, rank) * 0.01
Q = np.random.rand(num_items, rank) * 0.01

for iteration in range(maxi):

    # Update user factors
    for u in range(num_users):

        # Get only items user 'u' actually rated
        user_mask = mask[u, :]
        Q_u = Q[user_mask, :]
        R_u = R[u, user_mask]

        if Q_u.shape[0] > 0:
            P[u, :] = np.linalg.solve(np.dot(Q_u.T, Q_u) + reg * np.eye(rank), np.dot(Q_u.T, R_u))

    # Update item factors
    for i in range(num_items):

        # Get only users who actually rated item 'i'
        item_mask = mask[:, i]
        P_i = P[item_mask, :]
        R_i = R[item_mask, i]

        if P_i.shape[0] > 0:
            Q[i, :] = np.linalg.solve(np.dot(P_i.T, P_i) + reg * np.eye(rank), np.dot(P_i.T, R_i))

R_hat = P @ Q.T

SystemDS allows us to execute the same logic using high-level script-like functions that are internally optimized. It offers a wide variety of built-in algorithms Algorithms, including ALS. First, we import our algorithm.

from systemds.operator.algorithm import als

Then, we initialize the SystemDS context:

with SystemDSContext() as sds:

To tune the model for our specific dataset, we configure the following hyperparameters:

• rank = 20 The number of latent factors (hidden features) used to describe users and movies. A higher rank allows for more complexity but increases the risk of overfitting.

• reg = 1.0 The regularization parameter. This prevents the model from becoming too complex by penalizing large weights, helping it generalize better to unseen data.

• maxi = 20 The maximum number of iterations. ALS is an iterative process.

Then we can do the computation.

# Load data into SystemDS
R = sds.from_numpy(pivot_df.fillna(0).values)

# Approximate factorization of R into two matrices P and Q using ALS
P, Q = als(R, rank=20, reg=1.0, maxi=20).compute()

R_hat = P @ Q

To test how well our models generalize to new data, we performed an 80/20 split, using the first 80,000 ratings for training and the remainder for testing. We compared both approaches based on execution time and Root Mean Squared Error (RMSE).

Both implementations are mathematically consistent. SystemDS achieved a slightly better Test RMSE.

Linear Regression

Unlike Matrix Factorization, which relies purely on interaction patterns, Linear Regression allows us to incorporate “side information” about users and items. By using features like user demographics and movie genres, we can build a Content-Based Filtering model that predicts ratings based on specific attributes.

ratings_df = pd.read_csv('movie_data/ml-100k/u.data', sep='\t', names=['user_id', 'item_id', 'rating', 'timestamp'])
user_df = pd.read_csv('movie_data/ml-100k/u.user', sep='|', names=['user_id', 'age', 'gender', 'occupation', 'zip_code'])
item_df = pd.read_csv('movie_data/ml-100k/u.item', sep='|', names=[
    'item_id', 'title', 'release_date', 'video_release_date', 'IMDb_URL', 'unknown', 'Action', 'Adventure', 'Animation',
    "Children's", 'Comedy', 'Crime', 'Documentary', 'Drama', 'Fantasy', 'Film-Noir', 'Horror', 'Musical', 'Mystery',
    'Romance', 'Sci-Fi', 'Thriller', 'War', 'Western'], encoding='latin-1')

# Turn categorical data into numerical
user_df['gender'] = user_df['gender'].apply(lambda x: 0 if x == 'F' else 1)
user_df = pd.get_dummies(user_df, columns=['occupation'])
item_df['release_date'] = pd.to_datetime(item_df['release_date'], errors='raise', format='%d-%b-%Y')
item_df['release_year'] = item_df['release_date'].dt.year

# Normalize data
user_df['age'] = (user_df['age'] - user_df['age'].min()) / (user_df['age'].max() - user_df['age'].min())
item_df['release_year'] = (item_df['release_year'] - item_df['release_year'].min()) / (item_df['release_year'].max() - item_df['release_year'].min())

X_train_tensor = torch.from_numpy(X_train).float()
y_train_tensor = torch.from_numpy(y_train).float().reshape(-1, 1)
X_test_tensor = torch.from_numpy(X_test).float()
y_test_tensor = torch.from_numpy(y_test).float().reshape(-1, 1)

n_features = X_train.shape[1]

class linearRegression(torch.nn.Module):
    def __init__(self):
        super(linearRegression, self).__init__()
        # input size: n_features, output size: 1
        self.linear = torch.nn.Linear(n_features, 1)

    def forward(self, x):
        out = self.linear(x)
        return out

lr_model = linearRegression()
criterion = torch.nn.MSELoss()
optimizer = torch.optim.SGD(lr_model.parameters(), lr = 0.01)

for epoch in range(1000):

    # Forward pass and loss
    pred_y = lr_model(X_train_tensor)
    loss = criterion(pred_y, y_train_tensor)

    # Backward pass and optimization
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

lr_model.eval()
with torch.no_grad():
    y_pred_test = lr_model(X_test_tensor)

y_pred = y_pred_test.numpy().flatten()
y_true = y_test_tensor.numpy().flatten()
rmse = np.sqrt(np.mean((y_pred - y_true) ** 2))

from systemds.operator.algorithm import lm, lmPredict

X_ds = sds.from_numpy(X_train)
y_ds = sds.from_numpy(y_train)
X_test_ds = sds.from_numpy(X_test)

model = lm(X=X_ds, y=y_ds)

predictions = lmPredict(X_test_ds, model).compute()