Factorization Machines in Go

Welcome to this lesson on factorization machines, an important model in the realm of recommendation systems. Factorization machines (FM) excel at capturing interactions between variables, making them a powerful tool for both regression and classification tasks. For instance, they can predict a rating (regression) or calculate the likelihood of a recommendation (classification).

Before we delve into the implementation of a factorization machine, let's briefly revisit the dataset preparation process from the previous lesson. In the last lesson, you learned how to load JSON files and create a user-item interaction matrix using dummy variables in Go. You also enriched the dataset with additional features, such as user preferences and genre similarity. These steps resulted in a data matrix where each row represents a user-item interaction, and columns represent features such as user and item dummy variables, user features, item features, and the rating. This structured data is essential for training a factorization machine.

Factorization machines leverage interactions between variables by decomposing them into simpler, latent factors. Mathematically, the prediction for a factorization machine can be expressed as:

$\hat{y}(\mathbf{x}) = w_0 + \sum_{i=1}^{n} w_i x_i + \sum_{i=1}^{n} \sum_{j=i+1}^{n} \langle \mathbf{v}_i, \mathbf{v}_j \rangle x_i x_j$

Here's what each component represents:

• $\hat{y}(\mathbf{x})$: The predicted value.
• $w_0$: The global bias term.
• $w_i$: The weight associated with the feature $x_i$.
• $x_i$: The individual features of the input vector $\mathbf{x}$.
• $\langle \mathbf{v}_i, \mathbf{v}_j \rangle$: The dot product between the latent vectors of two features, capturing their pairwise interaction.

Latent vectors are fundamental components in factorization machines used to capture complex pairwise interactions between features. Each column in the dataset is represented by a latent vector, and the interaction between different columns is determined by the dot product of these latent vectors.

For example, suppose your dataset has the following columns:

• user1, user2, user3, item1, item2, item3 (one-hot encoded user/item features)
• uf1, uf2 (user features)
• if1, if2 (item features)
• r (the rating)

Each feature column has an associated latent vector, which is initialized randomly and learned during training. For demonstration, let's consider a latent factor size of 2 for simplicity.

For a row with active features user2, item3, uf1, uf2, if1, and if2, the interaction between any two features (e.g., user2 and item3) is captured by the dot product of their latent vectors.

For example, if $\mathbf{v}_{u2} = [0.5, 0.3]$ and $\mathbf{v}_{i3} = [0.4, 0.7]$, then:

$\langle \mathbf{v}_{u2}, \mathbf{v}_{i3} \rangle = 0.5 \cdot 0.4 + 0.3 \cdot 0.7 = 0.2 + 0.21 = 0.41$

Latent vectors "encode" features, allowing compact representations of complex relationships beyond direct correlations. By updating these latent vectors during training, the factorization machine learns how different features interact, improving its predictive performance.

Let's move on to the implementation of the factorization machine model in Go. We'll break this into parts to ensure clarity.

package main

import (
    "math"
    "math/rand"
)

// SimpleFactorizationMachine represents the FM model.
type SimpleFactorizationMachine struct {
    NFactors     int         // Number of latent factors
    NFeatures    int         // Number of features
    LearningRate float64     // Learning rate
    Epochs       int         // Number of training epochs
    Reg          float64     // Regularization parameter

    W0 float64                // Global bias
    W  []float64              // Linear coefficients (size: NFeatures)
    V  [][]float64            // Latent factors (size: NFeatures x NFactors)
}

// NewSimpleFactorizationMachine initializes a new FM model.
// Pass a *rand.Rand with a fixed seed for reproducibility.
func NewSimpleFactorizationMachine(nFactors, nFeatures int, learningRate float64, epochs int, reg float64, rng *rand.Rand) *SimpleFactorizationMachine {
    V := make([][]float64, nFeatures)
    for i := range V {
        V[i] = make([]float64, nFactors)
        for f := range V[i] {
            V[i][f] = rng.NormFloat64() * 0.1 // Small random values
        }
    }
    return &SimpleFactorizationMachine{
        NFactors:     nFactors,
        NFeatures:    nFeatures,
        LearningRate: learningRate,
        Epochs:       epochs,
        Reg:          reg,
        W0:           0,
        W:            make([]float64, nFeatures),
        V:            V,
    }
}

In this code, we define a struct to hold all model parameters. The NewSimpleFactorizationMachine function initializes the model, including the latent factor matrix V with small random values. Slices are used to represent arrays and matrices.

To keep our code concise, consistent, and free from duplication, we introduce a helper method called predictRow. This method is responsible for computing the prediction for a single input vector, encapsulating both the linear and interaction terms as defined by the factorization machine (FM) formula. By centralizing this logic, we ensure that both training and prediction use exactly the same computation, which reduces the risk of errors and makes the code easier to maintain.

