# Efficient Sparse Flow Decomposition Methods for RNA Multi-Assembly

Mathieu Besançon

► To cite this version:

Mathieu Besançon. Efficient Sparse Flow Decomposition Methods for RNA Multi-Assembly. ECAI 2025 - 28th European Conference on Artificial Intelligence, Oct 2025, Bologna, Italy. pp.1-8. hal-05173160

HAL Id: hal-05173160

https://hal.science/hal-05173160v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient Sparse Flow Decomposition Methods for RNA Multi-Assembly

Mathieu Besançon<sup>a,*</sup>

<sup>a</sup>Université Grenoble Alpes, Inria, LIG, CNRS
ORCID (Mathieu Besançon): https://orcid.org/0000-0002-6284-3033

**Abstract.** Decomposing a flow on a Directed Acyclic Graph (DAG) into a weighted sum of a small number of paths is an essential task in operations research and bioinformatics. This problem, referred to as Sparse Flow Decomposition (SFD), has gained significant interest, in particular for its application in RNA transcript multi-assembly, the identification of the multiple transcripts corresponding to a given gene and their relative abundance. Several recent approaches cast SFD variants as integer optimization problems, motivated by the NP-hardness of the formulations they consider. We propose an alternative formulation of SFD as a data fitting problem on the conic hull of the flow polytope. By reformulating the problem on the flow polytope for compactness and solving it using specific variants of the Frank-Wolfe algorithm, we obtain a method converging rapidly to the minimizer of the chosen loss function while producing a parsimonious decomposition. Our approach subsumes previous formulations of SFD with exact and inexact flows and can model different priors on the error distributions. Computational experiments show that our method outperforms recent integer optimization approaches in runtime, but is also highly competitive in terms of reconstruction of the underlying transcripts, despite not explicitly minimizing the solution cardinality.

## 1 Introduction

We consider the problem of decomposing flows into a weighted sum of paths in a Directed Acyclic Graph (DAG) $G = (V, E)$, specifically seeking a small support for the weights. We refer to the problem as the Sparse Flow Decomposition (SFD) problem and will state various mathematical formulations of SFD from the literature. We define a *pseudo-flow* as a function $E \to \mathbb{R}_+$ on the edges of the graph and equivalently, as a vector $\mathbf{r} \in \mathbb{R}_+^{|E|}$. A flow is a pseudo-flow respecting a conservation constraint on all nodes except the source $s$ and target $t$, i.e., the sum of flows coming into a node is equal to the sum of the flow going out of it. This problem has been studied intensively in the last years, in particular thanks to the key application of multi-assembly for RNA transcripts in bioinformatics. In this application, the DAG corresponds to a given gene *splice graph* in which nodes represent exons (RNA sections encoding information on the gene) to which an artificial source and sink are added, and edges between two exons correspond to reads with one exon following the other. The goal is to identify *transcripts* which are sequences of exons corresponding to one modality of expression of the gene. A transcript is equivalent to a path in the graph, and the weight associated to that

transcript corresponds to the relative abundance of that transcript to express the gene. We illustrate the problem setup in Figure 1.

<figure><img src="image_2.png" /><figcaption>(a) Example splice graph with exons *a* to *d* and associated flow.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) A decomposition of the flow into two paths $s-a-c-d-t$ (orange) and $s-b-c-t$ (purple) of weight one and two respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) A decomposition of the flow into three paths $s-a-c-t$ (green), $s-b-c-t$ (blue), and $s-b-c-d-t$ (red), all of weight one.</figcaption></figure>

**Figure 1:** Illustration of the problem setup. The problem input illustrated in Figure 1a is provided as a directed acyclic graph with flows assigned to edges. All non-terminal nodes *a...d* correspond to exons. Figure 1b and 1c show two flow decompositions of the input, producing two disjoint sets of paths.

The problem has received scrutiny from the bioinformatics, algorithm design, and optimization community in the last decade due to the high relevance for transcriptomics and the growing availability of corresponding data. We highlight the lines of research connected to our work and more generally SFD and that are compatible with inexact flow models, i.e., when the input pseudo-flow data are contaminated with errors and may not form a flow.

## Related work on Flow Decomposition

Most early work focused on sparse fitting approaches and designing

<footnote>* Corresponding Author. Email: mathieu.besancon@inria.fr.</footnote>
efficient two-step algorithms. In a first stream of work, convex optimization models designed for sparse regression were leveraged e.g., in [28, 27] to fit sparse predictive models on the weights associated with paths. One major drawback of such an approach is the need to first enumerate the exponential number of paths before fitting the sparse model. In order to avoid explicitly working on the exponential number of paths, [36] design a two-step approach first using a minimum cost flow model to produce a flow fitting the data under a given loss, followed by a heuristic to decompose that flow into a weighted sum of paths. In [37], the authors consider the same problem with the restriction that only a few paths can be used to decompose the flow. Indeed, the sparsity in the number of paths corresponds to a property observed on real splice graphs and corresponding transcript abundances. They also show that several formulations of sparse flow decomposition asking for a given upper bound on the number of allowed paths are NP-hard. Similarly, the method proposed in [1] is based on a sparse statistical model for the path weights and encodes the problem as a network flow problem producing a flow which is by design in the convex hull of paths (since the graph is acyclic). Unlike previous approaches, they model the flow values on the arcs as Poisson random variables and optimize the corresponding log-likelihood instead of the least-square loss. However, they then require a greedy heuristic to decompose their flow into the corresponding weighted set of paths, which for exact flows is as hard as the original problem, greatly hindering the practicality of the approach. By casting their formulation and the least-square minimization from [36] in a form amenable to FW algorithms, we revisit these formulations and show their merit when combined with a sparsity-inducing algorithmic approach instead of a modified formulation.

In another approach to handle errors in the RNA-Seq data, the authors of [40] propose a formulation relaxing the constraint that the weighted sum of paths exactly matches a flow and instead construct an interval of flow values for each edge, and then design a custom heuristic to handle these flow intervals. One major drawback of this approach is already requiring the production of these intervals, and then seeking a feasible solution in this interval. Even if the interval contains the underlying flow, it is not given that the minimum-cardinality solution respecting these bounds will fit this ground truth flow.

More recently, the continuous progress in mixed-integer optimization methods and solvers [25] allowed considering explicit mixed-integer formulations of sparse flow decomposition problems for instance scales that would not have been tractable a decade ago. The first integer optimization model was proposed in [15], also accommodating the inexact version of the problem from [40] in which the weighted sum of paths must lie within some distance of the input flow values. The formulation is based on a quadratic number of variables encoding paths using flow conservation constraints, the flow expressed as the sum of the paths and the weight variables. The authors also adapt the least-square formulation from [37], resulting in a mixed-integer quadratic (convex) optimization problem which is theoretically hard and computationally harder than mixed-integer linear optimization problems. As an alternative formulation handling inexact flows, [13] proposed to lift the uncertainty handling from edges to paths in a mixed-integer formulation. The rationale for handling errors at the path level is that edges appearing in multiple paths are more prone to read errors in the flow value.

