机器学习之矩阵分解：日撸Java三百行day69-70

1、什么是矩阵分解

矩阵分解（Matrix Factorization，MF）是推荐系统领域里的一种经典且应用广泛的算法。矩阵分解最初的想法是从奇异值分解（Singular Value Decomposition，SVD）借鉴而来。与其说是借鉴，不如直接称其为“伪奇异值分解”。在基于用户行为的推荐算法中，矩阵分解算法算的上是效果出众的方法之一，在推荐系统中发挥着重要作用。

从名字我们就可以了解到，该算法就是进行分解：将一个矩阵分解为两个或多个维数较低的矩阵。如：将一个规模为$m\ast n$的矩阵$R$分解为两个规模为$m\ast k$的矩阵$U$与规模为$n\ast k$的矩阵$V$。分解后的两个矩阵，使得两个矩阵内积后得到的值尽可能逼近矩阵$R$。在这里的内积，有点类似于M-distance推荐的处理方式，相似但又不相同。如果矩阵$R$表示电影评分，那么矩阵$U$代表用户对电影中某些环节的重视程度，而矩阵$V$代表电影对某些环节的满足程度。

2、矩阵分解算法思想

矩阵分解算法由SVD演变而来。传统的奇异值分解算法只能对数据稠密的矩阵进行分解，而评分矩阵是十分稀疏的。倘若要对缺失值进行填充，那么将导致算法复杂度上升、数据失真等问题。因此，为了能够对已有评分进行分解，出现了诸如BasicSVD、FunkSVD等矩阵分解方法。

·BasicSVD

BasicSVD是最简单的矩阵分解方法，它将评分矩阵分解为两个低阶矩阵，分别代表用户的隐含特征和物品的隐含特征，数学表达为：
$$\hat{R}=p_{uk}\cdot q_{ik}^T\text{\hspace{1em}(1)}$$
其中，$u$为用户数，$i$为物品数，$k$为隐含特征数。通过计算$R$与$\hat{R}$的误差平方和来来进行评分预测：
$$\sum _{u,i}(R-\hat{R}{)}^2=\sum _{u,i}(R-p_{uk}\cdot q_{ik}^T{)}^2\text{\hspace{1em}(2)}$$

·FunkSVD

FunkSVD是在传统SVD面临计算效率等问题时提出来的。FunkSVD在BasicSVD的基础上添加了正则化项，防止过拟合。同时，采用了线性回归的思想，使得用户的评分与项目满足度之间的成绩得到的评分残差尽可能的小。矩阵分解的数学表达为：
$$R_{m\ast n}=P_{m\ast k}^T{Q}_{k\ast n}\text{\hspace{1em}(3)}$$
对于某一个用户评分$m_{ij}$，若用FunkSVD进行矩阵分解，则对应的表示为$q_j^T{p}_i$，将方差作为损失函数，其表达式为：
$$min\sum _{i,j}(r_{ij}-q_j^T{p}_i{)}^2\text{\hspace{1em}(4)}$$
只要找到能使上面的式子取最小值对应的$u_i$和$v_j$，就能够得到矩阵$U$与$V$。

机器学习的本质是“猜”，在这里也不例外。首先我们通过随机初始化矩阵$U$与矩阵$V$，通过两个矩阵按照（3）式进行内积得到预测值。然后通过（4）式计算残差来调整矩阵中的每一个向量，直到将差值降低到某个程度。
我们可以通过RSME或MAE来评价本次预测的好坏。

·MSE

均方误差（Mean Square Error，MSE）是反映估计量与被估计量之间差异程度的一种度量。该值越小，预测效果越好。当该值为0时，说明预测值与真实值完全吻合。其数学表达式为：
$$MSE=\frac{1}{N}\sum _{i=1}^n(\widehat{y_i}-y_i{)}^2\text{\hspace{1em}(5)}$$
其中，$\widehat{y_i}$代表预测值，$y_i$代表真实值。

·RMSE

认识了MSE以后，RMSE就更好理解了。均方根误差（Root Mean Square Error，RMSE）就是MSE加了根号。可以更加直观的呈现预测值与实际值之间的差距。其数学表达式为：
$$RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^n(\widehat{y_i}-y_i{)}^2}\text{\hspace{1em}(6)}$$

·MAE

