Matrix Factorization for Recommendations

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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

import numpy as np
 
class BasicMatrixFactorization:
    """
    Basic matrix factorization for recommendation.
    
    The rating matrix R ≈ P @ Q.T where:
    - P[u, :] is the latent vector for user u
    - Q[i, :] is the latent vector for item i
    - Predicted rating r_hat[u, i] = P[u, :] @ Q[i, :]
    """
    
    def __init__(self, n_users: int, n_items: int, n_factors: int = 50):
        self.n_factors = n_factors
        
        # Initialize latent factor matrices with small random values
        # Using Xavier initialization for better convergence
        scale = 1.0 / np.sqrt(n_factors)
        self.P = np.random.normal(0, scale, (n_users, n_factors))
        self.Q = np.random.normal(0, scale, (n_items, n_factors))
    
    def predict(self, user_id: int, item_id: int) -> float:
        """Predict rating as dot product of latent vectors."""
        return np.dot(self.P[user_id], self.Q[item_id])
    
    def predict_all(self) -> np.ndarray:
        """Reconstruct the full rating matrix (for analysis)."""
        return self.P @ self.Q.T
    
    def get_user_vector(self, user_id: int) -> np.ndarray:
        """Get the latent representation for a user."""
        return self.P[user_id]
    
    def get_item_vector(self, item_id: int) -> np.ndarray:
        """Get the latent representation for an item."""
        return self.Q[item_id]
    
    def similar_items(self, item_id: int, top_k: int = 10) -> list:
        """Find most similar items based on latent space distance."""
        target = self.Q[item_id]
        # Compute cosine similarity to all items
        norms = np.linalg.norm(self.Q, axis=1)
        similarities = (self.Q @ target) / (norms * np.linalg.norm(target) + 1e-8)
        # Get top-k (excluding the item itself)
        most_similar = np.argsort(similarities)[::-1]
        return [idx for idx in most_similar if idx != item_id][:top_k]
 
 
