SVD 梯度下降

花絮:

         非常抱歉,之前我没有严格的分清梯度下降和随机行走的区别,导致前一篇SVD的博客部分说法上有点小问题。之前说要实现svd的稀疏矩阵做法,这里算是实现了,总体来说,只要调整一下块的大小,该源码距离实际应用来做大矩阵的分解应该已经很接近了。

SVD的python源代码:

近期看了一篇小文 ,很有收获,故而实现之。由于该文中的中间有详细的说明算法流程和原理,这儿不介绍。大家仔细看看其中的算法说明部分,尤其是那几个公式。和原作相比这里有几个改变值得一提:

  1. 稀疏块
  2. 步长控制
  3. 块操作
  4. 循环终止控制
这里先上源码。
# coding=utf-8
import numpy as np
import time
import math

__author__ = 'axuanwu'
# 2015年 9 月 25 日
class XMatrix():
    def __init__(self, m=1000, n=100, step=4):
        self.zero_like = 0.01  # 伪零: 差异小于该值认为无差异
        self.step = int(step)  # 数据块的切分方法
        self.m = m  # 原矩阵行数
        self.n = n  # 原矩阵列数
        self.splits = int(math.ceil(1.0 * self.n / self.step))  # 每行记录的分块数
        self.res = self.n % self.step
        self.Memory = []  # 存储数据块的实体
        self.dict_rowid = {}  # 记录数据块的位置字典
        self.Memory_max_rowid = -1
        self.UUU = np.random.random((self.m, 2))
        self.VVV = np.ones((2, 1))  # 本来用做归一化UUU 和 MMM后的 系数矩阵,目前未处理
        self.MMM = np.random.random((2, self.n))

    def set_data(self, tezhenshu):
        self.tezhenshu = tezhenshu
        self.UUU = np.random.random((self.m, tezhenshu))
        self.VVV = np.zeros((tezhenshu, tezhenshu))
        for i in xrange(0, tezhenshu):
            self.VVV[i, i] = 1
        self.MMM = np.random.random((tezhenshu, self.n))

    def intoMatrix(self, i, j, data):
        # 矩阵赋值
        row_id = (int(j / self.step) + i * self.splits)
        i_temp = self.dict_rowid.get(row_id, -1)
        if i_temp == -1:
            self.Memory_max_rowid += 1
            i_temp = self.Memory_max_rowid
            self.dict_rowid[row_id] = i_temp
            self.Memory.append(np.array([0.0] * self.step))  # 增加一块 数据块
        self.Memory[i_temp][j % self.step] = data

    def getitem(self, i, j=-1):
        # 读取稀疏矩阵
        if j == -1:
            # 返回一数据块
            i_temp = self.dict_rowid.get(i, -1)
            if i_temp == -1:
                temp = np.array([0.0] * self.step)
            else:
                temp = self.Memory[i_temp]
            if self.splits == (i % self.splits) + 1:
                return temp[0:self.res]
            else:
                return temp
        else:
            # 返回元素
            i_temp = self.dict_rowid.get((int(j / self.step) + i * self.splits), -1)
            if i_temp == -1:
                return 0
            return self.Memory[i_temp][j % self.step]

    def error_sum(self, k=0):
        error_sum = 0
        if k != 0:  # 计算合成矩阵 与 与矩阵 的 最大差异
            for i in xrange(0, self.m):
                for j in xrange(0, self.splits):
                    start = j * self.step
                    end = min(start + self.step, self.n)
                    temp = np.dot(self.UUU[i, :], self.MMM[:, start:end])
                    i_temp = self.dict_rowid.get((j + i * self.splits), -1)
                    if i_temp == -1:
                        error_sum = max(error_sum, max(np.abs(temp)))
                    else:
                        error_sum = max(error_sum, max(np.abs(temp - self.Memory[i_temp][0:(end - start)])))
            return error_sum
        else:   # 计算合成矩阵 与 与矩阵 所有元素差异的平方和
            for i in xrange(0, self.m):
                for j in xrange(0, self.splits):
                    start = j * self.step
                    end = min(start + self.step, self.n)
                    temp = np.dot(self.UUU[i, :], self.MMM[:, start:end])
                    i_temp = self.dict_rowid.get((j + i * self.splits), -1)
                    if i_temp == -1:
                        error_sum += sum(np.power(temp, 2))
                    else:
                        error_sum += sum(np.power(temp - self.Memory[i_temp][0:(end - start)], 2))
            # print (np.power(self.test-np.dot(self.UUU[:, :], self.MMM[:, :]),2)).sum(0).sum(0)-error_sum
            return error_sum

    def SVD(self):
        step = 1
        iter = 0
        t_start = time.time()
        pre_error = self.error_sum()
        print "-- ", iter, " --", time.time() - t_start, pre_error, self.error_sum(1)
        while (1):
            du = np.ones((self.m, self.tezhenshu))
            dm = np.zeros((self.tezhenshu, self.n))
            # UUU 阵
            for i in xrange(0, self.m):
                t = np.zeros((self.tezhenshu, self.n))
                # t0 = np.zeros((self.tezhenshu, self.n))
                for j in xrange(0, self.splits):  # 按块遍历
                    js = j * self.step
                    je = min(js + self.step, self.n)
                    t[:, js:je] = (self.getitem(i * self.splits + j) - np.dot(self.UUU[i, :],
                                                                              self.MMM[:, js:je])) * self.MMM[:, js:je]
                du[i, :] = t.sum(1)  # - ku * self.UUU[i,:]
            # MMM 阵
            for j in xrange(0, self.splits):  # 按块遍历
                js = j * self.step
                je = min(js + self.step, self.n)
                t = np.zeros((self.tezhenshu, je - js))
                for i in xrange(0, self.m):
                    t[:, 0: je - js] += np.mat(self.UUU[i, :]).T * np.mat(
                        self.getitem(i * self.splits + j) - np.dot(self.UUU[i, :], self.MMM[:, js:je]))
                dm[:, js:je] = t  # - ku * self.MMM[:,js:je]
            u = step / max((np.abs(dm)).max(), (np.abs(du)).max())
            self.MMM += u * dm
            self.UUU += u * du
            iter += 1
            error_sum = self.error_sum()
            if error_sum > pre_error:  # 到达当前极值
                step *= 0.3  # 步长减少
                print step, u, iter, " --", time.time() - t_start, error_sum, self.error_sum(1)
            pre_error = error_sum
            # print "-- ", iter, " --", time.localtime(), error_sum,self.error_sum(1)
            if (self.error_sum(1) < self.zero_like) | (step < 10 ** (-8)):
                print "-- ", iter, " --", time.time() - t_start, error_sum, self.error_sum(1)
                break


