双线性模型

在推荐系统领域常常使用
含义是二元函数固定一个自变量时，函数关于另外一个自变量是线性的。

矩阵分解模型MF

如果我们需要对一个用户推荐电影我们应该怎么推荐？
结合用户以及用户对电影的评分，可以猜测用户的偏好，如果给某部电影分值很高，这个电影是爱情剧、战争、悬疑、哪个导演、什么时间等等都是特征，将每部电影的特征表示成向量，向量的每一个维度代表应该特征，其值越大就表示电影在某个特征程度高，还可以构建用户画像特征，因此我们就把用户和电影拆分成两个矩阵用户矩阵和电影矩阵。两个矩阵的乘积表示用户对电影的评分。

import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm # 进度条工具data = np.loadtxt('./data/movielens_100k.csv', delimiter=',', dtype=int)
print('数据集大小：', len(data))
# 用户和电影都是从1开始编号的，我们将其转化为从0开始
data[:, :2] = data[:, :2] - 1# 计算用户和电影数量
users = set()
items = set()
for i, j, k in data:users.add(i)items.add(j)
user_num = len(users)
item_num = len(items)
print(f'用户数：{user_num}，电影数：{item_num}')# 设置随机种子，划分训练集与测试集
np.random.seed(0)ratio = 0.8
split = int(len(data) * ratio)
np.random.shuffle(data)
train = data[:split]
test = data[split:]# 统计训练集中每个用户和电影出现的数量，作为正则化的权重
user_cnt = np.bincount(train[:, 0], minlength=user_num)
item_cnt = np.bincount(train[:, 1], minlength=item_num)
print(user_cnt[:10])
print(item_cnt[:10])# 用户和电影的编号要作为下标，必须保存为整数
user_train, user_test = train[:, 0], test[:, 0]
item_train, item_test = train[:, 1], test[:, 1]
y_train, y_test = train[:, 2], test[:, 2]
#%%
class MF:def __init__(self, N, M, d):# N是用户数量，M是电影数量，d是特征维度# 定义模型参数self.user_params = np.ones((N, d))self.item_params = np.ones((M, d))def pred(self, user_id, item_id):# 预测用户user_id对电影item_id的打分# 获得用户偏好和电影特征user_param = self.user_params[user_id]item_param = self.item_params[item_id]# 返回预测的评分rating_pred = np.sum(user_param * item_param, axis=1)return rating_preddef update(self, user_grad, item_grad, lr):# 根据参数的梯度更新参数self.user_params -= lr * user_gradself.item_params -= lr * item_grad
#%%
def train(model, learning_rate, lbd, max_training_step, batch_size):train_losses = []test_losses = []batch_num = int(np.ceil(len(user_train) / batch_size))with tqdm(range(max_training_step * batch_num)) as pbar:for epoch in range(max_training_step):# 随机梯度下降train_rmse = 0for i in range(batch_num):# 获取当前批量st = i * batch_sizeed = min(len(user_train), st + batch_size)user_batch = user_train[st: ed]item_batch = item_train[st: ed]y_batch = y_train[st: ed]# 计算模型预测y_pred = model.pred(user_batch, item_batch)# 计算梯度P = model.user_paramsQ = model.item_paramserrs = y_batch - y_predP_grad = np.zeros_like(P)Q_grad = np.zeros_like(Q)for user, item, err in zip(user_batch, item_batch, errs):P_grad[user] = P_grad[user] - err * Q[item] + lbd * P[user]Q_grad[item] = Q_grad[item] - err * P[user] + lbd * Q[item]model.update(P_grad / len(user_batch), Q_grad / len(user_batch), learning_rate)train_rmse += np.mean(errs ** 2)# 更新进度条pbar.set_postfix({'Epoch': epoch,'Train RMSE': f'{np.sqrt(train_rmse / (i + 1)):.4f}','Test RMSE': f'{test_losses[-1]:.4f}' if test_losses else None})pbar.update(1)# 计算 RMSE 损失train_rmse = np.sqrt(train_rmse / len(user_train))train_losses.append(train_rmse)y_test_pred = model.pred(user_test, item_test)test_rmse = np.sqrt(np.mean((y_test - y_test_pred) ** 2))test_losses.append(test_rmse)return train_losses, test_losses
#%%
# 超参数
feature_num = 16 # 特征数
learning_rate = 0.1 # 学习率
lbd = 1e-4 # 正则化强度
max_training_step = 30
batch_size = 64 # 批量大小# 建立模型
model = MF(user_num, item_num, feature_num)
# 训练部分
train_losses, test_losses = train(model, learning_rate, lbd, max_training_step, batch_size)plt.figure()
x = np.arange(max_training_step) + 1
plt.plot(x, train_losses, color='blue', label='train loss')
plt.plot(x, test_losses, color='red', ls='--', label='test loss')
plt.xlabel('Epoch')
plt.ylabel('RMSE')
plt.legend()
plt.show()

