矩阵的模怎么求_推荐系统玩家 之 随机梯度下降（Stochastic gradient descent）求解矩阵分解...

前言

在上一篇的文章中，我们讲到了推荐系统中矩阵分解的三种方法。而这三种基本方法中，Funk-SVD由于其对稀疏数据的处理能力好以及空间复杂度低，是最合适推荐系统情景的，（Funk-SVD只是这三个基本方法里最好的，不代表就是推荐系统中最好的，还有更多衍生出来的优秀的方法，未来会给大家介绍）我们这篇文章就以Funk-SVD为基础，为大家介绍下如何求解矩阵分解时运用的梯度下降法以及其具体代码的实现。（梯度下降法的思想不仅运用于求解矩阵分解，同时也是神经网络中反向传播的基础，因此，希望大家多多留意。）

这是上一篇矩阵分解的基本方法和其优缺点，有兴趣的朋友可以看一看：

推荐系统玩家：推荐系统玩家 之 矩阵分解(Matrix Factorization)的基本方法 及其 优缺点​zhuanlan.zhihu.com

矩阵分解的认识

那么首先，我们来回顾一下矩阵分解的过程：

矩阵分解算法是将

维的矩阵

分解为

维的用户矩阵

和

的物品矩阵

相乘的形式。其中，

为用户的数量，

为物品的数量，

为隐向量的维度。基于分解后的两个矩阵

和

，通过物品行向量和对应矩阵行向量的内积，我们就可以得到任何用户

对物品

的预估评分，即预测值：

其中，

是用户

在矩阵

中对应的行向量，

是物品

在物品矩阵

中对应的行向量。从这里也可以很清楚看出矩阵分解在推荐系统中的意义：一个稀疏的矩阵通过分解成为两个低维的矩阵，再通过这两个矩阵的行向量之间的内积，就可以得到一个满秩的预测评分矩阵。

那么问题来了，我们如何将矩阵

分解成矩阵用户矩阵

和物品矩阵

的呢？

首先，从上面的式子我们可以知道，

是矩阵分解后，通过用户向量

和物品向量

的内积得到的预测评分。我们的目的是让预测评分

和真实评分

的差尽可能的小，即

最小。这样内积后的新矩阵才会更好拟合原矩阵，保存原始的信息的同时，对未知评分有准确的预测。

基于以上原理，我们采用均方差做为损失函数可以得到：

其中，

是所有用户评分样本的集合。因为为了减少拟合，我们在上式的后面加入正则化项。那么我们最终的

损失函数就变为：

在这里，我们需要说明下，损失函数的种类也有很多，Simon Funk 在这里运用均方根误差，也就是我们所熟悉的RMSE作为损失函数，同时这也就是说Funk-SVD引入了线性回归思想的体现。同时正则化的种类也不只以上这一种。因此，我们在后续的文章中还会专门介绍损失函数的种类以及什么是过拟合现象和正则化。

那么当我们已经有了损失函数，下一步就是要求损失函数的最小值，这就是我们常说的最小化损失函数。那么，如何最小化损失函数呢？

梯度下降法

至此，我们来到了本篇文章的重点。说实话，在学习这部分的时候，我始终处于懵的状态，因为无论是书籍还是知乎上，都用了大面积的数学公式来讲解导数，偏导数，多元函数求偏导等等，导致我很难理解，那么在这里我会用简单的思路来解释。但希望大家看后有不一样的理解。其中，导数，微分的定义与推导的基本知识我默认大家都有，想要复习的朋友可以看下这篇：

初中生能看懂的梯度下降算法​mp.weixin.qq.com

那么，我在这里主要会解决这三个问题：

• 什么是梯度？
• 为什么梯度反方向是函数值局部下降最快的方向？为什么沿着梯度的反方向更新参数？
• 梯度下降如何运用于推荐系统中求解矩阵分解

1.什么是梯度

首先明确几个概念，一元函数曲线上的切线斜率叫做导数， 也就是函数在该点的变化率。那么多元函数是个曲面，曲面上的切线可以有无数条。那么，我们找曲面上沿

轴方向的切线斜率和与沿着

轴方向的切线斜率就是我们常说的偏导数。因此