An experimental assessment of the minimality assumption was conducted in [24], showing that even though most ground truth solutions are of minimum support (i.e. use the smallest possible number of $s-t$ paths), this is not the case for all of them. This observation

naturally leads to seeking sparsity of the solution in terms of number of paths, rather than its minimality. Furthermore, the arguments for computing a set of paths of minimum cardinality do not necessarily hold in the realistic case where the input data are contaminated with errors. Finally, we note a recent line of work detecting paths that are required in optimal solutions of the decomposition [19, 32] and that can be exploited in integer optimization formulations.

## Frank-Wolfe algorithms

Frank-Wolfe (FW) or conditional gradient algorithms [17, 26] optimize differentiable functions over compact convex sets. They have benefitted from a strong interest in the last decade, in particular thanks to their advantages for large-scale machine learning applications [22], including their low cost per iteration and possible exploitation of the structure of the constraint set. At its core, FW produces iterates as convex combinations of a small number of extreme points of the feasible set, while only requiring that the function is differentiable (and Lipschitz-smooth in typical cases) and equipped with zeroth and first-order oracles, and that the feasible region can be accessed through a linear minimization oracle (LMO), i.e., an algorithm which, given a direction, computes an extreme point of the feasible region minimizing its inner product with the direction. On a generic polytope, this LMO can be implemented through linear optimization but on many structured sets, specialized algorithms can implement the LMO without forming the linear problem constraints explicitly. Finally, we highlight that FW algorithms produce a so-called Frank-Wolfe gap as a by-product at every iteration, which upper-bounds the unknown primal gap. We refer interested readers to the recent surveys [8, 9] for applications and important results on FW algorithms.

## Contributions

Our contributions are the following.

1. We propose a formulation of SFD on the flow polytope optimizing either the least-square error or Poisson log-likelihood and a corresponding solution approach based on Frank-Wolfe algorithms. Importantly, the resulting optimization problems are convex and formulated over the convex hull of $s-t$ paths. We provide theoretical justification and computational evidence for sparsity of the solutions obtained by our method.

2. We establish the convergence rate of our method and evaluate its cost per iteration, showing linear convergence under both the least-square and Poisson log-likelihood losses and despite the lack of strong convexity in both cases.

3. We evaluate our method compared to recent integer optimization approaches on reference multi-assembly datasets with and without error contamination in the flow data. The results show that our approach not only dominates the integer formulations in runtime, but also produces high-quality solutions in terms of reconstruction error, path identification, and solution sparsity.

## Notation and terminology

Vectors are denoted with bold small letters, scalars with standard small letters, and matrices with capital bold letters. For $n$, $[n] := \{1...n\}$. We use $|a|$ for a rounding of $a$ to the closest integer. $\Delta_n$ denotes the standard simplex. For $u$ a node in the graph, $\delta_u^{in}$, $\delta_u^{out}$ will denote the set of edges with that node as destination, origin respectively. When unspecified, the default norm $\|\cdot\|$ is the Euclidean norm in the appropriate vector space. The function $\log(\cdot)$ denotes the natural logarithm. For a function $f$, we denote with $D^k f[\mathbf{u}_1...\mathbf{u}_k]$ the $k$-th directional derivative of $f$ along directions $\mathbf{u}_1...\mathbf{u}_k$.
## 2 Sparse Flow Decompositions via Frank-Wolfe Approaches

In this section, we formulate SFD under the least-square and Poisson models, transforming both into constrained minimization problems over polytopes.

### 2.1 Least-square Problem Formulation

The flow decomposition problem can be viewed as finding a sparse approximation of a given (pseudo-)flow $\mathbf{r}$ on a DAG given by a conic combination of weighted paths:

$$ \begin{array}{ll} \min_{\mathbf{x},\mathbf{v},\mathbf{w}} & \frac{1}{2} \|\mathbf{x} - \mathbf{r}\|_2^2 \\ \text{s.t.} & \sum_{s=1}^{|E|-|V|+2} w_s \mathbf{v}_s = \mathbf{x} \\ & \mathbf{v}_s \in \mathcal{X}, \|\mathbf{w}\|_0 \le k, w_s \in \mathbb{Z}_+ \forall s \in \llbracket k \rrbracket, \end{array} $$

where $\mathcal{X}$ is the set of $s-t$ paths, and $k$ is an upper bound on the number of paths to use. If $\mathbf{r}$ is a flow, then the optimal solution has a zero objective value. Furthermore, the flow $\mathbf{x}$ can be decomposed as the weighted sum of at most $|E| - |V| + 2$ paths [39], i.e., the cyclomatic number of the DAG plus one, hence the upper bound on the number of paths. The path weights $w$ represents the abundance of each transcript, the integrality constraints on $\mathbf{w}$ are a modeling choice depending on prior knowledge on the input data. The formulation thus seeks the flow $\mathbf{x}$ that is the closest to $\mathbf{r}$ in the Euclidean sense and that can be formed as an integer conic combination of $k$ paths. From a statistical perspective, this modeling choice follows naturally, e.g., from a Gaussian assumption on the errors polluting the flow data. The formulation captures several previous approaches [37, 40, 15] and unifies both exact and inexact flow decompositions. An aspect of importance for RNA reconstruction and other applications of SFD is producing a sparse decomposition, i.e., a decomposition using a small number of paths. This consideration motivated several lines of work to formulate the problem objective solely on the solution sparsity, i.e. minimizing the weight support subject to fitness to the data. An objective function minimizing the number of paths leads to NP-hard versions of the problem, for which the natural solution method is based on mixed-integer formulations. Instead, we propose an algorithmic approach to sparsity, leveraging methods based on the Frank-Wolfe algorithm which naturally produce iterates as convex combinations of a small number of vertices. We reformulate and relax the problem to:

$$ \begin{array}{ll} \min_{\mathbf{x},\mathbf{v},\mathbf{w},\tau} & \frac{1}{2} \|\mathbf{x} - \tau\mathbf{r}\|_2^2 \\ \text{s.t.} & \sum_{k=1}^s w_s \mathbf{v}_k = \mathbf{x} \\ & \mathbf{v}_s \in \mathcal{X}, w_s \in \Delta_n \forall s \in \llbracket k \rrbracket \\ & \tau \ge 0. \end{array} $$

The additional $\tau$ variable scales the flow to the appropriate magnitude to be contained in the convex hull of paths instead of its conic hull. We can then expand its expression by minimizing the objective w.r.t. $\tau$:

$$ \tau_{\min} \in \underset{\tau}{\operatorname{argmin}} \|\mathbf{x} - \tau\mathbf{r}\|^2 \Leftrightarrow \tau_{\min} = \frac{\langle \mathbf{r}, \mathbf{x} \rangle}{\|\mathbf{r}\|^2}, $$

leading to our final least-square formulation:

$$ \min_{\mathbf{x} \in \operatorname{conv}(\mathcal{X})} \frac{1}{2} \left\| \mathbf{x} - \frac{\langle \mathbf{x}, \mathbf{r} \rangle}{\|\mathbf{r}\|^2} \mathbf{r} \right\|^2. $$

