Keywords: Flow decomposition; Integer Linear Programming; safety; RNA-seq; RNA transcript assembly; isoform

1 Introduction↩︎

1.0.0.1 Motivation.

Minimum Flow Decomposition (MFD) is the problem of finding a minimum number of weighted paths to decompose a flow on a directed graph, such that the sum of weights of the paths crossing an edge is equal to the flow on that edge. Hartman et al. [1] proved that MFD is NP-hard, even on directed acyclic graphs (DAGs), and even if the flow values are only in \(\{1,2,4\}\).

MFD is a crucial tool in many applications, for example, in networking, where the paths represent traffic going through a sequence of routers. Minimizing the number of these paths reduces the amount of maintenance work, among other additional benefits [1]. In RNA sequence reconstruction problems, MFD is used by e.g. [2]–[7] to decompose a splice graph, in which every node is a gene exon, and every weighted path is an RNA sequence that needs to be reconstructed. This splice graph can, for example, be constructed by aligning RNA-seq reads to a reference genome. MFD can also used to reconstruct viral quasispecies [8]. The input graphs in these applications are DAGs, and all previous work on MFD that we cite below is also on DAGs. As such, in this paper, all graphs are DAGs.

Moreover, also MFD variants are of practical interest. For example, one practical strategy is to incorporate longer reads as subpaths potentially spanning more than two nodes [2], [9]–[11], which must appear in some solution path (i.e. reconstructed transcript or viral genome). Moreover, the edge weights may not form an exact flow, due to sequencing errors and read mapping artifacts. An MFD variant adapted to this case is to consider intervals of edge weights instead, called inexact flows in [7], [12].

1.0.0.2 Related work.

Given the complexity of the problem, it is common to use heuristic methods [13], for example, by using a greedy approach, iteratively removing the path with the largest currently available flow until the flow is fully decomposed. Cáceres et al. [14] showed that for DAGs this approach performs well, with an approximation ratio of \(O(\log \mathrm{\small val}(f))\) (where \(\mathrm{\small val}(f)\) is the total flow of the graph), only if the width of the graph (the minimum number of paths to cover all edges) does not increase in the process; otherwise it gives an exponentially worse result than the optimal with an approximation ratio of \(\Omega(m/\log m)\).

Mumey et al. [15] gave an approximation algorithm on DAGs with ratio \(\lambda^{\log\lVert f\rVert}\log\lVert f\rVert\), where \(\lambda\) is the longest source-to-sink path length and \(\lVert f\rVert\) is the maximum-norm on the flow \(f\) (i.e., the largest flow value of the graph). Polynomial-time approximation algorithms with sublinear ratio in the size of the graph are not known, neither whether the problem is in APX, i.e. admits a constant-factor approximation algorithm.

Kloster et al. [16] proposed the first algorithm solving MFD optimally, which is linear-time fixed-parameter tractable (FPT), parameterized on the MFD size \(k\). While it performs well for small \(k\) values, it quickly becomes unfeasible for larger \(k\) values, since the parameterized runtime grows exponentially on \(k\) with a degree \(2\) polynomial in the exponent. This is a limiting factor to the usability of such algorithms in practical applications.

In a recent work, Dias et al. [17] formulated an Integer Linear Problem (ILP) model for MFD on DAGs and showed that it performs faster on graphs with larger solution sizes than the FPT algorithm of Kloster et al. [16]. Due to this better performance of an ILP, it is a natural to work on further optimizations. An additional strength of ILP formulations is their extensibility, since ILPs can easily be modified to handle other aspects of the applications, such as subpath constraints and inexact flows cited above, see [17].

1.0.0.3 Contributions.

In this paper, we significantly optimize the ILP for MFD using the notion of safe paths for all flow decompositions in DAGs, first studied by Khan et al. [18]. More specifically, a path \(P\) is said to be safe if all flow decompositions (of any size) have at least one path containing \(P\) as a subpath. Clearly, safe paths must appear in any MFD. We show how they can be used to reduce the ILP search space by using them to fix many binary variables of the ILP. At a high level, this approach is along the same lines as for other NP-hard problems, for example, Bumpus et al. [19] analyzed \(c\)-essential vertices, as those vertices belonging to any \(c\)-approximate solution, for some specific graph problems whose solutions are sets of vertices. Clearly, \(c\)-essential vertices also belong to all optimal solutions, and can be simply removed from the graph when running e.g. an FPT algorithm (see [19] for details).

In our case, safe solutions are sets of paths (not single vertices), and it is non-trivial how to use them to simplify the input graph, or how to integrate them into a combinatorial algorithm (e.g., we cannot simply remove them, they could overlap, etc.). However, incorporating safe paths into an ILP is easily supported since they can be modelled by additional constraints, or by fixing some ILP variables. Moreover, in order to use more than one safe path, we observe that pairwise unreachable paths (they are not subpaths of a superpath) must be present in different MFD solution paths. As a fast heuristic to select a set of such paths of large total length, we use a reduction to a maximum weight antichain problem [20]. In 10 we show that a maximum-length set of pairwise unreachable safe paths can also be found in polynomial time, but this procedure is overall computationally more involved, and thus may be too expensive as a pre-processing step.

