SVD奇异值分解

http://www.cnblogs.com/jerrylead/archive/2011/04/18/2020209.html  今天看了这篇博客博主的一些文章，感觉很有收获。突然想到自己也写过svd的东西，就来看看自己的，才发现被之前不小心删除了。理论方面的东西就不多说了，这里会贴上两个不同实现方法的一段源码，然后分析一下异同。

svd的源码1：

      这一份源码是根据《recommender system handbook》里面提到的隐语意模型做的。代码非常短，只实现了从矩阵中分解出一个主要成分的功能，为了能够方便对比，一开始的矩阵也是有两个向量相乘构建的。

__author__ = 'axuanwu'
# coding=utf-8

import numpy as np
# 构建一个只包含一个主成分的矩阵
aaa = np.zeros((4, 5))
(dim0, dim1) = aaa.shape
x0 = [1, 2, 3, 4, 5, 6, 7]
x1 = [7, 6, 5, 4, 3, 2, 1]
for i in xrange(0, dim0):
    for j in xrange(0, dim1):
        aaa[i, j] = x0[i]*x1[j]
# 以下为矩阵分解部分
lambda_1 = 0.02  # 经验参数 用于调整步长
y0 = np.matrix([1.0]*dim0)
y1 = np.matrix([1.0]*dim1)
first = True
i_step = 0
while i_step < 20:  # 经验参数  循环次数
    i_step += 1
    #  调整y0
    for i in xrange(0, dim0):
        error_sum = float((aaa[i, :] - y0[0,i] * y1).sum(1))
        y0[0, i] += error_sum * lambda_1
    #  调整y1
    for i in xrange(0, dim1):
        error_sum = (aaa[:, i] - y1[0, i] * y0).sum(1)[0,0]
        y1[0, i] += error_sum * lambda_1
print np.dot(y0.T, y1)  # 分解之后的两个向量的叉乘结果
print aaa  # 原始的矩阵

        如果看上面的源码，就烦请问自己几个问题：目标函数是什么？梯度下降的步长受什么影响？函数的目标是所凑的矩阵的每一列以及行的求和与目标矩阵每一列以及行相同，其实这里的目标与数学上的svd的目标是不同的，但是，管他呢，能工作而且效率高。这里的梯度下降虽然步长可变，前期表现很好，但是完全不受控制。下面给出另外一个方法的源码，在目标函数和梯度下降方法上都不同，其中梯度下降的步长是控制的。

svd的源码2：

__author__ = 'axuanwu'
# coding=utf-8
import numpy as np