# Example: Conceptual demonstration
def demonstrate_factorization():
    """Show how matrix factorization captures hidden structure."""
    
    # Original rating matrix (5 users, 4 movies)
    # Notice the implicit structure: users 0,2,4 prefer action; users 1,3 prefer romance
    R = np.array([
        [5, 2, 5, 1],  # Action fan
        [2, 5, 1, 5],  # Romance fan
        [5, 1, 5, 2],  # Action fan
        [1, 5, 2, 5],  # Romance fan
        [5, 3, 4, 2],  # Action fan (slightly mixed)
    ], dtype=float)
    
    # Perfect factorization with k=2 factors
    # Factor 0: "Action affinity", Factor 1: "Romance affinity"
    P_true = np.array([
        [1.0, 0.2],   # User 0: loves action
        [0.2, 1.0],   # User 1: loves romance
        [1.0, 0.1],   # User 2: loves action
        [0.1, 1.0],   # User 3: loves romance
        [0.9, 0.4],   # User 4: mostly action
    ])
    
    Q_true = np.array([
        [5.0, 1.0],   # Inception: action movie
        [1.0, 5.0],   # Titanic: romance movie
        [5.0, 1.0],   # Matrix: action movie
        [1.0, 5.0],   # Notebook: romance movie
    ])
    
    R_reconstructed = P_true @ Q_true.T
    
    print("Original R:")
    print(R)
    print("
Reconstructed (P @ Q^T):")
    print(R_reconstructed)
    print("
Reconstruction Error (Frobenius norm):")
    print(f"{np.linalg.norm(R - R_reconstructed):.4f}")
 
 
if __name__ == "__main__":
    demonstrate_factorization()

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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

import numpy as np
from typing import List, Tuple, Dict
 
def analyze_latent_factors(
    Q: np.ndarray,
    item_names: List[str],
    top_k: int = 5
) -> Dict[int, Tuple[List[str], List[str]]]:
    """
    Analyze what each latent factor captures by examining
    items with highest and lowest values for that factor.
    
    Args:
        Q: Item latent factor matrix (n_items x n_factors)
        item_names: Names of items for interpretation
        top_k: Number of items to show per factor
    
    Returns:
        Dictionary mapping factor index to (high_items, low_items)
    """
    n_factors = Q.shape[1]
    factor_analysis = {}
    
    for f in range(n_factors):
        # Get factor f values for all items
        factor_values = Q[:, f]
        
        # Find items with highest and lowest values
        sorted_indices = np.argsort(factor_values)
        high_indices = sorted_indices[-top_k:][::-1]
        low_indices = sorted_indices[:top_k]
        
        high_items = [item_names[i] for i in high_indices]
        low_items = [item_names[i] for i in low_indices]
        
        factor_analysis[f] = (high_items, low_items)
        
        print(f"
Factor {f}:")
        print(f"  High: {', '.join(high_items)}")
        print(f"  Low:  {', '.join(low_items)}")
    
    return factor_analysis
 
 
def find_factor_correlations(
    Q: np.ndarray,
    item_metadata: Dict[int, Dict[str, any]]
) -> Dict[int, str]:
    """
    Attempt to correlate latent factors with known metadata.
    
    This helps interpret what abstract factors might represent
    by correlating them with explicit attributes (genre, year, etc.)
    """
    correlations = {}
    
    for f in range(Q.shape[1]):
        factor_values = Q[:, f]
        
        # Correlate with various metadata fields
        best_correlation = 0
        best_attribute = "Unknown"
        
        for attr in ['genre_action', 'genre_romance', 'year', 'budget']:
            if attr in item_metadata.get(0, {}):
                attr_values = np.array([
                    item_metadata[i].get(attr, 0) 
                    for i in range(len(factor_values))
                ])
                corr = np.corrcoef(factor_values, attr_values)[0, 1]
                if abs(corr) > abs(best_correlation):
                    best_correlation = corr
                    best_attribute = f"{attr} ({'positive' if corr > 0 else 'negative'})"
        
        correlations[f] = f"{best_attribute} (r={best_correlation:.2f})"
    
    return correlations
 
 
def visualize_user_in_factor_space(
    P: np.ndarray,
    user_ids: List[int],
    factor_interpretations: List[str]
):
    """
    Print a user's position in the latent factor space
    with human-readable factor interpretations.
    """
    for user_id in user_ids:
        print(f"
User {user_id} preference profile:")
        user_vector = P[user_id]
        
        # Sort factors by absolute magnitude (most defining preferences first)
        sorted_factors = np.argsort(np.abs(user_vector))[::-1]
        
        for f in sorted_factors[:5]:  # Top 5 defining factors
            value = user_vector[f]
            interpretation = factor_interpretations[f] if f < len(factor_interpretations) else f"Factor {f}"
            direction = "strongly positive" if value > 1 else "positive" if value > 0 else "negative" if value > -1 else "strongly negative"
            print(f"  {interpretation}: {direction} ({value:.2f})")

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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

import numpy as np
from typing import List, Tuple
from dataclasses import dataclass
 
@dataclass
class Rating:
    """A single observed rating."""
    user_id: int
    item_id: int
    rating: float
 
def compute_mf_loss(
    P: np.ndarray,
    Q: np.ndarray,
    ratings: List[Rating],
    lambda_reg: float = 0.1
) -> Tuple[float, float, float]:
    """
    Compute the regularized matrix factorization loss.
    
    L = Σ (r_ui - p_u · q_i)² + λ(||P||² + ||Q||²)
    
    Returns:
        (total_loss, squared_error, regularization_term)
    """
    # Squared error over observed ratings
    squared_error = 0.0
    for rating in ratings:
        predicted = np.dot(P[rating.user_id], Q[rating.item_id])
        error = rating.rating - predicted
        squared_error += error ** 2
    
    # L2 regularization
    reg_P = lambda_reg * np.sum(P ** 2)
    reg_Q = lambda_reg * np.sum(Q ** 2)
    regularization = reg_P + reg_Q
    
    total_loss = squared_error + regularization
    
    return total_loss, squared_error, regularization
 
 
def compute_gradients(
    P: np.ndarray,
    Q: np.ndarray,
    user_id: int,
    item_id: int,
    rating: float,
    lambda_reg: float = 0.1
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Compute gradients of loss w.r.t. p_u and q_i for a single rating.
    
    For the loss term: (r_ui - p_u · q_i)² + λ(||p_u||² + ||q_i||²)
    
    ∂L/∂p_u = -2(r_ui - p_u · q_i) * q_i + 2λp_u
    ∂L/∂q_i = -2(r_ui - p_u · q_i) * p_u + 2λq_i
    """
    p_u = P[user_id]
    q_i = Q[item_id]
    
    error = rating - np.dot(p_u, q_i)
    
    # Gradients (without the factor of 2 for computational convenience)
    grad_p = -error * q_i + lambda_reg * p_u
    grad_q = -error * p_u + lambda_reg * q_i
    
    return grad_p, grad_q
 
 
def train_mf_sgd(
    n_users: int,
    n_items: int,
    ratings: List[Rating],
    n_factors: int = 50,
    learning_rate: float = 0.01,
    lambda_reg: float = 0.1,
    n_epochs: int = 100,
    verbose: bool = True
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Train matrix factorization model using stochastic gradient descent.
    
    This is the core algorithm used in production systems, with various
    optimizations (momentum, adaptive learning rates, etc.)
    """
    # Initialize latent factors
    scale = 1.0 / np.sqrt(n_factors)
    P = np.random.normal(0, scale, (n_users, n_factors))
    Q = np.random.normal(0, scale, (n_items, n_factors))
    
    for epoch in range(n_epochs):
        # Shuffle ratings for SGD
        np.random.shuffle(ratings)
        
        for rating in ratings:
            # Compute gradients for this rating
            grad_p, grad_q = compute_gradients(
                P, Q, rating.user_id, rating.item_id, rating.rating, lambda_reg
            )
            
            # Update latent vectors
            P[rating.user_id] -= learning_rate * grad_p
            Q[rating.item_id] -= learning_rate * grad_q
        
        # Compute and report loss periodically
        if verbose and (epoch + 1) % 10 == 0:
            loss, sq_err, reg = compute_mf_loss(P, Q, ratings, lambda_reg)
            rmse = np.sqrt(sq_err / len(ratings))
            print(f"Epoch {epoch + 1}: Loss = {loss:.4f}, RMSE = {rmse:.4f}")
    
    return P, Q

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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

import numpy as np
from typing import List, Tuple
from dataclasses import dataclass
 
@dataclass
class Rating:
    user_id: int
    item_id: int
    rating: float
 
class BiasedMatrixFactorization:
    """
    Matrix factorization with user and item biases.
    
    Prediction: r_hat = μ + b_u + b_i + p_u · q_i
    
    This model captures:
    - Global baseline (μ)
    - User-specific rating tendencies (b_u)
    - Item-specific quality (b_i)
    - User-item interaction (p_u · q_i)
    """
    
    def __init__(
        self,
        n_users: int,
        n_items: int,
        n_factors: int = 50,
        learning_rate: float = 0.01,
        reg_bias: float = 0.02,
        reg_factors: float = 0.02
    ):
        self.n_factors = n_factors
        self.lr = learning_rate
        self.reg_bias = reg_bias
        self.reg_factors = reg_factors
        
        # Biases
        self.global_mean = 0.0
        self.user_bias = np.zeros(n_users)
        self.item_bias = np.zeros(n_items)
        
        # Latent factors
        scale = 0.1
        self.P = np.random.normal(0, scale, (n_users, n_factors))
        self.Q = np.random.normal(0, scale, (n_items, n_factors))
    
    def predict(self, user_id: int, item_id: int) -> float:
        """Predict rating with biases and latent factors."""
        prediction = (
            self.global_mean +
            self.user_bias[user_id] +
            self.item_bias[item_id] +
            np.dot(self.P[user_id], self.Q[item_id])
        )
        # Clip to valid rating range
        return np.clip(prediction, 1.0, 5.0)
    
    def fit(self, ratings: List[Rating], n_epochs: int = 100):
        """Train the model using SGD."""
        
        # Compute global mean
        self.global_mean = np.mean([r.rating for r in ratings])
        
        for epoch in range(n_epochs):
            np.random.shuffle(ratings)
            
            for r in ratings:
                # Current prediction and error
                pred = self.predict(r.user_id, r.item_id)
                error = r.rating - pred
                
                # Update biases
                self.user_bias[r.user_id] += self.lr * (
                    error - self.reg_bias * self.user_bias[r.user_id]
                )
                self.item_bias[r.item_id] += self.lr * (
                    error - self.reg_bias * self.item_bias[r.item_id]
                )
                
                # Update latent factors
                p_u = self.P[r.user_id].copy()
                q_i = self.Q[r.item_id].copy()
                
                self.P[r.user_id] += self.lr * (
                    error * q_i - self.reg_factors * p_u
                )
                self.Q[r.item_id] += self.lr * (
                    error * p_u - self.reg_factors * q_i
                )
    
    def explain_prediction(self, user_id: int, item_id: int) -> dict:
        """
        Decompose a prediction into its components.
        
        Useful for debugging and interpretation.
        """
        return {
            "global_mean": self.global_mean,
            "user_bias": self.user_bias[user_id],
            "item_bias": self.item_bias[item_id],
            "interaction": np.dot(self.P[user_id], self.Q[item_id]),
            "final_prediction": self.predict(user_id, item_id)
        }
    
    def analyze_biases(
        self, 
        user_names: List[str] = None, 
        item_names: List[str] = None,
        top_k: int = 5
    ):
        """Analyze learned biases to understand rating patterns."""
        
        print("
=== Highest User Biases (Lenient Raters) ===")
        top_users = np.argsort(self.user_bias)[-top_k:][::-1]
        for u in top_users:
            name = user_names[u] if user_names else f"User {u}"
            print(f"  {name}: {self.user_bias[u]:+.2f}")
        
        print("
=== Lowest User Biases (Harsh Critics) ===")
        bottom_users = np.argsort(self.user_bias)[:top_k]
        for u in bottom_users:
            name = user_names[u] if user_names else f"User {u}"
            print(f"  {name}: {self.user_bias[u]:+.2f}")
        
        print("
=== Highest Item Biases (Universally Loved) ===")
        top_items = np.argsort(self.item_bias)[-top_k:][::-1]
        for i in top_items:
            name = item_names[i] if item_names else f"Item {i}"
            print(f"  {name}: {self.item_bias[i]:+.2f}")
        
        print("
=== Lowest Item Biases (Universally Disliked) ===")
        bottom_items = np.argsort(self.item_bias)[:top_k]
        for i in bottom_items:
            name = item_names[i] if item_names else f"Item {i}"
            print(f"  {name}: {self.item_bias[i]:+.2f}")