平均绝对误差（Mean Absolute Error）是所有单个预测值与实际值之差的绝对值的均值。MAE能够避免误差相互抵消的情况，更加准确地反映实际预测误差大小。其数学表达式为：
$$MAE=\frac{1}{N}\sum _{i=1}^n\left|\widehat{y_i}-y_i\right|\text{\hspace{1em}(7)}$$

介绍完预测评估标准后，再介绍一下如何通过梯度下降法进行优化。

·梯度下降

梯度下降法（Gradient descent）是一个一阶最优化算法，是求解无约束优化问题最简单、最经典的方法之一。要使用梯度下降法找到一个函数的局部最小值，必须向函数上当前点对应梯度的反方向的规定步长距离点进行迭代搜索。
梯度下降的基本过程与下山十分类似。当旅行者无法确定下山路径时，便可以通过梯度下降法来寻找路径。旅行者首先以他当前所在位置为基准，寻找这个位置最陡峭的地方，然后朝该方向走一步。然后，再以当前位置为基准，寻找最陡峭的地方，朝该方向走一步。以此类推，就能够下山了。

为了更加直观的理解梯度下降的思想，可以通过一个简单示意图理解：每一次前进都是通过寻找一个局部最优值来实现的。每一次前进，都能够朝着全局最优值逼近。当值不减反增时，就说明已经到达了最优值，此时必须停止优化。

小黑点想要在最短时间内逃出蓝色大椭圆，就需要先找到当前的“局部最优值”，先逃出最里层的白色椭圆。然后，就需要找到灰色圆中的“局部最优值”，逃出灰色圆。最后再在蓝色大椭圆中找到“局部最优值”。每一次前进，都在朝着全局最优值逼近。

其中，整个过程的最关键在于如何找到下降方向，也就是梯度的反向。之所以选择梯度的反方向主要是因为：梯度的方向是函数在指定点上升最快的方向，梯度的反方向就是下降最快的方向了。
对（4）式进行求偏导，来得到损失函数的负梯度：
$$\frac{\partial }{\partial \overrightarrow{P_u}}△r^2=-2(r_{u,v}-\overrightarrow{p_u}\overrightarrow{q_v})\overrightarrow{q_v}=-2△r\overrightarrow{q_v}\text{\hspace{1em}(8)}$$
$$\frac{\partial }{\partial \overrightarrow{q_u}}△r^2=-2(r_{u,v}-\overrightarrow{p_u}\overrightarrow{q_v})\overrightarrow{p_u}=-2△r\overrightarrow{p_u}\text{\hspace{1em}(9)}$$
由于为梯度的反向，所以加上了一个负号。同时，为了约束每一次下降的距离处于一个适中的值，需要引入学习率$\alpha$。学习率也被称为迭代的步长，梯度下降方向是不断变化的，因此需要该变量来控制步长。调节函数就变成了$-2\alpha △r\overrightarrow{q_v}$与$-2\alpha △r\overrightarrow{p_u}$。

三、算法的基本流程及操作

本实验数据集为电影评分表movielens，下载地址在此。数据量为100000。数据集共三种属性：user、movie、score。如：[0,64,4]表示用户0观看了64号电影，给出的评价为4分。

User	Movie	Score
0	64	4
0	65	4
…	…	…
942	1329	3

①初始化全局变量，读取数据文件，以及存储数据的三元组。

Random rand=new Random();//随机数
	int numUsers;//用户数量
	int numItems;//项目数量
	int numRatings;//评分数量
	Triple[] dataset;//三元组来存放数据
	double alpha;//学习参数
	double lamba;//控制学习速度
	int rank;//方阵排序
	double[][] userSubspace;//用户子空间矩阵
	double[][] itemSubspace;//项目子空间矩阵
	double ratingLowerBound;//评估值的下界
	double ratingUpperBound;//评估值的上界
	public class Triple{//存放数据的三元组
		public int user;//存放用户
		public int item;//存放项目
		public double rating;//存放评分
		public Triple() {//初始化
			user=-1;
			item=-1;
			rating=-1;
		}
		public Triple(int paraUser, int paraItem,double paraRating) {//初始化
			user=paraUser;
			item=paraItem;
			rating=paraRating;
		}
		public String toString() {
			return ""+user+", "+item+", "+rating;
		}
	}
		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);//读入数据
		}catch (Exception ee) {
			System.out.println("File " + paraFilename + " cannot be read! " + ee);
			System.exit(0);
			// TODO: handle exception
		}
	}