Importantly, the problem can now be tackled efficiently in a FW setting, since the constraint set admits an LMO implementable through a shortest path computation on the DAG. One requirement missing from this formulation is that $\mathbf{x}$ is formed as a sparse convex combination of paths. Instead of enforcing this in the problem formulation, we will ensure sparsity of the iterates and final solution through the algorithm leveraged for the solution process, namely Frank-Wolfe methods and in particular active set-based methods. The sparsity in the number of vertices used to construct the solutions $\mathbf{x}_t$ at any iteration $t$ has been a primary motivation for the strong interest in Frank-Wolfe methods in the last decade [22, 9]. For most practical cases, the empirical sparsity is much better than the upper bounds that could be obtained (at most one vertex per iteration for most Frank-Wolfe variants). More recently, results on active set identification of Frank-Wolfe FW methods provided theoretical evidence and guarantees on the sparsity of some FW methods. The active set identification property introduced in [7] ensures that when using an optimal step size (e.g. through line search), the Away-step Frank-Wolfe algorithm reaches the optimal face of the problem in a finite number of iterations and never leaves it afterwards. This property was extended to the Blended Pairwise Conditional Gradient (BPCG) [38] in [42], in both cases, without the need to assume strong convexity. The result on active set identification comes with a second result which is a non-trivial upper bound on the number of vertices that are required to represent the current iterate. After the finite number of iterations needed to reach the optimal face $\mathcal{F}^*$, [42] establishes a bound of $\dim(\mathcal{F}^*) + 1$ vertices required to form the iterate. Not only can this bound be reached, but the corresponding convex decomposition can be obtained as the basic solution of an auxiliary linear optimization problem. Furthermore, the computational evidence pointed out that BPCG was already sparse enough not to require that additional step and provided sufficient sparsity matching this bound. We present the BPCG algorithm applied to our setting in Algorithm 1. The quantity $g_t$ corresponds to the FW gap, $\mathrm{lm}_{\mathrm{G}}(\cdot)$ computes the shortest path incidence vector with the edge weights as argument. The function $\mathrm{weight}_S(\mathbf{v})$ returns the weight of the given vertex in the decomposition of the iterate.

### 2.2 Optimal Step Size

Numerous step-size strategies have been proposed for Frank-Wolfe algorithms, from function-agnostic step sizes based on the iteration count to algorithms efficiently approximating a line search. For our least-square formulation however, the optimal step size can easily be computed from the problem data as shown in Proposition 1.

**Proposition 1.** At any iteration of the BPCG algorithm applied to the least-square problem, given the current iterate $\mathbf{x}$ and the current direction $\mathbf{d}$, the optimal step size $\gamma^*$ is given by:

$$ \gamma^* = \min \left\{ \frac{\|\mathbf{r}\|^2 \langle \mathbf{x}, \mathbf{d} \rangle - \langle \mathbf{d}, \mathbf{r} \rangle \langle \mathbf{x}, \mathbf{r} \rangle}{\|\mathbf{d}\|^2 \|\mathbf{r}\|^2 - \langle \mathbf{d}, \mathbf{r} \rangle^2}, \gamma_{\max} \right\}. \quad (1) $$

*Proof.* The step size can be derived from the expression of the ob-
**Algorithm 1** Blended pairwise conditional gradient for flow decomposition

**Input:** reference flow $\mathbf{r}$, DAG $G$, loss function $f$

**Output:** Paths and weights $\{\mathbf{v}_k\}_{k \in \llbracket S_T \rrbracket}$, $\{\lambda_k\}_{k \in \llbracket S_T \rrbracket}$.

1: $\mathbf{x}_0 \leftarrow p_G(\mathbf{1})$, $g_0 \leftarrow +\infty$

2: **for** $t \in \llbracket T \rrbracket$ **do**

3: $\mathbf{v}_t \leftarrow \text{lm}_O_G(\nabla f(\mathbf{x}_t))$

4: $\mathbf{a}_t \leftarrow \text{argmax}_{\mathbf{v} \in S_t} \langle \nabla f(\mathbf{x}_t), \mathbf{v} \rangle$

5: $\mathbf{s}_t \leftarrow \text{argmin}_{\mathbf{v} \in S_t} \langle \nabla f(\mathbf{x}_t), \mathbf{v} \rangle$

6: $g_t \leftarrow \langle \nabla f(\mathbf{x}_t), \mathbf{x}_t - \mathbf{v}_t \rangle$

7: **if** $\langle \nabla f(\mathbf{x}_t), \mathbf{a}_t - \mathbf{s}_t \rangle \ge g_t$ **then**

8: $\mathbf{d}_t \leftarrow \mathbf{a}_t - \mathbf{s}_t$

9: $\gamma_{\max} \leftarrow \text{weight}_{S_t}(\mathbf{a}_t)$

10: $\gamma_t \leftarrow \text{argmin}_{\gamma \in [0, \gamma_{\max}]} f(\mathbf{x}_t - \gamma \mathbf{d}_t)$

11: **if** $\gamma_t < \gamma_{\max}$ **then**

12: $S_{t+1} \leftarrow S_t$

13: **else**

14: $S_{t+1} \leftarrow S_t \setminus \{\mathbf{a}_t\}$

15: **end if**

16: **else**

17: $\mathbf{d}_t \leftarrow \mathbf{x}_t - \mathbf{v}_t$

18: $\gamma_t \leftarrow \text{argmin}_{\gamma \in [0,1]} f(\mathbf{x}_t - \gamma \mathbf{d}_t)$

19: **if** $\gamma_t < 1$ **then**

20: $S_{t+1} \leftarrow S_t \cup \{\mathbf{v}_t\}$

21: **else**

22: $S_{t+1} \leftarrow \{\mathbf{v}_t\}$

23: **end if**

24: **end if**

25: $\mathbf{x}_{t+1} \leftarrow \mathbf{x}_t - \gamma_t \mathbf{d}_t$

26: **end for**

27: **return** $(\mathbf{x}_T, S_T)$

jective evaluated at $\mathbf{x} - \gamma\mathbf{d}$:

$$ 
\begin{aligned}
f(\mathbf{x} - \gamma \mathbf{d}) &= \frac{1}{2} \left\| \mathbf{x} - \gamma \mathbf{d} - \frac{\langle \mathbf{x} - \gamma \mathbf{d}, \mathbf{r} \rangle}{\|\mathbf{r}\|^2} \mathbf{r} \right\|^2 \\
&= \frac{1}{2} \left\| \mathbf{x} - \frac{\langle \mathbf{x}, \mathbf{r} \rangle}{\|\mathbf{r}\|^2} \mathbf{r} - \gamma \left( \mathbf{d} - \frac{\langle \mathbf{d}, \mathbf{r} \rangle}{\|\mathbf{r}\|^2} \mathbf{r} \right) \right\|^2 \\
&= \frac{1}{2} \|\mathbf{a} - \gamma \mathbf{b}\|^2,
\end{aligned}
$$