Furthermore, as observed by Cáceres et al. [14], the size of any edge antichain (like ours) is a lower bound on the MFD size. As such, we can also use this lower bound to check if a heuristic MFD solution size attains this lower bound, in which case it is optimal. To the best of our knowledge, this is the first time such a (non-trivial) lower bound is used in a solver.Lastly, before running the ILP solver, we also simplify the input graph using the Y-to-V graph compaction, also used in the implementation of the FPT algorithm of Kloster et al. [16], but not used in the one of the ILP by Dias et al. [17]. Using all these optimizations, we experimentally show that we obtain significant speedups over the original ILP, of up to 55-90\(\times\) on the harder inputs.

In addition, we also adapt these optimizations to the MFD variants considering inexact flows and subpath constraints, proving that such adaptations can easily be done also for MFD variants of practical interest. Furthermore, we also show (in 4) how the Y-to-V graph compaction can be modified for these two variants, an issue which has not been considered before in the literature. With these optimizations, on the hardest graphs among the instances with inexact flows, we obtain an impressive speedup of 123\(\times\) compared to the original ILP for MFD with subpath constraints. As a last contribution, we tackle an MFD variant where the weights of the solution paths must belong to a given set. We show that this problem on DAGs admits a simpler and thus faster ILP than the one for MFD. When the set of given weights consists of powers of 2 (up to the maximum flow value of any edge), this ILP provides a \(\log||f||\)-factor approximation ratio for the original MFD problem. For practical applications, we show that if we also add the flow values of all edges to the given set of weights, then this ILP returns an MFD solution with a number of paths that is significantly closer to the minimum one than the state-of-the-art heuristic MFD solver \(\mathsf{Catfish}\) [13], on instances with large MFD size. These are of particular interest, since our optimized solver runs in under 2 seconds on instances of solution size at most \(10\). The ILP for this problem variant is up to a further order of magnitude faster than our optimized ILP, and running in under 1 second on average.

2 Preliminaries↩︎

2.0.0.1 Definitions.

We let \(G=(V,E)\) be a graph and assume throughout the paper that \(G\) is a DAG (i.e. directed and acyclic). For simplicity we also assume that \(G\) has no parallel edges and that it has a single source node \(s\) and a single sink node \(t\). We denote by \(\mathbb{N}\) the natural numbers including \(0\) and by \(\mathbb{Z}^+ = \mathbb{N}\setminus\{0\}\) the positive integers. A \(u\)-\(v\) path \(P\) is a sequence of edges going from the node \(u\) to \(v\) and its length \(|P|\) is defined as its number of edges. Additionally, we identify paths \(P\) with functions \(P:E\to\{0,1\}\) where \(P(e) = 1\) iff \(e \in P\). We call two paths \(P_1, P_2\) (and resp. edges) independent, if there exists no path \(P^*\) such that \(P_1\) and \(P_2\) are subpaths of \(P^*\). For a path \(P=(e_1,\dots,e_{|P|})\) we denote by \(P[\ell..r]\mathrel{\vcenter{:}}=(e_\ell,\dots,e_r)\) a subpath of length \(r-\ell+1\). A flow network is the tuple \(G = (V,E,f)\), where \(f:E\to\mathbb{Z}^+\) is a flow, i.e. a function that preserves the flow conservation \(\sum_{(u,v) \in E} f_{uv} = \sum_{(v,w) \in E} f_{vw}\), \(\forall v \in V \setminus \{s,t\}\).

By a \(k\)-flow decomposition of a flow \(f\) (see 1) we denote a family of \(s\)-\(t\) paths \(\mathcal{P}= (P_1,\dots,P_k)\) with their associated weights \(w=(w_1,\dots,w_k)\in(\mathbb{Z}^+)^k\), such that: \[\label{eqn:superposition95of95flow} \sum_{i \in \{1,\dots,k\}} P_i(u,v)w_i = f_{uv}, \quad \forall (u,v) \in E.\tag{1}\] In this paper, we focus on positive integer flow values, which is motivated by its application in RNA sequence reconstruction. Vatinlen et al. [21] proved that there is always a \((|E|-|V|+2)\)-flow decomposition. The Minimum Flow Decomposition (MFD) problem asks for a \(k\)-flow decomposition with minimum \(k\).

In case of imperfect data, Williams et al. [7] considered the problem variant Minimum Inexact Flow Decomposition (MIFD) where for each edge \((u,v)\) we have a lower bound \(L_{uv}\) and an upper bound \(R_{uv}\), and constraint 1 is changed to require that \(\sum_{i \in \{1,\dots,k\}} P_i(u,v)w_i \in [L_{uv},R_{uv}]\), \(\forall(u,v) \in E\). Another MFD variant of practical interest is Minimum Flow Decomposition with Subpath Constraints (MFDSC) [10]. The only difference with respect to MFD is that in the input we also have a set of paths \(\mathcal{R}\), called subpath constraints, and it is required that every path in \(\mathcal{R}\) is a subpath of at least one path of the decomposition.