②生成数据集，并设置参数。

	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));//缓冲区读取内容
		
		//分配空间
		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]),Integer.parseInt(tempStringArray[2]));
			//分别代表：用户数量、项目数量、评分
		}
	}
	public void setParameters(int paraRank,double paraAlpha,double paraLamba) {
		rank=paraRank;
		alpha=paraAlpha;
	    lamba=paraLamba;
	}

③随机初始化子空间。

	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
	}

④预测，并且通过预测值与实际值之间的参加来进行梯度下降调整。

	public double predict(int paraUser,int paraItem) {//预测
		double resultValue=0;
		for(int i=0;i<rank;i++) {
			resultValue+=userSubspace[paraUser][i]*itemSubspace[paraItem][i];//用户期望值*项目满足值=预测值
		}
		return resultValue;//返回预测结果
	}
		private void updateNoRegular() {//调整
		for(int i=0;i<numRatings;i++) {
			int tempUserId=dataset[i].user;//读取用户id
			int tempItemId=dataset[i].item;//读取项目id
			double tempRate=dataset[i].rating;//读取评分
			
			double tempResidual=tempRate-predict(tempUserId, tempItemId);//计算残差
			
			//调整用户子空间
			double tempValue=0;
			for(int j=0;j<rank;j++) {
				tempValue=2*tempResidual*itemSubspace[tempItemId][j];//梯度下降，调整
				userSubspace[tempUserId][j]+=alpha*tempValue;//更新
			}//of for j
			
			//调整项目子空间
			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 updateNoRegualr

⑤评估本次预测水平，评估标准有：RSME、MAE

	public double rsme() {
		double resultRsme=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;//将最低阈值作为预测值
			}else if(tempPrediction>ratingUpperBound) {//若高于最高预测值
				tempPrediction=ratingUpperBound;//将最高阈值作为预测值
			}//of if
			
			double tempError=tempRate-tempPrediction;//差错
			resultRsme+=tempError*tempError;//将差错值累加
			tempTestCount++;//测试数量+1
		}//of for i
		return Math.sqrt(resultRsme/tempTestCount);
	}
	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;//将最低阈值作为预测值
			}else if(tempPrediction>ratingUpperBound) {//若高于最高预测值
				tempPrediction=ratingUpperBound;//将最高阈值作为预测值
			}//of if
			
			double tempError=tempRate-tempPrediction;//差错
			resultMae+=Math.abs(tempError);//取差错的绝对值
			tempTestCount++;//测试数量+1
		}//of for i
		return (resultMae/tempTestCount);
	}

⑥构建测试模型，输出最终的预测效果。

	public void train(int paraRounds) {//训练
		initializeSubspaces();//初始化子空间
		for(int i=0;i<paraRounds;i++) {
			updateNoRegular();//梯度下降调整
			if(i%50==0) {
				System.out.println("Round " + i);//输出当前轮数
				System.out.println("MAE: " + mae());
			}
		}
	}
		public static void testTrainingTesting(String paraFilename,int paraNumUsers,int paraNumItems,int paraNumRatings,double paraRatingLowerBound,double paraRatingUpperBound,int paraBounds) {
		try {
			//第一步：读入训练数据与测试数据
			MatrixFactorization tempMF=new MatrixFactorization(paraFilename, paraNumUsers, paraNumItems, paraNumRatings, paraRatingLowerBound, paraRatingUpperBound);//创建矩阵分类对象
			tempMF.setParameters(5, 0.0001, 0.005);//设置参数
			//第二步 调整+预测
			System.out.println("Begin Training !!!");
			tempMF.train(paraBounds);//进行训练
			
			double tempMAE=tempMF.mae();//计算MAE值
			double tempRSME=tempMF.rsme();//计算RSME值
			System.out.println("Finally, MAE = "+ tempMAE + ", RSME = " + tempRSME);
		} catch (Exception e) {
			e.printStackTrace();
			// TODO: handle exception
		}
	}
	public static void main(String args[]) {
		testTrainingTesting("C:/Users/11989/Desktop/sampledata-main/sampledata-main/movielens-943u1682m.txt", 943, 1682, 10000, 1, 5, 2000);//开始测试
	}