where $\mathbf{a}, \mathbf{b}$ are the appropriate expressions substituted here for con-
ciseness. We can differentiate the loss with respect to $\gamma$, resulting in
the unconstrained minimum $\gamma^* = \frac{\langle \mathbf{a}, \mathbf{b} \rangle}{\langle \mathbf{b}, \mathbf{b} \rangle}$. By expanding the terms $\mathbf{a}$
and $\mathbf{b}$, we have:

$$ \langle \mathbf{a}, \mathbf{b} \rangle = \langle \mathbf{x}, \mathbf{d} \rangle - 2 \frac{\langle \mathbf{d}, \mathbf{r} \rangle \langle \mathbf{x}, \mathbf{r} \rangle}{\|\mathbf{r}\|^2} + \frac{\langle \mathbf{x}, \mathbf{r} \rangle \langle \mathbf{d}, \mathbf{r} \rangle \langle \mathbf{r}, \mathbf{r} \rangle}{\|\mathbf{r}\|^4} $$

$$ \langle \mathbf{b}, \mathbf{b} \rangle = \|\mathbf{d}\|^2 - 2 \frac{\langle \mathbf{d}, \mathbf{r} \rangle^2}{\|\mathbf{r}\|^2} + \frac{\langle \mathbf{d}, \mathbf{r} \rangle^2}{\|\mathbf{r}\|^2} $$

$$ \gamma^* = \frac{\langle \mathbf{a}, \mathbf{b} \rangle}{\langle \mathbf{b}, \mathbf{b} \rangle} = \frac{\|\mathbf{r}\|^2 \langle \mathbf{x}, \mathbf{d} \rangle - 2 \langle \mathbf{d}, \mathbf{r} \rangle \langle \mathbf{x}, \mathbf{r} \rangle + \langle \mathbf{x}, \mathbf{r} \rangle \langle \mathbf{d}, \mathbf{r} \rangle}{\|\mathbf{d}\|^2 \|\mathbf{r}\|^2 - \langle \mathbf{d}, \mathbf{r} \rangle^2} $$

resulting in Equation (1) with the appropriate upper bound $\gamma_{\max}$. $\square$

## 2.3 Poisson Regression Formulation

In this section, we revisit the Poisson regression model developed in
[1]. In this formulation, the flow passing through each exon is mod-
eled as a Poisson random variable with a mean given by the sum of

flows passing through that node in the input data. The log-likelihood
minimization problem is expressed as:

$$ \min_{\mathbf{x} \in \text{cone}(\mathcal{X})} \sum_{u \in V} \left[ \sum_{e \in \delta_u^{\text{in}}} x_e - \left( \sum_{e \in \delta_u^{\text{in}}} r_e \right) \log \left( \sum_{e \in \delta_u^{\text{in}}} x_e \right) \right]. $$

This formulation computes the flow of maximum likelihood, which is a sum of weighted paths with nonnegative weights by the flow decomposition theorem applied to a DAG. The optimal rescaling technique from the least-square formulation cannot be applied here without losing convexity of the loss function. In order to formulate an equivalent compact set, we therefore replace cone($\mathcal{X}$) with

$$ \bar{\mathcal{X}} = \{ \|\mathbf{r}\|_{\infty} \mathbf{x} : \mathbf{x} \in \mathcal{X} \cup \mathbf{0} \}, $$

the convex hull of paths that can be scaled up to at most the maximum flow on the data. For a given direction $\mathbf{g}$, the corresponding LMO corresponds to 1) computing the shortest path $\mathbf{v}$ with edge lengths given by $\mathbf{g}$ and 2) if $\langle \mathbf{v}, \mathbf{g} \rangle \le 0$, returning that vertex, otherwise returning the origin as the minimizer. This ensures that the formulation is applicable to FW.

## 3 Convergence Analysis

The least-square objective function is not strongly nor strictly convex, its Hessian has one zero eigenvalue. We will however show it presents a quadratic growth property [23]:

**Definition 1** (Quadratic growth property). Let $f$ be a closed proper convex function defined over a compact $\mathcal{X}$, $\text{dist}_*(\cdot)$ the distance to its set of minimizers and $f^*$ its minimal value. Then it is said to satisfy a quadratic growth property with constant $\mu > 0$ if

$$ f(\mathbf{x}) - f^* \geq \mu \text{ dist}_*(\mathbf{x}) \quad \forall \mathbf{x} \in \mathcal{X}. $$

**Proposition 2.** The least-square objective function $f$ respects the quadratic growth condition with $\mu = \frac{1}{2}$.

**Proof.** From the expression of the gradient, we can deduce that unconstrained optima are of the form:

$$ X^* := \{ \mathbf{x}^* \in \mathbb{R}_+^{|E|} : \mathbf{x}^* = \alpha \mathbf{r}, \alpha \in \mathbb{R}_+ \}, $$

which implies that the projection $\text{proj}_*(\mathbf{x})$ of a point $\mathbf{x}$ onto the set of optimizers, and the distance of that point to the set of optimizers $\text{dist}_*(\mathbf{x})$ are respectively given by:

$$ \text{proj}_*(\mathbf{x}) = \frac{\langle \mathbf{r}, \mathbf{x} \rangle}{\|\mathbf{r}\|^2} \mathbf{r}, \quad \text{dist}_*(\mathbf{x}) = \left\| \mathbf{x} - \frac{\langle \mathbf{r}, \mathbf{x} \rangle}{\|\mathbf{r}\|^2} \mathbf{r} \right\|. $$

The objective function is precisely half of the squared distance, resulting in $\mu = \frac{1}{2}$ with the notation of Definition 1. $\square$

The BPCG algorithm thus converges linearly for the least-square problem, based on [41, Theorem 3.6] and Proposition 2, since the quadratic growth condition is a special case of sharpness, and our function is Lipschitz-smooth.

Unlike the least-square formulation, the Poisson objective function is not Lipschitz-smooth, the classic FW convergence results hence do not apply. However, it is three times differentiable on its domain and self-concordant. FW with some step-size strategies have been shown to converge on (generalized) self-concordant functions in [16] and [11] at the usual $O(1/t)$ rate and at a linear rate in specific settings
(including strongly convex functions) which do not match ours. Together with the fact that the feasible set is a polytope, [44] establishes linear convergence of Away-step Frank-Wolfe without requiring strong convexity if the objective is the composition of a logarithmically homogeneous self-concordant barrier (LHSCB) function (see Definition 2) with an affine map. Their result is extended in [20] to the BPCG algorithm we apply here, which typically produces sparser iterate than Away-step FW. We present below the definition of a LHSCB function and apply it to our objective function.

**Definition 2.** A convex function $g: \mathbb{R}^n \to \mathbb{R}$ that is three-times differentiable on its domain is a $\theta$-logarithmically homogeneous self-concordant barrier for a proper (closed, convex, pointed) cone $K$ iff:

1. It is self-concordant: $|D^3g(\mathbf{z})[\mathbf{u}, \mathbf{u}, \mathbf{u}]| \le 2(D^2g(\mathbf{z})[\mathbf{u}, \mathbf{u}])^{\frac{3}{2}}$.

2. It is a barrier for $K$: for any sequence $\{\mathbf{z}_k\}_{k \ge 1}$ such that $\mathbf{z}_k \to \text{boundary}(K)$, $g(\mathbf{z}_k) \to \infty$.