2.0.0.2 ILP formulations.

We now review the ILP formulation for MFD given by Dias et al. [17]. To model each path \(P_i\), \(i \in \{1,\dots,k\}\), binary variables \(x_{uvi}\) are introduced to represent each edge \((u,v) \in E\). We set \(x_{uvi}=P_i(u,v)\), i.e. \(x_{uvi} = 1\) if \((u,v)\in P_i\) and otherwise \(x_{uvi} = 0\). At the same time, each path is required to start from the source and end in the sink. For all the intermediary nodes of a path, a unit in-degree and a unit out-degree is required. Those requirements can be modeled by the following constraints. For all \(v \in V\) and for all \(i \in \{1,\dots,k\}\):

\[\label{eqn:path95definition} \sum_{(u,v) \in E} x_{uvi} - \sum_{(v,u) \in E} x_{vui} = \begin{cases} 0, & \text{if v \in V \setminus \{s,t\}}, \\ 1, & \text{if v = t}, \\ -1, & \text{if v = s}, \end{cases} \quad\quad\tag{2}\]

Constraint 1 is then modeled as: \[\label{eq:ILP-flow-superposition} f_{uv} - \sum_{i \in \{1,\dots,k\}} x_{uvi}w_i = 0,\quad\forall (u,v) \in E.\tag{3}\]

The above constraint contains the non-linear term \(x_{uvi}w_i\) that standard linear solvers cannot solve. However, using basic linearization techniques, each such constraint can be replaced by three linear constraints, see [17] for details. For completeness, we give the full formulation in 9.1.

Dias et al. [17] used the flexibility of this ILP formulation to easily model also the two variants MIFD and MFDSC. The former problem can be trivially modeled in ILP by changing the (linearized version of) constraint 3 to state that \(\sum\limits_{i \in \{1,\dots,k\}} x_{uvi}w_i\) belongs to the interval \([L_{uv},R_{uv}]\). The latter problem can be modeled by introducing additional indicator variables \(r_{ij} \in \{0, 1\}\) for every \(i \in \{1,\dots,k\}\), and every path \(R_j \in \mathcal{R}\), modeling if the subpath \(j\) is contained in the \(i\)-th path of the flow decomposition, via the additional constraints: \[\begin{align} \sum_{(u,v)\in R_j} x_{uvi} \geq |R_j|r_{ij},\quad &\forall i\in\{1,\dots,k\}, \forall R_j\in\mathcal{R},\tag{4} \\ \sum_{i\in\{1,\dots,k\}} r_{ij} \geq 1,\quad & \forall R_j\in\mathcal{R}. \tag{5} \end{align}\]

3 Reducing the ILP search space via safe paths↩︎

In this section, we optimize the ILP model for MFD presented above based on the notion of solution safety (see e.g. Tomescu and Medvedev [22]). Then, in [rem:safety-for-variants,rem:safety-for-subpaths], we discuss how this method can also be applied to the MIFD and MFDSC problems.

We call a path safe if it is part of every (not necessarily optimal) flow decomposition; that is, if it is a subpath of some path in every decomposition (not necessarily of minimum size). Khan et al. [18] characterized all safe paths via the notion of excess flow, which can be thought of as the flow value of the first edge of the path minus the flow values of the non-path edges out-going from the internal nodes of the path (i.e., the flow “leaking” from the path). Formally:

Since we assume the flow to be positive, every edge (as a path of length \(1\)) is safe. Khan et al. [18] gave a simple algorithm for finding all safe paths, as follows. We can find all safe subpaths of each \(s\)-\(t\) path \(P=(e_1,\dots,e_{|P|})\) of an arbitrary flow decomposition with a two-pointer algorithm. Starting with \(\ell = 1\) and \(r = 2\) as an inclusive interval (i.e., including two edges), we increase \(r\) as long as \(f_{P[\ell..r]} > 0\), and increase \(\ell\) as long as \(f_{P[\ell..r]} \leq 0\). That is, for every index \(r=2,\dots,|P|\), we find the minimum \(\ell<r\) such that the subpath \(P[\ell..r]\) is safe. The runtime of this procedure is \(O(\text{out-degree}(P))\), where \(\text{out-degree}(P)\mathrel{\vcenter{:}}=\sum_{(u,v)\in P}\text{out-degree}(u)\), since flows of edges across \(P\) are at most added and subtracted twice throughout the algorithm. To find the set \(\mathcal{S}\) of all safe paths of \(G\), we can find any flow decomposition \(\mathcal{P}\) and run the two-pointer algorithm on all paths in \(\mathcal{P}\). For quickly obtaining a flow decomposition, one could e.g. use the greedy approach of removing paths of largest currently available flow in runtime \(O(|\mathcal{P}|\cdot(n+m))\) and use it as upper bound to the optimal solution in the ILP [21]. In our implementations for the three problems MFD, MIFD, MFDSC we use flow decompositons obtained from the fastest heuristic solvers for these problems.

