《An Algorithm for Subgraph Isomorphism》论文总结

人生第一次看论文，用时5天，这里记录一下论文大体内容和自己的理解与对算法的修改。

论文链接http://theory.stanford.edu/~virgi/cs267/lecture1.pdf

子图同构的概念请自行了解，可以百度可以看wiki。通俗的讲两个图同构的意思就是对其中一张图的结点重新编号以后，两张图长的一样。也可以说是一种对应关系，只要点与点的对应关系合适，那么两张图的邻接矩阵一样。子图同构就是有A,B两个图，看能不能在B图中找到一个子图，这个子图和A同构。

子图同构的问题是NP完全问题，我们没有办法经过一定的步骤推导出这个合适的映射关系，只有办法验证枚举出的映射矩阵是否可以满足子图同构的定义。验证的方法也很简单，论文里也说了，就是

定理1： 
如果图A关于映射f子图同构与图B,令 

MC=M′(M′MB)T

其中T表示的是矩阵的转置。

则有: 

所以我们很自然的想法就是通过DFS把所有可能的映射矩阵找出来，然后放到定理一的式子里挨个判断。这也是论文第二部分讲的。今天的第一个重点就在这里了，第二部分的伪代码有问题。

第二部分的伪代码在第5步的时候，也就是第4步执行完或者第7步执行完的时候，会在当前深度重新找一个k。但是在重新找k之前，第7步没有还原F[k] = 0，第4步没有还原M，所以是枚举不完全的。第一次看论文不太懂这是什么情况，但是我按照自己的想法还原F和M之后，枚举的结果就完全了，树就会search完全。下面两份代码，第一份是我自己写的DFS代码，递归搜索。第二份是我根据伪代码实现出来后再修改第5步的C++代码。另外第二份代码还在第3步以及第7步对H进行了实时的修改，这个无关紧要，只是为了变量跟得上该表示的状态，不改也可以。

论文的第三部分就是论文的核心的，一种采用refine过程的求子图同构映射矩阵的方法。refine的依据是论文里的定理2。一种同构的可能under M的意思就是说，M对应的是树的非终端节点，也就是没映射完的，而这个同构是M的子孙，那么这个同构就under M.如果M中的m[i][j] == 1，那么对于i,j，定理二就必须满足。所以说，如果对于特定的i和j，定理二没有满足，那么m[i][j]一定应该为0.这是对(i,j)的refine，每次refine都是对M整体的，只要i和j范围对就可以了。

定理2：

这个m[i][j] -> 定理2通俗来讲的意思就是说，如果m[i][j] == 1了，那么任何一个在A图中和i点（M下，在B的映射点是j点）邻接的x点，必然可以也在M下映射到一个B图中和j点映射的y点上。所以说，如果y存在，那么定理满足，m[i][j] = 1也满足。如果在M下找不到这样的y，那么M一定是有问题的，或者说m[i][j]应该为0，i怎么着也不应该映射到j。每次对M的一行，也就是i不变，j遍历完之后都要检查一下这一行还有没有1，如果这一行没有1了，那么说明i点找不到可以映射的，则个M下面的所有子孙都不必尝试，直接退出即可。

在一次对M的refine中，如果所有的1都可以存在，那么这次refine就结束，进行下个A中点的安排就可以。如果这次refine之后有1被改成0，说明这个M是不可用的，那就再继续对它refine，直到没有1被改变。

这里就有一个问题，万一对同一个点的两个邻接点的安排（找到的y）是同一个点呢，那这个M不是非法的吗？不是这样的，首先如果m[i][j] == 1，那么j的度一定是大于i的，另外最终的M里，同一个列不会有两个1.refine的过程只能把1变成0，不会把0变成1的。

按照论文里的伪代码，每次选出新的k之后，都refine一下，最终能走到根结点的M，一定都是可用的。

算法的原理就是这样。正确性在论文有所说明。如果refine对M不起作用，那么说明M保持了一种A到B的映射，使得所有在A中存在的边在B中也会存在，A中邻接了的点，在B中也一定邻接。这个提升还是很大的，因为不仅对深度搜索树进行了剪枝，而且第3步找k的判断也变少了。和第二部分介绍的算法相比，提高很多。