3. It is $\theta$-logarithmically homogeneous: $g(\alpha\mathbf{z}) = g(\mathbf{z}) - \theta \log(\alpha)$.

We show that their framework is applicable to the Poisson objective in Proposition 3.

**Proposition 3.** The Poisson objective function can be written as:

$$f(\mathbf{x}) = h(A\mathbf{x}) + \langle \mathbf{b}, \mathbf{x} \rangle, \quad (2)$$

where $A: \mathbb{R}^{|E|} \to \mathbb{R}^{|V|}$ is a linear operator, $h: \mathbb{R}^{|V|} \to \mathbb{R}$ is a logarithmically homogeneous self-concordant barrier function and $\mathbf{b} \in \mathbb{R}^{|E|}$.

*Proof.* The expression of the objective as Equation (2) follows from the following elements

$$\begin{align*}
(A\mathbf{x})_u &= \sum_{e \in \delta_u^{\text{in}}} x_e \\
h(\mathbf{z}) &= -\sum_{u \in V} \left( \sum_{e \in \delta_u^{\text{in}}} r_e \right) \log z_u \\
\langle \mathbf{b}, \mathbf{x} \rangle &= \sum_{u \in V} \sum_{e \in \delta_u^{\text{in}}} x_e.
\end{align*}$$

The function $h$ is self-concordant as the weighted sum of self-concordant functions with nonnegative weights [35, Proposition 1]. It is also $\theta$-logarithmically-homogeneous with

$$\theta = \left( \sum_{e \in \delta_u^{\text{in}}} r_e \right) \text{ since } h(\alpha\mathbf{z}) = h(\mathbf{z}) - \sum_{u \in V} \left( \sum_{e \in \delta_u^{\text{in}}} r_e \right) \log(\alpha).$$

Finally, $h$ is a log-barrier for $\mathbb{R}_+^{|V|}$ since for any point $\bar{z} \in \text{boundary}(\mathbb{R}_+^{|V|})$, a sequence $\{\mathbf{z}_k\}_{k \ge 1}$ of componentwise positive vector such that $\lim_{k \to \infty} \mathbf{z}_k = \bar{z}$ results in $\lim_{k \to +\infty} h(\mathbf{z}_k) = +\infty$. $\square$

Equipped with the result of Proposition 3, we can use the convergence guarantees of [44, 20] to ensure that the BPCG algorithm we are using achieves a linear convergence rate

$$f(\mathbf{x}_t) - f^* \le \mathcal{O}(\exp(-ct)),$$

matching the empirical rate of our computational experiments.

## 3.1 Early Termination

First-order methods are known for their good scalability due to a low cost per iteration, although they can be hindered by a high number of iterations compared, e.g., to interior points. We propose an early stopping criterion for the least-square loss to reduce the number of iterations when the current decomposition is converging towards the scaled flow. For a given solution $\mathbf{x}$, we can compute its optimal scaling $\alpha$ by minimizing over $\alpha \ge 0$ the least-square error: $\|\alpha\mathbf{x} - \mathbf{r}\|^2$. Note that this scales up the solution $\mathbf{x}$ instead of scaling down the flow “r” in order to obtain an integer solution in the conic hull. We did not optimize this function directly since it would result in a non-convex objective. We can derive the optimal $\alpha^* = \langle \mathbf{x}, \mathbf{r} \rangle / \|\mathbf{x}\|^2$ and compute conic weights rounded to the closest integer from the current active set weights $S: \{\mu_k\}_{k \in 1...|S|} = \{|\alpha^*\lambda_k|\}_{k \in 1...|S|}$. At any iteration, a simple test can be performed to check whether the resulting decomposition exactly matches the original flow $\mathbf{r}$, in which case we can stop the algorithm.

## 3.2 Iteration Cost

We break down the cost of iterations of Algorithm 1 from individual components. A crude upper bound on the cost of the linear minimization oracle is $\mathcal{O}(|V||E|)$ provided by the Bellman-Ford algorithm performing a single-source shortest path with negative edge lengths. A finer bound is provided by the Goldberg-Radzik algorithm [18] which obtains a $\mathcal{O}(|V| + |E|)$ runtime with a modification proposed in [12]. Function and gradient evaluations, as well as the exact step size computation can all be performed in $\mathcal{O}(|E|)$. The last operation that could dominate the iteration cost is the inner product search over the active set performed in Line 4 and 5 of Algorithm 1. Indeed, in the worst case, the algorithm would add one vertex to the active set per iteration, resulting in a cost of the active set search $\mathcal{O}(T|E|)$ for the last iterations. However, this represents a worst-case that is not tight in the light of the active set identification and associated bounds from [7, 42] presented in Subsection 2.1, with the number of vertices after a finite number of iterations being bounded by the dimension of the optimal face plus one. For practical purposes of the SFD application, the cost of the active set search has not been limiting in the computational experiments. Future work could consider applying recent advancements in inner-product data structures such as the one presented in [43, 34] to derive a cost of the search almost linear in the number of vertices and dimension.

# 4 Computational Experiments

In this section, we evaluate the performance of our proposed FW approach compared to recent integer-based methods for multi-assembly problems with and without error. The source code for all experiments is available at [5] and archived at [4].

## Experimental Setup

We perform all our computations in Julia 1.11. The mixed-integer problems are modeled with JuMP 1.23 [29] and solved with SCIP 9.2 [6]. We load the DAGs into Graphs.jl and use its implementation of Bellman-Ford for shortest paths. We compare the solution quality to that of the integer optimization formulations from [15] and [13]. We use the BPCG implementation present in FrankWolfe.jl 0.4 [2, 3], noted FW, the same algorithm with early termination from Subsection 3.1, noted FW-C, the Poisson loss optimized with FW noted as FW-P, the integer optimization model from [15] noted IP, and the robust integer optimization version noted IP-R. FW algorithms
are limited to 5000 iterations. All methods are restricted to 1800 seconds as a time limit. Experiments are performed on a cluster with all nodes equipped with Intel Xeon Gold 6338 2GHz CPUs and 512GB of RAM.

## Instance Data

We use the dataset compiled in [14] which contains data for four species summarized in Table 1. These datasets notably include the ground truth transcripts and their abundance, i.e., the paths and corresponding weights and can thus be used to assess the ability of the various methods to recover the underlying structure.

<table><caption>**Table 1:** Statistics on the inexact dataset. Species are indicated by their first letter, human, zebrafish, salmon, mouse. Non-trivial graphs refer to graphs with an upper bound $k = |E|-|V|+2$ on the number of paths being greater than one.</caption><thead><tr><th>Species</th><th>h</th><th>z</th><th>s</th><th>m</th><th>total</th></tr></thead><tbody><tr><td># DAGs</td><td>11783</td><td>15664</td><td>40870</td><td>13122</td><td>81439</td></tr><tr><td>% non-triv.</td><td>0.45</td><td>0.29</td><td>0.36</td><td>0.36</td><td>0.36</td></tr></tbody></table>

## 4.1 Reconstruction of Inexact Flows