An (edge) antichain is a subset \(Q\subseteq E\) such that every pair of edges \(e_1, e_2\in Q\) is independent. The size of any edge antichain is also a lower bound on the size of an MFD, as noted by [14], since each edge in the antichain must be traversed by different paths in any flow decomposition. The idea is to use such edges (and extensions of them via safe paths, as we discuss next) to fix some ILP variables.

For every edge \(e \in E\), denote by \(a(e)\) the length of the longest safe path traversing \(e\), i.e. \(a(e) \mathrel{\vcenter{:}}= \max_{P \in \mathcal{S}, e \in P} |P|\); this can be obtained from the computation of safe paths described above. Let \(\mathcal{P}(Q) = \{P_1,\dots,P_{|Q|}\}\) be the corresponding longest safe paths passing through the edges in \(Q\); see 2 in 8 for an illustration. If \(Q \subseteq E\) is an antichain, then also the corresponding longest safe paths for the edges of \(Q\) must be fully traversed by different paths in any flow decomposition (otherwise, if some \(P_i\) and \(P_j\) were traversed by the same \(s\)-\(t\) path in some flow decomposition, then also the corresponding edges in \(Q\) could be traversed by the same \(s\)-\(t\), contradicting the fact that \(Q\) is an antichain).

As such, without loss of generality, we can fix in the ILP the \(i\)-th path of the MFD to contain the \(i\)-th safe path in \(\mathcal{P}(Q)\): \[\label{eq:safety-constraints} x_{uvi} = 1,\quad \forall (u,v)\in P_i, \forall P_i\in \mathcal{P}(Q).\tag{6}\]

Additionally, in order to further optimize the choice of the antichain \(Q\), we consider the following notion. Given a weight function \(a:E\to\mathbb{N}\), a maximum weight antichain is an antichain \(Q\) maximizing \(\sum_{e\in Q}a(e)\) [20]. Maximum weight antichains can be found by a reduction to minimum flows, with demands on the edges given by their weights, in runtime \(O(m\cdot\sum_{e\in E}a(e))\subseteq O(n\cdot m^2)\), followed by a depth-first search through the graph. Maximum weight antichains are cut-sets with edges whose minimum flow is equal to their demand [20]. In order to heuristically fix as many \(x_{uvi}\) variables as possible, we set the weight \(a(e)\) of each edge \(e\) to be the length of the longest safe path traversing \(e\). The overall runtime of this preprocessing is thus \(O(m\cdot\sum_{e\in E}a(e))\). Note that this approach does not yield the maximum total length of independent safe paths, we describe an approach on how to find them in 10.

As further optimization, note that symmetries are well known to slow down ILP solvers. In the previous ILP formulations, the paths can be arbitrarily permuted, which we can mitigate with the following constraints: \(w_{i+1} \leq w_i, \forall i \in \{|Q|+1,\dots,k-1\}.\)

Note that we exclude the first \(|Q|\) weights as they have already been fixed to use the longest safe paths crossing the antichain \(Q\). This approach does not introduce more variables but adds constraints. Moreover, it can be effective in case the paths are forced to have widely different weights, which is the case in practical inputs for the RNA reconstruction problem [23]. We show the complete optimized ILP in 9.2.

Moreover, in our optimized ILPs we have an additional check on whether the lower bound \(|Q|\) equals the size of the heuristic flow decomposition (which we need in any case to initially compute the safe paths, as mentioned above). If this holds, then the heuristic flow decomposition has minimum size, and we directly report it, without running the ILP.

4 Y-to-V reduction for subpath constraints and inexact flows↩︎

Koster et al. [16] used an optimization to reduce the graph size in a pre-processing step by suitably contracting all edges entering nodes with in-degree 1, and respectively, exiting from nodes with out-degree 1, as we review below. This reduction has sometimes been used, under the name “Y-to-V”, to simplify the input graphs to other problems, see e.g. [22]. Note that the ILP solver by Dias et al. [17] did not include pre-processing step. In this section, we show how to extent this pre-processing step for the two problem variants MIFD and MFDSC (missing proofs are in 8).

