首先要知道什么是割(cut)。割是把圖的節點划分成兩個集合S和T,那么有一些邊的端點是分別處於S和T中的。所謂最小割就是使這種邊的數目最少的划分。
理論分析
Karger算法是隨機算法,它的描述很簡單:
每次隨機選擇一條邊,把邊的兩個端點合二為一。原來與這兩個點鄰接的點,現在把邊連到合並后的節點去,把原來的點和邊刪除。
合並后可能會有平行邊,在鄰接矩陣里記錄邊的數目。把形成的自環刪除。可以看到,合並前的兩個點與“外界”的連接邊數、合並后的點與“外界”的連接邊數(即合並點的度數),這兩者是一樣的。所以當合並至只剩兩個點時,相當於將原圖的點划分成兩個集合。這兩個點之間的連邊數就是形成的割的邊數。
現在分析下算法正確的概率。因為一個圖可能有多個最小割,我們來計算算法得到某個特定最小割的概率。設最小割交叉邊的數目為c。那么圖中每個點的度數至少為c(因為如果某個點的度數小於c,可以把這個點和其他點分開,形成的最小割邊數小於c,與假設矛盾。注意這個結論在合並點過程中仍然成立,因為前面說過了,合並得到的點的度數等於它代表的集合與外界的連接數)。如果圖有n個節點,那么至少有n*c/2條邊。我們不能選特定的c條邊,否則就不是特定的割了。不選這c條邊的概率是(1-c/邊數)>= (1-2/n). 每合並一次,點的數目減一。那么在整個過程中都不選那c條邊的概率是>=(1-2/n)*(1-2/(n-1))*...*(1-2/3) = 2/(n*(n-1)) >= 1/n2.
這個正確概率看着好像不很高。不過我們運行算法多次。在每次運行中,沒得到特定最小割的概率<=1-1/n2。那么運行m次沒得到的概率<=(1-1/n2)m。由於1-1/n2<=e-1/n^2(因為1+x<=ex),當m=n2,失敗的概率小於1/e。當m=n2ln(n)時,失敗的概率更是不到1/n,考慮到一般圖的節點數量,可以忽略不計。
插播:n個節點的圖不同最小割的數量最多多少?
答案:n*(n-1)/2. 前面我們看到Karger算法找到特定最小割的概率>=2/(n*(n-1)). 設有r個最小割。因為得到這些割是互斥事件,所以得到任意一個最小割的概率>=2/(n*(n-1))*r. 因為概率最大為1,r最大為n*(n-1)/2。這是可以取到的。比如n個點的環,任意刪除兩條邊都得到一個獨特的最小割(最小割邊數為2),因此為n*(n-1)/2.
算法實現
1.
雖然這個算法的描述簡潔易懂,但要實現好還是讓我頗費思量。先來拋個磚,給個O(n2)的實現,n是結點數。這個實現按照算法的字面描述,對於隨機選定的一條邊的兩個結點v1,v2(v1<v2),把和v2相關的邊合並到v1中去。具體的操作就是更新相關的點的度數和鄰接矩陣的項。程序沒有儲存邊,當隨機生成一個邊編號時,先用度數確定第一個端點,再根據這個端點的鄰接矩陣項找到另一個端點。這兩個過程都是O(n)的。
#include <bits/stdc++.h> using namespace std; class Graph { private: //For i = 1 and larger, adjList[i][0] = i. //For j = 1 and larger, adjList[i][j] is adjacent to i. vector<vector<int> > adjList; int n; vector<vector<int> > adjMatrix; vector<int> deg; private: void getEdgeEnds(int edgeIndex, int& row, int& col); int initForContract(); int contract(); public: int minCut(int nIter); int getNodeNum(); friend void buildGraph(char *fname, Graph& g); }; int Graph::getNodeNum() { return n; } void Graph::getEdgeEnds(int edgeIndex, int& row, int& col) { int num = 0; for(row=1;row<=n;row++) { if(num + deg[row] >= edgeIndex) break; num += deg[row]; } for(col=1;col<=n;col++) { num += adjMatrix[row][col]; if(num >= edgeIndex) break; } } int Graph::initForContract() { if(deg.size() == 0) { deg.assign(n+1, 0); for(int i=0; i<=n; i++) { vector<int> tmp(n+1, 0); adjMatrix.push_back(tmp); } } else { for(int i=1; i<=n; i++) for(int j=1; j<=n; j++) adjMatrix[i][j] = 0; } int edgeNum = 0; for(int i=1; i<=n; i++) { deg[i] = 0; for(int j=1;j<adjList[i].size();j++) { int node = adjList[i][j]; adjMatrix[i][node] = 1; deg[i]++; } edgeNum += deg[i]; } return edgeNum; } int Graph::contract() { int edgeNum = initForContract(); for(int iter = 0;iter < n-2; iter++) { int edgeIndex = rand()%edgeNum+1; int v1, v2; getEdgeEnds(edgeIndex, v1, v2); if(v1 > v2) swap(v1, v2); for(int i=1; i <= n; i++) { if(i == v1) { deg[v1] -= adjMatrix[v1][v2]; edgeNum -= 2*adjMatrix[v1][v2]; adjMatrix[v1][v2] = adjMatrix[v2][v1] = 0; } else if(adjMatrix[v2][i]) { adjMatrix[v1][i] += adjMatrix[v2][i]; adjMatrix[i][v1] = adjMatrix[v1][i]; deg[v1] += adjMatrix[v2][i]; adjMatrix[v2][i] = adjMatrix[i][v2] = 0; } } deg[v2] = 0; } return edgeNum/2; } int Graph::minCut(int nIter) { srand(time(NULL)); int mincut = n*n; for(int i=0; i<nIter; i++) { mincut = min(mincut, contract()); } return mincut; } void parse(char str[], vector<vector<int> > &adjList) { int len = strlen(str); int offset=0; int node; sscanf(str,"%d%n",&node,&offset); vector<int> lst(1, node); while(offset < len) { int neibor, count; sscanf(str+offset, "%d%n", &neibor, &count); offset += count; lst.push_back(neibor); if(offset == len-1) break; } adjList.push_back(lst); } void buildGraph(char *fname, Graph &g) { freopen(fname, "r", stdin); vector<int> tmp(1, 0); g.adjList.push_back(tmp); char str[200]; while(gets(str)!=NULL) parse(str, g.adjList); g.n = g.adjList.size()-1; } int main() { time_t st = clock(); Graph g; buildGraph("kargerMincut.txt", g); int n = g.getNodeNum(); int nIter = n*n*int(log(n*1.0)); int mincut = g.minCut(nIter); time_t ed = clock(); printf("%f\n",(double)(ed-st)/CLOCKS_PER_SEC); printf("%d\n",mincut); return 0; }
2.
另外的方法是事先打亂邊的順序,然后遍歷邊,把邊的兩端點連起來(而不是刪邊了),直到只剩兩個連通分量為止,剩下的邊里端點在不同連通分量的邊的數目就是割的大小。
查詢邊的端點是否處於不同的連通分量,(是的話)把兩個連通分量連起來可以用並查集(此時總連通分量數減一)。
#include <bits/stdc++.h> using namespace std; class UnionFind { protected: int componentNum; vector<int> parent, rank; void MakeSet(int n); int Find(int x); void Union(int x, int y); }; void UnionFind::MakeSet(int n) { componentNum = n; parent.assign(n+1, 0); for(int i=1; i<=n; i++) parent[i] = i; rank.assign(n+1, 0); } int UnionFind::Find(int x) { if(x != parent[x]) parent[x] = Find(parent[x]); return parent[x]; } void UnionFind::Union(int x, int y) { int rootX = Find(x), rootY = Find(y); if(rootX != rootY) { if(rank[rootX] < rank[rootY]) parent[rootX] = rootY; else if(rank[rootY] < rank[rootX]) parent[rootY] = rootX; else { parent[rootY] = rootX; rank[rootX]++; } componentNum--; } } class Graph:public UnionFind { private: //For i = 1 and larger, adjList[i][0] = i. //For j = 1 and larger, adjList[i][j] is adjacent to i. vector<vector<int> > adjList; int n; vector<pair<int, int> > edgeList; private: int contract(); void makeEdgeList(); public: int minCut(int nIter); int getNodeNum(); friend void buildGraph(char *fname, Graph& g); }; void Graph::makeEdgeList() { for(int i=1; i<=n; i++) for(int j=1; j<adjList[i].size(); j++) { int node = adjList[i][j]; if(i < node) { pair<int, int> edge(i, node); edgeList.push_back(edge); } } } int Graph::getNodeNum() { return n; } int Graph::contract() { random_shuffle(edgeList.begin(), edgeList.end()); int i; for(i = 0; componentNum > 2; i++) { int v1 = edgeList[i].first, v2 = edgeList[i].second; Union(v1, v2); } int cut = 0; for(; i<edgeList.size(); i++) { int v1 = edgeList[i].first, v2 = edgeList[i].second; int root1 = Find(v1), root2 = Find(v2); if(root1 != root2) cut++; } return cut; } int Graph::minCut(int nIter) { makeEdgeList(); srand(time(NULL)); int mincut = n*n; for(int i=0; i<nIter; i++) { MakeSet(n); mincut = min(mincut, contract()); } return mincut; } void parse(char str[], vector<vector<int> > &adjList) { int len = strlen(str); int offset=0; int node; sscanf(str,"%d%n",&node,&offset); vector<int> lst(1, node); while(offset < len) { int neibor, count; sscanf(str+offset, "%d%n", &neibor, &count); offset += count; lst.push_back(neibor); if(offset == len-1) break; } adjList.push_back(lst); } void buildGraph(char *fname, Graph &g) { freopen(fname, "r", stdin); vector<int> tmp(1, 0); g.adjList.push_back(tmp); char str[200]; while(gets(str)!=NULL) parse(str, g.adjList); g.n = g.adjList.size()-1; } int main() { time_t st = clock(); Graph g; buildGraph("kargerMincut.txt", g); int n = g.getNodeNum(); int nIter = n*n*int(log(1.0*n)); int mincut = g.minCut(nIter); time_t ed = clock(); printf("%f\n",(double)(ed-st)/CLOCKS_PER_SEC); printf("%d\n",mincut); return 0; }
3.
我們來看看Karger論文里的方法。我們需要確定的其實是使得只剩兩個連通分量的邊序列的最短前綴。知道了這個,在剩下的邊里求交叉邊的數量就容易了。注意到遍歷邊的過程中,連通分量的數目是單調減少的,因此我們可以二分。然而每次需要O(m)(m是邊數)來確定連通分量數目(用bfs之類),總復雜度O(m*log(m)),不夠好。
另一種方法是先選前m/2條邊,如果形成了一個連通分量,那么多了,考慮前m/4。如果形成多於兩個,再考慮后m/2的前一半。如果剛好是兩個,則可以直接計算割了。其實也是二分,關鍵是在已有的圖上加邊看連通分量的變化,還沒想好怎么寫。
參考資料:
Algorithms: Design and Analysis, Part 1
Karger, David (1993). "Global Min-cuts in RNC and Other Ramifications of a Simple Mincut Algorithm"
