矩阵分解(java)

矩阵分解（java）

时间：2022/6/20

1.奇异值分解

对于任意实矩阵$A\in R^{m\ast n}$都可以分解为

$$A=U\sum V^T\text{\hspace{1em}(1)}$$

其中，$U\in R^{m\ast m}$是满足$U^{T}U=I$的m阶酋矩阵；$V\in R^{n\ast n}$是满足$V^{T}V=I$的n阶酋矩阵；$\sum \in R^{m\ast n}$是mxn的的矩阵，其中$（\sum {）}_{ii}={\sigma }_i$且其他位置的元素均为0，${\sigma }_i为非负实数且{\sigma }_1\ge {\sigma }_2\ge ...\ge 0$.

（1）式的分解便是奇异值分解(singular Value Decomposition,SVD)，矩阵U的列向量称为A的左奇异向量，矩阵V的列向量，被称为A的右奇异向量，${\sigma }_i$则称为奇异值(singular value).矩阵A的秩r就等于非零奇异值的个数。

2.矩阵分解——Funk-SVD

由矩阵乘法可得，对于一个矩阵A，可分解成两个子矩阵的乘积，即

$$A_{m\ast n}=P_{m\ast k}\ast Q_{k\ast n}\text{\hspace{1em}(2)}$$

通过这种方法可以将一个高阶矩阵分解成低阶矩阵。这种方法和奇异值分解很类似，常用于推荐系统中。通过将原始数据进行分解得到两个子矩阵，使用子矩阵的值不断拟合已知数据，达到一定目标之后，便可以用子矩阵的乘积去预测原始数据矩阵中的未知项。以此来达到推荐的目的。

这种方法是借鉴的奇异值分解的方法。可以说是伪奇异值分解。

为什么不使用奇异值分解呢？原因是奇异值分解主要适用于稠密矩阵，不适用于稀疏矩阵。同时奇异值分解也不能处理缺失值。所以我们选择改进的FUNK—SVD.

对于分解后的两个子矩阵，便是两个隐空间。隐空间的维度k是需要人工调整的，不断调整，直至达到满意的结果。对于k值和隐空间维度的一种解释便是，以观众对电影的评分为例，k值便是影响电影评分的相关参数，包括导演水平，演员演技，电影剧情等等。矩阵P则为用户矩阵，P的每一行对应于一个科幻对于电影相关参数的需求。而矩阵Q则是电影矩阵，Q的每一列对应于一个电影对相关参数的满足程度。$P\ast Q$便是该用户对该电影的评分。当然，这种分析看起来很合理，但并不是所有的情况都能这样解释。因为隐空间的维度其实很难进行合理的解释，算法的可解释性是很差的。

这样我们就很容易得到算法的具体流程：

1. 首先使用随机数按照给定的k值，构造两个子矩阵。

2. 使用子矩阵通过式（2）对已知数据进行拟合。

3. 计算预测误差，调整子矩阵。

4. 重复（2）（3）两个步骤，使得结果满足一定条件结束。例如迭代多少轮或者误差小于某个给定的阈值。

这里有个问题，如何使得调整之后的矩阵快速到达我们满意的结果。这里可以采用梯度下降法。

3.梯度下降与损失函数

3.1.损失函数

对于FUNK-SVD算法，用户评分矩阵$R_{m\ast n}$被分解成$P_{m\ast k},Q_{k\ast n}$两个矩阵，则相应的预测评分则为：

$${\hat{r}}_{ij}=p_i\ast q_j^T\text{\hspace{1em}(3)}$$

则对于的平方损失函数则为：

$$f(u,i)=\Delta r^2=(r_{u,i}-{\hat{r}}_{u,i}{)}^2\text{\hspace{1em}(4)}$$

对于数据集S而言，最终的优化目标便是：

$$\underset{p_u,q_i}{\min }\sum _{(u,i)\in S}f(u,i)\text{\hspace{1em}(5)}$$

使得所有项的损失最小。

3.2.梯度下降

如何使得算法快速迭代至我们满意的情况呢？这里便要使用梯度下降法。

梯度下降法（gradient descent）是一种常用的一阶优化方法，是求解无约束优化问题最简单、最经典的方法之一。

对于无约束优化问题${\min }_{x}f(x)$，其中$f(x)$是连续可微函数。欲构造一个数列$x_0,x_1,x_2,...$满足：

$$f(x_{t+1})<f(x_t),t=0,1,2,...\text{\hspace{1em}(6)}$$

使得函数不断递减，直至趋于极小点，欲满足上式，则根据泰勒展开式有：