The Y-to-V reduction for MFD works the following way. Consider the graph \((V,E)\) and let \(v\in V\) be a node of in-degree \(1\), and let \(u\in V\) be the unique node with \((u,v) \in E\). Let \(v\) have out-degree \(\ell\) and let \(w_i \in V\), \(i=1,\dots,\ell\) be all nodes with \((v, w_i)\in E\). The reduction creates a new graph \((V,E')\) with edges \((u,w_i)\) and flow \(f':E'\to\mathbb{N}\) with \(f'(u,w_i) = f(v,w_i)\) for \(i=1,\dots,\ell\). This process can be repeated until no nodes of in-degree \(1\) exist anymore, while all other edges are copied to the new graph. Analogously this can be done for nodes of out-degree \(1\).

For the MFDSC problem (with subpath constraints), note that subpaths are given with respect to the original graph, and hence we need to modify them for the Y-to-V contracted graph. The induced subgraphs of nodes of in-degree (resp. out-degree) \(1\) define a forest of out-trees (resp. in-trees). Subpaths can potentially pass through, begin in or end in such trees.

If a path completely passes through a tree, we can merely merge the edges of the tree in the path to the contracted edge. In that case, the subpath which intersects with the tree gets translated a single new edge. If a path stops inside an out-tree or begins inside an in-tree, we can not do that because we do not know which leaf the path will cross in the decomposition. In that case, the subpath which intersects with the tree gets translated to several new, parallel edges.

We first contract all out-trees and describe how to reduce the path to the contracted graph. We can follow up by contracting all in-trees analogously. We can partition the graph uniquely in maximal-size out-trees, such that no adjacent out-trees can be merged to a new out-tree. An out-tree in this partition consists of a single edge if and only if this edge will not be contracted in the Y-to-V reduction. Let the subpath \(R_j\in \mathcal{R}\) intersect the out-trees \(T_1,\dots,T_{\ell_j}\), ending at node \(u\in T_{\ell_j}\). We can contract the path to obtain \(R'_j\) in the following way. Let the set of leaf nodes of \(T_{\ell_j}\) that can be reached from \(u\) be denoted by \(L(T_{\ell_j}, u)\). For every tree \(T_1,\dots,T_{\ell_j-1}\), we add an edge from the root of the tree to the leaf the path is crossing to \(R'_j\). For \(T_{\ell_j}\) we add the parallel edges \((\text{root}(T_{\ell_j}), v)\) for all \(v\in L(T_{\ell_j}, u)\). That means \(R'_j\) is a path of length \(\ell_j-1\) followed by some parallel edges. Let the obtained set of “subpaths” \(R'_j\) be \(\mathcal{R}'\). We change the constraints to: \[\begin{align} \sum_{(u,v)\in R'_j} x_{uvi} \geq |\ell_j|r_{ij},\quad &\forall i\in\{1,\dots,k\}, \forall R'_j\in\mathcal{R}', \tag{7} \\ \sum_{i\in\{1,\dots,k\}} r_{ij} \geq 1,\quad & \forall R'_j\in\mathcal{R}'. \tag{8} \end{align}\] Note that only the set \(\mathcal{R}\) and the integers \(|R_j|\) change compared to the constraints .

The flow conservation is a necessary property for the Y-to-V reduction to work, since because of it, a single out-edge carries no further information. All paths of a flow decomposition entering a node of out-degree \(1\) must continue through this edge, and will together decompose it (and similarly for nodes of in-degree \(1\)). In the case of MIFD, we do not require the input to have any conservation of flow, and unlike in MFD, valid inputs can also be infeasible to solve.

For nodes of in-degree \(1\) (resp. out-degree \(1\)) we generalize flow conservation to the property of inexact flow conservation, which nodes must have in order for its edges to be contracted: \[L_\text{in} \leq \sum_{(u,v)\in E} L_{uv} \leq \sum_{(u,v)\in E} R_{uv} \leq R_\text{in}, \label{eq:inexact-fc}\tag{9}\] where \([L_\text{in}, R_\text{in}]\) is the interval of the single in-edge of \(u\) (respectively for the sum of all outgoing edges and the interval \([L_\text{out}, R_\text{out}]\) of the single out-edge of \(u\)). Note that the second inequality in 9 is always fulfilled since \(L_{uv}\leq R_{uv}\). The Y-to-V reduction for MIFD works the same as the Y-to-V reduction for MFD, but defines intervals \([L'_{uw_i},R'_{uw_i}]\) (resp. \([L'_{w_iu},R'_{w_iu}]\)) instead of flow values \(f'(u,w_i)\) (resp. \(f'(w_i,u)\)) and is restricted on nodes of in-degree (resp. out-degree) \(1\) fulfilling the inexact flow conservation.

5 MFD with given weights↩︎

In this section we consider an MFD variant in which the path weights in the solution are restricted to belong to a given set \(W = \{w_1, \dots, w_\ell\}\), see also 3 in 8 for an example and note that not all weights in \(W\) must be used by the solution paths. This problem was also defined by Kloster et al. [16], and in this paper we call it MFD with given weights (MFDW). More formally, given a flow network \(G=(V,E,f)\) and a set of weights \(W = \{w_1,\dots,w_\ell\} \subseteq \mathbb{Z}^+\), find a minimum-sized set of \(s\)-\(t\) paths \(\mathcal{P}= (P_1,\ldots,P_k)\) with associated weights \((\tilde{w}_1,\ldots,\tilde{w}_k)\), with each \(\tilde{w}_i \in W\), such that \(\sum_{i \in \{1,\dots,k\}} P_i(u,v)\tilde{w}_i = f_{uv}\), \(\forall (u,v) \in E\).

