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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197

"""
Distributed Item-Item Similarity using PySpark
 
This implementation scales to billions of ratings across
hundreds of machines.
"""
 
from pyspark.sql import SparkSession
from pyspark.sql import functions as F
from pyspark.sql.window import Window
from pyspark.sql.types import StructType, StructField, IntegerType, FloatType
import numpy as np
 
def compute_item_similarities_distributed(
    ratings_df,  # DataFrame with columns: user_id, item_id, rating
    min_co_raters: int = 5,
    top_k: int = 50,
):
    """
    Compute top-k similar items for each item using distributed computation.
    
    Algorithm:
    1. Compute user means
    2. Create item pairs from co-ratings
    3. Compute adjusted cosine for each pair
    4. Keep top-k for each item
    """
    
    # Step 1: Compute user means
    user_stats = (
        ratings_df
        .groupBy("user_id")
        .agg(
            F.mean("rating").alias("user_mean"),
            F.stddev("rating").alias("user_std")
        )
    )
    
    # Join back to get deviations
    ratings_with_dev = (
        ratings_df
        .join(user_stats, "user_id")
        .withColumn("deviation", F.col("rating") - F.col("user_mean"))
        .select("user_id", "item_id", "rating", "deviation")
    )
    
    # Step 2: Self-join to create item pairs
    # For each user, join their ratings to themselves to get pairs
    ratings_left = ratings_with_dev.alias("left")
    ratings_right = ratings_with_dev.alias("right")
    
    item_pairs = (
        ratings_left
        .join(
            ratings_right,
            (F.col("left.user_id") == F.col("right.user_id")) &
            (F.col("left.item_id") < F.col("right.item_id"))  # Avoid duplicates
        )
        .select(
            F.col("left.item_id").alias("item_i"),
            F.col("right.item_id").alias("item_j"),
            F.col("left.deviation").alias("dev_i"),
            F.col("right.deviation").alias("dev_j"),
        )
    )
    
    # Step 3: Aggregate by item pair
    pair_stats = (
        item_pairs
        .groupBy("item_i", "item_j")
        .agg(
            F.sum(F.col("dev_i") * F.col("dev_j")).alias("dot_product"),
            F.sum(F.col("dev_i") ** 2).alias("norm_i_sq"),
            F.sum(F.col("dev_j") ** 2).alias("norm_j_sq"),
            F.count("*").alias("n_co_raters"),
        )
        .filter(F.col("n_co_raters") >= min_co_raters)
    )
    
    # Step 4: Compute adjusted cosine similarity
    similarities = (
        pair_stats
        .withColumn(
            "similarity",
            F.col("dot_product") / (
                F.sqrt(F.col("norm_i_sq")) * F.sqrt(F.col("norm_j_sq"))
            )
        )
        .filter(F.col("similarity") > 0)  # Keep positive similarities only
        .select("item_i", "item_j", "similarity", "n_co_raters")
    )
    
    # Step 5: Get top-k for each item
    # Need to consider both directions (i->j and j->i)
    sim_both_dirs = (
        similarities
        .select(
            F.col("item_i").alias("item"),
            F.col("item_j").alias("neighbor"),
            "similarity"
        )
        .union(
            similarities.select(
                F.col("item_j").alias("item"),
                F.col("item_i").alias("neighbor"),
                "similarity"
            )
        )
    )
    
    # Rank by similarity within each item
    window = Window.partitionBy("item").orderBy(F.desc("similarity"))
    
    top_k_neighbors = (
        sim_both_dirs
        .withColumn("rank", F.row_number().over(window))
        .filter(F.col("rank") <= top_k)
        .select("item", "neighbor", "similarity", "rank")
    )
    
    return top_k_neighbors
 
 