$$f(x+\Delta x)\approx f(x)+\Delta x^T\nabla f(x)\text{\hspace{1em}(7)}$$

故欲使得$f(x+\Delta x)<f(x)$可选择：

$$\Delta x=-\gamma \nabla f(x)\text{\hspace{1em}(8)}$$

其中$\gamma$是步长，由人为调整参数。通过选取合适的步长，可以快速的达到函数的极小点。

对于本次实验，对损失函数（4）按上式进行计算偏导数，可以得到：

$$\frac{\partial \Delta r^2}{\partial {\vec{p}}_u}=-2\Delta r\ast {\vec{q}}_i\text{\hspace{1em}(9)}$$
$$\frac{\partial \Delta r^2}{\partial {\vec{q}}_i}=-2\Delta r\ast {\vec{p}}_u\text{\hspace{1em}(10)}$$

故：

$$\begin{gathered} \Delta x_p=\gamma \ast 2\ast \Delta r\ast {\vec{q}}_i \\ \Delta x_q=\gamma \ast 2\ast \Delta r\ast {\vec{p}}_u\text{\hspace{1em}(11)} \end{gathered}$$

据此式可以调整分解后的两个矩阵$P,Q$.

4.算法实现

/**
 * MyMatrixFactorization.java
 *
 * @author zjy
 * @date 2022/6/19
 * @Description: 矩阵分解
 * @version V1.0
 */
package swpu.zjy.ML.MartixFactorization;

import java.io.BufferedReader;
import java.io.FileReader;
import java.util.Random;

public class MyMatrixFactorization {
    //数据集，按照user-item-rating存储
    Triple[] dataset;
    //用户数量
    int numUser;
    //电影数量
    int numItem;
    //评分数量
    int numRating;
    //随机数生成器
    Random random=new Random();
    //控制学习规则的参数
    double alpha;
    //用于控制学习速率的参数
    double lambda;
    //隐空间维度k
    int rank;
    //用户子空间
    double[][] userSubspace;
    //电影项目子空间
    double[][] itemSubspace;
    //预测评分最低边界
    double ratingLowerBound;
    //预测评分最高边界
    double ratingUpperBound;

    /**
     * 三元组类，存储数据元组。user-item-rating
     */
    public class Triple{
        int user;
        int item;
        double rating;
        public Triple(int user,int item,double rating){
            this.user=user;
            this.item=item;
            this.rating=rating;
        }

        @Override
        public String toString() {
            return "Triple{" +
                    "user=" + user +
                    ", item=" + item +
                    ", rating=" + rating +
                    '}';
        }
    }

    /**
     * 构造方法，初始化相应参数，读入数据集。
     * @param datasetFileName 数据集地址
     * @param numUser 用户数量
     * @param numItem 电影数量
     * @param numRating 评分数量
     * @param ratingLowerBound 最低评分
     * @param ratingUpperBound 最高评分
     */
    public MyMatrixFactorization(String datasetFileName,int numUser,
                                 int numItem,int numRating,double ratingLowerBound,double ratingUpperBound){
        //初始化
        this.numItem=numItem;
        this.numUser=numUser;
        this.numRating=numRating;
        this.ratingLowerBound=ratingLowerBound;
        this.ratingUpperBound=ratingUpperBound;
        //读入数据集
        try {
            FileReader fileReader=new FileReader(datasetFileName);
            BufferedReader bufferedReader=new BufferedReader(fileReader);
            dataset=new Triple[numRating];
            String tempString;
            String[] tempData;
            for (int i = 0; i < numRating; i++) {
                tempString=bufferedReader.readLine();
                tempData=tempString.split(",");
                dataset[i]=new Triple(Integer.parseInt(tempData[0]),Integer.parseInt(tempData[1]),
                        Double.parseDouble(tempData[2]));
            }

            bufferedReader.close();
        } catch (Exception e) {
            e.printStackTrace();
        }
    }

    /**
     * 设置参数
     * @param paraRank 隐空间维度k
     * @param paraAlpha 梯度下降参数
     * @param paraLambda
     */
    public void setParameters(int paraRank, double paraAlpha, double paraLambda) {
        rank = paraRank;
        alpha = paraAlpha;
        lambda = paraLambda;
    }