We apply the FW-based approaches, the integer optimization formulation and the robust integer formulation to the inexact flow instances from [14]. We quantify the reconstruction quality with two metrics, the *path error* and *flow error*. The path error is computed as the number of paths in the solution with a weight different from the ground truth decomposition. The flow error is the Euclidean distance between the true flow and the reconstructed one. The results are summarized in Table 2.

We observe that despite producing sparser solutions on average, the integer optimization model performs worse than the two FW-based approaches in terms of path and flow error. The robust optimization model yields the best path error performance at the cost of an increased runtime. Surprisingly, the robust integer model is not costlier than the original integer model despite the increased number of variables and constraints. We first analyze the solution qualities on the different metrics. The differences in runtime are analyzed in more details below. The FW-C method performs slightly better than FW on both error metrics, meaning our specific early termination criterion from Subsection 3.1 yields a better solution than reaching the least-square optimum. Note that this is however dependent on the assumption that the underlying flow is formed from integer weights of the individual paths. Importantly, we observe that despite producing sparser solutions, the IP model rarely results in the best reconstruction in terms of path and flow errors. This may imply that least-cardinality is not sufficient as a solution concept to reconstruct transcripts from data, and a statistical approach should be prefer to model the data-generating process, with sparsity being handled through the algorithmic process instead of the formulation. When considering all instances for which a primal solution was returned, the flow reconstructed by FW methods was superior to IP-R. The runtime distribution of the different methods are presented in Figure 2.

The FW methods applied to the least-square loss clearly outperform all others in runtime. There are 919 instances on which FW-P reaches the iteration limit due to numerical instabilities which we further analyze below. The methods IP, IP-R reach the time limit for 3829, resp. 119 instances. The FW-C method terminates faster for small instances than FW, showing that exploiting the weight integral-

<figure><img src="image_5.png" /><figcaption>**Figure 2:** Runtime distribution of the different methods on all inexact instances. IP is the integer optimization formulation with intervals, IP-R is the robust optimization formulation. The prefix S- indicates the time to the best found solution. FW is the blended pairwise Frank-Wolfe method, FW-C includes the early termination criterion based on the flow integrality assumption. FW-P is the Poisson regression problem.</figcaption></figure>

ity assumption can accelerate the solution process, in addition to the improved solution quality shown above in Table 2.

In Figure 3, we illustrate the numerical difficulty on an instance created from the inexact dataset with additional Poisson noise on the flow data. The step size used is the adaptive step size from [30]. We observe in particular that when using standard 64-bit precision, the primal value and FW gap can stall because of numerical errors. These numerical errors can be linked to the logarithm terms in the Poisson objective which induces a large range in the magnitude of the coefficients in the gradient. In addition, the non-smoothness yields challenging subproblems in the adaptive line search, potentially causing excessive estimates of the local Lipschitz constant. Such numerical challenges with the line search on self-concordant functions were already reported in [10, 11]. Other similar line searches such as the one introduced in [31] or [21] improve stability but at the cost of additional gradient evaluations in the line search procedure.

<figure><img src="image_6.png" /><figcaption>**Figure 3:** Convergence of the Poisson and least-square formulations on the instance salmon-895, with the primal and FW gap trajectories displayed for the least-square model “L”, the Poisson model using standard 64-bit floating points “PF” and the same loss optimized in extended precision using Julia’s `BigFloat` labeled “PE”.</figcaption></figure>
<table><caption>**Table 2:** Aggregate statistics on the inexact flow decomposition problem using all non-trivial instances. The path and flow errors and number of paths are computed on 29328 instances, excluding 33 instances for which either IP or IP-R did not find any primal solution in the time limit. For the path and flow errors and the number of paths, the first number is the arithmetic mean of the metric, the second is the shifted geometric mean (with shift 1), and the last is the number of instances on which the method performed the best on the given metric. The relative flow error is computed as the flow error divided by the number of edges of the graph. The flow error with the (o) label reports the mean (resp. geometric mean) for instances where both integer models were optimized to proven optimality, thus removing instances for which they computed a solution but could not close the gap. Only the shifted geometric mean of the runtime is displayed for conciseness.</caption><thead><tr><th>Metric</th><th>FW</th><th>FW-C</th><th>FW-P</th><th>IP</th><th>IP-R</th></tr></thead><tbody><tr><td>Path error</td><td>3.64 / 2.23 / 13078</td><td>3.62 / 2.22 / 13189</td><td>4.79 / 4.03 / 1846</td><td>4.38 / 3.68 / 4819</td><td>2.91 / 1.68 / 16584</td></tr><tr><td>Flow error</td><td>4.82e-01 / 1.15e-01 / 27118</td><td>4.81e-01 / 1.15e-01 / 27125</td><td>6.27e+01 / 5.55e+00 / 16171</td><td>4.38e+00 / 3.68e+00 / 104</td><td>2.91e+00 / 1.68e+00 / 10245</td></tr><tr><td>Flow err. (o)</td><td>3.80e-01 / 5.52e-02</td><td>3.80e-01 / 5.49e-02</td><td>3.93e+01 / 3.18e+00</td><td>3.49e+00 / 3.14e+00</td><td>2.06e+00 / 1.26e+00</td></tr><tr><td>Flow rel. err.</td><td>1.56e-02 / 6.79e-03</td><td>1.55e-02 / 6.77e-03</td><td>2.74e+00 / 9.55e-01</td><td>2.04e-01 / 2.00e-01</td><td>1.14e-01 / 1.09e-01</td></tr><tr><td># paths</td><td>4.69 / 3.96 / 10204</td><td>4.69 / 3.96 / 10204</td><td>4.79 / 4.03 / 9910</td><td>3.80 / 3.16 / 24659</td><td>4.19 / 3.54 / 17682</td></tr><tr><td>Time (s)</td><td>2.46e-04</td><td>7.91e-04</td><td>3.67e-01</td><td>3.43e+00</td><td>7.62e-01</td></tr></tbody></table>

## 4.2 Reconstruction with Error Distributions

In order to test the dependence of the different models' performance on the error distribution, we produce more instances from a subset of large splice graphs. We take all salmon instances for which the upper bound on the number of paths *k* is at least 16 and the number of edges at least 81, resulting in 109 instances. On these instances, we perturb the true flow $\bar{r}$ either with a Poisson distribution of parameter $\lambda_e = r_e$, or with a binomial distribution of parameters *p* = 0.5, *n* = 2$r_e$ so that in both cases, the pseudo-flow used by the methods remains nonnegative and of expectation $r_e$. The results are presented in Table 3; we removed FW-C since it performs similarly to FW on both groups of instances and IP given its poor performance on previous instances, and its need for an interval for the value of the flow on each edge.

As a first observation, the relative performance and behavior of all three models is similar for the two flow distributions, leading us to conclude that the models are robust beyond their initial modeling assumptions. We still note that the flow error of FW-P is lower under a Poisson distribution than under a binomial distribution, while the flow error of the least-square model is higher under a Poisson model than under the binomial one. The effect of the flow distribution on the flow error of IP-R does not follow as clear of a trend, the average relative flow error and the geometric mean of the flow error being higher under a Poisson distribution but the arithmetic mean of the absolute flow error being lower. The robust integer model IP-R performs the best in a majority of instances in terms of path error, followed by the Poisson regression model FW-P, followed by the least-square model FW. In terms of flow error, the least-square model outperforms the two other methods by far and achieves the best performance on almost all instances under both flow distributions.