偏导数就是指的是多元函数 沿坐标轴的变化率。那么，我们可以看出，偏导数是沿着坐标轴的变化率，当我要想求沿任意方向的变化率怎么办呢？ 这时候 方向导数的概念就出现了。因为数学公式对很多朋友都不太友好，我在这里就不大面积列举公式了。我们只需要知道方向导数可以用来表示曲面上任意方向的变化率。那么既然可以求任意方向的变化率了，是不是就可以找到一个变化率最大的方向呢？答案是肯定，而这个函数变化率最大的方向就是 梯度的方向，这个点函数的最大变化率就是就是 梯度。

这里有一个三维的视频非常好帮助大家理解的导数，偏导数以及梯度的意义：

推荐系统玩家 之 导数，偏导数，方向导数的三维视图_哔哩哔哩 (゜-゜)つロ 干杯~-bilibili​www.bilibili.com

2.为什么梯度反方向是函数值局部下降最快的方向？（梯度下降的推导）

高中的时候，我们学过一个公式叫做泰勒公式，在这里，我们以一阶泰勒展开式为例：

如上图，函数

在

上用黑色曲线表示。红色的线为函数在

上的切线，那么利用极限近似的思想

就可以用上式近似。 其中，

为函数的导数，也就是切线的斜率。(因为导数在物理上有瞬时速度的意义，我们可以抽象的将上式理解为纵坐标路程就等于横坐标时间乘以斜率速度。）

理解了一阶泰勒展开式我们继续向下看，公式中的

是一个微小的矢量，它代表着变化的距离，它的大小我们可以用

来表示，

为标量，即一个确定的数。

的单位向量则用

来表示:

因此上式就变为：

我们做梯度下降的意义就在于最小化损失函数，也就是最小化

，同时，梯度下降是个迭代的过程。也就是说，我们希望每次更新

后，都能够让

的值变小。因此

要小于0，则有：

其中，

为标量且一般为正值，可以忽略，因此：

在这里，我们需要仔细研究下这个简单的不等式。首先，

和

均为向量，且

为导数，在这里我们叫做它梯度，它的方向是函数变化率最大的方向。

就是我们下一步要前进的方向。从公式来看以上两个向量相乘需要小于0。向量相乘的公式是包含夹角的：

从上图可以看出，当

也就是第三种情况，两个向量的方向完全相反，就可以让向量相乘的积小于0且是最小的。也就是说，当

和

的方向相反时，

且最小。

从这里我们就可以看出，为什么梯度反方向是函数值局部下降最快的方向。在上面我们提到过梯度方向是指函数变化率最大的方向，也就是斜率越来越大的方向，曲线会越来越陡。那么梯度的反方向，就是函数值局部下降最快的方向，斜率是越来越小的，斜率最终会为0。这也就是为什么我们高中求解二元一次方程组的最小值的时候，求导后需要导数等于0。

当我们知道了

的方向是

的反方向后，我们就可以得到：

因为

是单位向量，所以除以

的模。将上式带入到

中可得：

因为向量的模是个标量，因此：

其中，

就是当前参数，

就代表

步长，也就是我们所说的 学习率，

为梯度的方向，

减号就代表沿着梯度的负方向。

至此，我们就得到了梯度下降法中参数

更新的表达式

。

说了这么多，我相信还有很多朋友很懵。因为我在学习梯度下降初期，也不明白为什要理解什么表达式，为什么要知道这么多导数微分的知识。所以，我用一个例子，来给大家解释下我们上面那么长篇幅所学到内容的意义所在。

首先我们来看一个最简单的例子，我们目标函数为：

。我们想最小化这个函数，即求这个函数的最小值。很多朋友会说，求这个函数的解很简单啊，不就求导数，然后另导数等于0不就行了么。是的，从初高中的知识来看，我们确实可以这么做，因为这个目标函数是一个一元二次方程。如果当目标函数是个多变量函数,例如：

。还能那种方法来求么？显然是不行的。因此，我们就需要运用近似的方法，通过梯度下降的迭代方法逐渐逼近最优解。

那么我来看下以

为目标目标函数做梯度下降的过程：

首先我们设

的初始值为想

，学习率

。迭代的次数为10次，也就是说

更新10次。那么，根据梯度下降参数更新的公式：

的导数为

。因此，

第一次更新:

，得

，将更新后的

带入继续更新，以此类推。在这里我们放上python的代码：

%matplotlib inline
import d2lzh as d2l
import math
from mxnet import nd
import numpy as np

def gd(eta):
    x = 10
    results = [x]
    for i in range(10):
        x -= eta * 2 * x  # f(x) = x * x的导数为f'(x) = 2 * x
        results.append(x)
    print('epoch 10, x:', x)
    return results

res = gd(0.2)

绘制

迭代轨迹图可以看到：

def show_trace(res):
    n = max(abs(min(res)), abs(max(res)), 10)
    f_line = np.arange(-n, n, 0.1)
    d2l.set_figsize()
    d2l.plt.plot(f_line, [x * x for x in f_line])
    d2l.plt.plot(res, [x * x for x in res], '-o')
    d2l.plt.xlabel('x')
    d2l.plt.ylabel('f(x)')

show_trace(res)

在最后的值为：0.0604661759999999，也就是最后的结果。我们都知道，原方程的解应该为

，那么为什么最后我们的解跟它有误差呢？因为梯度下降法是通过迭代的方法无线逼近于最小值，其结果和学习率以及迭代次数都有着很大的关系。因此，通过调整参数，就可以优化我们想要的结果。

我们在里稍微讲一下学习率

在梯度下降中起到的作用。学习率

是一个超参数，在实际的实验中需要我们人为设定，合适的学习率可以很大程度上提高样本训练的精度。这个值，不能过大也不能过小，学习率太小可能会导致无法到达最优解，迭代就已经完成了。过大会让解来回震荡，无法走到最小值，我们来看个例子：

同样是目标函数

。当我们调整学习率

时，其他参数不变，最终的

如下图：

当我们调整学习率为

时：

的值会在曲线两边来回震荡，无法达到最优解。

因此可以看出，

学习率的值是个非常关键的参数，可能直接影响模型的好坏。我们在实际运用过程中，需要经验以及不断的测试来找到最优参数。当然，如今也有很多研究在寻找一个自动调参数的方法，我们也希望有更好的方法来帮助我们找到合适的参数。

3.随机梯度下降如何运用于推荐系统中求解矩阵分解

能看到这里，我相信大家对于梯度下降应该已经有了一个很直观的理解。梯度下降的方法其实还分为批量梯度下降法BGD，随机梯度下降SGD以及min-batch 小批量梯度下降法MBGD，那么，在这里，我们会讲解随机梯度下降是如何求解推荐系统中的矩阵分解的。

我们以上的概念和例子都是基于最简单的一维函数的。那么对于多维函数，又该怎么办呢？道理其实是一样，对一元函数，我们求其导数。那么对于多元函数，我们就分别求其偏导数。

那么梯度下降参数更新的公式就变为：

还记得我们做矩阵分解时得到的损失函数么?

这其中

,

为我们需要优化的变量，

为正则化稀疏，和

学习率一样是一个参数，标量。

为用户对物品真实的打分。因此，我们的目的就是优化

和

找到损失函数的最小值，从而拟合出一个用户矩阵和一个物品矩阵，他们内积后得到一个与原评分矩阵误差最小的预测评分矩阵。

那么，我们的对损失函数中的

,

分别求偏导可得：

将以上两个偏导分别带入梯度下降参数更新公式

得：

在以上公式中，

即为当前真实评分和预测评分的误差

,因此化简可得：

至此，随机梯度下降法的原理就解释完毕了。当然，梯度下降的方法还分为批量梯度下降法BGD，随机梯度下降SGD以及min-batch 小批量梯度下降法MBGD，我们就不在这里一一具体讲解，有兴趣的朋友可以看看这篇文章：

忆臻：详解梯度下降法的三种形式BGD、SGD以及MBGD​zhuanlan.zhihu.com

当我们有了参数更新的公式，就可以付诸于代码了。以下便为随机梯度下降求解矩阵分解的python代码和Matlab代码，同时数据集可以用Movielens 100k 开源数据集进行测试：

MovieLens 100K Dataset​grouplens.org

Python：

def MF(train_list, test_list, N, M, K=10, learning_rate=0.001, lamda_regularizer=0.1, max_iteration=50):
    # train_list: train data 
    # test_list: test data
    # N:the number of user
    # M:the number of item
    # learning_rate: the learning rata
    # K: the number of latent factor
    # lamda_regularizer: regularization parameters
    # max_iteration: the max iteration
    
    P = np.random.normal(0, 0.1, (N, K))
    Q = np.random.normal(0, 0.1, (M, K))
    
    train_mat = sequence2mat(sequence = train_list, N = N, M = M)
    test_mat = sequence2mat(sequence = test_list, N = N, M = M)
    
    records_list = []
    for step in range(max_iteration):
        los=0.0
        for data in train_list:
            u,i,r = data
            P[u],Q[i],ls = update(P[u], Q[i], r=r, learning_rate=learning_rate, lamda_regularizer=lamda_regularizer)
            los += ls
        mae, rmse, recall, precision = evaluation(P, Q, train_mat, test_mat)
        records_list.append(np.array([los, mae, rmse, recall, precision]))
        
        if step % 10 ==0:
            print(' step:%d n loss:%.4f,mae:%.4f,rmse:%.4f,recall:%.4f,precision:%.4f'
                  %(step,los,mae,rmse,recall,precision))
        
    print(' end. n loss:%.4f,mae:%.4f,rmse:%.4f,recall:%.4f,precision:%.4f'
          %(records_list[-1][0],records_list[-1][1],records_list[-1][2],records_list[-1][3],records_list[-1][4]))
    return P,Q,np.array(records_list)


def update(P, Q, r, learning_rate=0.001, lamda_regularizer=0.1):
    error = r - np.dot(P, Q.T)            
    P = P + learning_rate*(error*Q - lamda_regularizer*P)
    Q = Q + learning_rate*(error*P - lamda_regularizer*Q)
    loss = 0.5 * (error**2 + lamda_regularizer*(np.square(P).sum() + np.square(Q).sum()))
    return P,Q,loss

Matlab：

function [F, G] = mySVDBiasedV2(X,testSet,option)

k = option.k;
lRate = option.lRate;
lambda1 = option.lambda1;
iteLimit = option.iteLimit;

[rowNum colNum] = size(X);

XInd = find(X);

a = -1;
b = 1;

F = (b-a).*rand(rowNum, k) + a;
G = (b-a).* rand(colNum, k) + a;

rowSum = (b-a) .* rand(rowNum, 1) + a;
colSum = (b-a) .* rand(colNum, 1) + a;

M = F * G';

error = [Inf];
errorRMSE = [Inf];

trainErr = [Inf];
trainRMSE = [Inf];

iteNum = 0;

while true
    for rInd = 1:length(XInd)
        curError = X(XInd(rInd)) - M(XInd(rInd));
        
        [curR, curC] = ind2sub([rowNum colNum], XInd(rInd));
        
        for kInd = 1:k
            F(curR, kInd) = F(curR, kInd) + lRate * curError * G(curC, kInd) - lambda1 * F(curR, kInd);
            G(curC, kInd) = G(curC, kInd) + lRate * curError * F(curR, kInd) - lambda1 * G(curC, kInd); 
        end
    end
    
    M = F*G';
    
    p = [];
    
    trainRMSE = [trainErr sqrt(sum((M(XInd) - X(XInd)).^2)/length(XInd))];
    trainErr = [trainErr sum(abs(M(XInd) - X(XInd)))/length(XInd)];
    
    if (~isempty(testSet))
        for testInd = 1:size(testSet, 1)
            p = [p M(testSet(testInd, 1), testSet(testInd, 2))];
        end
        errorRMSE = [error sqrt(sum((p - testSet(:, 3)').^2)/size(testSet, 1))];
        error = [error sum(abs(p - testSet(:, 3)'))/size(testSet, 1)];
    end
    
    iteNum = iteNum + 1;
    
    fprintf('iteNum: %d, (trainMAE = %f, testMAE = %f), (trainRMSE = %f, testRMSE = %f)n', iteNum, trainErr(end), error(end), trainRMSE(end), errorRMSE(end));
    if trainErr(end) > trainErr(end-1) || iteNum > iteLimit
        plot(trainErr, '-ob');
        if (~isempty(testSet))
            hold on
            plot(error, '-sr');
            hold off
        end
        return;
    end
end

上一篇：

推荐系统玩家：推荐系统玩家 之 矩阵分解(Matrix Factorization)的基本方法 及其 优缺点​zhuanlan.zhihu.com

下一篇：

推荐系统玩家：推荐系统玩家 之 推荐系统入门——推荐系统的发展历程（上）​zhuanlan.zhihu.com

参考文献

1.https://zh.d2l.ai/chapter_optimization/gd-sgd.html

2.为什么局部下降最快的方向就是梯度的负方向？_红色石头的专栏-CSDN博客_负梯度方向

3.https://zhuanlan.zhihu.com/p/30759733

4.https://mp.weixin.qq.com/s/jrMs8U_Mg82lPxZB3h4SbQ

5.https://zh.d2l.ai/chapter_optimization/gd-sgd.html

6.https://zhuanlan.zhihu.com/p/33321183

7.https://zhuanlan.zhihu.com/p/25765735