    /**
     * 初始化用户子空间和项目子空间，使用随机数填入
     */
    public void initalizeSubspace(){
        //初始化两个子空间
        userSubspace=new double[numUser][rank];
        itemSubspace=new double[numItem][rank];
        //用随机数填入
        for (int i = 0; i < numUser; i++) {
            for (int j = 0; j < rank; j++) {
                userSubspace[i][j]=random.nextDouble();
            }
        }

        for (int i = 0; i < numItem; i++) {
            for (int j = 0; j < rank; j++) {
                itemSubspace[i][j]=random.nextDouble();
            }
        }
    }

    /**
     * 使用子空间进行预测用户对项目的评分
     * @param userId 用户ID
     * @param itemId 项目ID
     * @return 预测评分
     */
    public double predict(int userId,int itemId){
        double result=0;
        //两个矩阵相乘
        for (int i = 0; i < rank; i++) {
            result+=userSubspace[userId][i]*itemSubspace[itemId][i];
        }
        return result;
    }

    /**
     * 按照梯度下降调整子空间
     */
    public void updateNoRegular(){
        for (int i = 0; i < numRating; i++) {
            int tempUserId=dataset[i].user;
            int tempItemId=dataset[i].item;
            double tempRating=dataset[i].rating;
            //计算误差
            double tempError=tempRating-predict(tempUserId,tempItemId);
            //更新矩阵
            //使用梯度下降法
            double tempValue=0;
            for (int j = 0; j < rank; j++) {
                tempValue=2*tempError*itemSubspace[tempItemId][j];
                userSubspace[tempUserId][j]+=tempValue*alpha;
            }

            for (int j = 0; j < rank; j++) {
                tempValue=2*tempError*userSubspace[tempUserId][j];
                itemSubspace[tempItemId][j]+=tempValue*alpha;
            }
        }
    }

    /**
     * 使用已有数据集进行训练
     * @param paraRounds 训练次数
     */
    public void train(int paraRounds){
        initalizeSubspace();
        for (int i = 0; i < paraRounds; i++) {
            updateNoRegular();
            //输出MAE
            if(i%50==0){
                System.out.println("Round="+i);
                System.out.println("MAE="+MAE());
            }
        }
    }

    /**
     * 计算均方根误差
     * @return 均方根误差
     */
    private double RSME(){
        double resultRSME=0;
        int tempTestCount=0;
        for (int i = 0; i < numRating; i++) {
            int tempUserId=dataset[i].user;
            int tempItemId=dataset[i].item;
            double tempRating=dataset[i].rating;
            double tempPredict=predict(tempUserId,tempItemId);
            //防止评分越界
            if(tempPredict<ratingLowerBound){
                tempPredict=ratingLowerBound;
            }else if(tempPredict>ratingUpperBound){
                tempPredict=ratingUpperBound;
            }
            double tempError=tempRating-tempPredict;
            resultRSME+=tempError*tempError;
            tempTestCount++;
        }
        return Math.sqrt(resultRSME/tempTestCount);
    }

    /**
     * 平均绝对误差
     * @return MAE
     */
    private double MAE() {
        double resultMAE=0;
        int tempTestCount=0;
        for (int i = 0; i < numRating; i++) {
            int tempUserId=dataset[i].user;
            int tempItemId=dataset[i].item;
            double tempRating=dataset[i].rating;
            double tempPredict=predict(tempUserId,tempItemId);

            if(tempPredict<ratingLowerBound){
                tempPredict=ratingLowerBound;
            }else if(tempPredict>ratingUpperBound){
                tempPredict=ratingUpperBound;
            }
            double tempError=tempRating-tempPredict;
            resultMAE+=Math.abs(tempError);
            tempTestCount++;
        }
        return resultMAE/tempTestCount;
    }

    /**
     * 测试算法
     * @param args
     */
    public static void main(String[] args) {
        //计算开始时间
        long startTime = System.currentTimeMillis();
        try {
            // Step 1. 读取数据集
            MyMatrixFactorization mf=new MyMatrixFactorization(
                    "E:\\DataSet\\movielens-943u1682m.txt",
                    943, 1682, 10000,
                    1, 5);

            //step 2. 设置参数
            mf.setParameters(5, 0.0001, 0.005);

            // Step 3. 进行训练
            System.out.println("Begin Training ! ! !");
            mf.train(2000);

            //输出相关误差
            double tempMAE = mf.MAE();
            double tempRSME = mf.RSME();
            System.out.println("Finally, MAE = " + tempMAE + ", RSME = " + tempRSME);
            //计算运行时间
            long endTime = System.currentTimeMillis();
            System.out.println("运行时间："+(endTime-startTime));
        } catch (Exception e) {
            e.printStackTrace();
        } // Of try
    }
}

5.运行结果

运行2000轮，只用了662ms，java的速度还是很快的。