<table><caption>**Table 3:** Results on the 109 large salmon instances with binomial and Poisson distributions of the observed flow. The presented metrics are identical to Table 2.</caption><thead><tr><th>Dist.</th><th>Metric</th><th>FW</th><th>FW-P</th><th>IP-R</th></tr></thead><tbody><tr><td rowspan="5">Binomial</td><td>Path error</td><td>25.61 / 24.84 / 4</td><td>20.56 / 19.88 / 51</td><td>18.86 / 18.40 / 65</td></tr><tr><td>Flow error</td><td>33.40 / 31.22 / 106</td><td>408.30 / 291.30 / 0</td><td>201.71 / 96.57 / 3</td></tr><tr><td>Flow rel. err.</td><td>0.36 / 0.35</td><td>4.47 / 3.33</td><td>1.97 / 1.24</td></tr><tr><td># paths</td><td>25.53 / 24.76 / 2</td><td>20.56 / 19.88 / 38</td><td>17.32 / 16.98 / 77</td></tr><tr><td>Time (s)</td><td>1.21e-02</td><td>8.00e+00</td><td>1.80e+03</td></tr><tr><td rowspan="5">Poisson</td><td>Path error</td><td>26.13 / 25.29 / 3</td><td>20.54 / 19.90 / 48</td><td>18.46 / 17.96 / 71</td></tr><tr><td>Flow error</td><td>49.63 / 45.86 / 107</td><td>388.32 / 285.85 / 0</td><td>191.70 / 123.45 / 2</td></tr><tr><td>Flow rel. err.</td><td>0.53 / 0.51</td><td>4.23 / 3.23</td><td>2.00 / 1.48</td></tr><tr><td># paths</td><td>26.06 / 25.29 / 1</td><td>20.54 / 19.90 / 33</td><td>16.91 / 16.43 / 84</td></tr><tr><td>Time (s)</td><td>2.61e-02</td><td>7.84e+00</td><td>1.80e+03</td></tr></tbody></table>

## 5 Conclusion

In this paper, we presented novel formulations for the sparse flow decomposition problem, formulating it as regression problems over the flow polytope and proposed an algorithmic framework producing sparse iterates and corresponding convex decompositions. The proposed methods are efficient and converge to the optimum of the corresponding formulation at a linear rate, ensuring their applicability to transcriptomics for genes with large splice graphs. The derived algorithm performs at each iteration a computation of the loss gradient and computes a shortest path on the DAG, thus leveraging specialized low-cost algorithms. Furthermore, they offer a high reconstruction performance, on par with the best integer optimization model from the literature, showing that even though previous formulations were NP-hard, they are not necessarily the (unique) way to approach transcript multi-assembly. In particular, these results show that even though the underlying decompositions are sparse, they are not necessarily the sparsest and seeking a minimum decomposition does not provide the best reconstruction error or path recovery. Due to the generic nature of the flow fitting optimization problem, our framework is applicable under a variety of loss functions, leaving the possibility to test more distributional assumptions in future work. Beyond the problem tackled in this paper, future work will also consider extending our proposed approach to other problems in which one seeks a solution constructed as convex or conic combinations of elementary combinatorial objects, including for instance flow decomposition of non-acyclic graphs [33].

## Acknowledgements

We thank Miroslav Kratochvíl, St. Elmo Wilken and Laurène Pfajfer for insights into the nature of transcriptomics problems and Panayotis Mertikopoulos for fruitful discussions on the objective function of our formulation. This work benefitted from the support of the FMJH Program PGMO.
References

[1] E. Bernard, L. Jacob, J. Mairal, and J.-P. Vert. Efficient RNA iso-
form identification and quantification from RNA-Seq data with network
flows. *Bioinformatics*, 30(17):2447–2455, 2014.

[2] M. Besançon, A. Carderera, and S. Pokutta. FrankWolfe.jl: A high-
performance and flexible toolbox for Frank-Wolfe algorithms and con-
ditional gradients. *INFORMS Journal on Computing*, 34(5):2611–2620,
2022.

[3] M. Besançon, S. Designolle, J. Halbey, D. Hendrych, D. Kuzinowicz, S. Pokutta, H. Troppens, D. V. Herrmannsdoerfer, and E. Wirth. Improved algorithms and novel applications of the FrankWolfe.jl library. *arXiv preprint arXiv:2501.14613*, 2025.

[4] M. Besançon. Code for efficient sparse flow decomposition methods for RNA multi-assembly, July 2025. URL https://doi.org/10.5281/zenodo.16033781.

[5] M. Besançon. Repository for efficient sparse flow decomposition methods for RNA multi-assembly, July 2025. https://github.com/inria-ghost/sparse-flow-decomposition.

[6] S. Bolusani, M. Besançon, K. Bestuzheva, A. Chmiela, J. Dionf-sio, T. Donkiewicz, J. van Doornmalen, L. Eifler, M. Ghannam, A. Gleixner, et al. The SCIP optimization suite 9.0. *arXiv preprint arXiv:2402.17702*, 2024.

[7] I. M. Bomze, F. Rinaldi, and D. Zeffiro. Active set complexity of the away-step Frank-Wolfe algorithm. *SIAM Journal on Optimization*, 30 (3):2470–2500, 2020.

[8] I. M. Bomze, F. Rinaldi, and D. Zeffiro. Frank-Wolfe and friends: a journey into projection-free first-order optimization methods. *4OR*, 19: 313–345, 2021.

[9] G. Braun, A. Carderera, C. W. Combettes, H. Hassani, A. Karbasi, A. Mokhtari, and S. Pokutta. Conditional gradient methods. *arXiv preprint arXiv:2211.14103*, 2022.

[10] A. Carderera, M. Besançon, and S. Pokutta. Simple steps are all you need: Frank-Wolfe and generalized self-concordant functions. *Advances in Neural Information Processing Systems*, 34:5390–5401, 2021.

[11] A. Carderera, M. Besançon, and S. Pokutta. Scalable Frank-Wolfe on generalized self-concordant functions via simple steps. *SIAM Journal on Optimization*, 34(3):2231–2258, 2024.

[12] B. V. Cherkassky, A. V. Goldberg, and T. Radzik. Shortest paths algo-
rithms: Theory and experimental evaluation. *Mathematical programming*, 73(2):129–174, 1996.

[13] F. H. Dias and A. I. Tomescu. Accurate flow decomposition via robust integer linear programming. *IEEE/ACM Transactions on Computational Biology and Bioinformatics*, 2024.

[14] F. H. Dias and A. I. Tomescu. Dataset for: Accurate flow decom-
position via robust integer linear programming, Mar. 2024. Zenodo,
https://doi.org/10.5281/zenodo.10775004.

