Question: I have a digraph consisting of a strongly connected component (blue) and a set of nodes (orange) that are the inputs to it. The challenge is to break as many cycles as possible with a minimum of removed edges. In addition, there must be a path from each orange node to each blue node.  I solve the problem with a brute force:  Removing the random edge Check for a path from every orange node to every blue one. If everything is ok, I add an edge to the list and count the number of cycles. I return the edge to the graph and go to step 1 until I iterate over all the edges Next, from the resulting list (of length n) I generate combinations C (n, k) where k = {2 ... n} I perform operations 1, 2, 3 for all combinations of edges  The core of the code looks like this:  for level in range(2, len(edges)): stop = True edges2 = combinations(edges,level) for i, e in enumerate(edges2): g.remove_edges_from(e) test = True for node in orange_nodes: d = nx.algorithms.descendants(g, node) test = blue_nodes == d if not test: break if test: stop = False cycles_count = len(list(nx.simple_cycles(g))) print(f'{i}\t{level}\t{cycles_count}\t{e}') g.add_edges_from(e) if stop: break  Questions:  Is it possible to somehow optimize the code (nx.algorithms.descendants() and nx.simple_cycles() are dramatically slow)? Is it possible to rewrite code using Spanning tree or Feedback arc set? Maybe there is a fast search algorithm for not the best solution, but a good one?  Additionally: I rewrote the code as it is using the graph-tool, which gave a ~20x...50x speed boost. But this still does not allow us to approach the set practical task =(

Answer: The problem as stated is NP-Hard. Not sure if it is in NP either. In order to verify NP-hardness of the problem, consider graphs such that every blue node has an incoming edge from an orange node. For such graphs, what we need is that the graph after removing edges continues to be strongly connected. We also assume that maximum number of cycles need to be removed. Now, in order to break as many cycles as possible with a minimum of removed edges, assume that the maximum number of cycles that can be removed for a graph G while continuing to be strongly connected be removable(G) = k. This is a well-defined quantity for any graph G. Thus we need a graph G' that is a subgraph of G with number of cycles being cycles(G)-k. Now maximizing k is equivalent to minimizing the number of cycles that survive in G'. This is what makes the problem hard. Consider the Hamiltonian Cycle problem that is known to be NP-hard. Assume we have a program breakCycles(G) that computes a graph G' as a subgraph of G with maximum number of cycles removed (with minimal number of edges removed) or cycles(G') = cycles(G) - k. Then, it is straightforward to see that the Hamiltonian cycle problem can also be solved using breakCycles(G) by just providing input graph G to breakCycles to obtain the graph G' and return true iff G' is a simple cycle involving all vertices (of G).  Update : In order to obtain a practical solution, let's look at obtaining a graph with minimal cycles, that is a subgraph of the blue nodes such that removing any edge will result in loss of connectivity for those nodes that have an orange node incident to it. Solving the above problem is much easier and should work well in practice def getRemovableEdges(G, edgeLst, initNodes): import networkx as nx removableEdgeLst = [] for (u,v) in edgeLst: G.remove_edge(u, v) f = nx.floyd_warshall(G) addEdge = True for s in initNodes: if 'inf' in list(map(str, f[s].values())): G.add_edge(u,v) addEdge = False break if addEdge: removableEdgeLst.append((u,v)) return removableEdgeLst  To try it on the example provided, we need to first initialize the graph DG = nx.DiGraph() DG.add_nodes_from(range(1,8)) DG.add_edges_from([(1,2), (2,3), (3,4), (3,5), (4,5), (5,1), (5,4), (5,7), (6,4), (7,6)]);  With our graph initialized above, we execute the function as below...  In [5]: %time eL = getRemovableEdges(DG, list(DG.edges()), [2, 5]) CPU times: user 791 µs, sys: 141 µs, total: 932 µs Wall time: 936 µs In [6]: DG.remove_edges_from(eL); In [7]: plt.subplot(121) ...: nx.draw(DG, with_labels=True, font_weight='bold'); ...: plt.show();  We get the graph as below,