class Svd_0():
    #  把 svd 的 vd 乘积 作为 d；s 矩阵 任然是 单位矩阵
    def __init__(self):
        self.aaa = np.zeros((4, 5))
        (dim0, dim1) = self.aaa.shape
        self.num_dim = 2 #分解出来的主要 因子的数量
        x0 = [1, 2, 3, 4, 5, 6, 7]
        x1 = [7, 6, 5, 4, 3, 2, 1]
        for i in xrange(0, dim0):
            for j in xrange(0, dim1):
                self.aaa[i, j] = float(x0[i]*x1[j] + x0[j]*x1[j])
        self.remain = self.aaa - 0
        self.s_matrix = np.ones((dim0,self.num_dim))
        self.d_matrix = np.ones((self.num_dim,dim1))
        self.pre_error0 = np.array([0.0] * dim0)
        self.pre_error1 = np.array([0.0] * dim1)
        self.sign0 = np.array([1] * dim0) #  前进方向
        self.sign1 = np.array([1] * dim1) #  前进方向
        self.step_my = 0.0  # 当前处理步长
        self.step_start = 0.5  # 起始步长 设置为 0.5
        self.step_end = 10 ** -4  # 终止精度

    def set_pre_error(self, chenfen_i):
        (dim0, dim1) = self.aaa.shape
        # dim 0
        for i in xrange(0, dim0):
            self.pre_error0[i] = sum(np.power(np.dot(self.s_matrix[i, chenfen_i], self.d_matrix[chenfen_i, :]) - self.remain[i, :], 2)) 
        for j in xrange(0, dim1):
            self.pre_error1[j] = sum(np.power(np.dot(self.s_matrix[:, chenfen_i], self.d_matrix[chenfen_i, j]) - self.remain[:, j], 2)) 

    def accept_change(self, i_index, chenfen_i, axi_num):
        #  chenfen_i 主成分编号
        #  axi_num 维度编号
        #  元素下标 i_index
        (dim0, dim1) = self.aaa.shape
        if axi_num == 0:
            temp_i = self.s_matrix[i_index, chenfen_i] + self.sign0[i_index] * self.step_my
            temp_e = sum(np.power(np.dot(temp_i, self.d_matrix[chenfen_i, :]) - self.remain[i_index, :], 2))
            if temp_e < self.pre_error0[i_index]:
                self.s_matrix[i_index, chenfen_i] = temp_i
                self.pre_error0[i_index] = temp_e
                return True
            else:
                self.sign0[i_index] *= -1
                return False
        else:
            temp_i = self.d_matrix[chenfen_i, i_index] + self.sign1[i_index] * self.step_my
            temp_e = sum(np.power(np.dot(self.s_matrix[:, chenfen_i], temp_i) - self.remain[:, i_index], 2))
            if temp_e < self.pre_error1[i_index]:
                self.d_matrix[chenfen_i, i_index] = temp_i
                self.pre_error1[i_index] = temp_e
                return True
            else:
                self.sign1[i_index] *= -1
                return False

    def svd_fun(self):
        (dim0, dim1) = self.aaa.shape
        for i_chenfen in xrange(0, self.num_dim):
            self.step_my = self.step_start
            mark_pre = True  # 刷新 标志
            self.set_pre_error(i_chenfen)
            while self.step_my > self.step_end:
                mark = False  # 刷新 标志
                for i in xrange(0, dim0):
                    mark |= self.accept_change(i, i_chenfen, 0)
                for i in xrange(0, dim1):
                    mark |= self.accept_change(i,i_chenfen, 1)
                if not (mark | mark_pre):  # 连续2次没刷新，达到当前最优解
                    self.step_my *= 0.5
                    mark_pre = True
                else:
                    mark_pre = mark
            self.remain = self.aaa - np.dot(self.s_matrix[:,0:(i_chenfen+1)], self.d_matrix[0:(i_chenfen+1), :])
        #  归一化
        for i_chenfen in xrange(0, self.num_dim):
            normal = float(sum(self.s_matrix[:, i_chenfen] * self.s_matrix[:, i_chenfen])) ** 0.5  # 归一化常数
            if normal == 0:
                continue
            self.s_matrix[:, i_chenfen] = self.s_matrix[:, i_chenfen] / normal
            self.d_matrix[i_chenfen, :] = self.d_matrix[i_chenfen, :] * normal

if __name__ == "__main__":
    a = Svd_0()
    a.svd_fun()
    print np.dot(a.s_matrix, a.d_matrix)
    print "----------------------------------------"
    print a.aaa
    print a.s_matrix
    print a.d_matrix

        这一段源码里面的优化目标是 构建出来的矩阵与目标矩阵每个元素差的平方和最小，这里和线性拟合的最小二乘类似，但是与数学上的svd的目标有出入。在这里面首先分解出一个主成分，让后将剩下的部分取出再分解，凑出另外一个主成分。不过这里的梯度下降算法却是比较有意思的，大家可以仔细看看。其中的思想是先找到大步长的最优解后，再将步长缩短，寻找小步长的最优解，是不是有点类似于动态规划的先找父问题的最优解在从这个最优解上寻找子问题的最优解。虽然这个方法看似很巧妙，但是前期收敛速度没有第一个代码快，不过最后的拟合效果要好一些。这一原因是当凑出来的矩阵和目标矩阵相近时，第一个算法的步长变得太小，几乎不能修正结果，不过这时也基本到了人可以接受的误差范围了。