因子分解机FM

FM和MF有区别，FM是希望通过特征和用户点击物品的记录，预测用户点击其他物品的概率，点击率CTR，被点击和未点击是一个二分类问题，CTR的预测可以用逻辑斯蒂回归，但是逻辑斯蒂回归是线性化的，输入的特征之间是独立的，但是在现实中可能会有不独立的情况，如元宵和元宵节之间，灯笼和春节、对联之间，因此我们可以对逻辑斯蒂回归进行改进，引入双线性的部分“线性+双线性”。

用向量表示某一事物的离散特征时，用到独热编码，就是形如这样子的 一个向量
（x1，x2，x3，c1，c2,c3,c4,v1,v2,b1,b2,b3)
x类表示城市有3个城市，c表示月份有四个月份，v表示类型，b表示等级，有这个特征的值是1，没有的是0。
如一个物品来自城市x1，月份是c3，那么对应的值是1，其余是0。这就是一个独热编码，其稀疏性很高。

class FM:def __init__(self, feature_num, vector_dim):# vector_dim代表公式中的k，为向量v的维度self.theta0 = 0.0 # 常数项self.theta = np.zeros(feature_num) # 线性参数self.v = np.random.normal(size=(feature_num, vector_dim)) # 双线性参数self.eps = 1e-6 # 精度参数def _logistic(self, x):# 工具函数，用于将预测转化为概率return 1 / (1 + np.exp(-x))def pred(self, x):# 线性部分linear_term = self.theta0 + x @ self.theta# 双线性部分square_of_sum = np.square(x @ self.v)sum_of_square = np.square(x) @ np.square(self.v)# 最终预测y_pred = self._logistic(linear_term \+ 0.5 * np.sum(square_of_sum - sum_of_square, axis=1))# 为了防止后续梯度过大，对预测值进行裁剪，将其限制在某一范围内y_pred = np.clip(y_pred, self.eps, 1 - self.eps)return y_preddef update(self, grad0, grad_theta, grad_v, lr):self.theta0 -= lr * grad0self.theta -= lr * grad_thetaself.v -= lr * grad_v
#%%
# 超参数设置，包括学习率、训练轮数等
vector_dim = 16
learning_rate = 0.01
lbd = 0.05
max_training_step = 200
batch_size = 32# 初始化模型
np.random.seed(0)
model = FM(feature_num, vector_dim)train_acc = []
test_acc = []
train_auc = []
test_auc = []with tqdm(range(max_training_step)) as pbar:for epoch in pbar:st = 0while st < len(x_train):ed = min(st + batch_size, len(x_train))X = x_train[st: ed]Y = y_train[st: ed]st += batch_size# 计算模型预测y_pred = model.pred(X)# 计算交叉熵损失cross_entropy = -Y * np.log(y_pred) \- (1 - Y) * np.log(1 - y_pred)loss = np.sum(cross_entropy)# 计算损失函数对y的梯度，再根据链式法则得到总梯度grad_y = (y_pred - Y).reshape(-1, 1)# 计算y对参数的梯度# 常数项grad0 = np.sum(grad_y * (1 / len(X) + lbd))# 线性项grad_theta = np.sum(grad_y * (X / len(X) \+ lbd * model.theta), axis=0)# 双线性项grad_v = np.zeros((feature_num, vector_dim))for i, x in enumerate(X):# 先计算sum(x_i * v_i)xv = x @ model.vgrad_vi = np.zeros((feature_num, vector_dim))for s in range(feature_num):grad_vi[s] += x[s] * xv - (x[s] ** 2) * model.v[s]grad_v += grad_y[i] * grad_vigrad_v = grad_v / len(X) + lbd * model.vmodel.update(grad0, grad_theta, grad_v, learning_rate)pbar.set_postfix({'训练轮数': epoch,'训练损失': f'{loss:.4f}','训练集准确率': train_acc[-1] if train_acc else None,'测试集准确率': test_acc[-1] if test_acc else None})# 计算模型预测的准确率和AUC# 预测准确率，阈值设置为0.5y_train_pred = (model.pred(x_train) >= 0.5)acc = np.mean(y_train_pred == y_train)train_acc.append(acc)auc = metrics.roc_auc_score(y_train, y_train_pred) # sklearn中的AUC函数train_auc.append(auc)y_test_pred = (model.pred(x_test) >= 0.5)acc = np.mean(y_test_pred == y_test)test_acc.append(acc)auc = metrics.roc_auc_score(y_test, y_test_pred) test_auc.append(auc)print(f'测试集准确率：{test_acc[-1]}，\t测试集AUC：{test_auc[-1]}')