This problem has a smaller search space than MFD, because the weights are restricted to \(W\), and thus can potentially admit faster solvers. Moreover, if the weights of an optimal solution are already known, then an optimal solution to MFDW is also an optimal solution to MFD, and this problem can be used as potentially faster solver. For instance, Kloster et al. [16] have observed that path weights can be found in the flow in graphs used for RNA sequence reconstruction.

In the rest of this section we show that knowing such a set \(W\), we are able to formulate a model with substantially fewer variables than the ILP model for MFD, potentially decreasing the runtime of the ILP solver by a substantial amount. This formulation also has the potential to be an alternative efficient heuristic algorithm for MFD, if \(W\) is well chosen.

Given the set of possible weights, our ILP model will able to answer the following question: In an optimal solution (that minimizes the number of used paths), how many paths of weight \(w_i\) pass through an edge \(e\)? In other words, to decompose a flow \(f\) using the weights in \(W\), we aim to find flows \(X_i\) such that \(f = \sum_{i=1}^{|W|} w_iX_i\), where each \(X_i(e)\) answers that question. The ILP model for MFDW is the following: \[\begin{align} \allowdisplaybreaks &\emph{Minimize} \quad k\mathrel{\vcenter{:}}=\sum_{(s,u)\in E} \sum_{i=1}^{|W|} x_{sui} \nonumber \\ &\emph{Subject to:} \nonumber \\ & f_{uv} = \sum_{i=1}^{|W|} w_i x_{uvi}, \quad && \forall (u,v)\in E,\\ & \sum_{(u,v)\in E} x_{uvi} - \sum_{(v,u)\in E} x_{vui} = 0, \quad &&\forall v\in V\setminus \{s,t\}, \forall i \in \{1,\dots,|W|\},\\ & x_{uvi} \in \mathbb{N}, \quad && \forall (u,v)\in E, \forall i \in \{1,\dots,|W|\}. \end{align}\]

Note that the flow \(X_i(u,v)=x_{uvi}\) can be decomposed trivially into weight \(1\) paths. In addition, the product \(x_{uvi}w_i\) does not require linearization due to \(w_i\) being an input in this problem.

In 8 we further argue that MFDW is strongly NP-hard and that solving MFDW of \(G\) with weights \(\{2^i \mid i=0,\dots,\ceil{\log ||f||}\}\) optimally is a \(\ceil{\log ||f||}+1\)-approximation of MFD. Depending on the application, weight sets other than powers-of-two can give a better approximation or faster runtime. We can for example use powers of \(2^i\) for any \(i \in \mathbb{Z}^+\), which yields an approximation scheme with approximation factor \((2^i-1)\ceil{(\log||f||+1)/i}\). Increasing \(i\), we expect a decreased runtime due to the reduced dimension of the ILP model, but the size of the decomposition increases.

6 Experimental results↩︎

6.0.0.1 Solvers and datasets.

Our implementation of the ILPs uses the GUROBI Python API under default settings and is available at https://github.com/algbio/optimized-fd. We compare the ILPs for MFD, MFDSC and MIFD also to the heuristic algorithms for them by [13] (\(\mathsf{Catfish}\)), by [10] (\(\mathsf{Coaster Heuristic}\)) and by [7] (\(\mathsf{IFDSolver}\)), respectively. For MFDW we used the weight set \(W = \{ 2^i \mid i \leq \log ||f|| \} \cup \{ f(e) \mid e\in E \}\).

We also experimented with the implementation from the \(\mathsf{Catfish}\) solver of the standard greedy algorithm for MFD [21], but since we observed that it performs worse than the \(\mathsf{Catfish}\) heuristic algorithm, we did not include it in the results. We also did not include in the results the \(\mathsf{Toboggan}\) implementation by Kloster et al. [16] of the FPT algorithm for MFD, nor the \(\mathsf{Coaster}\) implementation by Williams et al. [10] of the FPT algorithm for MFDSC since it was already observed in [17] that they do not scale for minimum flow decomposition sizes larger than \(6\) (which we also confirmed experimentally).

The runtimes of all ILPs (except MFDW) include a linear scan in increasing order to find the smallest \(k\) for which there is a flow decomposition in \(k\) paths. As discussed in 3, the size of the maximum weight antichain \(Q\) is a lower bound on \(k\), thus the linear scans for MFD-optimized, MFDSC-optimized and MIFD-optimized start at \(|Q|\). We set up a time limit of up to 30 minutes for each input and each method. Our experiments were performed in a server running Linux with one AMD Ryzen Threadripper PRO 3975WX 32-core CPU with 512 GB RAM.