// predictRow computes the FM prediction for a single input vector.
func (fm *SimpleFactorizationMachine) predictRow(xi []float64) float64 {
    n := fm.NFeatures
    linearTerms := fm.W0
    for j := 0; j < n; j++ {
        linearTerms += fm.W[j] * xi[j]
    }
    interactionTerm := 0.0
    for f := 0; f < fm.NFactors; f++ {
        sumVx, sumVx2 := 0.0, 0.0
        for j := 0; j < n; j++ {
            vjf := fm.V[j][f]
            xij := xi[j]
            sumVx += vjf * xij
            sumVx2 += vjf * vjf * xij * xij
        }
        interactionTerm += (sumVx*sumVx - sumVx2) / 2.0
    }
    return linearTerms + interactionTerm
}

Let's break down how this method works:

• Linear Terms:
  The first part of the prediction is the linear component. This is calculated as the sum of the global bias (fm.W0) and the dot product of the feature coefficients (fm.W) with the input vector (xi). This captures the individual contribution of each feature to the prediction, similar to a standard linear regression model.

• Interaction Terms:
  The second part of the prediction captures the pairwise interactions between features using latent factors. For each latent factor (dimension), the method computes two quantities:

  • sumVx: The sum of the products of each feature value and its corresponding latent factor for the current dimension.
  • sumVx2: The sum of the squared products for each feature and its latent factor. The interaction term for each latent factor is then calculated as (sumVx*sumVx - sumVx2) / 2.0, which efficiently computes the sum of all pairwise interactions for that factor. The total interaction term is the sum over all latent factors.

