布尔矩阵分解 代码实现（BMF）--MEBF论文阅读

关于布尔矩阵分解

最近在做有关布尔矩阵分解方面的研究，因为自己的方向需要，找了有关布尔矩阵分解的一干论文看了看。
关于布尔矩阵分解，实际就是对于一个布尔矩阵$A_{m\ast n}$（其中元素非0即1）来分解成$B_{m\ast k}$与$C_{k\ast n}$，其中B与C也同样是布尔矩阵（其中元素非0即1），也就是将一个大型的布尔矩阵给分解为两个低秩矩阵的乘积。
其英文为Boolean Matrix Factorization，或者Boolean Matrix decomposition也行，前者用的更多更普遍。其有着很重要的应用，对于一些数据分析方面有着很大的帮助。
不过布尔矩阵分解本身就是一个np难的问题，它不像普通矩阵分解，可以用浮点数，其中方法比如SVD，能够很快的得到可以说是非常准确的分解结果。布尔矩阵分解成布尔矩阵的乘积，势必会有误差，难以准确的分解，同时，很多算法运行起来效率也是非常的低。
可以说，现在应该是没有如同SVD在非布尔矩阵分解方向如此高效准确的布尔矩阵分解方法的，布尔矩阵分解方向的方法可能说只是尽量去提高分解速度以及降低分解误差，根本就达不到说有一个公认的完美的方法。同时，对于稀疏程度不同，规模不同的布尔矩阵的分解，不同方法又会显示出不同的效果，所以该方向可以说还在努力的发展中。
关于布尔矩阵的分解方法，现在提出的有不少，比如ASSO，PANDA以及其plus版本，还有NMF等，本文中介绍的是2019年普渡大学出的一篇文章《Fast and Efficient Boolean Matrix Factorization by Geometric Segmentation》中提出的MEBF方法。

布尔矩阵分解的思路

下面就直接来介绍该方法的思想及算法流程。
首先所谓布尔矩阵分解，实际上就可以等价理解为去在原布尔矩阵上去确定有哪些1稠密的子矩阵。
我们举下面一个例子：

对于这样一个布尔矩阵，我们找它的1稠密子矩阵，很明显，这个矩阵比较规律也比较简单，我们可以进行一个近似，将它看做是两个矩形的叠加。

如上图所示，是红矩形与黑矩形的叠加，我们可以将该矩阵给分解成两个子矩阵的和。

我们将这两个矩阵分别分解为列向量与行向量的乘积（这个很好做到，因为矩形嘛，很均匀），所以便得到了原布尔矩阵的一个分解结果。
结果如下图所示。

如此便是布尔矩阵分解的思路，当然这里的例子十分的简单，所以可以这样简单的理解，实际上，并不可能这样的做。

MEBF思路

介绍完了布尔矩阵分解的大致思路，下面来具体讲讲论文中提出的MEBF。
MEBF的全称是Median Expansion for Boolean Factorization，姑且翻译为中位数扩展balabala，等你看玩整个算法，基本就理解了中位数扩展的含义。

其大致基于这样一些引理：

• 首先我们的问题可以刻画为这样一个表达式： a r g m i n A ∈ { 0 , 1 } n ∗ k , B ∈ { 0 , 1 } k ∗ m ∣ X − A ∗ B ∣ arg min_{A\in\{0,1\}^{n*k},B\in \{0,1\}^{k*m}}|X-A*B| argminA∈{0,1}n∗k,B∈{0,1}k∗m​∣X−A∗B∣这样一个问题，其中A有k列，而B有k行，A中的列向量与B中相对于的行向量的乘积便是一个子矩阵$X_{I_j}$，这k个子矩阵累加起来便是与原矩阵一个最接近的近似。当然，这里的子矩阵因为布尔乘积的问题，一定是矩形。所以根据此引理，我们的任务就是找子矩阵们。
• 然后本文中一个比较重要的引理，就是UTLmatrix with direct SC1P。
  首先说说什么是C1P，C1P就是一个矩阵的的所有行向量中的1都是连续的。而SC1P则是C1P的升级版，不仅对行向量如此，对列向量也是一样。然后UTL是说，一个矩阵，它的列和和行和满足从左到右列和逐渐增加，从上到下，行和逐渐减小的的一个规律。
  引理是当一个矩阵满足UTL和directSC1P且没有全零行或者全零列时，则该矩阵中全1的的子区域由该矩阵的中位数列或中位数行来确定。

基于上面的引理，下面就进行介绍MEBF算法：

我们针对论文中的例子来进行讲解。

对于一个矩阵，我们需要尽可能的去找其median列或者行来去找到（覆盖）最大的稠密的矩形区域。（这是我们的思想）

给定一个Original矩阵$A_{11\ast 11}$。目标将其分解成$B_{11\ast 3}$和$C_{3\ast 11}$的乘积。

根据引理，我们将该矩阵构建成UTL矩阵，即对行与列进行排序，按照行和与列和来排序，将1s集中到矩阵的右上角去。

这其中的黄色列便是此时的median列，明显，这一列有6个1。
然后用这一列去覆盖其他列（即找到长为median列，宽尽可能大的，1尽量多的矩形）
很清晰，这个很好找。

这里整个红色与黄色部分就是我们找到的矩形（可能不连续，因为实际矩阵重排序之后也不一定是连续的，这个没办法，但是总体上是矩形，这就可以了，毕竟是重排序过的，原始矩阵也基本不是连续的，这个不影响），黄色的列是median列。
所以这就是我们找到的第一个子矩阵了（可以说是subMatrix，也可以说是patterns）。
这个子矩阵是一个矩形，那我们就可以得到B的第一列和C的第一行。
明显B的第一列就是我们的median列，而C的第一行就是这个矩形覆盖了哪些列，覆盖了为1，未覆盖为0。