论文第4部分讲了实验，对算法的性能进行测试。对算法稍加修改，就可以解决找最大团问题（每个K都要消除序列，不找自同构）、图同构问题和有向图同构的问题。每个实验还与其他人的结果进行了对比。论文的最后提出了一种实现此算法的硬件设计。个人觉得这部分不是很难，边翻译边看就行了。

下面这份代码是refine版本的ULL算法，其中refine是我自己根据理解写的实现代码，refineX是根据论文附件1写的实现代码。拿出来大家一起交流。

第一次看论文，可能写的不是太好太清楚，希望大家斧正，给我提出意见。

#include<vector>#include<iostream>#include<time.h>using namespace std;#define DEBUG#define MATRIX vector<vector<int>>vector<int> getDegree(MATRIX m){int p = m.size();vector<int> degree(p, 0);for (int i = 0; i < p; i++){for (int j = 0; j < p; j++){degree[i] += m[i][j];}}return degree;}MATRIX multi(MATRIX& M1,MATRIX& M2){int rowN = M1.size();int columnN = M2[0].size();MATRIX result;for (int i = 0; i < rowN; i++){result.push_back(vector<int>(columnN, 0));for (int j = 0; j < columnN; j++){for (int k = 0; k < M2.size(); k++){result[i][j] += M1[i][k] * M2[k][j];}}}return result;}MATRIX transposition(MATRIX& M){int rowN = M.size();int columnN = M[0].size();MATRIX result;for (int j = 0; j < columnN; j++){result.push_back(vector<int>(rowN));for (int i = 0; i < rowN; i++){result[j][i] = M[i][j];}}return result;}void printMatrix(MATRIX& M){for (MATRIX::iterator it = M.begin(); it != M.end(); it++){for (vector<int>::iterator itN = (*it).begin(); itN != (*it).end(); itN++){cout << *itN << " ";}cout << endl;}}MATRIX getMatrix(int** a,int m, int n){MATRIX result;for (int i = 0; i < m; i++){result.push_back(vector<int>((int*)a + n*i, (int*)a + n*i + n));}return result;}bool refine(MATRIX& M,MATRIX& ma, MATRIX& mb){int elim = 0;for (int i = 0; i < ma.size(); i++){for (int j = 0; j < mb.size(); j++){//修改一行的Mif (M[i][j] == 1){bool flag = true;for (int x = 0; x < ma.size(); x++){if (ma[i][x] == 1){bool existY = false;for (int y = 0; y < mb.size(); y++){if (M[x][y] & mb[y][j] == 1){existY = true;break;}}if (existY){continue;}else{flag = false;}}}if (!flag){M[i][j] = 0;//这个映射不可用elim++;}}}bool one = false;for (int j = 0; j < mb.size(); j++){if (M[i][j] == 1){one = true;}}if (!one){return false;}}if (elim != 0)refine(M,ma,mb);return true;}//the super version of refinebool refineX(MATRIX& M, MATRIX& ma, MATRIX& mb){int pa = ma.size();int pb = mb.size();int elim = 0, i = 0,j = 0,k = 0,h = 0,x = 0;vector<int> lst;bool allZero;vector<int> degreeA = getDegree(ma);step1:elim = 0;i = 0;lst.clear();step2:k = 0;vector<int> sc(pa, 0);sc[0] = 1;h = 0;step3:allZero = true;for (int n = 0; n < ma[0].size(); n++){if (sc[n] & ma[i][n] == 1){allZero = false;break;}}if (allZero){goto step4;}else{lst.push_back(h);k = k + 1;}step4://sc = sc * 1/2for (int n = 0; n < sc.size(); n++){if (sc[n] == 1){sc[n] = 0;if (n + 1 <= sc.size() - 1){sc[n + 1] = 1;}break;}}h = h + 1;if (k != degreeA[i]){goto step3;}step5:j = 0;sc.clear();sc.assign(pb, 0);sc[0] = 1;step6:allZero = true;for (int n = 0; n < M[0].size(); n++){if (M[i][n] & sc[n] == 1){allZero = false;}}if (allZero){goto step9;}h = 0;step7:x = lst[h];allZero = true;for (int n = 0; n < M[0].size(); n++){if (M[x][n] & mb[n][j] == 1){allZero = false;break;}}if (allZero){goto step8;}h = h + 1;if (h != degreeA[i]){goto step7;}else {goto step9;}step8: for (int n = 0; n < M[0].size(); n++){M[i][n] = M[i][n] & (~sc[n]);}elim = elim + 1;h = h + 1;step9://sc = sc * 1/2for (int n = 0; n < sc.size(); n++){if (sc[n] == 1){sc[n] = 0;if (n + 1 <= sc.size() - 1){sc[n + 1] = 1;}break;}}j = j + 1;if (j != pb)goto step6;step10:allZero = true;for (int n = 0; n < M[0].size(); n++){if (M[i][n] == 1){allZero = false;break;}}if (allZero){goto FAIL;}i = i + 1;if (i != pa)goto step2;if (elim != 0){goto step1;}goto SUCCESS;FAIL:return false;SUCCESS:return true;}void refinULL(MATRIX &ma,MATRIX &mb){bool(*refineFUNC)(MATRIX& M, MATRIX& ma, MATRIX& mb) = &refine;int pa = ma.size(), qa = ma[0].size();int pb = mb.size(), qb = mb[0].size();vector<int> degreeA = getDegree(ma);vector<int> degreeB = getDegree(mb);vector<vector<int>> M0;for (int i = 0; i < pa; i++){M0.push_back(vector<int>(pb));for (int j = 0; j < pb; j++){M0[i][j] = degreeB[j] >= degreeA[i] ? 1 : 0;}}vector<int> F(pb, 0);vector<int> H(pa, -1);vector<MATRIX> M(pa);int k = -1;int sum = 0;step1:MATRIX CM = M0;int d = 0;if (!refineFUNC(CM, ma, mb)){goto END;}step2:bool existJ = false;for (int j = 0; j < pb; j++){if (CM[d][j] == 1 && F[j] == 0){existJ = true;break;}}if (existJ){M[d] = CM;k = d == 0 ? H[0] : -1;}else{goto step7;}step3:k = k + 1;if (CM[d][k] == 0 || F[k] == 1)goto step3;H[d] = k;//record it by timefor (int j = 0; j < pb; j++){if (j != k){CM[d][j] = 0;}}if (!refineFUNC(CM, ma, mb)){goto step5;}step4:if (d < pa - 1){goto step6;}else{//get a M'//cout << "---#" << ++sum << "---" << endl;//printMatrix(CM);#ifdef DEBUGcout << "---#" << ++sum << "---" << endl;for (int n = 0; n < H.size(); n++){cout << H[n] + 1 << " ";}cout << endl;#endif}step5:F[k] = 0;CM = M[d];//复原existJ = false;for (int j = k + 1; j < pb; j++){if (CM[d][j] == 1 && F[j] == 0){existJ = true;break;}}if (existJ){goto step3;}else{goto step7;}step6:H[d] = k;F[k] = 1;d = d + 1;goto step2;step7:if (d == 0){goto END;}else{H[d] = -1;F[k] = 0;d = d - 1;//CM = M[d];k = H[d];goto step5;}END:#ifdef DEBUGcout << "---end---" << endl;#endifreturn;}int main(){int a[3][3] = { { 0, 1, 1 }, { 1, 0, 0 }, { 1, 0, 0 } };int b[4][4] = { { 0, 1,0, 0 }, { 1, 0, 1, 1 }, { 0, 1, 0, 1 }, { 0, 1, 1, 0 } };MATRIX ma = getMatrix((int **)a, 3, 3);MATRIX mb = getMatrix((int **)b, 4, 4);clock_t t1 = clock();refinULL(ma, mb);clock_t t2 = clock();cout << "with refine time is : " << t2 - t1 << endl;return 0;}