• Centralization and Reusability:
  By encapsulating the prediction logic in this helper method, we ensure that both the training process (when updating parameters) and the prediction process (when making predictions on new data) use the exact same computation. This follows the DRY (Don't Repeat Yourself) principle, making the codebase easier to maintain and less prone to bugs.

This approach not only streamlines the code but also guarantees consistency and correctness throughout the model's lifecycle.

Before implementing the training method, let's briefly discuss gradient descent. Gradient descent is an optimization algorithm that minimizes a function by iteratively moving in the direction opposite to the gradient. In the context of our factorization machine, we adjust model parameters to minimize the error between predicted and actual outputs.

The update rule for a parameter $\theta$ is:

$\theta := \theta - \alpha \frac{\partial}{\partial \theta} J(\theta)$

where $\alpha$ is the learning rate, and $J(\theta)$ is the cost function.

To prevent overfitting, we add a regularization term to the cost function for the linear coefficients ($w$) and latent factors ($V$). Regularization discourages large parameter values, helping the model generalize better to unseen data. In our implementation, we do not regularize the global bias ($w_0$), as it simply captures the overall mean of the target variable and does not contribute to overfitting in the same way as the other parameters. Regularizing $w_0$ can unnecessarily restrict the model's ability to fit the data's mean, so it is typically left unregularized.

The regularization parameter (Reg in the code, often denoted as $\lambda$) controls the strength of regularization. Setting it too high can cause underfitting, while setting it too low can lead to overfitting. The optimal value depends on your dataset and should be determined empirically.

Next, let's implement the training logic for our factorization machine using Go slices and explicit loops. The Fit method will update the model parameters using gradient descent. Notice how we use the predictRow helper to compute predictions, which keeps the method concise and consistent.

// Fit trains the factorization machine on the given data.
func (fm *SimpleFactorizationMachine) Fit(X [][]float64, y []float64) {
    m, n := len(X), fm.NFeatures
    for epoch := 0; epoch < fm.Epochs; epoch++ {
        for i := 0; i < m; i++ {
            xi := X[i]
            pred := fm.predictRow(xi)
            err := pred - y[i]

            // Update global bias (w0) - no regularization
            fm.W0 -= fm.LearningRate * err

            // Update linear coefficients (w) with regularization
            for j := 0; j < n; j++ {
                gradW := err*xi[j] + fm.Reg*fm.W[j]
                fm.W[j] -= fm.LearningRate * gradW
            }

            // Update latent factors (V) with regularization
            for f := 0; f < fm.NFactors; f++ {
                for j := 0; j < n; j++ {
                    vjf, xij := fm.V[j][f], xi[j]
                    sumVx := 0.0
                    for jp := 0; jp < n; jp++ {
                        if jp != j {
                            sumVx += fm.V[jp][f] * xi[jp]
                        }
                    }
                    gradV := err*xij*sumVx + fm.Reg*vjf
                    fm.V[j][f] -= fm.LearningRate * gradV
                }
            }
        }
    }
}

• Prediction:
  For each training instance, we use the predictRow helper method to compute the current prediction based on the model's parameters. This ensures that the prediction logic is consistent and centralized.

• Error Calculation:
  The error (err) is calculated as the difference between the predicted value and the actual target value for the current instance.

• Global Bias Update (w0):
  The global bias term is updated by subtracting the product of the learning rate and the error. This term helps the model adjust for the overall average of the target variable.

• Linear Coefficient Updates (w):
  For each feature, we compute the gradient of the loss with respect to the linear coefficient. The gradient consists of two parts:

  • The error term scaled by the feature value (err * xi[j])
  • The regularization term (fm.Reg * fm.W[j]), which discourages large weights and helps prevent overfitting. The coefficient is then updated by moving in the direction opposite to the gradient.

• Latent Factor Updates (V):
  For each latent factor and each feature, we update the corresponding entry in the latent factor matrix. The gradient for each latent factor is:

  • The error term, scaled by the feature value and the sum of the products of the other features' latent factors and their values (err * xij * sumVx)
  • The regularization term (fm.Reg * vjf) The update ensures that the model learns how each feature interacts with every other feature through the latent factors.

• Efficiency and Maintainability:
  By centralizing the prediction logic in the predictRow helper, we avoid code duplication and reduce the risk of inconsistencies. This makes the code easier to maintain and less error-prone.

• Regularization:
  Regularization terms (fm.Reg * ...) are included in both the linear and latent factor updates to help prevent overfitting, especially when the number of features or latent factors is large.

This method iteratively updates all model parameters over multiple epochs, gradually reducing the prediction error on the training data. The explicit use of slices and loops makes the implementation clear and idiomatic in Go.

Now, let's implement the prediction logic as a method on the struct. This method will use the trained parameters to make predictions for new data. Again, we use the predictRow helper for each input row.

// Predict returns predictions for the given data.
func (fm *SimpleFactorizationMachine) Predict(X [][]float64) []float64 {
    m := len(X)
    yPred := make([]float64, m)
    for i := 0; i < m; i++ {
        yPred[i] = fm.predictRow(X[i])
    }
    return yPred
}

The code defines a Predict method for the SimpleFactorizationMachine struct, which generates predictions for a batch of input data. Here's what the code does, step by step:

• It takes a 2D slice X as input, where each row represents a feature vector for a user-item interaction or data point.
• It creates a slice yPred to store the predicted values, with the same length as the number of input rows.
• It loops over each row in X, and for each row, it calls the predictRow helper method to compute the prediction using the model's current parameters.
• The predicted value for each row is stored in the corresponding position in yPred.
• After processing all input rows, it returns the slice yPred containing the predictions for the entire dataset.

This method allows you to efficiently generate predictions for multiple data points at once, using the trained factorization machine model.

Before training machine learning models, it's important to standardize features to ensure they are on a similar scale. Standardization transforms features to have zero mean and unit variance, which helps gradient descent converge faster and prevents features with larger magnitudes from dominating the learning process.

The standardization formula for a feature $x$ is:

$z = \frac{x - \mu}{\sigma}$

where $\mu$ is the mean and $\sigma$ is the standard deviation of the feature.

Let's implement a Standardizer struct that can fit on training data and transform both training and test data:

type Standardizer struct {
    Means []float64
    Stds  []float64
}

func (s *Standardizer) Fit(features [][]float64) {
    if len(features) == 0 {
        return
    }
    numFeatures := len(features[0])
    numSamples := len(features)
    s.Means = make([]float64, numFeatures)
    s.Stds = make([]float64, numFeatures)
    
    for j := 0; j < numFeatures; j++ {
        sum := 0.0
        for i := 0; i < numSamples; i++ {
            sum += features[i][j]
        }
        s.Means[j] = sum / float64(numSamples)
    }
    
    for j := 0; j < numFeatures; j++ {
        sumSq := 0.0
        for i := 0; i < numSamples; i++ {
            diff := features[i][j] - s.Means[j]
            sumSq += diff * diff
        }
        variance := sumSq / float64(numSamples)
        s.Stds[j] = math.Sqrt(variance)
        if s.Stds[j] == 0 {
            s.Stds[j] = 1
        }
    }
}

func (s *Standardizer) Transform(features [][]float64) {
    if len(features) == 0 {
        return
    }
    numFeatures := len(features[0])
    numSamples := len(features)
    
    for j := 0; j < numFeatures; j++ {
        for i := 0; i < numSamples; i++ {
            features[i][j] = (features[i][j] - s.Means[j]) / s.Stds[j]
        }
    }
}

The Standardizer works in two steps:

• Fit: Computes the mean and standard deviation for each feature from the training data. If a feature has zero standard deviation (constant value), we set it to 1 to avoid division by zero.
• Transform: Applies the standardization transformation using the computed statistics. This method modifies the input data in-place.

Important: Always fit the standardizer on the training data only, then use those same statistics to transform both training and test data. This prevents data leakage and ensures the test set remains truly unseen.

Advanced note: In the examples we standardize all columns for brevity. In practice, you typically do not standardize one-hot (dummy) columns (e.g., userX, trackY). Keep them 0/1 to preserve sparsity and semantics. Standardize only continuous features (for example: track_likes, user_listening_avg, genre_similarity). Standardizing dummies turns many zeros into non-zero values (−μ/σ), which densifies the matrix and can hurt both efficiency and learning dynamics. If you want to try selective standardization here, standardize only the last 3 columns of the feature matrix produced by ExtractDataset (these are the continuous features).

Let's see how to use the factorization machine for training and evaluation in Go. We'll split the data into training and test sets, standardize the features, train the model, make predictions, and compute the Mean Absolute Error (MAE).

package main

import (
    "fmt"
    "math"
    "math/rand"
)

// splitData splits the dataset into training and test sets.
func splitData(X [][]float64, y []float64, testRatio float64, rng *rand.Rand) ([][]float64, [][]float64, []float64, []float64) {
    n := len(X)
    indices := rng.Perm(n)
    testSize := int(float64(n) * testRatio)

    XTest := make([][]float64, testSize)
    yTest := make([]float64, testSize)
    XTrain := make([][]float64, n-testSize)
    yTrain := make([]float64, n-testSize)

    for i, idx := range indices {
        if i < testSize {
            XTest[i] = X[idx]
            yTest[i] = y[idx]
        } else {
            XTrain[i-testSize] = X[idx]
            yTrain[i-testSize] = y[idx]
        }
    }
    return XTrain, XTest, yTrain, yTest
}

// meanAbsoluteError computes the MAE between true and predicted values.
func meanAbsoluteError(yTrue, yPred []float64) float64 {
    sum := 0.0
    for i := range yTrue {
        sum += math.Abs(yTrue[i] - yPred[i])
    }
    return sum / float64(len(yTrue))
}

func main() {
    seed := int64(42)
    rng := rand.New(rand.NewSource(seed))

    df := ExtractDataset()
    X := df.Features
    y := df.Labels

    // Split dataset into train and test sets
    XTrain, XTest, yTrain, yTest := splitData(X, y, 0.2, rng)

    // Standardize features
    std := &Standardizer{}
    std.Fit(XTrain)
    std.Transform(XTrain)
    std.Transform(XTest)

    // Initialize and train the factorization machine
    nFeatures := len(XTrain[0])
    fmModel := NewSimpleFactorizationMachine(3, nFeatures, 0.001, 1000, 0.01, rng)
    fmModel.Fit(XTrain, yTrain)

    // Make predictions on the test set
    yPred := fmModel.Predict(XTest)

    // Evaluate model using Mean Absolute Error
    mae := meanAbsoluteError(yTest, yPred)
    fmt.Printf("Mean Absolute Error: %.4f\n", mae)
}

The Mean Absolute Error (MAE) is defined as:

$\text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i|$

where $n$ is the number of observations, $y_i$ is the actual value, and $\hat{y}_i$ is the predicted value.

The code demonstrates a complete workflow for training and evaluating a factorization machine model:

• It first loads the dataset using ExtractDataset() (from the previous lesson).
• The splitData function randomly divides the data into training and test sets based on a specified ratio.
• The Standardizer is fitted on the training data and then used to transform both training and test sets, ensuring features are on the same scale.
• The factorization machine model is initialized with chosen hyperparameters and trained on the standardized training set using the Fit method.
• The trained model then predicts ratings for the test set using the Predict method.
• Finally, the meanAbsoluteError function calculates the average absolute difference between the predicted and actual test ratings, providing a measure of the model's performance.

Note: While MAE is a common metric for regression tasks and is useful for illustrative purposes here, in real-world recommender system evaluation, ranking metrics such as Recall@K, NDCG, or MAP are often more appropriate. These metrics better reflect the quality of the ranked recommendations presented to users, especially when the goal is to recommend top items rather than predict exact ratings.

In this lesson, we successfully implemented and evaluated a factorization machine model for recommendation systems in Go. We covered parameter initialization, feature standardization, training with gradient descent, making predictions, and evaluating model performance. This concludes our exploration of factorization machines and marks the end of this course module.

Congratulations on completing the course! The skills you've acquired here form a strong foundation for building and understanding recommendation systems. Continue exploring other models and refine your expertise in this dynamic field. Well done!