回顾MEBF思路
矩阵分解
N
N
N:对象数,
M
M
M:属性数,
L
L
L:标签数
X
:
N
∗
M
X: N*M
X:N∗M
Y
:
N
∗
L
Y: N*L
Y:N∗L
线性模型:
W
:
M
∗
L
W: M*L
W:M∗L
希望:
Y
=
X
∗
W
Y=X*W
Y=X∗W
问题:
- 线性模型假设每个标签之间是没有关系的. 写一个算例:
N
=
4
,
M
=
3
,
L
=
2.
N = 4, M = 3, L = 2.
N=4,M=3,L=2.
线性模型可以写成
y = f ( x ) = W x y = f(x) = Wx y=f(x)=Wx
x x x是 1 ∗ M 1 * M 1∗M - 一般的神经网络模型也不能很好地体现标签相关性. 但总是非线性模型.
网络模型写成: y = f ( x ) y = f(x) y=f(x)
f f f是一个非线性的函数
解决方案:
使用布尔矩阵分解, 把标签相关性提取出来.
U
:
M
∗
K
U: M * K
U:M∗K
V
:
K
∗
L
V: K * L
V:K∗L
K
K
K: latent factors/vector length, 隐含标签 (属性) 数
Y = U ∗ V Y = U*V Y=U∗V 布尔矩阵分解, 不一定是完全的拟合, K K K 可以自己确定, 当然要小于 L L L
f ( X ) = U f(X) = U f(X)=U
训练神经网络,只需要拟合 K < L K < L K<L 个标签. 预测的时候, 将这 K K K 个标签用 V V V 映射到 L L L 个.
先不管 missing label, 进行一般的布尔矩阵分解.
布尔矩阵分解的实例及代码实现(MEBF)
首先输入是Original矩阵X,然后t是用于是该列(行)是否覆盖的一个标准,还有一个字母是收敛的一个判定(这个在代码中比较好理解)
输出自然不用说,两个布尔矩阵。
进入循环,当不收敛时一直继续,
首先是使用Bidirectional_growth算法来获取subMatrix的分解结果列向量a和行向量b(结合上面讲的MEBF的思路,应该就知道这里a是什么,b是什么了)。(Bidirectional_growth就是Median expansion)
然后下面就是算一下这里的提取出来的这个patterns是否能让矩阵的误差减小,如果减小的话,便更新,如果不能减小的话,则说明,目前矩阵的信号很小了,需要用weak signal detection来处理。
之后对矩阵进行一个更新。
这便是MEBF整体的一个逻辑。
import numpy as np
def MEBF(Thres,MAT,DIM=1000,COVER=0.99):
'''
:param Thres: 论文中的t,作为衡量是否能够覆盖的阈值 高了则会覆盖过少,低了则会覆盖过多
:param MAT: 输入的矩阵
:param DIM: patterns的个数 也是希望得到子矩阵的个数
:param COVER: 作为收敛的判断条件
:return:
'''
if min(np.shape(MAT))<DIM: #目标patterns的个数一定要小于矩阵的维度,这个显然
DIM=min(np.shape(MAT))
M1 = MAT #M1作为residual matrix
SUM = np.sum(MAT) #计算一下原矩阵的1的个数
MAT_B = np.empty([np.shape(MAT)[0],0]) #论文中的A*
MAT_C = np.empty([0, np.shape(MAT)[1]])#给定shape
while np.sum(M1)>(1-COVER)*SUM and min(MAT_B.shape)<DIM:#COVER值只在这里出现
#循环条件,M1的和是原矩阵的一个1-cover的倍数,这个cover感觉是一个收敛的程度
#另一个循环条件是当MAT_B的k达到了目标的维数
#这里的循环条件意为:因为是布尔矩阵,所以其1的个数其实就是算法进行的一个衡量标准
#当1的个数只有原来的0.01以下时,说明表示的差不多了
e1= np.sum(M1) # 我改成e1了 来作为一个对比 M1的值
B1 = np.zeros(np.shape(M1)[0]) #对于MAT_B的列向量 初始化 全0
B1_use = B1
B2 = np.zeros(np.shape(M1)[1]) #对于MAT_C的行向量 初始化 全0
B2_use = B2
COL = np.sum(M1, axis=0) #列和 是一个向量
ROW = np.sum(M1, axis=1) #行和 也是一个向量
### start with column
if np.median(COL[COL>0])>1: #这里是非0值们的中位 COL[COL>0] 是选取其中大于0的元素作为一个切片 这里相当于是将全0的列给剔除了
TEMP=np.argwhere(COL==min(COL[COL>=np.median(COL[COL>0])])) #返回满足条件的索引 是个二维数组,第二维对于向量是一个元素
#条件为:col= col中大于均值的的最小值(略大于均值)
#意为:找到最接近中位数的数的位置
if np.shape(TEMP)[0]==1: #如果这个条件的只有1个数 这很正常
B1=M1[:,TEMP[0,0]] #找到那一列 ,找到列和接近均值的那一列 也就是median column 列向量 B1就是pattrens
#这里就是求各列向量与median column的相似性
#B1==1是直接找到在B1为1的位置上,1的sum大于阈值的列们
#然后B2就是论文中的e,也就是patterns的覆盖的情况
#expansion
B2[np.sum(M1[B1==1,:],axis=0)>=min(Thres*sum(B1)+1,sum(B1))]=1 #这里B1==1的结果是[false false true false true] 即找到为1的那些索引
#所以M1[B1==1,:]是找到其中为1的那些行 ,然后axis=0,又是得到列和向量
#然后sum(B1)就是B1中有多少1,然后min这个应该是Thes有关,大概率是第一个,第二个的话,太少了,最少只有B1那一个类是的
#找到大于这个值的列们,然后对应位置1
e2 = error(M1,B1,B2)
if e1>e2: #如果B1,B2的1的数量更多,那么就更新B1_use和B2_use
B1_use=B1
B2_use=B2
e1=e2 #这个标准是逐渐提高的
B1 = np.zeros(np.shape(M1)[0])
B2 = np.zeros(np.shape(M1)[1]) #这里相当于初始化,来重新计算行的
else: #如果最接近中位数的数不值一个
for j in range(np.shape(TEMP)[0]): #对每列进行遍历
B1=M1[:,TEMP[j,0]] #TEMP[j,0]是第j个数的行
B2[np.sum(M1[B1 == 1, :], axis=0) >= min(Thres * sum(B1) + 1, sum(B1))] = 1 #这个和上面一样
e2 = error(M1,B1,B2)
if e1 > e2:
B1_use = B1
B2_use = B2
e1 = e2
B1 = np.zeros(np.shape(M1)[0])
B2 = np.zeros(np.shape(M1)[1])
#如此得到的是B1_use和B2_use ,然后C1也更新了
### start with row 这里的逻辑应该一样,和上面差不多 只是这里用的是ROW
if np.median(ROW[ROW>0])>1:
TEMP = np.argwhere(ROW == min(ROW[ROW >= np.median(ROW[ROW > 0])]))
if np.shape(TEMP)[0] == 1:
B2 = M1[TEMP[0,0],:]
B1[np.sum(M1[:,B2 == 1], axis=1) >= min(Thres * sum(B2) + 1, sum(B2))] = 1
e2 = error(M1,B1,B2)
if e1 > e2:
B1_use = B1
B2_use = B2
e1 = e2
B1 = np.zeros(np.shape(M1)[0])
B2 = np.zeros(np.shape(M1)[1])
else:
for j in range(np.shape(TEMP)[0]):
B2 = M1[TEMP[j, 0], :]
B1[np.sum(M1[:, B2 == 1], axis=1) >= min(Thres * sum(B2) + 1, sum(B2))] = 1
e2 = error(M1,B1,B2)
if e1 > e2:
B1_use = B1
B2_use = B2
e1 = e2
B1 = np.zeros(np.shape(M1)[0])
B2 = np.zeros(np.shape(M1)[1])
#weak signal detection algorithm
if(True):#我让Bidirectional growth算法和weak signal detection 并行走了,直接看谁的效果好就用谁
#实际上,局部最优不一定是全局最优,这样处理实际不是特别合适的其实
COL_order=np.argsort(COL)[::-1] #返回 按照sum和最大来排序的索引
ROW_order=np.argsort(ROW)[::-1]
B1 = np.zeros(np.shape(M1)[0])
B2 = np.zeros(np.shape(M1)[1])
### start from COL
B1[(M1[:,COL_order[0]]+M1[:,COL_order[1]]==2)]=1 #把1最稠密的两列给叠加起来,然后来构成新的pattern
B2[np.sum(M1[B1 == 1, :], axis=0) >= min(Thres * sum(B1) + 1, sum(B1))] = 1 #expansion
e2 = error(M1,B1,B2)
if e1>e2:
B1_use = B1
B2_use = B2
e1 = e2
### start from ROW
B1 = np.zeros(np.shape(M1)[0])
B2 = np.zeros(np.shape(M1)[1])
B2[(M1[ROW_order[0],] + M1[ROW_order[1],] == 2)] = 1
B1[np.sum(M1[:, B2 == 1], axis=1) >= min(Thres * sum(B2) + 1, sum(B2))] = 1
e2 = error(M1,B1,B2)
if e1> e2:
B1_use = B1
B2_use = B2
e1 = e2
# 如果处理之后还等于0, 则直接brek
if e1==np.sum(M1):
break
#确实更新了,所以这里进行修正矩阵
else:
M1 = M1 - np.outer(B1_use,B2_use) #减去subMatrix
M1[M1<0] = 0 #减法可能会造成-1,那么置零
B1_use=B1_use.reshape(B1_use.shape[0],1)
B2_use=B2_use.reshape(1,B2_use.shape[0])
MAT_B = np.concatenate((MAT_B, B1_use), axis=1)#矩阵的拼接 所以前面是empty(0)就可以接受了
MAT_C = np.concatenate((MAT_C, B2_use), axis=0)
print("本轮误差是:"+str(e2))
return MAT_B,MAT_C
def error(residuMatrix, columnVector,rowVector):
'''
计算误差,原残差矩阵-向量之积,然后计算矩阵之和(非绝对值)
:param residuMatrix: 残差矩阵,不能修改
:param columnVector: 列向量,patterns
:param rowVector: 行向量 覆盖情况
:return: int型的误差值
'''
e_matrix = residuMatrix - np.outer(columnVector,rowVector)
e_matrix[e_matrix<0] =0 #将负数置0,因为我上面的减法就是这样操作的,需要和那个一样,不然就没有意义了
e = np.sum(e_matrix)
return e
#论文中实例的a1例子
matrix = np.array([[1,1,0,1,1,0,0,0,0,0,0],[0,1,1,1,1,0,0,0,0,0,0],[1,1,0,1,0,1,1,0,0,0,0],[1,1,1,1,1,1,1,1,1,1,0],[0,0,0,1,0,0,0,0,0,0,0],[1,1,1,1,1,1,1,0,1,0,0,],[1,1,1,1,1,1,1,1,1,0,0,],[0,1,0,0,0,0,0,0,0,0,0],[1,1,0,1,1,1,1,0,1,0,1],[1,1,1,1,1,1,1,1,1,1,1],[1,1,0,1,1,1,0,1,0,0,0]])
print("原本矩阵的能量:"+str(np.sum(matrix)))
B,C = MEBF.MEBF(0.85,matrix,DIM=6)
reconstruction_matrix = np.dot(B,C)
reconstruction_matrix[reconstruction_matrix>1]=1
error_matrix = matrix-reconstruction_matrix
error_matrix[error_matrix!=0]=1
print("最终误差是:"+str(np.sum(error_matrix)))
print("done")