For the MFD problem, we experimented with the datasets created by Shao and Kingsford [13]. This dataset was created from human transcriptome using the quantification tool Salmon [24] and also contains datasets that were created using Flux-Simulator [25] from human and mouse transcriptomes. First, a sample file from their dataset referred as SRR020730 Salmon (corresponding to the file rnaseq/sparse_quant_SRR020730.graph), followed by the entire archive present in the directories rnaseq/human/ and rnaseq/mouse/, referred as Catfish Human and Catfish Mouse, respectively. From those two directories, only the graphs that had 50 nodes (prior to Y-to-V) were considered. For MFDSC, two datasets were used: an adaptation of the SRR020730 Salmon, where subpaths were generated and another dataset SRR30790 StringTie, which was created by Khan et al. [18] from human RNA-seq reads SRR307903 assembled using the StringTie tool [2]). In both datasets, a limited of four subpaths per test instance were considered (due to performance limitations regarding previous tools used as benchmarks). Finally, for MIFD, we simulate interval flows similar to what was done in [7] and, as also described in [7], it is possible that infeasible instances are generated in the process. In this case, we repeated the process until only feasible instances were generated.

6.0.0.2 Experimental setup.

To show the behaviour of the solvers for graphs of increasingly large MFD, we group the input graphs in ranges based on their MFD size, computed with MFD-optimized. If MFD-optimized does not finish within the time limit on an instance, we exclude it also from all other solvers – note that only the heuristic solvers are faster than MFD-optimized, which are not optimal in general. For each solver, we report the average runtime per range of graphs (column Avg.), and the total runtime in that range on the graphs (column Total), both in seconds, and only for the graphs on which that solver finished within the time limit. Thus, these numbers are an underestimation of the time one would need to run the solver in practice, since it would run for at least 30 minutes on unsolved instances. In column #Solved we list the number of instances on which the solver finished. We captured the runtime with the GNU time tool by separating the graphs into individual files and running the tools separately on each instance.

Column shows the approximation accuracy of the heuristic methods as follows: for each instance, we compute the difference between the number of paths reported by each formulation and the minimum number of paths in an MFD (computed with MFD-optimized). In each table cell, corresponding to a specific range of inputs, we list the sum of these differences, and in parentheses their averages.

6.0.0.3 Discussion.

[tab:FD95standard] illustrates the performance of the solvers for the MFD problem. MFD-optimized improves over MFD-original in all ranges of MFD size. We generally observe that the larger the MFD size, the more significant the runtime improvement. For example, for Catfish Human (\(\min k \geq 16\)), MFD-optimized is 49\(\times\) faster than MFD-original on average and for Catfish Mouse (\(\min k \geq 21\)), MFD-optimized is 70\(\times\) faster than MFD-original on average.

As mentioned above, in the reported runtimes of the solvers we are not including instances which did not finish within the time limit. Thus, in practice MFD-optimized has even larger speedups compared to MFD-original, since MFD-optimized solves more instances within the time limit than MFD-original, for Catfish Human and Mouse this happening already for \(\min k \geq 6\). In SRR020730 Salmon MFD-optimized solves all instances, and on the other two datasets MFD-optimized solves 22 and 46 more instances, respectively.

Among the two heuristic solvers, \(\mathsf{Catfish}\) finishes on all instances in a few seconds in total, whereas MFDW is slower than \(\mathsf{Catfish}\), but still running in under one second on average per graph and faster than MFD-optimized in almost all cases. When comparing the approximation accuracy, note that the instances of practical interest are those for which the MFD size is more than 10, since for smaller sizes MFD-optimized already returns an optimal solution seconds in under 2 on average. On such instances, MFDW is more accurate than \(\mathsf{Catfish}\), and interestingly, the larger the MFD size, the more accurate MFD-optimized becomes (for \(\min k \geq 21\) being fully accurate on all datasets). One reason might be that in more complex graphs, the weights of an optimal MFD appear among the flow values of the edges, which we added to the set \(W\).

In [tab:MFDSC95results], we show the results for the MFDSC problem. We observe a similar and even more pronounced trend as for the MFD solvers, where the larger the MFDSC size, the larger the improvement of our optimized ILP. For example, for dataset SRR020730 Salmon, MFDSC-optimized is 11\(\times\), 58\(\times\) and 54\(\times\) faster than MFDSC-original for \(\min k\) in ranges 11-15, 16-20, 21-max, respectively. For both datasets, MFDSC-optimized manages to solve one more instance than MFDSC-original, and in SRR020730 Salmon it solves all instances. Moreover, the total running time of \(\mathsf{Coaster Heuristic}\) is more than half of the running time of MFDSC-optimized, while not giving optimal solutions (which are also less accurate on average than the heuristic solvers for MFD).

In [tab:MIFD95results], we show the results for MIFD problem and again observe a significant speedup compared to MIFD-original on complex instances. For example, for \(\min k\) in 11-15, MIFD-optimized is up to 15\(\times\) faster on SRR020730 Salmon, up to 123\(\times\) faster on Catfish Human and up to 47\(\times\) faster on Catfish Mouse. On the latter two datasets, MIFD-optimized drastically runs in less than \(1\) second on average. In addition, although \(\mathsf{IFDSolver}\) has, on average, a better runtime, its difference compared to MIFD-optimized is negligible.

7 Conclusion↩︎