# Memory-efficient variant using blocking
def compute_similarities_blocked(
    ratings_df,
    n_blocks: int = 100,
    min_co_raters: int = 5,
):
    """
    Block-based similarity computation for very large item catalogs.
    
    Splits items into blocks and processes block pairs sequentially.
    Reduces memory from O(n²) to O(n²/n_blocks²) per worker.
    """
    
    # Assign items to blocks based on hash
    blocked_ratings = (
        ratings_df
        .withColumn(
            "item_block",
            F.abs(F.hash("item_id")) % n_blocks
        )
    )
    
    all_similarities = None
    
    # Process each block pair
    for block_i in range(n_blocks):
        for block_j in range(block_i, n_blocks):  # Upper triangle
            
            left = blocked_ratings.filter(F.col("item_block") == block_i)
            right = blocked_ratings.filter(F.col("item_block") == block_j)
            
            # Same logic as before but only for this block pair
            block_pairs = (
                left.alias("l")
                .join(
                    right.alias("r"),
                    (F.col("l.user_id") == F.col("r.user_id")) &
                    (
                        (block_i < block_j) |  # Different blocks: all pairs
                        (F.col("l.item_id") < F.col("r.item_id"))  # Same block: upper triangle
                    )
                )
                # ... compute similarities for this block pair
            )
            
            # Accumulate results
            if all_similarities is None:
                all_similarities = block_pairs
            else:
                all_similarities = all_similarities.union(block_pairs)
    
    return all_similarities
 
 
# Example cluster configuration
print("Example Spark Configuration for Item Similarity:")
print("=" * 50)
print("""
spark-submit \\
  --master yarn \\
  --deploy-mode cluster \\
  --num-executors 100 \\
  --executor-cores 4 \\
  --executor-memory 16g \\
  --driver-memory 8g \\
  --conf spark.sql.shuffle.partitions=1000 \\
  --conf spark.dynamicAllocation.enabled=true \\
  item_similarity.py
 
Expected performance (100M users, 10M items, 10B ratings):
- Time: 2-4 hours
- Cluster: 400 cores, 1.6TB RAM
- Output: ~500M item pairs with similarity > 0.1
""")

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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181

"""
Approximate Nearest Neighbors for Fast Similarity Lookup
 
This example shows how to use FAISS for item-based CF at scale.
"""
 
import numpy as np
import faiss
from typing import List, Tuple
import time
 
class FastItemNeighbors:
    """
    Fast approximate nearest neighbor lookup using FAISS.
    
    Suitable for:
    - 1M+ items
    - Real-time queries (<10ms)
    - 95%+ recall on true nearest neighbors
    """
    
    def __init__(
        self,
        embedding_dim: int = 128,
        n_lists: int = 1000,  # IVF clusters
        n_probe: int = 50,    # Clusters to search at query time
        use_gpu: bool = False,
    ):
        self.dim = embedding_dim
        self.n_lists = n_lists
        self.n_probe = n_probe
        self.use_gpu = use_gpu
        
        self.index = None
        self.item_ids = None
    
    def build_index(
        self, 
        item_embeddings: np.ndarray,
        item_ids: np.ndarray,
    ):
        """
        Build FAISS index from item embeddings.
        
        Args:
            item_embeddings: (n_items, embedding_dim) float32 array
            item_ids: (n_items,) array of item IDs
        """
        n_items = len(item_embeddings)
        self.item_ids = item_ids
        
        print(f"Building index for {n_items:,} items...")
        start = time.time()
        
        # Normalize for cosine similarity
        embeddings = item_embeddings.copy()
        faiss.normalize_L2(embeddings)
        
        # Create IVF index with flat (exact) quantizer
        quantizer = faiss.IndexFlatIP(self.dim)  # Inner product = cosine on normalized
        self.index = faiss.IndexIVFFlat(
            quantizer,
            self.dim,
            self.n_lists,
            faiss.METRIC_INNER_PRODUCT
        )
        
        # Train on subset if large
        if n_items > 100000:
            train_sample = embeddings[np.random.choice(n_items, 100000, replace=False)]
        else:
            train_sample = embeddings
        
        print("Training index...")
        self.index.train(train_sample)
        
        print("Adding vectors...")
        self.index.add(embeddings)
        
        # Set search parameters
        self.index.nprobe = self.n_probe
        
        # Move to GPU if available
        if self.use_gpu:
            res = faiss.StandardGpuResources()
            self.index = faiss.index_cpu_to_gpu(res, 0, self.index)
        
        print(f"Index built in {time.time() - start:.1f}s")
    
    def find_neighbors(
        self,
        item_ids: List[int],
        k: int = 50,
    ) -> List[List[Tuple[int, float]]]:
        """
        Find k nearest neighbors for given items.
        
        Returns:
            List of [(neighbor_id, similarity), ...] for each input item
        """
        # Get indices for requested items
        id_to_idx = {iid: idx for idx, iid in enumerate(self.item_ids)}
        indices = [id_to_idx[iid] for iid in item_ids if iid in id_to_idx]
        
        if not indices:
            return [[] for _ in item_ids]
        
        # Get embeddings and search
        # Note: In production, you'd store embeddings separately
        # Here we assume they're retrievable from the index
        
        # For IVF indices, we can reconstruct
        embeddings = np.zeros((len(indices), self.dim), dtype=np.float32)
        for i, idx in enumerate(indices):
            self.index.reconstruct(int(idx), embeddings[i])
        
        # Search
        faiss.normalize_L2(embeddings)
        similarities, neighbor_indices = self.index.search(embeddings, k + 1)
        
        # Convert to item IDs and exclude self
        results = []
        for i, (sims, nbrs) in enumerate(zip(similarities, neighbor_indices)):
            self_id = item_ids[i]
            neighbors = [
                (int(self.item_ids[nbr_idx]), float(sim))
                for nbr_idx, sim in zip(nbrs, sims)
                if nbr_idx >= 0 and self.item_ids[nbr_idx] != self_id
            ][:k]
            results.append(neighbors)
        
        return results
    
    def benchmark(self, n_queries: int = 1000):
        """Benchmark query performance."""
        if self.item_ids is None:
            print("Build index first!")
            return
        
        # Random query items
        query_ids = np.random.choice(self.item_ids, n_queries, replace=False)
        
        start = time.time()
        results = self.find_neighbors(list(query_ids), k=50)
        elapsed = time.time() - start
        
        print(f"\nBenchmark Results:")
        print(f"  Queries: {n_queries}")
        print(f"  Total time: {elapsed:.2f}s")
        print(f"  Per query: {1000 * elapsed / n_queries:.2f}ms")
        print(f"  QPS: {n_queries / elapsed:.0f}")
 
 