[15] F. H. Dias, L. Williams, B. Mumey, and A. I. Tomescu. Efficient min-
imum flow decomposition via integer linear programming. *Journal of
Computational Biology*, 29(11):1252–1267, 2022.

[16] P. Dvurechensky, K. Safin, S. Shtern, and M. Staudigl. Generalized self-concordant analysis of Frank-Wolfe algorithms. *Mathematical Programming*, 198(1):255–323, 2023.

[17] M. Frank, P. Wolfe, et al. An algorithm for quadratic programming. *Naval research logistics quarterly*, 3(1-2):95–110, 1956.

[18] A. V. Goldberg and T. Radzik. *A heuristic improvement of the Bellman-Ford algorithm*. Citeseer, 1993.

[19] A. Grigorjew, F. H. C. Dias, A. Cracco, R. Rizzi, and A. I. Tomescu. Accelerating ILP solvers for minimum flow decompositions through search space and dimensionality reductions. In L. Liberti, editor, *22nd International Symposium on Experimental Algorithms (SEA 2024)*, volume 301 of *Leibniz International Proceedings in Informatics (LIPICS)*, pages 14:1–14:19, Dagstuhl, Germany, 2024. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. ISBN 978-3-95977-325-6. doi: 10.4230/LIPICS.SEA.2024.14. URL https://drops.dagstuhl.de/ entities/document/10.4230/LIPICS.SEA.2024.14.

[20] D. Hendrych, M. Besançon, and S. Pokutta. Solving the optimal ex-
periment design problem with mixed-integer convex methods. *arXiv
preprint arXiv:2312.11200*, 2025.

[21] D. Hendrych, S. Pokutta, M. Besançon, and D. Martínez-Rubio. Secant line search for Frank-Wolfe algorithms. In *Forty-second International Conference on Machine Learning*, 2025.

[22] M. Jaggi. Revisiting Frank-Wolfe: Projection-free sparse convex optimization. In *International conference on machine learning*, pages 427–435. PMLR, 2013.

[23] H. Karimi, J. Nutini, and M. Schmidt. Linear convergence of gradient and proximal-gradient methods under the Polyak-lojasiewicz condition.

In *Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings. Part I 16*, pages 795–811. Springer, 2016.

[24] K. Kloster, P. Kuinke, M. P. O'Brien, F. Reidl, F. S. Villaamil, B. D. Sullivan, and A. van der Poel. A practical FPT algorithm for flow decomposition and transcript assembly. In *2018 Proceedings of the Twentieth Workshop on Algorithm Engineering and Experiments (ALENEX)*, pages 75–86. SIAM, 2018.

[25] T. Koch, T. Berthold, J. Pedersen, and C. Vanaret. Progress in math-
ematical programming solvers from 2001 to 2020. EURO Journal on
Computational Optimization, 10:100031, 2022.

[26] E. S. Levitin and B. T. Polyak. Constrained minimization methods. *USSR Computational mathematics and mathematical physics*, 6(5):1–50, 1966.

[27] J. J. Li, C.-R. Jiang, J. B. Brown, H. Huang, and P. J. Bickel. Sparse linear modeling of next-generation mRNA sequencing (RNA-Seq) data for isoform discovery and abundance estimation. *Proceedings of the National Academy of Sciences*, 108(50):19867–19872, 2011.

[28] W. Li, J. Feng, and T. Jiang. Isolasso: a LASSO regression approach to RNA-Seq based transcriptome assembly. In *Research in Computational Molecular Biology: 15th Annual International Conference, RECOMB 2011, Vancouver, BC, Canada, March 28-31, 2011. Proceedings 15*, pages 168–188. Springer, 2011.

[29] M. Lubin, O. Dowson, J. D. Garcia, J. Huchette, B. Legat, and J. P. Vielma. JuMP 1.0: Recent improvements to a modeling language for mathematical optimization. *Mathematical Programming Computation*, 15(3):581–589, 2023.

[30] F. Pedregosa, G. Negiar, A. Askari, and M. Jaggi. Linearly convergent Frank-Wolfe with backtracking line-search. In *International conference on artificial intelligence and statistics*, pages 1–10. PMLR, 2020.

[31] S. Pokutta. The Frank-Wolfe algorithm: a short introduction. *Jahresbericht der Deutschen Mathematiker-Vereinigung*, 126(1):3–35, 2024.

[32] F. Sena and A. I. Tomescu. Safe paths and sequences for scalable ILPs in RNA transcript assembly problems. *arXiv preprint arXiv:2411.03871*, 2024.

[33] F. Sena, E. Ingervo, S. Khan, A. Prjibelski, S. Schmidt, and A. Tomescu. Flowtigs: safety in flow decompositions for assembly graphs. *iScience*, 27(12), 2024.

[34] Z. Song, Z. Xu, Y. Yang, and L. Zhang. Accelerating Frank-Wolfe algo-
rithm using low-dimensional and adaptive data structures, 2022. URL
https://arxiv.org/abs/2207.09002.

[35] T. Sun and Q. Tran-Dinh. Generalized self-concordant functions: a recipe for Newton-type methods. *Mathematical Programming*, 178(1): 145–213, 2019.

[36] A. I. Tomescu, A. Kuosmanen, R. Rizzi, and V. Mäkinen. A novel min-
cost flow method for estimating transcript expression with RNA-Seq.
In *BMC bioinformatics*, volume 14, pages 1–10. Springer, 2013.

[37] A. I. Tomescu, T. Gagie, A. Popa, R. Rizzi, A. Kuosmanen, and V. Mäkinen. Explaining a weighted DAG with few paths for solving genome-guided multi-assembly. *IEEE/ACM transactions on computational biology and bioinformatics*, 12(6):1345–1354, 2015.

[38] K. K. Tsuji, K. Tanaka, and S. Pokutta. Pairwise conditional gradients without swap steps and sparser kernel herding. In *International Conference on Machine Learning*, pages 21864–21883. PMLR, 2022.

[39] B. Vatinlen, F. Chauvet, P. Chrétienne, and P. Mahey. Simple bounds and greedy algorithms for decomposing a flow into a minimal set of paths. *European Journal of Operational Research*, 185(3):1390–1401, 2008.

[40] L. Williams, G. Reynolds, and B. Mumey. RNA transcript assembly using inexact flows. In *2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM)*, pages 1907–1914. IEEE, 2019.

[41] E. Wirth, J. Pena, and S. Pokutta. Fast convergence of Frank-Wolfe algorithms on polytopes. *arXiv preprint arXiv:2406.18789*, 2024.

[42] E. Wirth, M. Besançon, and S. Pokutta. The pivoting framework: Frank-Wolfe algorithms with active set size control. In *The 28th International Conference on Artificial Intelligence and Statistics*, 2025.

[43] Z. Xu, Z. Song, and A. Shrivastava. Breaking the linear iteration cost barrier for some well-known conditional gradient methods using MaxIP data-structures. *Advances in Neural Information Processing Systems*, 34:5576–5589, 2021.

[44] R. Zhao. New analysis of an away-step Frank-Wolfe method for minimizing log-homogeneous barriers. *Mathematics of Operations Research*, 2025.