if __name__ == "__main__":
    m = 10  # 需要分解的矩阵行数
    n = 8  # 需要分解的矩阵列数
    a = np.zeros((m, n))
    aaa = XMatrix(m, n, 6)
    aaa.set_data(2)  # 3 个主要成分
    for i in xrange(0, m):
        for j in xrange(0, n):
            a[i, j] = i + j % 4 + (i % 5) * (5 - j % 4)
            if a[i,j] == 0:
                continue
            aaa.intoMatrix(i, j, a[i, j])  # 将需要分解的矩阵的元素放入对象 aaa 中
    aaa.SVD()
    print aaa.error_sum(1)


稀疏块:

          源码中矩阵的存储形式并不是用一个m行n列的矩阵来存储数据,而是将每一行元素按照 self.step 个元素做为一块,存储在self.Memory列表中,如果整个块取值都为0,那么这个数据块将不被存储。因此这个方法与稀疏矩阵相比:块的个数少于元素的个数,因此需要记录位置减少;而与原本的实矩阵相比,部分零元素是不用记录的。这种方法介于稀疏矩阵与实矩阵之间。以商品商品推荐为例,如果按照热度排行来排列商品序号,稀疏块同样能够极大地减少记录矩阵需要的存储空间。

块操作:

          源码中在计算du,dm以及误差时均直接以数据块为单位进行矩阵运算,这样for循环的次数可以大约减少x倍(x为稀疏块的元素个数),而使用矩阵操作代替循环。就像R中的apply函数可以大幅增加计算性能一样,块操作对性能有很大提升。


步长控制:

         很多梯度下降算法都喜欢用经验参数,换个数据集就可能不适应,各种无穷大,精度损失,这还真是有点坑。而本源码中步长调节参数u通过step步长间接控制,在找到当前step的最小值后step会减小,从而进行更精细的调整。能够确保误差稳定快速的收敛。


循环控制:

         循环终止的条件有两个,1)精度达到要求,亦即任意元素与原矩阵对应元素的差的绝对值小于self.zero_like,2)矩阵稳定不再改变,亦即步长极小(亿分之一)。程序不会浪费太多时间做无用功。


目前的测试下,程序基本实现矩阵分解的功能,并且稳定不跳出异常。不过分解之后的矩阵各向量不正交。看了原文的同学可能华发现,与原来的相比,du,dm中我去掉了正则化项的影响,由于正则化的目的是减少过拟合,作为矩阵分解而言,在这里并不是有必要的。程序仍有两个不足:稀疏块是一个条形块,而不是一个方块,概率上来说,元素个数相同时方块更容易出现全0的块;目前并没有充分使用并行化技术来优化(块操作应该也属于并行运算的范畴)。

外话:

其实分解之前,我们不知道最少需要几个维度才能分解完矩阵,并且到精度要求。因此我也改了个一个维度一个维度分解的,由于未完成块操作部分的编码,里面仍然使用的是按元素操作。故而没有放这一个的源码。该方法的逻辑是将第一个维度分解到稳定后,第一个维度不变,分解下一个维度,每当前面的维度稳定就不再改变,如果精度未达到要求就增加一个维度,直到精度达标。 可惜的是在测试中,仍然不能满足正交。那么问题就来了:怎样通过梯度下降算法等价的实现奇异值分解?是最优函数的问题?还是……


©️2020 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客 返回首页
实付0元
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值