# Demo with synthetic data
if __name__ == "__main__":
    np.random.seed(42)
    
    # Simulate 1M items with 128-dim embeddings
    n_items = 100000  # Use 100K for demo
    dim = 128
    
    print(f"Generating {n_items:,} item embeddings...")
    embeddings = np.random.randn(n_items, dim).astype(np.float32)
    item_ids = np.arange(n_items)
    
    # Build index
    ann = FastItemNeighbors(
        embedding_dim=dim,
        n_lists=1000,
        n_probe=50,
    )
    ann.build_index(embeddings, item_ids)
    
    # Benchmark
    ann.benchmark(n_queries=1000)
    
    # Example query
    print("\nExample: Neighbors for item 0:")
    neighbors = ann.find_neighbors([0], k=5)[0]
    for item_id, sim in neighbors:
        print(f"  Item {item_id}: similarity = {sim:.3f}")

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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

"""
Scalable Item Similarity via Dimensionality Reduction
 
This approach can handle 10M+ items by:
1. Reducing to d=100-500 dimensions
2. Computing similarities in reduced space
3. Using batch operations for efficiency
"""
 
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import svds
from sklearn.metrics.pairwise import cosine_similarity
import time
 
class ScalableItemSimilarity:
    """
    Two-phase similarity computation:
    1. Offline: SVD reduction
    2. Similarity: Compute in reduced space
    """
    
    def __init__(
        self,
        n_factors: int = 100,
        chunk_size: int = 10000,
    ):
        self.d = n_factors
        self.chunk_size = chunk_size
        
        self.item_factors = None
        self.item_ids = None
    
    def fit_transform(
        self,
        ratings_matrix: csr_matrix,
        item_ids: np.ndarray,
    ):
        """
        Learn item factors via SVD.
        
        Args:
            ratings_matrix: n_items × n_users sparse matrix
            item_ids: Array of item IDs
        """
        self.item_ids = item_ids
        n_items = ratings_matrix.shape[0]
        
        print(f"Computing SVD with d={self.d} on {n_items:,} items...")
        start = time.time()
        
        # Truncated SVD
        U, sigma, Vt = svds(ratings_matrix.astype(float), k=self.d)
        
        # Item factors are U * sqrt(Sigma) for balanced representation
        self.item_factors = U * np.sqrt(sigma)
        
        # Normalize for cosine similarity
        norms = np.linalg.norm(self.item_factors, axis=1, keepdims=True)
        norms = np.where(norms > 0, norms, 1)  # Avoid division by zero
        self.item_factors /= norms
        
        print(f"SVD completed in {time.time() - start:.1f}s")
        print(f"Item factors shape: {self.item_factors.shape}")
        
        return self
    
    def compute_all_neighbors(
        self,
        k: int = 50,
    ) -> dict:
        """
        Compute top-k neighbors for all items using blocked computation.
        
        Memory-efficient: processes items in chunks rather than
        computing full n×n similarity matrix.
        """
        n_items = len(self.item_ids)
        neighbors = {}
        
        print(f"Computing {k}-nearest neighbors for {n_items:,} items...")
        start = time.time()
        
        # Process in chunks
        for chunk_start in range(0, n_items, self.chunk_size):
            chunk_end = min(chunk_start + self.chunk_size, n_items)
            
            # Compute similarities of this chunk against all items
            chunk_factors = self.item_factors[chunk_start:chunk_end]
            
            # Chunk × All similarity matrix
            sims = cosine_similarity(chunk_factors, self.item_factors)
            
            # Extract top-k for each item in chunk
            for i, row_sims in enumerate(sims):
                item_idx = chunk_start + i
                item_id = self.item_ids[item_idx]
                
                # Get top-k+1 (excluding self)
                top_indices = np.argpartition(row_sims, -k-1)[-k-1:]
                top_sims = row_sims[top_indices]
                
                # Sort by similarity
                sorted_order = np.argsort(-top_sims)
                top_indices = top_indices[sorted_order]
                top_sims = top_sims[sorted_order]
                
                # Exclude self and store
                neighbors[item_id] = [
                    (int(self.item_ids[idx]), float(sim))
                    for idx, sim in zip(top_indices, top_sims)
                    if idx != item_idx
                ][:k]
            
            if (chunk_start // self.chunk_size) % 10 == 0:
                print(f"  Processed {chunk_end:,} / {n_items:,} items")
        
        print(f"Neighbors computed in {time.time() - start:.1f}s")
        return neighbors
 
 
def compute_similarity_matrix_gpu(
    item_factors: np.ndarray,
    top_k: int = 50,
):
    """
    GPU-accelerated similarity computation using CuPy.
    
    Can process 1M+ items in seconds on modern GPUs.
    """
    try:
        import cupy as cp
        
        # Move to GPU
        factors_gpu = cp.asarray(item_factors)
        
        # Compute full similarity matrix in batches
        n_items = len(item_factors)
        batch_size = 10000
        
        all_neighbors = []
        
        for i in range(0, n_items, batch_size):
            batch = factors_gpu[i:i+batch_size]
            
            # Similarity: batch × all_items
            sims = cp.dot(batch, factors_gpu.T)  # Already normalized
            
            # Top-k per row
            top_k_indices = cp.argpartition(sims, -top_k, axis=1)[:, -top_k:]
            
            all_neighbors.append(cp.asnumpy(top_k_indices))
        
        return np.vstack(all_neighbors)
    
    except ImportError:
        print("CuPy not installed. Using CPU fallback.")
        return None
 
 
# Complexity comparison
print("Complexity Comparison:")
print("=" * 60)
print()
print("Direct similarity computation (n=10M items, m=100M users):")
print("  Time: O(n² · m) = 10^22 operations (infeasible)")
print()
print("With SVD reduction (d=100):")
print("  SVD: O(n · m · d) = 10^14 (hours)")
print("  Similarity: O(n² · d) = 10^16 (hours)")
print("  Total: Still challenging at this scale")
print()
print("With SVD + chunked computation + GPU:")
print("  SVD: Distributed (minutes)")
print("  Similarity: GPU-accelerated chunks (minutes)")
print("  Total: Practical for 10M+ items")