注意：本篇为50天后的Java自学笔记扩充，内容不再是基础数据结构内容而是机器学习中的各种经典算法。这部分博客更侧重于笔记以方便自己的理解，自我知识的输出明显减少，若有错误欢迎指正！
目录
前情提要
之前在54~55天介绍过老师和师姐发表的推荐算法论文——M-distance，当时就提到过推荐系统的设计，而今天的矩阵分解是也推荐系统的进一步运用与设计。但是矩阵分解采用的预测思想与M-distance具有较大差异：M-distance利用的相邻的潜在邻居的均值特征作为预测；而矩阵分解是基于随机矩阵的奇异值分解，有全局性与随机性，而且误差会更小。
一、算法简介
1.1 基本特性
所谓的矩阵分解其实就是将一个矩阵分解为两个或者多个低维矩阵的，这两个低维矩阵能够代表原矩阵特性并且预测原矩阵中未知的特性——在推荐系统矩阵中的描述就是：通过评估低维矩阵乘积来拟合评分矩阵。
沿用之前图，面对一个有m个用户与n个项目的稀疏的矩阵\(R_{m×n}\)，第i行表示第i个用户对于每个项目的评分，图中的问号部分表示这部分没有具体的评分；第j列表示某个项目不同用户给予的评分状况。
现在假设将此矩阵分解为两个或者多个矩阵的乘积，这里假设有两个低维矩阵\(P_{m×k}\)与\(Q_{k×n}\)，存在\(R_{m×n} = P_{m×k}Q_{k×n}\)（这里为了方便，我们将矩阵R中的那些未评分项目默认写为0）
这里采用的技术是 奇异值分解（Singular Value Decomposition，SVD）这个算法在降维、数据压缩、推荐系统中都有广泛的应用。但是！今天使用的SVD是一种基于机器学习的“ 伪 ”SVD，因为SVD本身是要求将原矩阵分解为可取的两个矩阵，但是这里的SVD我们是做近似分解(所以上图我用的约等于)，然后通过拟合让分解得到的矩阵的部分元素不断接近原矩阵的有效元素，最终使得原矩阵的无效元素(0元素)也能被拟合预测出来。这种SVD最早被Simon Funk在博客中发布，因此又叫做Funk-SVD。
具体的思路是先随机构建矩阵\(P\)与\(Q\)，然后通过矩阵\(P \times Q\)取得的每个数据同主要矩阵R的有效值进行差值计算，并通过得到的差值逐步调整\(P\)与\(Q\)的每个向量，如此反复操作直到这样的差值达到足够小的程度，或者说执行调整的次数满足某个上限后（比如MAE或者RSME在某个可接受的范围）
可以发现，在\(P\)矩阵中可以取出\(m\)个\(k\)维向量，在\(Q\)矩阵中可以取出\(n\)个\(k\)维向量。这俩矩阵所具有的这\(m+n\)个向量可以唯一表示R矩阵中的某个用户或者项目情况，而\(m \times n\)相互交叉可以得到对于某个用户-项目的评分的估计。于是不妨将那\(m\)个向量命名为用户向量，其余\(n\)个为项目向量。然后下面继续为这个定义建立数学模型：
取出稀疏矩阵中的实值定义为一个三元组数组，正如我们数据集的txt文本所示这样：
0,0,5
0,1,3
0,2,4
0,3,3
0,4,3
0,5,5
0,6,4
0,7,1
0,8,5
0,9,3
0,10,2
0,11,5
0,12,5
0,13,5
...分别表示：{用户，项目，评分}，这里的评分是确定的实值，不能为0（因为我们定义0为无实际评分），有集合\(S = \{(u,v)|r_{uv}≠\varnothing\}\)，这里\(u\)表示用户下标，\(v\)为项目下标，需要注意，集合\(S\)拥有的组合数目是要小于 用户数 * 项目数，集合\(S\)中的组合\((u,v)\)也并非全部\(u\)取值与\(v\)取值的枚举。
通过矩阵分解将用户向量与项目向量定义为：
- \(\vec p_{u} = (u_1,u_2,...,u_k)\)
- \(\vec q_{v} = (v_1,v_2,...,v_k)^{T}\)
因此不难得出对于用户\(u\)与物品\(v\)的预测评分\( \hat{r}_{uv} = \vec p_{u}\vec q_{v}\)。这里要注意一个细节，\(\vec p_{u} \)是面向矩阵\(P\)的向量，\(\vec q_{v} \)是面向\(Q\)的向量，而这俩矩阵本身是基于随机值生成的，因此得到的某个预测评分\(\hat{r}_{uv}\)有可能在原R矩阵中\(r_{uv} = 0\)，这个细节正是奇异值分解具有预测能力的原理。所以有等式\(R_{m \times n} \approx P_{m \times k} \times Q_{k \times n}=\hat{R}_{m \times n}\)，而算法正是通过\(\hat{R}_{m \times n}\)预测\(R_{m \times n}\)，预测准确度由\(r_{uv} ≠ 0\)时\(|\hat{r}_{uv} - r_{uv}|\)的大小保证。
那么这些向量是否具有一些实际含义呢，一种可以用来解释的含义就是：向量\(k\)维可以表示这个向量的语义对于\(k\)个隐含因子的重视程度。例如\(\vec p_{u} \)是面向用户的向量，假设\(k=3\)有\(\vec p_{u} = (u_1,u_2,u_3)\)，这里的\(u_1\)可以是用户觉得的电影情节是否精彩对于她来说的重要性，\(u_2\)可以是演员演技水平对于他来说的重要性......而对于\(\vec q_{v}\)，我们假设是电影，那么这里就可以是电影本身的精彩程度、电影演员水平......而得到\(\hat{r}_{uv}\)正式用户自己认为的重要性同电影本身这个属性的占比作乘积，然后各属性加权的结果。本身各维度在逻辑上是可以自由扩展的隐因子，不必解释也性，这里用语义去解释能更加方便理解吧。
另一个问题，特征空间\(k\)的大小应该设置为多大合适呢，这个并没有非常确定的设置，其本身也没有非常好的解释性，常常称之为" 隐因子向量 "，只用知道其存在即可，是用户向量\(\vec p_{u} \)与项目向量\(\vec q_{v} \)的基本维度。
二、损失函数与梯度下降(Gradient descent)
2.1 损失函数
通过刚刚分析，预测评分是\( \hat{r}_{uv} = \vec p_{u}\vec q_{v}\)，有预测误差\(\Delta r = \hat{r}_{uv} - r_{uv}((u,v) \in S)\)，随着预测与调整的循环，\(||\Delta r||\)会越来越小，自然对于某些\(\hat{r}_{uv}((u,v) \notin S) \)的预测会更接近于这个用户会实际给出的分数。
于是下面就是具体将应用问题向数学化模型的转换，我们将这种误差定义为平方损失函数\[ f(\vec p_{u},\vec q_{v})=\Delta r^{2} = \left(r_{u v}-\vec p_{u} \vec q_{v}\right)^{2} \tag{1}\] 最终我们需要求解的所有属于集合\(S\)的平方损失函数和的最小值：\[\min _{\vec p, \vec q} \sum_{(u, v) \in S} f(\vec p_{u},\vec q_{v})\] 也就是说求出所有观众真正评了分组合当中误差平方误差最小的那个作为评估标准值，这个标准值过大说明预测值与训练的实际值差别过大，即欠拟合。这里做平方是保证正取值。
2.2 梯度下降
确定了平方损失函数之后，沿着损失的方向利用梯度下降方法进行求解。
梯度下降法 是一个一阶最优化算法，通常也称为最速下降法。 要使用梯度下降法找到一个函数的局部最小值，必须向函数上当前点对应梯度（或者是近似梯度）的反方向的规定步长距离点进行迭代搜索。
梯度下降方法基于以下的观察：如果实值函数\(F(\mathbf{x})\)于点\(a\)处有可微定义，那么函数\(F(\mathbf{x})\)在\(a\)沿着梯度相反方向\(-\nabla F(\mathbf{a})\)下降最快。
因而，如果\[\mathbf{b}=\mathbf{a}-\gamma \nabla F(\mathbf{a})\] 对于\(\gamma>0\)为一个足够小的数时成立那么\(F(\mathbf{a}) \geq F(\mathbf{b})\)。考虑如此，可从函数\(F\)的局部最小值初始估计\(\mathbf{x}_0\)出发并考虑如下序列：\(\mathbf{x}_0,\mathbf{x}_1,\mathbf{x}_2...\)可得：\[\mathbf{x}_{n+1}=\mathbf{x}_{n}-\gamma_{n} \nabla F\left(\mathbf{x}_{n}\right), n \geq 0\]
简单来说，其实梯度下降法就是让变量沿着目标函数的负梯度方向进行下降直到局部最低点，因为我们当前这个问题本身就是求取损失函数的最小值，因此只要对损失函数求偏导然而确定一个两个向量（所代表\(P\)与\(Q\)向量）梯度下降的迭代式\(\vec p^{\prime}_{u}\)与\(\vec q^{\prime}_{v}\)即可。
就如同上图，从\(x_0\)得到逆梯度方向（正向）的变量位移移动\(x_1\)，但在\(x_1\)的到逆梯度方向（负向）的位移移动到\(x_2\)吗，如此反复直到一个稳定的局部最小。当然本文的案例中抽象出来是二维的，而并非此处一维的，此图于此只为说明此方法罢了。
首先对于平方损失函数分别求向量\(\vec p_{u}\)与\(\vec q_{v}\)方向的偏导，从而得到损失函数的负梯度：\[\frac{\partial}{\partial \vec p_{u}} \Delta r^{2}=-2\left(r_{u, v}-\vec p_{u} \vec q_{v}\right) \vec q_{v}=-2 \Delta r \vec q_{v} \tag{2}\] \[\frac{\partial}{\partial \vec q_{v}} \Delta r^{2}=-2\left(r_{u, v}-\vec p_{u} \vec q_{v}\right) \vec p_{u}=-2 \Delta r \vec p_{u} \tag{3}\] 然后依据负梯度方向更新\(P\)与\(Q\)的每个向量，同时自定义一个极小的参数\(\alpha > 0\)，这样便得到了基础的调节函数：\[\vec p^{\prime}_{u} = \vec p_{u} - \alpha \frac{\partial}{\partial \vec p_{u}} \Delta r^{2}=\vec p_{u}-2\alpha \Delta r \vec q_{v} \tag{4} \] \[\vec q^{\prime}_{v} = \vec q_{v} - \alpha \frac{\partial}{\partial \vec q_{v}} \Delta r^{2}=\vec q_{v}-2\alpha \Delta r \vec p_{u} \tag{5} \]
这里的\(\alpha\)其实有步长之功效，即是学习速率，也是梯度向量，需要调参确定。不同的步长对于地图下降的方式都会存在影响，当步长过小其下降速度会过慢：（下述动图来自：什么是梯度下降法？ - 知乎）
合适时，很快就能抵达谷底：
当这个值为1的时候，会在两侧同等移动，陷入僵局：
当这个值大于1时，反而会想希望到达的最小值的反方向走
2.3 加入正则项的损失函数
还有一种损失函数的评判标准，就是引入正则项\( \left\|\vec p_{u}\right\|^{2}+\left\|\vec q_{v}\right\|^{2}\)来进一步丰富损失函数的含义。之前的平方误差可以权衡模型的欠拟合程度，保证系统对于训练集的基本学习程度。而正则项就如同方差，能权衡模型的稳定性，如此一来：一旦损失函数\(f\)高也预示模型不稳定、对非训练集数据可能表现笨拙、过拟合。一般来说，在\(R_{n\times m}\)过于稀疏时会出现这种过拟合情况。
至于使用二范数限制而不是一范数限制似乎也是有一定说法(二范数：向量全维度的平方和，然后开根号；一范数：向量全维度绝对值和)，一范数会使得很多因子为0，削减模型大小，而二范数值会令因子接近0而不是为0。（参考论文：https://www.jstor.org/stable/2346178）
最终有损失函数：\[f_{2}(\vec p_{u},\vec q_{v}) = \left(r_{u v}-\vec p_{u} \vec q_{v}\right) ^{2}+\frac{\lambda}{2}\left(\left\|\vec p_{u}\right\|^{2}+\left\|\vec q_{v}\right\|^{2}\right) \tag{6}\] 其中映入了控制参数\(\frac{\lambda}{2}\)，正则项的获取数学细节可以参考老师的这篇博客：正则项的物理意义--概率矩阵分解的角度_闵帆的博客-CSDN博客。继续，模仿之前的梯度下降的方式，对于\(f_{2}\)求偏导可得损失函数的负梯度：\[\frac{\partial}{\partial \vec p_{u}} f_{2}=-2\left(r_{i, j}-\vec p_{u} \vec q_{v}\right) \vec q_{v}+\lambda \vec p_{u}=-2 \Delta r \vec q_{v}+\lambda \vec p_{u} \tag{7}\]\[\frac{\partial}{\partial \vec q_{v}} f_{2}=-2\left(r_{i, j}-\vec p_{u} \vec q_{v}\right) \vec p_{u}+\lambda \vec q_{v}=-2 \Delta r \vec p_{u}+\lambda \vec q_{v} \tag{8}\] 然后引入\(\alpha\)，构造原向量按照负梯度方向进行更新的调整函数：\[\vec p^{\prime}_{u} = \vec p^{\prime}_{u} - \alpha \frac{\partial}{\partial \vec p_{u}} f_{2}=\vec q^{\prime}_{u}-\alpha \left( 2\Delta r \vec q_{v}+\lambda \vec p_{u}\right) \tag{9}\]\[\vec q^{\prime}_{v} = \vec q^{\prime}_{v} - \alpha \frac{\partial}{\partial \vec q_{v}} f_{2}=\vec q^{\prime}_{v}-\alpha \left( 2\Delta r \vec p_{u}+\lambda \vec q_{v}\right) \tag{10}\]
如此便得到携带正则式的求解迭代式，这个迭代过程引入了方差的限制，控制了最终的\(P\)与\(Q\)矩阵预测的结果不会过拟合，但是实际测试结果来看，如果\(\lambda\)设置一般，这个限制对于最终结果似乎影响不大(此结果可待后续有不同意见时修改)
2.4 算法思路总结
- 确定k值之后，针对用户-项目矩阵\(R_{n\times m}\)生成矩阵\(P_{n \times k}\)与\(Q_{k \times m}\)，其中生成内容为随机数。
- 计算矩阵\(R_{n\times m}\)中的实际取值\(r_{uv}\)同\(P_{n \times k}Q_{k \times m}\)中得到的预测值\(\hat{r}_{uv}\)的差值
- 通过第二步取得的差值，利用公式4、5或者公式进行9、10去调整P矩阵的n个向量\(\vec p_{u}\)与\(Q\)矩阵的\(m\)个向量\(\vec q_{v}\)，得到全新的矩阵\(P\)与\(Q\)
- 重复2、3两步，直到最终MAE或者RSME收敛于一个预定义可靠的范围内。或者单纯循环若干次即可（我们选择后者）
三、代码描述
3.1 基本参数定义与构造函数、文件读取
矩阵分解的代码并不难，只要推导出梯度迭代的算式，在关键地方只需要将公式转换为代码即可。
/**
* Used to generate random numbers.
*/
Random rand = new Random();
/**
* Number of users.
*/
int numUsers;
/**
* Number of items.
*/
int numItems;
/**
* Number of ratings.
*/
int numRatings;
/**
* Training data.
*/
Triple[] dataset;除开随机对象，分别表示用户数目，项目数目，以及集合\(S\)的组合对象数目(已知评分个数)。最后一个是对于已知评分构成的三元组数组。
/**
* A parameter for controlling learning regular.
*/
double alpha;
/**
* A parameter for controlling the learning speed.
*/
double lambda;
/**
* The low rank of the small matrices.
*/
int rank;\(\alpha\)与\(\lambda\)是两个用于控制的参数，而rank表示\(P\)、\(Q\)矩阵的阶。因为\(K \ll \min \{M, N\}\)，故常称之为“ 低秩 ”
/**
* The user matrix P.
*/
double[][] userSubspace;
/**
* The item matrix Q.
*/
double[][] itemSubspace;
/**
* The lower bound of the rating value.
*/
double ratingLowerBound;
/**
* The upper bound of the rating value.
*/
double ratingUpperBound;分别表示奇异值分解时的分解矩阵\(P\)与矩阵\(Q\)，ratingLowerBound与ratingUpperBound是我们评分的下界和上界，当最终值超过这个范围，我们将默认将其设置为边界值（这个是针对问题做出的语义限制）
/**
************************
* The first constructor.
*
* @param paraFilename
* The data filename.
* @param paraNumUsers
* The number of users.
* @param paraNumItems
* The number of items.
* @param paraNumRatings
* The number of ratings.
************************
*/
public MatrixFactorization(String paraFilename, int paraNumUsers, int paraNumItems,
int paraNumRatings, double paraRatingLowerBound, double paraRatingUpperBound) {
numUsers = paraNumUsers;
numItems = paraNumItems;
numRatings = paraNumRatings;
ratingLowerBound = paraRatingLowerBound;
ratingUpperBound = paraRatingUpperBound;
try {
readData(paraFilename, paraNumUsers, paraNumItems, paraNumRatings);
// adjustUsingMeanRating();
} catch (Exception ee) {
System.out.println("File " + paraFilename + " cannot be read! " + ee);
System.exit(0);
} // Of try
}// Of the first constructor构造函数完成基本的初始化性赋值，其中涉及到读数据的方案：
/**
************************
* Read the data from the file.
*
* @param paraFilename
* The given file.
* @throws IOException
************************
*/
public void readData(String paraFilename, int paraNumUsers, int paraNumItems,
int paraNumRatings) throws IOException {
File tempFile = new File(paraFilename);
if (!tempFile.exists()) {
System.out.println("File " + paraFilename + " does not exists.");
System.exit(0);
} // Of if
BufferedReader tempBufferReader = new BufferedReader(new FileReader(tempFile));
// Allocate space.
dataset = new Triple[paraNumRatings];
String tempString;
String[] tempStringArray;
for (int i = 0; i < paraNumRatings; i++) {
tempString = tempBufferReader.readLine();
tempStringArray = tempString.split(",");
dataset[i] = new Triple(Integer.parseInt(tempStringArray[0]),
Integer.parseInt(tempStringArray[1]), Double.parseDouble(tempStringArray[2]));
} // Of for i
tempBufferReader.close();
}// Of readData任何代码的读数据的操作的关注都不在“读”这个本身，代码语言已经为我们封装好了，重点永远在于数据集到算法采用的数据结构的转换。这里的转换远不如M-distance麻烦，M-distance要动态处理稀疏矩阵与三元组之间的转化，需要使用一些信息辅助。而这里只需要读入到三元组即可，后续但凡有\(R_{n \times m}\)与\(\hat{R}_{n \times m}\)之间的验证都是全局性的，因此没必要一一取出稀疏矩阵某列或者某行信息。因此没必要设置繁杂的辅助数组。
这里为了方便，自定义了一个三元组：
public class Triple {
public int user;
public int item;
public double rating;
/**
*********************
* The constructor.
*********************
*/
public Triple() {
user = -1;
item = -1;
rating = -1;
}// Of the first constructor
/**
*********************
* The constructor.
*********************
*/
public Triple(int paraUser, int paraItem, double paraRating) {
user = paraUser;
item = paraItem;
rating = paraRating;
}// Of the first constructor
/**
*********************
* Show me.
*********************
*/
public String toString() {
return "" + user + ", " + item + ", " + rating;
}// Of toString
}// Of class Triple
为了多元化测试方便，并没有把\(alpha\)与\(lambda\)的声明放到构造函数中，而是单独设置了一个setter：
/**
************************
* Set parameters.
*
* @param paraRank
* The given rank.
* @throws IOException
************************
*/
public void setParameters(int paraRank, double paraAlpha, double paraLambda) {
rank = paraRank;
alpha = paraAlpha;
lambda = paraLambda;
}// Of setParameters3.2 分解矩阵定义与梯度调整
矩阵定义方面为了遍历方便，代码层面令矩阵\(Q\)转置。因为代码遍历是习惯从左上开始，行优先遍历的，这样是大多数语言的物理存储顺序。
/**
************************
* Initialize subspaces. Each value is in [0, 1].
************************
*/
void initializeSubspaces() {
userSubspace = new double[numUsers][rank];
for (int i = 0; i < numUsers; i++) {
for (int j = 0; j < rank; j++) {
userSubspace[i][j] = rand.nextDouble();
} // Of for j
} // Of for i
itemSubspace = new double[numItems][rank];
for (int i = 0; i < numItems; i++) {
for (int j = 0; j < rank; j++) {
itemSubspace[i][j] = rand.nextDouble();
} // Of for j
} // Of for i
}// Of initializeSubspaces
采用随机数生成这两个矩阵，\(k\)值是我们定义的分解矩阵秩的大小，往往来说，变量rank的大小是大概率小于numUsers与numItems的，因此这个rank大概率就是这个矩阵的秩大小。注意，rand.nextDouble()生成数是介于[0,1]的（大多数情况不会生成出边界值，跑几亿次都遇不到，因此你可以认为这是一个开区间）。
/**
************************
* Update sub-spaces using the training data.
************************
*/
public void updateNoRegular() {
for (int i = 0; i < numRatings; i++) {
int tempUserId = dataset[i].user;
int tempItemId = dataset[i].item;
double tempRate = dataset[i].rating;
double tempResidual = tempRate - predict(tempUserId, tempItemId); // Residual
// Update user subspace
double tempValue = 0;
for (int j = 0; j < rank; j++) {
tempValue = 2 * tempResidual * itemSubspace[tempItemId][j];
userSubspace[tempUserId][j] += alpha * tempValue;
} // Of for j
// Update item subspace
for (int j = 0; j < rank; j++) {
tempValue = 2 * tempResidual * userSubspace[tempUserId][j];
itemSubspace[tempItemId][j] += alpha * tempValue;
} // Of for j
} // Of for i
}// Of updateNoRegular首先完成了非无正则项的梯度方向的迭代代码。首先遍历全部的\(S\)集合（就是用户确确实实评了分的组合，即三元组），得到了任意基本评分\(\hat{r}_{uv}\)，并通过参数\(u\)，\(v\)确定了矩阵P中的向量\(\vec p_{u}\)与向量Q中的向量\(\vec q_{v}\)，从而得到\(r_{u, v}-\vec p_{u} \vec q_{v}\)。这就是第12行的含义，毫无疑问，predice就是计算\(\vec p_{u} \vec q_{v}\)的：
public double predict(int paraUser, int paraItem) {
double resultValue = 0;
for (int i = 0; i < rank; i++) {
// The row vector of an user and the column vector of an item
resultValue += userSubspace[paraUser][i] * itemSubspace[paraItem][i];
} // Of for i
return resultValue;
}// Of predict
继续刚刚updateNoRegular()的代码，后续的两个for循环就是对于公式4、5的代码翻译。updateNoRegular()代码就不赘述了，只需要将这两个for循环向公式9、10方向修改，引入变量lambda的使用（下述为附带正则项的代码）
// Update user subspace
double tempValue = 0;
for (int j = 0; j < rank; j++) {
tempValue = 2 * tempResidual * itemSubspace[tempItemId][j] - lambda * userSubspace[tempUserId][j];
userSubspace[tempUserId][j] += alpha * tempValue;
} // Of for j
// Update item subspace
for (int j = 0; j < rank; j++) {
tempValue = 2 * tempResidual * userSubspace[tempUserId][j] - lambda * itemSubspace[tempItemId][j];
itemSubspace[tempItemId][j] += alpha * tempValue;
} // Of for j3.3 循环训练与度量
所谓循环训练即我们判定收敛的手段是设置次数，这是最基本的判断收敛的方案，因为设置阈值需要经验指导，最初在不知道数据的情况下设置单纯的循环次数是省时省力的策略。
/**
************************
* Train.
*
* @param paraRounds
* The number of rounds.
************************
*/
public void train(int paraRounds) {
initializeSubspaces();
for (int i = 0; i < paraRounds; i++) {
updateNoRegular();
//updateWithRegular();
if (i % 50 == 0) {
// Show the process
System.out.println("Round " + i);
System.out.println("MAE: " + mae());
} // Of if
} // Of for i
}// Of train这里提到了MAE误差判断，后续给出代码
矩阵分解的主体代码：
/**
************************
* Compute the MAE.
*
* @return MAE of the current factorization.
************************
*/
public static void testTrainingTesting(String paraFilename, int paraNumUsers, int paraNumItems,
int paraNumRatings, double paraRatingLowerBound, double paraRatingUpperBound,
int paraRounds) {
try {
// Step 1. read the training and testing data
MatrixFactorization tempMF = new MatrixFactorization(paraFilename, paraNumUsers,
paraNumItems, paraNumRatings, paraRatingLowerBound, paraRatingUpperBound);
tempMF.setParameters(5, 0.0001, 0.005);
// Step 3. update and predict
System.out.println("Begin Training ! ! !");
tempMF.train(paraRounds);
double tempMAE = tempMF.mae();
double tempRSME = tempMF.rsme();
System.out.println("Finally, MAE = " + tempMAE + ", RSME = " + tempRSME);
} catch (Exception e) {
e.printStackTrace();
} // Of try
}// Of testTrainingTesting第13行是调用构造函数初始化，16行设置基本参数（\(\alpha = 0.0001, \lambda = 0.005, k = 5\)）；20行调用循环测试代码，最后22与24行输出度量信息MAE与RSME。
至于MAE与RSME的描述我在我之前M-distance推荐系统的博客中有详细阐述。简单回顾就是：平均绝对误差(Mean Abusolute Error: MAE) 对于每个\(\Delta r\)求和并且计算平均值，这里MAE的计算是针对实际的评分推荐系统而考虑的误差均值，因此算法预测的\(\hat{r}_{uv}\)要满足评分的上下限制(设置一个约束上下限的函数\(g\))，因此我们的\(\Delta r\)应该满足：\(\Delta r = r_{uv} - g(\hat{r}_{uv}) \)。在这样条件下有MAE表达式：\[\text{MAE}=\frac{\sum_{(u,v)\in S} \Delta r}{|S|} \tag{11}\] 而均方根误差(Root Mean Square Error: RMSE) 不过是另一种评价指标，但是其主要是用于反应偏差，其普遍大于MAE，有公式：\[\text{RMSE}= \sqrt{ \frac {\sum_{(u,v)\in S} \Delta r^{2}}{|S|}} \tag{12}\] 代码是对公式的复写。
/**
************************
* Compute the MAE.
*
* @return MAE of the current factorization.
************************
*/
public double mae() {
double resultMae = 0;
int tempTestCount = 0;
for (int i = 0; i < numRatings; i++) {
int tempUserIndex = dataset[i].user;
int tempItemIndex = dataset[i].item;
double tempRate = dataset[i].rating;
double tempPrediction = predict(tempUserIndex, tempItemIndex);
if (tempPrediction < ratingLowerBound) {
tempPrediction = ratingLowerBound;
} // Of if
if (tempPrediction > ratingUpperBound) {
tempPrediction = ratingUpperBound;
} // Of if
double tempError = tempRate - tempPrediction;
resultMae += Math.abs(tempError);
// System.out.println("resultMae: " + resultMae);
tempTestCount++;
} // Of for i
return (resultMae / tempTestCount);
}// Of mae
/**
************************
* Compute the MAE.
*
* @return MAE of the current factorization.
************************
*/
public static void testTrainingTesting(String paraFilename, int paraNumUsers, int paraNumItems,
int paraNumRatings, double paraRatingLowerBound, double paraRatingUpperBound,
int paraRounds) {
try {
// Step 1. read the training and testing data
MatrixFactorization tempMF = new MatrixFactorization(paraFilename, paraNumUsers,
paraNumItems, paraNumRatings, paraRatingLowerBound, paraRatingUpperBound);
tempMF.setParameters(5, 0.0001, 0.005);
// Step 3. update and predict
System.out.println("Begin Training ! ! !");
tempMF.train(paraRounds);
double tempMAE = tempMF.mae();
double tempRSME = tempMF.rsme();
System.out.println("Finally, MAE = " + tempMAE + ", RSME = " + tempRSME);
} catch (Exception e) {
e.printStackTrace();
} // Of try
}// Of testTrainingTesting四、数据测试
测试循环为2000次。
/**
************************
* @param args
************************
*/
public static void main(String args[]) {
testTrainingTesting("D:/Java DataSet/movielens-943u1682m.txt", 943, 1682, 10000, 1, 5, 2000);
}// Of main4.1 常规测试(包括\(\lambda\))
在令\(\alpha = 0.0001\)、\(k=5\)的情况下得到无正则项调整的结果：
令\(\alpha = 0.0001\)、\(\lambda = 0.005\)、\(k=5\)的情况下得到正则项调整的结果：
是否携带正则项似乎对最终结果并没有非常大的差异，但是通过对比之前M-distance的偏差结果，效果还是值得令人满意的（MAE: 0.77->0.50 RSME: 0.97->0.66）。(但是真的如此吗？具体见下面的4.3)
此外，为了更加明显观察到收敛的情况，我列出了 \(\alpha = 0.0001\)时两个度量偏差随着循环次数增加而收敛的情况：
可以发现，误差值在最开始都以非常快的速度收缩，而之后，RSME值在跌破1.0时开始减缓收敛，MAE则是在0.75后开始减缓。这样的信息似乎还不够，我们需要进一步了解\(k\)值与\(\alpha\)的变化对于训练结果可能的影响。
首先先观察梯度向量\(\alpha\)对于实验效果的影响：
4.2 步长\(\alpha\)测试
（为了方便观测，下面的测试不再测试RSME）固定\(k\)为5，可以发现随着循环增加，不同\(\alpha\)引导下的收敛明显存在差异。
随着 梯度向量\(\alpha\)的变大，收缩的速度也在加快，这印证了之前关于梯度下降的理论。同时因为我们采样样例有限（以50为间隔采样40次，总共循环2000次），对于那些收缩得比较急的曲线没能在短的间隔内统计足够多的采样点，可能显得不是那么平滑（比如\(\alpha = 0.001\)）
但是这个\(\alpha\)绝无可能设置无限大，一旦设置过大，对于大批量数据存在越界风险。经过测试，一旦设置过大，数据变为NaN的概率也增加了。
当然，好戏还在后面。最开始我也以为\(\alpha\)绝不能无限大是因为数据越界问题，但是，直到我大胆地尝试枚举\(\alpha>0.0001\)的所有情况时，奇妙的事情发生了，可见下图(我将数据为NaN的过滤了，所以右侧采样点显得有点稀疏)
效果非常amazing啊！我本来以为MAE会随着步长\(\alpha\)增加而越来越低，就如同刚刚的收敛图所示，但是结果却让我傻眼了。为了确保无误，削减随机影响，我测试了两次，但是效果是一致的。
我再度审视了一下梯度下降方法，对于这样的情况试做如下推测：正如2.2所示的梯度下降的诸多情况，其实对于凹函数的收敛过程中是存在某次收束“ 落入对称点 ”的情况。就像下面这个草图一样，我简单模拟了4次调整过程，其中橙色是\(\alpha\)较小的情况，红色是比较大的情况。通过4次下降可以发现，橙色的下落幅度是可能高于红色的，言下之意，在有限的下落过程，\(\alpha\)的大小与下落幅度并不是单调的关系。
其实不妨说，某一个\(\alpha\)取值得到的效果会在另一个\(\alpha\)取值中复现，比如下图的中假设的四种单次下落（图中从\(\alpha_0\)到\(\alpha_3\)的过程中，\(\alpha\)取值逐步变大），\(\alpha_0\)与\(\alpha_2\)就是单次下落过程中彼此的复现，\(\alpha_1\)是单次的最低下落，而\(\alpha_3\)的值相对接近1，当它再增大到1时，图中所示的下降线就会与横坐标平行，自然就无法下降；当其值大于1，下降变成上升。这与2.2中的理论不谋而合！
当然有了这个观点我们就可以通过枚举测试找到(固定k值后)下降最快的\(\alpha_1\)。那么显然本数据集最佳的\(\alpha\)应该是0.003。
4.3 k值测试与过拟合
见识了\(\alpha\)特性，接下来看看\(k\)值吧（代码中的rank变量，\(P\)与\(Q\)的秩）
为控制变量此处将\(\alpha\)固定为0.0001，然后枚举了各种k：
效果更加亮眼了，这个收敛效果更加明显！而且在\(k\)=100时，MAE最终竟然下降到几乎0.02左右！我猜测是对于每个向量携带的隐因子信息更多，对于事物的描述更加全面丰富，单次收敛能兼容更多的信息量。
那么...代价是什么？
1.开销
整个代码中唯一使用了稀疏矩阵的长宽是在初始化时，假设稀疏矩阵是一个\(n \times m\)的矩阵，那么初始化开销就是\(O(k(n+m))\)，当然因为初始化操作简单，这个过程开销还算可观。如果说循环总次数是\(M\)次，而稀疏矩阵中有效元素个数是\(N\)个的话，那么训练的开销就是\(O(kNM)\)。同初始化开销合并在一起可以发现总开销为\(O(k(n+m+N+M))\)，相信你已经发现了——\(k\)是一个全过程都存在的关键参数，以至于最后计算总开销时，甚至可以将其拎出来，那么自然随着\(k\)上升，整个算法的时间开销也为呈倍数增加。
从数据上能直观发现这个过程：
| k = 5 | 480ms |
| k = 10 | 633ms |
| k = 50 | 1770ms |
| k = 1000 | 3213ms |
所以\(k\)能带来误差值缩小，应该是因为它扩大了分解矩阵的信息量，以大量训练为代价实现了这种目的，是不同于\(\alpha\)那样基于数学算式实现的低复杂度优化。
2.过拟合
这个是当头一棒。刚刚提到了我的猜测，任务向量携带的隐因子信息更多，对于事物的描述更加全面丰富，单次收敛能兼容更多的信息量。但是反过来想，信息更多会不会限制更加严格呢？“若无所需，勿增实体”(奥卡姆剃刀原则)其实说的还是很有意义的，若把评分约束设置过多，那么突然出现违反常规的评分就会导致偏差过大，这种可能性是完全存在的。你可以想，假如说人类通过日复一日的生活经验而设计了各种条条框框的规则，这对于框架内的人自然不错，但是对于那些想一反常规去尝试自己生活的人来说何尝又不是禁锢？
于是改变分割训练集，使用经典的7:3的方式筛选出部分测试集(作为那些突破禁锢的人)，然后在训练的过程中阶段性地核查训练集的MAE与测试集的MAE，从而得到下面这样的图：
可以发现\(k\)变大后训练集的准确度疯狂提高其实暗藏了严重的过拟合，其中\(k\)为100时，最终MAE的极差尽然高达1.1左右，因此不要随便地提高k值！(这里通过测试，\(k\)为100的RSME极差会更大，接近1.4)
此外，请把视线移到常规的\(k\)为5的情况，最开始我们将\(k=5\)的Funk-SVD同M-distance做比较的时候感慨过Funk-SVD的MAE优秀结果，但是现在再来看呢？\(k\)为5时候的Funk-SVD测试集的MAE最低点大概在0.75~0.8左右，这个结果似乎与M-distance无异，而M-distance是采用leave-one-out得到的，是绝对公平预测。因此，你还能说Funk-SVD严格优于M-distance吗？
这种过拟合会随着拟合次数的增加而变得严重，并且无法通过调整\(\alpha\)实现补救，反而步长\(\alpha\)设置的越合适(刚刚通过测试，最佳的\(\alpha\)是0.003左右)，会导致训练集的拟合更快，而且测试的过拟合数据也会变快，体现到图像上就是测试集数据的MAE“ 翻脸 ”地更快！(曲线变化更加尖锐)
662
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
