枚举无向图中的所有循环
5.00/5 (14投票s)
找到一个基本的循环集，形成一个完整的基，以枚举给定无向图的所有循环。
1. 引言
图可以用于许多不同的应用，从描述电路的电子工程到描述分子网络的理论化学。有时有必要枚举图中的环，或者在图中查找满足特定条件的某些环。在本文中，我将解释原则上如何枚举图中的所有环，但我们会看到这个数量很容易增长，以至于不可能遍历所有环。然而，枚举所有可能环的能力允许使用蒙特卡洛或进化算法等启发式方法来回答有关图中环的特定问题（例如，查找最小或最大的环，或特定长度的环），而无需实际访问所有环。
这里，我将讨论无向无权图（参见图 1a 的示例），但该算法可以轻松地转换为有权图。
2. 基础知识
在本节中，将解释理解后续各节绝对必要的所有工具。
a) 图的表示
第一个主题是如何在程序代码中表示给定的图（例如，如图 1a 所示）。一种常见且实用的方法是邻接矩阵 (A)。它由 NxN 个元素组成，其中 N 是图中的节点数。如果两个节点 \(i\) 和 \(j\) 相连，则每个元素 \(A_{ij}\) 等于 1，否则等于零。图 1a 所示图的邻接矩阵如图 1b 所示。由于我们处理的是无向图，因此邻接矩阵是对称的，也就是说，只需要下半部分或上半部分即可完全描述该图，因为如果节点 A 连接到节点 B，则自动意味着 B 连接到 A。此外，对角线元素也被忽略，这些元素仅用于指示一个节点连接到自身。因此，这会自动成为整个图的一个基本节点，因为它无法进一步分割。
代码提供了一个名为 HalfAdjacencyMatrix 的类，用于表示图。如上所述，它仅存储矩阵的一半，并另外忽略对角线元素。
该类还可以用于存储图中的环、路径或任何类型的子结构。在这种情况下，可能存在不属于子结构的节点，因此没有边。邻接矩阵也可能包含两个或更多个不相交的子结构（见下文）。
b) 合并两个路径/环
为了确定一组基本环并稍后枚举图中的所有可能环，需要能够合并两个邻接矩阵（可能包含路径、环、图等）。
接下来，我们将通过对两个邻接矩阵的每个边应用逻辑 XOR 运算符来完成此操作。在接下来的两个示例中，将展示 XOR 运算符如何用于生成合并的路径和环。
在图 2a 中，XOR 运算符应用于从给定图中的根元素发出的两条路径。结果是闭环 B-C-D-B，其中根元素 A 被排除。此方案将用于根据第 3 节中所述的图的生成树生成基本环。
图 2b 中合并了两个环，生成了一个新环。此方案将在第 4 节中使用，从图的环基中形成新环。
XOR 运算符（operator^）的实现非常简单。该函数遍历矩阵中存在的每个位，并单独对每个位（边）应用 XOR。该类还提供了 operator^= 以方便使用。
// performs a xor operation on the two matrices and returns a new one.
inline HalfAdjacencyMatrix operator^(const HalfAdjacencyMatrix& rhs) const
{
if(m_nNodes != rhs.m_nNodes)
throw std::runtime_error("HalfAdjacencyMatrix::operator^():
The two matrices MUST be of the same size!");
// constructor initializes A_{ij} = 0
HalfAdjacencyMatrix result(m_nNodes);
for(size_t i = 0; i < m_aBits.size(); ++i)
{
// XOR for each bit: If the bit is true for any of the two matrices
// AND the bits in both matrices are not equal
// the bit is again true in the result matrix.
if((m_aBits[i] || rhs.m_aBits[i]) && (m_aBits[i] != rhs.m_aBits[i]))
{
result.m_aBits[i] = 1;
++result.m_nEdges;
}
}
return result;
}
3. 查找基本环集
a) 生成树
此处描述的算法遵循 Paton [1] 发表的算法。核心思想是从无向图中生成一个生成树。生成树构建完成后，我们需要查找图中存在但树中不存在的所有边。将其中一条缺失的边添加到树中将形成一个环，称为基本环。请注意，根据选择的根节点和构建树的方式，一个图可以有许多不同的生成树。因此，每棵生成树都构建了自己的基本环集。所有基本环构成一个环基，即图的环空间的基。由于基是完整的，因此使用哪种生成树来生成环基并不重要，每个基都同样适用于构建图中所有可能的环。图 1a 所示示例图的两种可能的生成树如图 3 所示，它们分别使用深度优先（a）和广度优先搜索（b）构建。红色的虚线显示了树中缺失但图中存在的边。
请注意，Paton 倾向于使用深度优先搜索而不是广度优先搜索，因为使用深度优先搜索时，每个节点仅与主分支相差一条边。这可以用于更有效地构造基本环。为了简单起见，我使用 XOR 运算符合并生成树的两条路径，因此深度优先和广度优先搜索的效率都一样。
b) 环基的大小
如上一节所述，环基中的基本环会因选择的生成树而异。但是，基本环的数量始终相同，并且可以轻松计算。
对于任何给定的无向图，具有 \(V\) 个节点和 \(E\) 条边，基本环的数量 \(N_{\text{FC}}\) 是
假设图最初是完全连接的 [2]。这个数字也称为“环秩”或“电路秩”[3]。
c) 伪代码
V: All Vertices
E: All Edges
pick start Vertex v in V
T = {v} // Spanning tree
Q = {v} // Process Queue
while not Q.empty:
i := Q.top
for j < |V|
if (i,j) in E
if j in T
// Fundamental Cycle Found!
pi := T.pathFromRootToElement(i)
pj := T.pathFromRootToElement(j)
cycle := merge pi and pj
else
Q.push(j)
insert j to T and set i as parent of j
上面的伪代码找到了给定图（由 V 和 E 描述）的一组基本环。
这里有几件事需要解决
- 如何选择起始顶点
v？
如前所述，生成树将取决于选择的起始顶点v，因此基本环集也将取决于。但是，由于环基是等价的，所以选择哪个顶点并不重要，因此在此实现中，只选择图中的第一个节点。 - 如何合并两条路径
pi和pj？
合并两条路径的最简单方法是 XOR 运算符。
这里，直接去掉了同时存在于两条路径中的所有边，即示例中的 A-C 和 C-F。请注意，边 E-D（红色虚线）不存在于树中，必须添加它才能得到基本环。
d) 实现
实现遵循标准的深度搜索算法。一旦找到一个已经访问过的节点，就找到了图的一个环。通过使用 XOR 运算符合并到当前节点和找到的节点的路径，即可获得由邻接矩阵表示的环，并将其存储在类中以备将来使用。
void Graph::computeFundamentalCycles()
{
// Lazy evaluation; save the fundamental cycles in the Graph class and
// just compute it ONCE when its needed
if(!m_aFundamentalCycles.empty())
return;
std::unique_ptr<TreeNode[]> aTree(new TreeNode[m_aNodes.size()]);
std::stack<size_t> nodeStack;
// Start arbitrarily with the first Node!
nodeStack.push(0);
// Copy the adjacency matrix as it will be necessary to remove edges!
HalfAdjacencyMatrix adjMat = m_adjMat;
// At the beginning, all tree nodes point to itself as parent!
// The tree is built on the fly
for(size_t i = 0; i < m_aNodes.size(); ++i)
{
aTree[i].parent = &aTree[i];
aTree[i].index = i;
}
// Loop until all nodes are removed from the stack!
while(nodeStack.size() > 0)
{
// Next node index:
size_t currentNodeIndex = nodeStack.top();
nodeStack.pop();
TreeNode& currentTreeNode = aTree[currentNodeIndex];
// Iterate though all edges connecting this node:
for(size_t j = 0; j < m_aNodes.size(); ++j)
{
// not connected?
if(!adjMat.isConnected(currentNodeIndex, j))
continue;
// Is the foreign node already in the tree?
// This is the case, if the parent element of the TreeNode does not point to itself!
if(aTree[j].parent != &aTree[j])
{
// Fundamental Cycle found!
// Get unique paths from both nodes within the spanning tree!
HalfAdjacencyMatrix pi(m_aNodes.size()), pj(m_aNodes.size());
unique_tree_path(&aTree[currentNodeIndex], pi);
unique_tree_path(&aTree[j], pj);
// also the connection between currentNodeIndex and j has to be inserted
// to ONE of the two paths (which one does not matter)
pi.connect(currentNodeIndex, j);
// combine the two matrices with XOR (^) to obtain the fundamental cycle.
m_aFundamentalCycles.push_back(pi ^ pj);
}
else
{
// The foreign node is not contained in the tree yet; add it now!
aTree[j].parent = ¤tTreeNode;
// add the node to the search stack!
nodeStack.push(j);
}
// Either way remove this connection!
adjMat.disconnect(currentNodeIndex, j);
}
}
}
4. 遍历所有可能的环
在最后一个部分，我们使用基本环集作为基础来生成图中的所有可能环。由于基本环集是完整的，因此保证能够获得所有可能的环。
a) 如何合并基本环？
要再次合并两个环，可以使用 XOR 运算符。假设我们的算法找到了图 1a 示例中以红虚线表示的三个基本环（A-B-E-F-C-A；B-D-E-B；D-E-F-D）作为完整基。
例如，使用 XOR 合并两个环 B-D-E-B 和 D-E-F-D 将删除边 D-E 并生成圆 B-D-F-E-B（蓝线）。
然而，仅仅合并成对的圆是不够的，因为这样就找不到包围的环（A-B-D-F-C-A），而这个环只有在结合所有三个基本环，删除边 B-E、D-E 和 E-F 时才能获得。总的来说，有必要遍历所有可能的基本环元组，从对开始，一直到 \(N_\text{FC}\) 元组（基本环的总数）。
b) 如何表示基本环的元组？
在代码中，基本环的元组可以直接用长度为 \(N_\text{FC}\) 的位字符串表示。对于示例图，位字符串的长度将为 3，产生以下三个基本环（FC）的可能组合：
| 位字符串 | XOR 组合 | 环 |
| 100 | FC1 | A-B-E-F-C-A |
| 010 | FC2 | B-D-E-B |
| 001 | FC3 | D-E-F-D |
| 110 | FC1 ^ FC2 | A-B-D-E-F-C-A |
| 101 | FC1 ^ FC3 | A-B-E-D-F-C-A |
| 011 | FC2 ^ FC3 | B-D-F-E-B |
| 111 | FC1 ^ FC2 ^ FC3 | A-B-D-F-C-A |
在位字符串表示中，所有可能的环都被枚举，即被访问，如果枚举了所有长度为 \(N_\text{FC}\) 的位字符串的所有排列，其中 \(k\) 是字符串中 1 的数量，并且 \(2 \le k \le N_\text{FC}\)。
例如，如果一个图有四个基本环，我们将需要遍历所有位字符串的排列，1100、1110 和 1111，总共是 11 次迭代。
c) 组合学
此时，让我们来谈谈数学，看看这种方法是如何扩展的。从对开始，我们需要知道在长度为 \(N_\text{FC}\) 的位字符串中，有多少个长度为 2 的 1 的排列是可能的。这个数字直接由 \(N_\text{FC}\) 选择 2 的二项式系数给出。一般来说，如果我们想知道长度为 \(N_\text{FC}\) 的位字符串中，有多少个长度为 \(k\) 的 1 的排列是可能的，这个数字由 \(N_\text{FC}\) 选择 \(k\) 的二项式系数给出。要获得基本环的总组合数，需要将从 \(k=2\) 到 \(k=N_\text{FC}\) 的二项式系数相加，得出以下方程
因此，代码的扩展性与图中的基本环数量成指数关系。指数扩展性始终是一个问题，因为迭代次数非常多，一旦 \(N\) 增大，就不可能遍历所有组合。为了对扩展性有一个初步了解，我们估计一次迭代需要 10 毫秒才能计算完成。那么 \(N=10\) 需要 10 秒，而 \(N=35\) 大约需要 11 年。人们可以轻易地看出，一旦 \(N\) 足够大，一次迭代所需的时间就会变得可以忽略不计，从而导致一个无法解决的问题。
请注意，这仅在人们真的想枚举每一个可能的环时才成立。然而，对于大多数问题来说，能够原则上访问每个环而无需实际这样做就足够了，例如启发式算法、蒙特卡洛或进化算法。总的来说，因此，如果确实需要枚举图中的所有可能环，重新思考提出的问题是一个好主意。
d) 环验证
现在我们知道如何组合不同的基本环，仍然有一个问题与 XOR 运算符相关：用 XOR 运算符组合两个不相交的环将再次得到两个不相交的环。因此，必须验证每个组合以确保生成一个连接的环。
让我们先看看如何检查一对基本环是否生成一个连接的环。这相当直接，因为我们只需要应用AND运算符并检查是否有属于两个环的边。此检查可以直接集成到 XOR 操作中：如果操作清除了一个或多个边，则这两个环至少有一条边是共有的，并生成一个新的有效环。
对于更高级的元组，验证并不那么简单：考虑合并三个环，那么在 XOR 操作过程中清除至少两条边是必要的。但是，此测试不是充分的，因为这三个环中的两个可能具有两条共有的边，而第三个环是不相交的。一种选择是跟踪所有对，并检查边是否在有效对和第三个环之间被清除，但这将导致两个主要缺点：
- 在计算三元组之前，必须计算所有可能的基本环对。一旦我们必须处理四元组、五元组或更高阶的元组，所有“较低”元组都必须在更高阶元组被评估之前计算。因此，无法再随机访问任何可能的位字符串。
- 回想一下，根据组合学，这种方法需要大量的内存来存储有效的组合。
因此，我将使用一种非常简单的方法，这可能不是最高效的方法：对于每个 \(k>2\) 的 \(k\)-元组组合，将使用深度搜索算法来检查 CycleMatrix（类型定义为 HalfAdjacencyMatrix）中合并的子结构是否完全连接。这很容易实现，因为只需要计算访问过的边数。如果这个数字等于边的总数，那么该元组就构成了一个连接的环。
bool validateCycleMatrix(const CycleMatrix& m)
{
// We will use our knowledge on the cycle matrices we are using:
// We know that all nodes in the matrix which belong to the cycle have exactly 2 edges.
// when we now start a deep search from any node in the matrix and counting the path length
// to the starting node this length must be equal to the
// total number of edges
// Again this is exhaustive but it is a very simple approach validating the cycles
size_t pathLength = 0;
// Find any edge in the matrix:
for (size_t i = 0; i < m_aNodes.size(); ++i)
{
for (size_t j = 0; j < m_aNodes.size(); ++j)
{
if (m.isConnected(i, j))
{
// Increment the pathLength and start the recursion
++pathLength;
std::set<size_t> aVisited;
aVisited.insert(i);
validateCycleMatrix_recursion(m, pathLength, j, i, aVisited);
// Version 3:
// - From the recursion, the path length will not account
// for the last edge connecting the starting node
// with the last node from the recursion.
return pathLength + 1 == m.getNumEdges();
}
}
}
// When we are here, the matrix does not contain any edges!
throw std::runtime_error("Graph::validateCycleMatrix():
Given Cycle Matrix does not contain any edges!");
}
方法 validateCycleMatrix 仅获取要验证的 CycleMatrix。然后它查找第一条存在的边，并使用 validateCycleMatrix_recursion 递归地启动深度搜索（这与用于确定生成树的算法相同）。如果深度搜索访问的边数等于 CycleMatrix 中的总边数，则该环有效。
// i: The node which has to be investigated in the current step
// previousNode: The node which was investigated before node i; necessary to avoid going backwards
// startNode: The node which was investigated first; necessary to determine
// when the recursion can be stopped
void validateCycleMatrix_recursion(const CycleMatrix& m, size_t& pathLength,
const size_t i, size_t previousNode, std::set<size_t>& aVisited) const
{
// The path length is also a measure for the recursion steps.
// If the recursion takes too long, we abort it and throw an error message.
// If you expect cycles which are longer than 500 edges, you have to increase this number.
// Also note that there is a limit of maximal recursion levels which cannot be exceeded.
// If your cycles exceed that maximum length,
// you will have to come up with another validation method.
if (pathLength > 500)
throw runtime_error
("Graph::validateCycleMatrix_recursion(): Maximum recursion level reached.");
// Find the next connection of the given node, not going back
for (size_t j = 0; j < m_aNodes.size(); ++j)
{
// Are the two elements connected? Ensure that we are not going backwards
if (m.isConnected(i, j) && j != previousNode)
{
// Was this node not visisted before?
auto ppVisited = aVisited.find(j);
if (ppVisited != aVisited.end())
{
// This node was already visited, therefore we are done here!
return;
}
// This node was not visited yet, increment the path length and
// insert this node to the visited list:
++pathLength;
aVisited.insert(i);
// Call the next recursion:
validateCycleMatrix_recursion(m, pathLength, j, i, aVisited);
return;
}
}
// When we are here, we have found a dead end!
throw std::runtime_error("Graph::validateCycleMatrix_recursion(): Found a dead end!");
}
e) 实现
下面，枚举图中所有环所需的所有步骤都汇总到一个函数中，该函数尝试在可能的情况下保存类中的所有环。代码还提供了一个迭代器（CycleIterator），该迭代器遵循 C++ 输入迭代器。请注意，此函数的目的主要是说明如何将前面各节中描述的所有部分结合起来，如果给定图的环秩足够大，它将花费非常长的时间。
// Exhaustive!!!
void computeAllCycles()
{
// if the fundamental cycles are not determined yet do it now!
if(m_aFundamentalCycles.empty())
computeFundamentalCycles();
// all fundamental cycles also are cycles...
m_aCycles = m_aFundamentalCycles;
// The bitstring
std::vector<bool> v(m_aFundamentalCycles.size());
// Combine each fundamental cycle with any other.
// attention: not only pairing (M_i ^ M_j) is relevant but also all other tuples
// (M_i ^ M_j ^ ... ^ M_N)! quite exhausting...
// we pick r cycles from all fundamental cycles; starting with 2 cycles (pairs)
for(size_t r = 2; r <= m_aFundamentalCycles.size(); ++r)
{
// Fill the bitstring with r times true and N-r times 0.
std::fill_n(v.begin(), r, 1);
std::fill_n(v.rbegin(), v.size() - r, 0);
// The following code in the original source caused an error and is
// therefore replaced by the line above
// if (r < 5)
// std::fill_n(v.begin() + r + 1, 5 - r - 1, 0);
// Iterate through all combinations how r elements can be picked from N total cycles
do
{
// Building the cycle matrix based on the current bitstring
CycleMatrix M(m_aNodes.size());
size_t nEdges = 0;
for (size_t i = 0; i < m_aFundamentalCycles.size(); ++i)
if (v[i])
{
M ^= m_aFundamentalCycles[i];
nEdges += m_aFundamentalCycles[i].getNumEdges();
}
// as long as pairs are merged the validation is straightforward.
if(r == 2)
{
// When at least one edge was deleted from the adjacency matrix
// then the two fundamental cycles form one connected cycle
// as they shared at least one edge.
if(nEdges > M.getNumEdges())
m_aCycles.push_back(M);
}
else
{
// Here we have combined more than two cycles and the
// matrix is validated via depth-first search
if(validateCycleMatrix(M))
m_aCycles.push_back(M);
}
}
// the bitstring is build up with 11...00, therefore prev_permutation
// has to be used instead of next_permutation.
while (std::prev_permutation(v.begin(), v.end()));
}
}
5. 前景
代码可以轻松扩展以处理每条边的权重，并且使用位字符串表示每个环允许直接使用遗传算法来查找满足特定约束的最长路径或最短路径，而无需实际访问所有可能的环。
6. 代码
提供的代码包含所有描述的类和函数。还有一个示例代码，用于枚举图 1a 中所有环。CreateRandomGraph 函数通过给定的边连接概率生成一个随机图。代码在 VC++ 2017（Windows）和 GCC 6.4.0（Linux）上进行了测试。请注意，代码使用了一些 C++11 功能，因此必须使用 -std=c++11 或更高版本（GCC）进行编译。
7. 历史
2018 年 6 月 - 版本 2
不幸的是，原始帖子中存在一个代码错误，调试代码仍留在上传的版本中。在函数“Graph::computeAllCycles()”和“Graph::CycleIterator::next()”中替换了以下代码行：
if (r < 5)
std::fill_n(v.begin() + r + 1, 5 - r - 1, 0);
// Was replaced by
std::fill_n(v.rbegin(), v.size() - r, 0);
2018 年 9 月 - 版本 3
我上传了一个补丁，用于修复 validateCycleMatrix 方法中的一个错误：在第 666 行，将以下行
return pathLength == m.getNumEdges();
替换为
return pathLength + 1 == m.getNumEdges();
此更改是必要的，因为用于验证 CycleMatrix 的深度搜索算法确定了环的长度，但没有考虑到连接最后一个访问节点与起始节点的、封闭环的最后一条边。因此，CycleMatrix 中的总边数必须等于深度搜索算法获得的路径长度加一。
添加了一个比图 1a 稍大的图进行测试，以测试补丁。
代码在文章和下载源中都进行了更改。
参考文献
- [1] Paton, Keith (1969), "An Algorithm for Finding a Fundamental Set of Cycles of a Graph", Scientific Applications 12:514
- [2] Berge, Claude (2001), "Cyclomatic number", The Theory of Graphs, Courier Dover Publications, pp. 27–30
- [3] Circuit rank. In Wikipedia. Retrieved January 07, 2018, from https://en.wikipedia.org/wiki/Circuit_rank