Python实战：用NumPy手撕奇异值分解（SVD）及其在推荐系统中的应用

Python实战：用NumPy手撕奇异值分解（SVD）及其在推荐系统中的应用

当你打开视频平台时，首页总能精准推荐你感兴趣的内容；当你在电商网站浏览商品时，"猜你喜欢"栏目仿佛读懂了你的心思。这些智能推荐背后，往往隐藏着一个强大的数学工具——奇异值分解（SVD）。本文将带你用Python和NumPy从零实现SVD，并揭秘它在推荐系统中的神奇应用。

1. 奇异值分解的数学本质

奇异值分解是线性代数中一种强大的矩阵分解技术，它能将任意实数矩阵分解为三个特殊矩阵的乘积。给定一个m×n的实数矩阵A，其SVD形式为：

A = U @ Σ @ V.T

其中：

• U是m×m的正交矩阵，列向量称为左奇异向量
• Σ是m×n的对角矩阵，对角线元素σ₁≥σ₂≥...≥σₚ≥0称为奇异值
• V是n×n的正交矩阵，列向量称为右奇异向量

提示：正交矩阵的列向量是标准正交的，即U.T @ U = I，V.T @ V = I

奇异值分解的几何意义可以理解为：

1. 旋转（V.T）：将输入空间旋转到标准正交基
2. 缩放（Σ）：沿坐标轴进行不同程度的缩放
3. 旋转（U）：将结果旋转到输出空间

2. 从零实现SVD算法

2.1 理论基础与实现步骤

SVD的核心计算步骤如下：

1. 计算A.T @ A的特征值和特征向量
2. 对特征值排序并取其平方根得到奇异值
3. 构造右奇异矩阵V
4. 计算左奇异矩阵U

让我们用NumPy实现这个过程：

import numpy as np def svd_manual(A): # 计算A^T A的特征值和特征向量 ATA = A.T @ A eigvals, V = np.linalg.eig(ATA) # 确保特征值按降序排列 idx = eigvals.argsort()[::-1] eigvals = eigvals[idx] V = V[:, idx] # 计算奇异值 sigma = np.sqrt(eigvals) rank = np.sum(sigma > 1e-10) # 数值精度阈值 # 构建Σ矩阵 m, n = A.shape Sigma = np.zeros((m, n)) for i in range(rank): Sigma[i, i] = sigma[i] # 计算U矩阵 U = np.zeros((m, m)) for i in range(rank): U[:, i] = (A @ V[:, i]) / sigma[i] # 补充U的正交基 if rank < m: from scipy.linalg import null_space U[:, rank:] = null_space(U[:, :rank].T) return U, Sigma, V.T

2.2 验证我们的实现

让我们用一个简单矩阵测试我们的实现：

A = np.array([[1, 2, 0], [2, 0, 2]]) U, S, VT = svd_manual(A) print("U矩阵:\n", U) print("奇异值:\n", np.diag(S)) print("V转置矩阵:\n", VT) # 重构原始矩阵 reconstructed = U @ S @ VT print("重构矩阵:\n", reconstructed)

输出结果应与原始矩阵A非常接近，验证了我们实现的正确性。

3. SVD在推荐系统中的应用

3.1 推荐系统基础

推荐系统的核心是用户-物品评分矩阵，其中行代表用户，列代表物品，元素代表评分。这个矩阵通常是稀疏的，因为大多数用户只对少数物品有评分。

SVD通过以下方式帮助推荐：

1. 降维：提取最重要的特征
2. 去噪：过滤掉不重要的变异
3. 预测：在低维空间计算用户和物品的相似度

3.2 实战：电影推荐系统

让我们用MovieLens数据集构建一个简单的推荐系统：

import pandas as pd from scipy.sparse.linalg import svds # 加载数据 ratings = pd.read_csv('ratings.csv') movies = pd.read_csv('movies.csv') # 创建用户-电影评分矩阵 user_movie_ratings = ratings.pivot( index='userId', columns='movieId', values='rating' ).fillna(0) # 中心化数据（减去每个用户的平均评分） mean_user_ratings = user_movie_ratings.mean(axis=1) ratings_centered = user_movie_ratings.sub(mean_user_ratings, axis=0) # 执行SVD U, sigma, VT = svds(ratings_centered, k=50) sigma = np.diag(sigma) # 预测评分 predicted_ratings = U @ sigma @ VT + mean_user_ratings.values.reshape(-1, 1) predicted_ratings_df = pd.DataFrame( predicted_ratings, columns=user_movie_ratings.columns, index=user_movie_ratings.index ) def recommend_movies(user_id, n=5): # 获取用户未评分的电影 user_ratings = user_movie_ratings.loc[user_id] unrated = user_ratings[user_ratings == 0].index # 预测这些电影的评分 predictions = predicted_ratings_df.loc[user_id, unrated] # 返回评分最高的n部电影 top_movie_ids = predictions.sort_values(ascending=False).head(n).index return movies[movies['movieId'].isin(top_movie_ids)] # 为用户1推荐电影 recommend_movies(1)

4. SVD的高级应用与优化

4.1 截断SVD与计算效率

在实际应用中，我们通常只需要前k个最大的奇异值及其对应的奇异向量，这称为截断SVD：

from sklearn.decomposition import TruncatedSVD # 使用截断SVD svd = TruncatedSVD(n_components=50) reduced = svd.fit_transform(ratings_centered) print("解释方差比例:", svd.explained_variance_ratio_.sum())

截断SVD的优势：

• 计算效率高：只需计算部分奇异值
• 内存占用少：存储更小的矩阵
• 去噪效果：自动过滤小奇异值对应的噪声

4.2 处理大规模数据的增量SVD

对于超大规模矩阵，我们可以使用增量SVD：

from sklearn.decomposition import IncrementalPCA # 分批处理数据 batch_size = 1000 ipca = IncrementalPCA(n_components=50) for i in range(0, len(ratings_centered), batch_size): batch = ratings_centered.iloc[i:i+batch_size] ipca.partial_fit(batch) # 转换整个数据集 reduced = ipca.transform(ratings_centered)

4.3 SVD与其他技术的结合

现代推荐系统通常结合多种技术：

# 示例：结合SVD和k近邻 from sklearn.neighbors import NearestNeighbors # 在SVD降维后的空间计算相似用户 nbrs = NearestNeighbors(n_neighbors=5, metric='cosine').fit(reduced) # 为用户寻找邻居 user_index = 0 distances, indices = nbrs.kneighbors([reduced[user_index]]) # 基于邻居的评分推荐 neighbor_ratings = user_movie_ratings.iloc[indices[0]] recommendations = neighbor_ratings.mean(axis=0).sort_values(ascending=False)[:5]

在实际项目中，SVD的实现和调优需要根据具体场景进行调整。比如在电商推荐中，可能需要考虑时间衰减因素；在新闻推荐中，可能需要结合实时点击数据。