Python实现NMF（非负矩阵分解）

出现背景

• 数据分布不均，维度大（导致计算效率低下）
• 难以解释负值的实际意义（e.g.图像数据中不可能有负值的像素点）

分解方法

将一个非负矩阵分解为两个非负矩阵的乘积。

倍增更新规则（2001）：
在欧氏距离下，随机初始化，然后不断迭代，直到W（系数矩阵） H（基矩阵）稳定：

NMP求得的是局部最优解local optima（对应global optima问题），对于局部最优算法，一般考虑以下几点：

• 是否收敛
• 收敛速度如何
• 迭代的终止条件是什么（误差函数，最小化损失函数）

应用

• 压缩。将原本的矩阵由n*m转换为了（n+m）*r（与聚类的关系：新基是类中心，系数表示归属度）
• 实际应用：图像识别（机场、太空垃圾、识别星体、机器人识别环境），文本处理（通过NMF发现局部相关性），语音特征提取

实例

对于一个交通流量矩阵V^d，行表示路段，列表示该路段在96个时段内的交通流量。
即对于V^d_ij 表示第d天路段i在时段j内的交通流量。
其中，V^d的维数 = 路段数 * 时段数。

分解形式：
V_{m * h} ≈ C_{m * r} P_{r * h}

tips:

• m:路段数量
• h:时段数量
• r:交通模式数量（想分几类就把r定为几，对应着PCA中选取了前几个特征向量）

分解后得到基本模式矩阵P，系数矩阵C。
对于基本模式矩阵P_{r * h} ：

• 每一行对应一个基本的交通模式（流量随时间变化的模式）
• P_kj:第k种交通模式时段j内的流量

基本模式矩阵刻画了每一种模式的情况。

对于系数矩阵C_{m * r}:

• 每一行表示该路段在基本模式上的投影系数，C_ik:路段i在基本模式k上的投影系数
• 系数矩阵可以代表实际道路交通模式

系数矩阵的列对应三种模式，可以理解为一种道路有三种特征，这三种特征在这个道路上的占比。

import numpy as np
import pandas as pd

data = pd.DataFrame
data = pd.read_csv('20111116.csv',header = None)
data = data.iloc[:,3:]
data = data.values


def train(V, r, k, e):
    m, n = np.shape(V)
    W = np.mat(np.random.random((m, r)))
    H = np.mat(np.random.random((r, n)))
    data = []

    for x in range(k):
        V_pre = W * H
        E = V - V_pre  # 误差
        err = 0.0
        for i in range(m):
            for j in range(n):
                err += E[i, j] * E[i, j]
        data.append(err)
        print(err)
        if err < e:  # e 阈值
            break

        a = W.T * V  # Hkj
        b = W.T * W * H

        for i_1 in range(r):
            for j_1 in range(n):
                if b[i_1, j_1] != 0:
                    H[i_1, j_1] = H[i_1, j_1] * a[i_1, j_1] / b[i_1, j_1]

        c = V * H.T
        d = W * H * H.T
        for i_2 in range(m):
            for j_2 in range(r):
                if d[i_2, j_2] != 0:
                    W[i_2, j_2] = W[i_2, j_2] * c[i_2, j_2] / d[i_2, j_2]

    return W, H, data

W, H ,error= train(data, 4, 100, 1e-5 )

print(W)
print('\n')
print(H)