在提取了第一个subMatrix或者说是patterns之后，我们要对矩阵进行一个更新，用原矩阵减去subMatrix，（这里因为我们得到的矩形是一个近似的嘛，可能将0的位置也覆盖到了，所以我们做减法的时候可能会得到负值，所以我们需要将负数置为0）。
然后就会得到一个新的矩阵，我们将它称为residual Matrix，我们就将其看做是新一轮迭代的original matrix。然后重新进行排序，形成UTL矩阵。



如此迭代下去，如果得到的patterns不能让整个矩阵的值变得更小，那么就需要进行weak signal detection，来发现patterns。（这里我们下面介绍具体算法的时候说）

最终的分解结果是这样，彩色的是分解正确覆盖到的，黑色的是未覆盖到的，可以说是误差。

然后分解的两个矩阵是：

MEBF算法流程及分析

下面介绍以下算法流程

首先输入是Original矩阵X，然后t是用于是该列（行）是否覆盖的一个标准，还有一个字母是收敛的一个判定（这个在代码中比较好理解）
输出自然不用说，两个布尔矩阵。

进入循环，当不收敛时一直继续，
首先是使用Bidirectional_growth算法来获取subMatrix的分解结果列向量a和行向量b（结合上面讲的MEBF的思路，应该就知道这里a是什么，b是什么了）。（Bidirectional_growth就是Median expansion）
然后下面就是算一下这里的提取出来的这个patterns是否能让矩阵的误差减小，如果减小的话，便更新，如果不能减小的话，则说明，目前矩阵的信号很小了，需要用weak signal detection来处理。
之后对矩阵进行一个更新。

这便是MEBF整体的一个逻辑。
下面来讲Bidirectional_growth算法。
这个算法就是median expansion操作

这里，很明显，首先对传入的矩阵进行一个重排序操作，使之UTL。然后我们这里分别找median row和median column来覆盖，看看哪个覆盖的效果更好就采用哪个（覆盖的意思就是，我们要用median列去构造一个矩形，尽可能的多包含原矩阵中的1，就相当于我们以median column为刷子，向右刷，尽可能多刷到1，这就是覆盖，也就是expansion）。
举个例子：对于median column，
这里是取引进UTL过的矩阵X的中位数列，即median column赋给d，然后以d来去覆盖其他列，怎么判断d能不能覆盖某一列呢，这里的操作是计算该列与d的相似性，如果相似性高于t，则便是能覆盖，否则便是不能覆盖，所以这里的t等于就是判断两列是否能够近似等价的一个阈值。
明显$t\in (0,1]$，t越大，覆盖的越少，但结果越准确，t越小，覆盖的越多，结果越粗糙。
如此e便是一个形如[0,0,0,0,0,1,1,1,1,0]的一个向量。
在判断完这两个向量的乘积是否都能让矩阵的误差便小，判断其与median行的误差谁小之后，边会输出一个误差最小对应的两个向量了。

当信号微弱，不足以使得误差下降时，就会启用weak signal detection算法

Bidirectional_growth不能保证整体的误差能够连续的变小，所以此时需要weak signal detection算法，to identify from a residual matrix，这是原话。
分别从行与列入手，看谁误差小。
以列举例，找出两个列和最大的列，让它们做与运算，则得出一个新的列d，用这个d来去expansion，从而得到e。
如此便是weak detection。

以上便是对于MEBF的介绍与分析。

代码实现

这里的代码是根据论文进行了些许的修改，对原文贴出的代码进行了大幅的修改，因为原文的代码跑不通，而且代码bug不少（可能是其主要实现是R语言，python版本的并未实现成功），将weak signal detection和Bidirectional growth算法结合在了一块。这里我还加上了详细的注释，便于理解。
具体代码如下

'''
author:pmy
refer-url:https://blog.csdn.net/qq_34687559/article/details/109109766
time:2020-10-16-9:40
last-edit:2020-10-16-9:40
note:
对原论文代码进行改进 
改进的方式很简单，就是e = np.sum(M1（residual matrix）- np.outer(B1_use,B2.use)) 这里只要e值大于0就行
然后比较一下这里的e1和e2谁大，然后取对应的行或者列就行，这里就很轻松
这里要把误差计算给封装成方法，不然使用起来，修改起来也很麻烦
paper: Fast and Efficient Boolean Matrix Factorization by Geometric Segmentation
progress: done
version:v1
'''



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


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

布尔矩阵分解测试

下面对该代码进行测试，我是直接用的论文中的Figure3 a1的original matrix(11*11)进行的测试。
测试代码如下：

import numpy as np

#论文中实例的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")
1
2
3
4
5
6
7
8
9
10
11
12

运行结果：

分解与重构的效果：
（原矩阵）

（分解因子重构矩阵）

结合误差，可以看到重构的误差还是可以的，当然，因为修改了代码逻辑，所以分解的情况是与原文不大一样的，然后我这里dim设置成了6，比原文要多一些。

我之后还拿自己的数据集测试了一下。
结果为：试了一个矩阵400*1174的，它的误差停在了3566，当然，原始矩阵误差（1s）在3w多，因为只给了100个因子，变成200个之后，能将误差讲到800。应该还是不错的了，因为800/3w的话在2%。
总体效果还可以，不过同样，有改进的余地。

参考文献

Wan C, Chang W, Zhao T, et al. Fast And Efficient Boolean Matrix Factorization By Geometric Segmentation[J]. arXiv, 2019: arXiv: 1909.03991.