In this paper, we proposed optimizing the ILP formulations for MFD using the notion of safe paths for all flow decompositions. Since safe paths cannot be simply be removed from the graph, we observed that we can use a set of independent safe path to suitably fix a large number of ILP variables corresponding to their edges. Combined with the Y-to-V reduction, and the first usage of an antichain lower bound in a solver to detect the optimality of a heuristic solution, this resulted in a significantly faster MFD solver (up to 91\(\times\) on the harder instances). We also developed an ILP that can work as a heuristic for MFD, running in under 1 second on average, and more accurate than the state-of-the-art heuristic \(\mathsf{Catfish}\) on graphs of practical interest. We also showed that these optimizations can be applied to two MFD variants of practical interest (by also adapting the Y-to-V reduction the first time for these variants), resulting in even bigger speedups.

Future research encompasses extending such improvements to the ILPs for MFD in general graphs with cycles [26]. Another point of improvement is a further search space reduction through further optimizations on the application of safety, for example, by answering the following question: Is it tractably possible to find subpaths that are part of all \(\alpha\)-approximations of MFD, where \(\alpha\) can either be a constant or depend on the input size? Such \(\alpha\)-safe paths (following the notion of \(c\)-essential vertices by Bumpus et al. [19]) would be longer the closer \(\alpha\) is to \(1\), while safe paths in this paper would refer to \(\infty\)-safe paths. For the problem variants MFDSC and MIFD, can we find all \(\infty\)-safe paths in polynomial time?

8 Additional proofs and figures↩︎

This section contains the proofs and figures missing from the main matter, and some additional results for the MFDW problem.

When dealing with instances in which there is no known small set \(W\) of which the path weights of an optimal solution of MFD come from, we can also use MFDW to approximate the solution of MFD.

9 ILP formulations↩︎

9.1 MFD-original↩︎

\[\label{eq:mod:ST} \begin{align} & {\forall i \in \{1, \ldots, k\}:} \nonumber \\ & \quad \sum_{(s,v) \in E} x_{svi} = 1 \\ & \quad \sum_{(u,t) \in E} x_{uti} = 1 \\ & \quad \sum_{(u,v) \in E} x_{uvi} - \sum_{(v,w) \in E} x_{vwi} = 0 && \forall v \in V \setminus \{s, t\},\\ & f_{uv} = \sum_{i \in \{1,\dots,k\}} \phi_{uvi} && \forall (u,v) \in E, \\ & {\forall (u,v) \in E, \forall i \in \{1,\dots,k\}:} \nonumber \\ & \quad \phi_{uvi} \leq f_\text{max} x_{uvi}\\ & \quad \phi_{uvi} \leq w_i \\ & \quad \phi_{uvi} \geq w_i - (1-x_{uvi})f_\text{max} \\ & w_i \in \mathbb{Z}^+ && \forall i \in \{1, \ldots, k\}, \\ & x_{uvi} \in \{0,1\} && \forall (u,v) \in E, \forall i \in \{1,\dots,k\},\\ & \phi_{uvi} \in \mathbb{N}&& \forall (u,v) \in E, \forall i \in \{1,\dots,k\}. \end{align}\tag{10}\]

9.2 MFD-optimized↩︎

\[\allowdisplaybreaks \label{mod:OP} \begin{align} x_{uvi} = 1, \quad &&\forall (u,v)\in P_i, \forall P_i\in \mathcal{P}(Q) \\ w_{i+1} \leq w_i, \quad && \forall i \in \{ |Q|+1,\dots,k-1\} \end{align}\tag{11}\]

9.3 MFDW↩︎

\[\allowdisplaybreaks \label{mod:IA} \begin{align} & \emph{Minimize}\quad \sum_{(s,u)\in E} \sum_{i=1}^{|W|} x_{sui}\nonumber\\ & \emph{Subject to:} \nonumber \\ & f_{uv} = \sum_{i=1}^{|W|} w_i x_{uvi}, \quad && \forall (u,v)\in E\\ & \sum_{(u,v)\in E} x_{uvi} - \sum_{(v,u)\in E} x_{vui} = 0, \quad && \forall v\in V\setminus \{s, t\}, \forall i \in \{1,\dots,|W|\} \\ & x_{uvi} \in \mathbb{N}, \quad && \forall (u,v)\in E, \forall i \in \{1,\dots,|W|\} \end{align}\tag{12}\]

10 Maximum length independent safe paths↩︎

The approach used in 3 to calculate independent safe paths of a flow network \((G,f)\) does not return them of maximum length, but works only as a heuristic. In order to find independent safe paths of maximum length in polynomial time, one can, as in the heuristic approach, reduce to maximum weight antichains, on the following dependency graph \(\mathcal{D}=(V(\mathcal{D}),E(\mathcal{D}))\): nodes represent safe paths, with given weights of their length. A directed edge \((u,v)\) is added if there exists a path in \(G\) which first traverses both safe paths represented by \(u\) and \(v\), entering the safe path of \(u\) first and then entering the safe path of \(v\).

References↩︎