Michael A Burr
email: burr2@clemson.edu

Shuhong Gao
email: sgao@clemson.edu

Elias Tsigaridas
email: elias.tsigaridas@inria.fr

The Complexity of an Adaptive Subdivision Method for Approximating Real Curves *

We present the first complexity analysis of the algorithm by Plantinga and Vegter for approximating real implicit curves and surfaces. This approximation algorithm certifies the topological correctness of the output using both subdivision and interval arithmetic. In practice, it has been seen to be quite efficient; our goal is to quantify this efficiency.

We focus on the subdivision step (and not the approximation step) of the Plantinga and Vegter algorithm. We begin by extending the subdivision step to arbitrary dimensions. We provide a priori worst-case bounds on the complexity of this algorithm both in terms of the number of subregions constructed and the bit complexity for the construction. Then, we use continuous amortization to derive adaptive bounds on the complexity of the subdivided region. We also provide examples showing our bounds are tight.

INTRODUCTION

Subdivision-based algorithms are one of the most commonly used algorithms in many fields, from computational geometry and graphics to solving of polynomials and mathematical programming [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF][START_REF] Heiden | Fast generation of molecular surfaces from 3d data fields with an enhanced "marching cube" algorithm[END_REF][START_REF] Schröder | Subdivision as a fundamental building block of digital geometry processing algorithms[END_REF][START_REF] Babuvška | Error estimates for adaptive finite element computations[END_REF][START_REF] Mantzaflaris | On continued fraction expansion of real roots of polynomial systems, complexity and condition numbers[END_REF][START_REF] Yap | Towards exact numerical Voronoi diagrams[END_REF][START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF][START_REF] Allgower | On the complexity of exclusion algorithms for optimization[END_REF].

Subdivision-based algorithms are intrinsically adaptive. Starting with the domain of interest, usually an axis-aligned box, they recursively split it into sub-domains, eliminating those that do not contain a solution or interesting features of the problem at hand. At the end, the algorithm produces a union of sub-domains (boxes) which lie in the initial domain.

The main advantages of subdivision-based algorithms are their great flexibility and their local nature. Because of their recursive character, they are easy to implement, and this makes them popular among practitioners. Moreover, they are often efficient in practice because they only perform additional subdivisions near difficult features. It is exactly because of these advantages, however, that the complexity analysis of subdivision-based algorithms is particularly challenging. To analyze them we need to understand, in depth, the local complexity of the input instance and how the inclusion/exclusion predicates behave on them.

For univariate problems, the analysis of subdivision algorithms is well-understood, and there are several results, especially for the case of approximating the roots of polynomials, e.g., [START_REF] Yap | Analytic root clustering: A complete algorithm using soft zero tests[END_REF][START_REF] Sagraloff | A simple but exact and efficient algorithm for complex root isolation and its complexity analysis[END_REF]. However, for higher dimensions, very little is known. For example, there are no explicit complexity results for pure subdivision-based algorithms for approximating curves and surfaces.

In this paper we fill this gap. We consider the following problem: Plantinga and Vegter [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] presented a subdivisionbased algorithm for correctly approximating curves and surfaces, see Figure 1. We call this the PV algorithm. It takes, as input, a polynomial f ∈ R[x, y] or R[x, y, z], whose real zero set is bounded and smooth, and a region I ⊆ R 2 or R 3 . From this input data, it constructs a piecewise-linear approximation to the zero set of f in I. The approximation has the correct topology in the sense that there is an ambient isotopy between the approximation and the zero set; additionally, by further subdivisions, the Hausdorff distance between the approximation and the zero set can be made as small as desired. The authors of [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] claim that the PV algorithm is efficient in practice, but, to our knowledge, there is no prior complexity analysis of the PV algorithm. The approximation produced by the PV algorithm [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] as well as the boxes constructed by the algorithm.

We provide the first complexity analysis for approximating curves and surfaces using the PV algorithm. Our bounds exploit local point estimates, called local size bounds, to quantify the amount of work needed by the algorithm at each point of the input domain. For the case of curves, we prove a worst-case bit complexity bound of 2 O(d 3 τ ) . Additionally, we extend the predicates for the PV algorithm to all dimensions and analyze the complexity of these tests in two and higher dimensions. Moreover, we provide an adaptive bound on the size of the subdivision using continuous amortization, first developed in [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF] and extended to higher dimensions in [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]. Our results are the first application of continuous amortization to a pure high dimensional problem. Furthermore, we prove that our bounds are tight in Lemma 6.1.

Related Work

The design of efficient subdivision-based algorithms that are output-sensitive, precision-sensitive, certified, and exploit the underlying structure of the problem is a great challenge and an active area of research. A big step in this direction is the introduction of soft tests [START_REF] Wang | On soft predicates in subdivision motion planning[END_REF][START_REF] Yap | Analytic root clustering: A complete algorithm using soft zero tests[END_REF] that, roughly speaking, replace hard exact tests, usually comparisons with zero, with approximate computations, and they are exact in the limit. They introduce a new notion of correctness called resolution-exactness. In this context, it is exactly the continuous amortization tool [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF][START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF] that captures the complexity of the soft predicates. Therefore, continuous amortization is a key tool for the analysis of such algorithms.

The previous work on subdivision methods and exclusion/inclusion predicates is quite extensive and so we only scratch its surface. For works that focus on classical exclusion/exclusion algorithms but without bit complexity bounds we refer the reader to [START_REF] Giusti | On location and approximation of clusters of zeros of analytic functions[END_REF][START_REF] Yakoubsohn | Approximating the zeros of analytic functions by the exclusion algorithm[END_REF][START_REF] Dedieu | Computing the real roots of a polynomial by the exclusion algorithm[END_REF]. For other approaches for approximating curves and surfaces we refer the reader to [START_REF] Cheng | Sampling and meshing a surface with guaranteed topology and geometry[END_REF][START_REF] Boissonnat | Isotopic implicit surface meshing[END_REF][START_REF] Boissonnat | Meshing of surfaces[END_REF][START_REF] Cheng | Certified rational parametric approximation of real algebraic space curves with local generic position method[END_REF] and the references therein. Recently, there is an extension to the case of analytic functions [START_REF] Giusti | On location and approximation of clusters of zeros of analytic functions[END_REF][START_REF] Yap | Analytic root clustering: A complete algorithm using soft zero tests[END_REF]. For the problem of isolating the roots of polynomials we refer the reader to [START_REF] Mantzaflaris | On continued fraction expansion of real roots of polynomial systems, complexity and condition numbers[END_REF][START_REF] Krawczyk | Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken[END_REF][START_REF] Henrici | Methods of search for solving polynomial equations[END_REF][START_REF] Collins | Polynomial real root isolation using Descarte's rule of signs[END_REF][START_REF] Yakoubsohn | Numerical analysis of a bisection-exclusion method to find zeros of univeriate analytic functions[END_REF][START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF][START_REF] Burr | Sqfreeeval: An (almost) optimal real-root isolation algorithm[END_REF][START_REF] Sagraloff | A simple but exact and efficient algorithm for complex root isolation and its complexity analysis[END_REF][START_REF] Cheng | Root isolation of zero-dimensional polynomial systems with linear univariate representation[END_REF] and the references therein. There are also approaches [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF] that achieve locally quadratic convergence towards the simple roots of polynomial systems and there very efficient in practice. Another interesting direction of subdivision algorithms, of more geometric nature, concerns the approximation of algebraic varieties [START_REF] Snyder | Interval analysis for computer graphics[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF][START_REF] Burr | Complete subdivision algorithms, ii: Isotopic meshing of singular algebraic curves[END_REF][START_REF] Sharma | Isotopic arrangement of simple curves: an exact numerical approach based on subdivision[END_REF][START_REF] Yap | Adaptive isotopic approximation of nonsingular curves: the parameterizability and nonlocal isotopy approach[END_REF][START_REF] Lin | Non-local isotopic approximation of nonsingular surfaces[END_REF], and the computation of the approximate Voronoi diagrams [START_REF] Yap | Towards exact numerical Voronoi diagrams[END_REF]. There are also quite important applications of these algorithms to the problem of robot motion planning [START_REF] Wang | On soft predicates in subdivision motion planning[END_REF].

Notation

We use dist C (x, f ) for the distance, in C, of the point x to the hypersurface defined by the polynomial f . The point x might be real or complex. We also use Bdista(V, W ) to indicate that we are interested in the minimum distance between V, W ⊂ C n inside a hyperbox of C n whose corners have bit-size at most a.

For an interval or a (hyper-)box J we denote the midpoint as m = m(J) and side length as w = w(J). The notation f (J) denotes the interval (over-)approximation of f applied to the region J. In this paper, because of its nice properties, we use the centered interval form for f which is based on the Taylor expansion of f at the midpoint of J, see Remark 2.2 or [START_REF] Ratschek | Computer methods for the range of functions[END_REF] for further details.

We use O(•) and OB(•) to denote the arithmetic complexity and bit complexity, respectively. The soft-O notation, O(•) and OB(•), mean that we are ignoring logarithmic factors. By V C (f ), V R (f ), or VI (f ) we denote the zero set of a polynomial f over the complex numbers, real numbers, or the region I, respectively.

Outline

The rest of the paper is organized as follows: In Section 2, we recall the PV algorithms and the inclusion and inclusion predicates that it uses. We unify and generalize these predicates to make them applicable to higher dimensions. In Section 3, we derive the local size bounds for the inclusion/exclusion predicates for the PV algorithms. In Section 4, we use the local size bounds to provide worst-case estimates on the number of boxes and the bit complexity of the PV algorithms in terms of the size of the input. In Section 5, we present an introduction to and apply continuous amortization to present adaptive complexity bounds for the PV algorithms. Finally, in Section 6 we present example to demonstrate the tightness of our bounds.

THE MODIFIED PLANTINGA AND VEG-TER ALGORITHM

In this section, we present the subdivision part of the PV algorithm. We give details on how PV exploits interval arithmetic to achieve certified computations. We slightly modify the algorithm to unify the tests and to make the subdivisions applicable in arbitrary dimensions.

2.1

The PV Algorithm [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] Let f ∈ R[x, y] or R[x, y, z] be a square-free polynomial such that its real zero set V R (f ) is smooth and bounded. The PV algorithm recursively subdivides an initial bounding square or cube I for the variety with a quad-tree or oct-tree data structure until at least one of the following two tests holds on each subregion J. In the literature, these tests are often referred to as C0 and C1: C0(J) := 0 ∈ f (J) C1(J) := 0 ∈ ∇f (J), ∇f (J) .

When C0(J) holds, the variety does not enter the region J and so J can be discarded. On the other hand, when C1(J) holds, the curve or surface does not bend much within the region J. More precisely, given f and I, the PV algorithms construct a partition P of I so that C0 or C1 is true on each region in P . Initially, P = {I}:

Algorithm 2.

Main subdivision of PV algorithm

Repeatedly subdivide (into 4 or 8 children) each J ∈ P until one of the following conditions hold:

C0(J) is True or C1(J) is True.
After every sub-region J satisfies C0(J) or C1(J), the authors of [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] perform post-processing steps, which include balancing the tree, evaluating the sign of f on the corners of each J in P , and using sign changes along the sides of regions J to detect and approximate the curve or surface. This approximation is topologically correct as there is an ambient isotopy between the approximation and the variety V R (f ). Additionally, by further subdivision, the isotopy can be made sufficiently small so that the Hausdorff distance between the approximation and the variety is as small as desired. We note that it is possible to extend the PV algorithm in the plane to provide an approximation even when V R (f ) is unbounded, V R (f ) is singular, and I is not a bounding box, see [START_REF] Burr | Complete subdivision algorithms, ii: Isotopic meshing of singular algebraic curves[END_REF]. In this paper, however, we focus on the original PV algorithm without the restriction of a bounded curve.

The Subdivisions of the PV Algorithm

Our main focus is on computing the number of regions that the PV algorithms construct, see Algorithm 2.1, and not on the approximation of the curve or surface, per se. Therefore, we focus exclusively on the C0 and C1 tests and apply them in arbitrary dimensions. More precisely, let

f ∈ R[x1, • • • , xn] be such that its real zero set V R (f ) is smooth. Let I ⊆ R n be a n-dimensional real cube.
Then, we can generalize the tests C0 and C1, along with Algorithm 2.1, to n dimensions, where the subdivision splits an n-cube into 2 n children. However, we mention that, in this case, we no longer use the output of the algorithm to construct an approximation to V R (f ).

Extending the C 1 test

The predicate C1(J) has the following two consequences that are fundamental in the proof of correctness of the PV algorithm in [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]: (1) If a region J satisfies the conditions, then, in J, there cannot be any pair of gradient vectors which are orthogonal to each other. [START_REF] Babuvška | Error estimates for adaptive finite element computations[END_REF] The variety V R (f ) is parametrizable in the direction of at least one of the coordinate axes. Fact ( 2) is a direct consequence of Fact ( 1), but it is used so frequently in the proofs in [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF], that it is worthwhile to mention it explicitly.

We now modify and extend test C1 so that C0 and C1 have the same form. Let the function g : R n ×R n → R, defined as

g(x1, • • • , xn, y1, • • • , yn) = ∇f (x1, • • • , xn), ∇f (y1, • • • , yn) .
It follows that for a region J, if 0 ∈ g(J × J), then there is no pair of gradient vectors in J which are orthogonal to each other. We use this reformulation because it implies the same two consequences, Facts (1) and ( 2), as the original C1 test, but the application of interval arithmetic appears as the last step as opposed to an intermediate step. In particular, both tests C0 and C1 are of the same type, i.e., they consist of testing where 0 appears in the interval formulation of a function applied to a region. For the rest of the paper, all references to the C1 test refer to this new C1 test.

Remark 2.2. We use the centered interval form for interval arithmetic, see [START_REF] Ratschek | Computer methods for the range of functions[END_REF], which is based on the Taylor expansion of the polynomial at the midpoint of the region. In particular, let J be the n-cube with midpoint m = m(J) and side length w = w(J). Then, test C0 simplifies to

(C0) |f (m)| > |α|≥1 |∂ α f (m)| α! w 2 |α| where α = (α1, • • • , αn) ∈ N n is a multi-index. In par- ticular, |α| = αi is the sum of the entries in α, ∂ α = ∂ α 1 1 • • • ∂ αn n
is the operator of partial derivatives applied with multiplicity determined by the entries of α, and α! = αi! is the product of the factorials of the entries in α.

Since the C1 test is based on the function g whose domain is 2n-dimensional and the square J × J has midpoint (m, m) and the same side length as J, the C1 test becomes

(C1) |∇f (m)| 2 > |α|+|β|≥1 n i=1 ∂ α+e i f (m) • ∂ β+e i f (m) α! • β! w 2 |α|+|β|
where ei is the i-th standard basis vector, and α, β ∈ N n are multi-indices.

LOCAL SIZE BOUND FOR THE PV AL-GORITHM

The key to the bounds in this paper is a function, called a local size bound, on R n which locally describes the maximum amount of work that is required at each point in R n (cf [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]):

Definition 3.1. Let C be a predicate on n-dimensional cubes. A local size bound for C is a function F : R n → R ≥0 with the property that F (x) ≤ inf J, n-dim cube J x C(J)=False

Voln(J).

In other words, F (x) is a lower bound on the n-dimensional volume of a n-dimensional cube which contains x, but fails the stopping criterion test.

We develop a local size bound for the PV algorithm. For the PV algorithm, the local size bound relates the value of the partial derivatives at a point in R n to the distance from that point to the variety. The goal is to derive point estimates for the C0 and C1 predicates, Corollaries 3.6 and 3.7. As our first step, we reduce a higher-dimensional problem to a collection of one-dimensional problems as follows:

p v α 3 α 1 α 2 Figure 2: For a polynomial f ∈ R[x,
y] and a point p ∈ R 2 , we consider the roots of f , α1, α2, α3, in the direction of a unit vector v.

Definition 3.2. Let f ∈ R[x1, • • • , xn], p ∈ R n
, and v ∈ S n-1 . Define fv(t) to be the univariate polynomial passing through p and in the direction v, i.e., fv(t) = f (p + tv), see Figure 2. Then, define Σ fv be the sum of the reciprocals of the complex roots of fv, i.e.,

Σ fv (p) = s∈V (fv ) 1 |s| .
As in [START_REF] Burr | Sqfreeeval: An (almost) optimal real-root isolation algorithm[END_REF], Σ fv (p) links the Taylor coefficients of f to the geometry of the zero set of f . This relationship is explicitly explored in the following two results.

Lemma 3.3. Let f ∈ R[x1, • • • , xn], p ∈ R n , and v ∈ S n-1 . Then 1 f (p) • d k f (p + tv) dt k t=0 ≤ (Σ fv (p)) k ≤ deg(f ) dist C (p, f ) k .
Proof. The claim is trivial when k = 0. Since fv is a univariate polynomial, the first inequality follows directly from [7, Lemma 2.1]. The second inequality follows because, in the sum for Σ fv (p), there are at most deg(f ) terms and each element of the sum is the inverse of the distance between p and a point on V C (f ), which is bounded above by the inverse of the distance to the closest point.

The interval approximations in the PV predicates are based on the centered form, which, in turn is based on the Taylor expansion of the input polynomial, see Remark 2.2. We, therefore, use the geometric bounds from Lemma 3.3 to bound the Taylor coefficients in terms of the geometry of the zero set of f . The proof of Proposition 3.4 and Corollary 3.5 as well as other results that we need to effectively bound the local size bound appear in the full version of the paper.

Proposition 3.4. Let f ∈ R[x1, • • • , xn] and p ∈ R n . Then, for all multi-indices α ∈ N n , 1 f (p) k α ∂ |α| f ∂x α (p) ≤ 2 (n-1)(|α|+1) deg(f ) dist C (p, f ) |α| where k α = k! α!
is the multinomial coefficient. Each of the predicates for the PV algorithm test whether 0 is included in an interval approximation to the value of an appropriate function on a region. In practice, this test is performed by comparing the sizes of the Taylor coefficients and the size of the input region to the value of the polynomial at a given point, see Remark 2.2. In particular, if the inequality appearing in the conclusion of the following corollary holds, then 0 is not included in the interval approximation, for additional details, see [START_REF] Burr | Sqfreeeval: An (almost) optimal real-root isolation algorithm[END_REF].

Corollary 3.5. Let f ∈ R[x1, • • • , xn] and p ∈ R n . Sup- pose that 0 < w ≤ dist C (p,f ) ln(1+2 2-2n ) 2 n-1 deg(f ) . Then deg(f ) k=1 |α|=k 1 k! k α 1 f (p) • ∂ |α| f ∂x α (p) w 2 k < 1.
We now develop bounds on the size of J which guarantee the success of the PV predicates which can be applied to any point within J, not merely the midpoint.

Corollary 3.6. Let f ∈ R[x1, • • • , xn] and J ⊆ R n . As- sume that there is a point x ∈ J such that w(J) ≤ 2 ln 1 + 2 2-2n dist C (x, f ) 2 n deg(f ) + √ n ln (1 + 2 2-2n
) .

Then, C0(J) is true.

Proof. This result follows from Corollary 3.6 and [7, Section 3].

Corollary 3.7. Let f ∈ R[x1, • • • , xn] and g ∈ [x1, • • • , xn, y1, • • • , yn] where g(x1, • • • , xn, y1, • • • , yn) = ∇f (x1, • • • , xn), ∇f (y1, • • • , yn) . Let J ⊆ R n , and suppose that there is a point (a, b) ∈ J × J such that w(J) ≤ 2 ln 1 + 2 2-4n dist C ((a, b), g) 2 2n+1 (deg(f ) -1) + √ 2n ln (1 + 2 2-4n
) .

Then, C1(J) is true.

Proof. This result follows from Corollary 3.6 and [7, Section 3].

WORST-CASE BOUNDS

Observe that Corollary 3.6 applies to a polynomial in n variables and Corollary 3.7 uses a polynomial in 2n variables. Therefore, the varieties V C (f ) and V C (g) are in different dimensional spaces. In the arguments that follow, it is easier to study both of these varieties in the same dimensional space. Therefore, we consider V C (f ) as a subset of the diagonal in 2n-dimensional space. In particular, let the variables of C 2n be {x1, • • • , xn, y1, • • • , yn}; the diagonal ∆ consists of the points of the form xi = yi. Then, ∆ is n-dimensional and we identify C n with ∆, in particular, we write V ∆ C (f ) for the collection of all points (x1,

• • • , xn, y1, • • • , yn) such that f (x1, • • • , xn) = 0 and xi = yi for all i.

The Global Depth of the Subdivision Tree

We prove a global bound on the depth of the subdivision tree.

Proposition 4.1. Let f ∈ R[x1, • • • , xn], and let g ∈ R[x1, • • • , xn, y1, • • • , yn] where g(x1, • • • , xn, y1, • • • , yn) = ∇f (x1, • • • , xn), ∇f (y1, • • • , yn) . Suppose that I ⊆ R n and VI (f ) is smooth, i.e.
, or all x ∈ I, f (x) and g(x, x) are not both zero. In particular, let

0 < δ ≤ min x∈I max{dist C (x, f ), dist C ((x, x), g)}. Define D = max 2 ln 1 + 2 2-2n 2 n deg(f ) + √ n ln (1 + 2 2-2n ) , 2 ln 1 + 2 2-4n 2 2n+1 (deg(f ) -1) + √ 2n ln (1 + 2 2-4n
) .

The PV algorithm performs at most 2w(I)

Proof. 1 2 Dδ is a lower bound on the width of a region which occurs as a leaf in the subdivision tree. Therefore, we can divide the total n-dimensional volume into regions of this size.

Let f ∈ R[x1, • • • , xn] define a smooth variety, and let g ∈ R[x1, • • • , xn, y1, • • • , yn] where g(x1, • • • , xn, y1, • • • , yn) = ∇f (x1, • • • , xn), ∇f (y1, • • • , yn)
. Since f is smooth as a complex variety, there do not exist points x ∈ R n such that x ∈ V C (f ) and (x, x) ∈ V C (g). Suppose that ε is a separation bound between V ∆ C (f ) and V C (g) and between R n and V C (f, ∇f • ∇f ). Then, for all x ∈ R n , either the distance

dist C (x, f ) is at least ε 2 √ 2 or the distance dist C ((x, x), g) is at least ε 2 . Therefore, we may let δ = ε 2 √ 2 in Proposition 4.1.

Bounds on the distance between varieties

Suppose that the input polynomial f has integer coefficients. Then, we can express the bound in Proposition 4.1 in terms of the degree and the maximum bitsize of the coefficients of the input polynomial. For this, we need to bound ε, that appears in the previous section, from below. In other words, we need to bound the distance between V ∆ C (f ) and V C (g). We restrict our attention to the interior of an ndimensional cube such that a is an upper bound on the bitsize of the corners of the box. Following [START_REF] Sombra | Estimaciones para el teorema de ceros de Hilbert[END_REF], we define the following distance function, for a > 0,

Bdista(V, W ) = inf{ p -q : p ∈ V, q ∈ W, lg( p ), lg( q ) ≤ a}
We observe that Bdista(V, W ) is based on the distance between complex varieties, and it is not restricted to their real portions. The following theorem presents a lower bound on the distance of two varieties, when at least one of them is defined as a complete intersection, and is due to Martin Sombra, see also [31, §2.3.3], [START_REF] D'andrea | Heights of varieties in multiprojective spaces and arithmetic nullstellensätze[END_REF].

Theorem 4.2. Let V be of pure dimension r and W be of dimension s and defined by f1, . . . , fn-s ∈ Z[x1, . . . , xn]. Let a > 0, d = maxi{deg(fi)} and the maximum logarithmic height of the polynomials fi be h = maxi{h(fi)}. We can bound the distance between V and W within a ball B(0, 2 a ) as

-lg Bdista(V, W ) ≤ d n-s (h(V ) + deg(V )) (n -s) h d + (4r + 10) lg(n + 2) + 2 ,
where deg(V ) and h(V ) denote the degree and the canonical height of V , respectively. Furthermore, if we assume that V is a complete intersection defined by polynomials g1, . . . , gn-r, such that deg(gi) ≤ d2 and h(gi) ≤ h2, then

h(V ) ≤ (n -r) d n-r 2 h2 d2 + lg(n + 1) . ( 1
)
In this case the distance between V and W becomes:

lg ∆ ≥ -c(n)d n-s d n-r 2 h d + h2 d2 + a + 1 , (2)
where c(n) is a constant depending only on n.

To use Theorem 4.2 and Equation ( 2) to bound the distance between V ∆ C (f ) and V C (g), we consider the hypersurface V defined by the following equation:

g(x1, • • • , xn, y1, • • • , yn) = ∇f (x1, • • • , xn)•∇f (y1, • • • , yn).
Moreover, W is defined by the system of equations

f (x1, • • • , xn) and {yi = xi}.
Since g is a single equation, V is a r = (2n -1)-dimensional complex hypersurface. On the other hand, f defines an s = (n -1)-dimensional complex variety. Both varieties are complete intersections. We assumed that the input polynomial f has integer coefficients, that is

f ∈ Z[x1, • • • , xn] such that d = deg(f )
is the degree of f and that the maximum bitsize of the coefficients, aka logarithmic height, is h(f

) = h2 = τ . Then, deg(W ) = d, deg(∇f ) ≤ d -1, and h(∇f ) = τ + lg(d).
From this, it follows that h1 = h(g) = O(τ + lg(nd)) and

δ1 = deg(V ) = deg(g) ≤ 2d -2.
We can use the bound of Eq. ( 1) to bound the (logarithmic) height of a variety, and in our case V , as

h(V ) ≤ 2τ + (4n + 2) lg d = O(τ + n lg(nd)).
Combining these inequalities with Theorem 4.2, see also Equation ( 2), we obtain the following bound for the distance between V ∆ C (f ) and V C (g) within an n-dimensional cube with coordinates having bitsize bounded by a:

-lg(Bdista(V ∆ C (f ), V C (g)) ≤ d n+1 (6nτ + 40nd lg(nd) + 2da) = O(d n+1 (nτ + nd lg(nd) + da)).
The previous bound is used to bound ε and thus δ in Proposition 4.1. In particular, we apply the previous bound to a covering of I by small boxes the size of the separation bound for V C (f, ∇f • ∇f ) by translating and scaling these boxes. Therefore, we obtain the following explicit bound for the number of steps, or in other words the number of boxes, of the PV algorithm.

Theorem 4.3. The number of steps that PV algorithm performs is at most

2 O(nd 10n+1 (nτ +nd lg(nd)+9n+d)a) .
where a > 0 is the maximum bitsize of the coordinates of the corners of I.

Overall complexity bound

The bound on the number of steps that we have computed in Theorem 4.3 is an important quantity needed for the goal of estimating the overall bit complexity of the PV algorithm. PV is a subdivision algorithm that, roughly speaking, mimic the process of binary search. Therefore, to bound its Boolean complexity we have to multiply the number steps (subdivisions) that it performs with the worst case Boolean complexity of each step. Theorem 4.3 provides the number of steps. It remains to estimate the complexity of each step. A closer look at the predicates needed for the realization of the algorithm reveals that each step of the PV algorithm consists of a multivariate Taylor shift. Given a polynomial F ∈ Z[x1, . . . , xn] and numbers a1, . . . , an, we need to compute the coefficients of F (x1 + a1, . . . , xn + an).

We perform this operation recursively. First, we consider the bivariate case. We are given a polynomial F ∈ Z[x1, x2] of total degree d and maximum coefficient bitsize τ , and integers a1, and a2 of bitsize . We want to estimate the complexity of computing the polynomial F (x1 + a1, x2 + a2).

We consider the polynomial as

F ∈ (Z[x1])[x2]
, that is as a univariate polynomial in x2 with coefficients in x1; we write it as

F = d i=0 Fi(x1) x i 2 . The Taylor shift is F (x1 + a1, x2 + a2) = d i=0 Fi(x1 + a1) (x2 + a2) i .
There are d + 1 univariate Taylor shifts, Fi(x1+a1), that we have to perform. Each costs OB(d 2 + dτ ) [START_REF] Zur Gathen | Fast algorithms for taylor shifts and certain difference equations[END_REF]; hence the cost of all of them is OB(d 3 + d 2 τ ). We also have to compute (expand) the d+1 polynomials (x2 +a2) i . Each "expansion" costs at most OB(d 2 ), as it corresponds, in the worst case, to O(lg d) multiplications of two polynomials of degree at most d and bitsize O(d ). We perform all of them in OB(d 3 ).

Finally, to get To obtain the worst case complexity of each step of the PV algorithm, we replace in Corollary 4.5 with the (exponent of the) bound of Theorem 4.3. To see this, notice that at each step of the subdivision we increase by 1 the number of bits of the number we have to perform Taylor shifts with. Hence, in the worst case, we have to perform a Taylor shift with number have as bitsize the number of subdivisions. For the overall complexity of the algorithm we multiply this bound with the number of steps of Theorem 4.3. In the worst case, since we construct a polynomial of degree 2d -2 in 2n variables, each multivariate Taylor shift costs OB(n 2 d 2n+3 + nd 2n+2 (τ + da)) and the overall complexity of the algorithm is

F (x1+a1, x2+a2) = d i=0 Fi(x1 + a1) (x2+ a2) i ,
2 d 10n+1 (6n 2 τ +42n 2 d lg(nd)+2(d+n)a)+n OB(n 2 d 2 (τ + d(1 + a))),
where a is the maximum bitsize of the coordinates of the corners of I.

For the 2D case, n = 2, that commonly appears in applications, we obtain the following bound for applying PV is Theorem 4.6. The bit complexity of the PV algorithm for curves, is OB(2 d 21 (24τ +168 lg(d)+2da+1) ).

ADAPTIVE BOUNDS

Continuous amortization was introduced in [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF] as a way to adaptively analyze the complexity of subdivision-based algorithms. The theory of continuous amortization was extended to higher dimensions in [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]. We recall this technique in the special case for computing the number of regions created by a subdivision algorithm in R n which recursively subdivides n-dimensional cubes in R n into 2 n identical smaller n-dimensional cubes.

Theorem 5.1 ( [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF][START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]). Let F be a local size bound (see Definition 3.1) for a stopping criterion C and I be an ndimensional cube. The number of regions formed by a subdivision algorithm which recursively subdivides n-dimensional cubes in R n into 2 n identical smaller n-dimensional cubes is at most max 1,

I 2 n dVn F (x)
where dVn is the n-dimensional volume form. If the algorithm does not terminate, then the integral is infinite.

We can use continuous amortization to express the complexity of the PV algorithm. We express the PV subdivision as a subdivision of 2n-dimensional real space, but restrict our attention to the diagonal ∆. Since the real part of the diagonal ∆ can be identified with R n via projection onto the first n coordinates, we can interpret all of our subdivisions as subdivisions of R n . In particular, for J × J ∈ R 2n , under a standard subdivision, this region would be split into 2 2n children. For our case, however, we only need to consider the subregions which intersect the diagonal; in other words, we only consider subregions of the form L × L. By via the identification to R n this subdivision corresponds to the standard subdivision in R n where L is a child of J, and J is 2 n times bigger than L.

Proposition 5.2. Let f ∈ R[x1, • • • , xn], and let g ∈ R[x1, • • • , xn, y1, • • • , yn] where g(x1, • • • , xn, y1, • • • , yn) = ∇f (x1, • • • , xn), ∇f (y1, • • • , yn) . Suppose that I ⊆ R n .
The number of regions after the subdivision performed by the PV algorithm (before balancing) is bounded above by the maximum of 1 and

2 n I min 2 n deg(f ) + √ n ln 1 + 2 2-2n 2 ln (1 + 2 2-2n ) dist C (x, f ) n , 2 2n+1 (deg(f ) -1) + √ 2n ln 1 + 2 2-4n 2 ln (1 + 2 2-4n ) dist C ((x, x), g) n dVn
where dVn is the n-dimensional volume form.

Proof. These results follow almost directly from a standard application of continuous amortization, see [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]. The only difference is that the argmax in the proof of continuous amortization must be taken over the diagonal, (J × J) ∩ ∆ = {(a, a) : a ∈ J} instead of all of J × J.

This integral provides a more adaptive and accurate estimate on the complexity than the worst-case a priori bounds based on the size of the input. Moreover, this integral can be evaluated even when the input polynomial has complex (but not real) singularities. Examples of the evaluation of this integral appear in the remainder of this paper.

EXAMPLES

Both bounds are exponential with respect to the degree of the polynomial f and the number of variables. They remain exponential even if we assume that the number of variables is constant. In [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF], the authors show that for several examples the computation time is efficient in practice. The following lemma, which is also illustrated in Figure 3(a), proves that this exponential behavior is optimal, up to constants in the exponents.

(a) (b) The width of boxes between vertical lines is at most 2a -n 2 -1 and they extend the entire length of the initial region. The number of regions is bounded by: Ω(w(I)a n 2 +1 ) where w(I) is the width of he initial region. (b) The approximation of f (x, y) = x 2 + y 2 + ε 2 . The number of regions is bounded by: O(lg(w(I)) -lg(ε)). Proof. Following the construction in [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF], consider the Mignotte polynomial P (x) = x d -2(ax -1) 2 and the related polynomial P2(x) = x d -(ax -1) 2 where a is a sufficiently large positive integer. The product P (x)P2(x) is of degree 2d and the largest coefficient is of size 2a 4 . In [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF], it is shown that the product P (x)P2(x) has (at least) three roots in the interval (a -1 -h, a -1 + h) where h = a -d/2-1 . Treating P (x)P2(x) as a polynomial in n variables, we see that the PV algorithm to approximate the variety in an n-dimensional cube I of side length w(I) requires subdividing to regions of side length at most 2h to separate the three vertical hyperplanes in the interval (a -1 -h, a -1 + h). Since this occurs along an entire hyperplane of the input region, the number of small boxes is, at least, w(I) n-1 2h = 1 2 w(I) n-1 a d/2+1 , which is exponential in both the size of the input of the box and the size of the coefficients of the polynomial.

Even though our bounds are optimal, in practice, these are quite pessimistic, as the actual separation bounds do not follow the worst case behavior, see Figure 3 Moreover, for each of these examples, the minimum distance between V ∆ C (f ) and V C (g) is at most ε. Therefore, a bound coming from Proposition 4.1 would be much larger than the bound continuous amortization provides.

It remains an open question to deduce the bit complexity bounds for the PV algorithms in the general case from the adaptive bounds of Proposition 5.2. Since the complexity of the algorithm can be exponential in the inputs, the integral must be described in terms of additional geometric and intrinsic parameters.

Figure 1 :

 1 Figure 1: (a) The graph of the real trace of the curve f = 3y 3 +3xy 2 -2x 3 -3y 2 +xy +3x 2 -3y +3x+2. (b) The approximation produced by the PV algorithm [25] as well as the boxes constructed by the algorithm.

  we need to perform the d + 1 multiplications between polynomials Fi(x1 + a1) and (x2 + a2) i to obtain each term in the sum, and d + 1 additions. The multiplications dominate the complexity. As each multiplication corresponds, in the worst case, to the multiplication of two polynomials of degree d and bitsize O(τ + d ) and O(d ), it costs OB(dτ + d 2 ). All of them cost OB(d 3 + d 2 τ ). Lemma 4.4. The Taylor shift for a bivariate polynomial of total degree d and maximum coefficient bitsize τ costs OB(d 3 + d 2 τ ). Using induction we obtain the following result: Corollary 4.5. The Taylor shift for a multivariate polynomial of total degree d in n variables and maximum coefficient bitsize τ costs OB(d n+1 + d n τ ).

Figure 3 :

 3 Figure 3: (a) The output of the PV algorithm for f (x, y) = (x n -2(ax -1) 2 )(x n -(ax -1) 2 ). The solutions to f (x, y) = 0 are close vertical lines (n = a = 3).The width of boxes between vertical lines is at most 2a -n 2 -1 and they extend the entire length of the initial region. The number of regions is bounded by: Ω(w(I)a

Lemma 6 . 1 .

 61 The bound of Theorem 4.3 is asymptotically tight.

Example 6 . 2 .+ x 2 2 + ε 2 . 2 + ε 2 2 1 + x 2 2 .Example 6 . 3 .+ x 2 2 -4ε 2 and

 622222226322 (b). This is illustrated in the following two examples: Fix ε > 0 and consider f (x1, x2) = x 2 1 Then, dist C ((x1, x2), f ) = and dist C ((x1, x2, x1, x2), g) = x For any square I, the number of regions constructed by the PV algorithm isO(lg(w(I)) -lg(ε)). Fix ε > 0 and consider f (x1, x2) = x 2 1 ε 2 . Then, dist C ((x1, x2), f ) = dist C ((x1, x2, x1, x2), g) = x 2 1 + x 2 2 .For any square I, the number of regions constructed by the PV algorithm is O(lg(w(I)) -lg(ε)).

Acknowledgments. The authors are greatful to Martin Sombra for his comments and for providing them the lower bound on the distance between varieties.

* This work was significantly advanced while the first and third authors were supported at the Fields Workshop on Algebra, Geometry, and Proof in Symbolic Computation† Partially supported by a grant from the Simons Foundation (#282399 to Michael Burr) and National Science Foundation Grant CCF-1527193. ‡ Partially supported by the National Science Foundation under Grants CCF-1407623, DMS-1403062 and DMS-1547399. § Partially supported by HPAC (ANR ANR-11-BS02-013) and an FP7 Marie Curie Career Integration Grant.


M. Dyer 1 and L. Stougie 2 Hanasusanto et al. [HKW] pointed out a mistake in our paper. To be precise, the error occurs in the first formula on page 429. We gratefully acknowledge their careful reading of our paper. We repair the error here in this corrigendum. We have chosen to completely rewrite Section 3.2 of the paper. We have done this so that the correction can be read almost independently of the rest of the paper, and to take the opportunity to improve the exposition.

Continuous distributions

For two-stage stochastic programming problems with continuously distributed parameters, P-hardness of obtaining the optimal solution can be established under even the mildest conditions on the distributions. For the proof we use a reduction from the problem of computing the volume of a knapsack polytope, proved P-hard in [START_REF] Dyer | On the complexity of computing the volume of a polyhedron[END_REF].

Let Z + denote the nonnegative integers, and Q + the rational numbers. Definition 3.2. Let P = {w ∈ [0, 1] n | n j=1 α j w j ≤ β} be the knapsack polytope, where α j ∈ Z + (j = 1, . . . , n) and β ∈ Q + . We consider the following two computational problems. (ii) Volume of knapsack polytope. Compute Vol(P ), the volume of P . These problems are P-hard. In [START_REF] Dyer | On the complexity of computing the volume of a polyhedron[END_REF], Volume of knapsack polytope was proved P-hard by reduction from the known P-hard problem Counting knapsack solutions.

Claim 1. Counting knapsack solutions remains P-hard if α

n = β + 1. Proof. If α n = β +1, let P = {w n-1 ∈ [0, 1] n-1 | n-1 j=1 α j w j ≤ β }|.
Then N (P ) = N (P ), since α j ∈ Z + (j = 1, . . . , n -1), and w ∈ {0, 1} n ∩ P implies w n = 0. Clearly, computing N (P ) is P-hard. Hence computing N (P ) is P-hard when α n = β + 1.

Claim 2. Volume of knapsack polytope remains P-hard if α n > β.

1 School of Computing, University of Leeds, Leeds, LS2 9JT, UK; e-mail: dyer@comp.leeds.ac.uk 2 Department of Econometrics and Operations Research, VU University, Amsterdam and CWI, P.O. Box 94079, NL-1090 GB Amsterdam, The Netherlands; e-mail: stougie@cwi.nl, Corresponding author.

Proof. For P as in Claim 1, let P i = {w ∈ {0, 1} n | n j=1 α j w j ≤ β + i/n}, where (0 ≤ i < n). So N (P i ) = N (P ) for all 0 ≤ i < n, since α j ∈ Z + (j = 1, . . . , n), and β ≤ β + i/n < β + 1. Thus computing N (P i ) is P-hard.

The reduction on p. 970 of [START_REF] Dyer | On the complexity of computing the volume of a polyhedron[END_REF] implies that, if computing N (P ) is P-hard, then computing Vol(P k ) is P-hard, for some 0 ≤ k < n. It follows that computing Vol(P ) is P-hard if

β = β + k/n < β + 1 = α n .
Theorem 3.2. Determining the optimal solution of a two-stage stochastic programming problem with continuously distributed parameters is P-hard, even if all the stochastic parameters have the uniform [0, 1] distribution. Proof. Define i.i.d. random variables q 1 , . . . , q n-1 , each uniformly distributed on [0, 1], and write q n-1 = (q 1 , . . . , q n-1 ). Now, given an instance of computing the volume of a knapsack polytope as in Claim 2, consider the following two-stage stochastic programming problem with continuously distributed parameters:

max{-cx + Q(x) | 0 ≤ x ≤ 1}, where Q(x) = E q n-1 max{βy -n-1 i=1 q j y j | 0 ≤ y ≤ x, y j ≥ α j y (j = 1, . . . , n)} .
For any realisation q n-1 of q n-1 , the optimal solution of the second stage problem is y j = α j y ≥ 0, (j = 1, . . . , n -1), so

Q(x) = E q n-1 max{(β -n-1 i=1 q j α j )y | 0 ≤ y ≤ x} . Thus y = x if n-1
j=1 q j α j ≤ β, and y = 0 if n-1 j=1 q j α j > β. Notice that, since α n > β, the constraint w n ≤ 1 is redundant in P in Claim 2. Define ϕ(w n-1 ) = (β -n-1 j=1 α j w j )/β. Thus ϕ(w n-1 ) ≤ 1 in P as defined in the proof of Claim 1. In P this implies that α n w n /β ≤ ϕ(w n-1 ). Now, the solution value of the two-stage problem can be written as max βΦ -c x | 0 ≤ x ≤ 1 , where Φ = E q n-1 max(ϕ(q n-1 ), 0) . Thus x = 1, with optimal objective value βΦ -c, if Φ > c/β, and x = 0, with optimal objective value 0, otherwise. Therefore, computing the optimal solution to the stochastic program requires determining whether or not Φ = E q n-1 max(ϕ(q n-1 ), 0) > c/β. Claim 3. Computing Φ is P-hard.

Proof. Since, from Claim 2 above, computing Vol(P ) is P-hard, this claim follows from Φ = E q n-1 max(ϕ(q n-1 ), 0) = P ϕ(w n-1 ) dw n-1 , since 0 ≤ ϕ(w n-1 ) ≤ 1 in P ,

= P ϕ(w n-1 ) 0 1 d(α n w n /β) dw n-1 , = (α n /β) P 1 dw n dw n-1 , since 0 ≤ α n w n /β ≤ ϕ(w n-1 ) in P, = (α n /β)Vol(P ).
Now, from [START_REF] Dyer | On the complexity of computing the volume of a polyhedron[END_REF], it follows that Vol(P ) is a rational number with known denominator, when P is a knapsack polytope. Hence we can compute Vol(P ) exactly if we can compute it to a close enough approximation. Hence there must exist values of c for which deciding Φ > c/β is P-hard. Otherwise we could use bisection to compute the P-hard quantity Vol(P ) exactly, in a polynomial number of iterations, a contradiction.

Showing that this problem is in P P would require additional conditions on the input distributions. We note that a result of Lawrence [START_REF] Lawrence | Polytope volume computation[END_REF] implies that exact computation may not even be in pspace.

( i )

 i Counting knapsack solutions. Compute N (P ) = |{0, 1} n ∩ P |.


Arne C Reimers
email: arne.c.reimers@gmail.com

Frank J Bruggeman

Brett G Olivier

Leen Stougie

Fast Flux Module Detection using Matroid Theory

 , that the space of optimal-yield fluxes decomposes into flux modules. These decompositions allow a much easier but still comprehensive analysis of the optimal-yield flux space.

Using the mathematical definition of module introduced by Müller and Bockmayr (2013b) , we discovered useful connections to matroid theory, through which efficient algorithms enable us to compute the decomposition into modules in a few seconds for genome-scale networks. Using that every module can be represented by one reaction that represents its function, in this paper, we also present a method that uses this decomposition to visualize the interplay of modules. We expect the new method to replace flux variability analysis in the pipelines for metabolic networks.

Introduction

The metabolic capabilities and behaviors of biological cells are often modeled using metabolic networks. A metabolic network is constituted of a set of chemical compounds and a set of reactions describing the possible transformations of compounds. In the last years it became possible to reconstruct such networks on the genome-scale. This means that on one hand nearly all the reactions that can happen in a biological cell are included. On the other hand such networks consist of thousands of reactions.

Constraint based methods have proven to be very successful in the analysis of metabolic networks [START_REF] Papin | Comparison of network-based pathway analysis methods[END_REF][START_REF] Price | Genome-scale models of microbial cells: evaluating the consequences of constraints[END_REF]. In constraint based methods no detailed information on reaction kinetics is needed. Often, the knowledge of reaction stoichiometries is sufficient. Rows of the stoichiometric matrix correspond to metabolites and columns to reactions: s ij the i, j-th entry of S is the number of molecules of compound i consumed (s ij < 0), produced (s ij > 0), or not involved (s ij = 0) in reaction j. In steady state this results in linear constraints that express flow conservation on all internal metabolites, with possibly lower-and upper bound on fluxes, yielding a polyhedron of feasible flux-vectors P = {v : Sv = 0, ℓ ≤ v ≤ u}.

Among the most prominent analysis methods is flux balance analysis (FBA) [START_REF] Varma | Metabolic flux balancing: Basic concepts, scientific and practical use[END_REF][START_REF] Orth | What is flux balance analysis[END_REF][START_REF] Santos | A practical guide to genome-scale metabolic models and their analysis[END_REF]. It is, for example, used to compute the optimal biomass yield that can be achieved by a cell under some growth medium [START_REF] Feist | The biomass objective function[END_REF]. This amounts to solving a linear programming problem of the form max{cv : v ∈ P }, where cv is the linear function expressing the (weighted) amount of biomass. In general such optimal flows are not unique [START_REF] Mahadevan | The effects of alternate optimal solutions in constraintbased genome-scale metabolic models[END_REF]. If this is ignored, it can lead to wrong predictions of by-product flux rates [START_REF] Khannapho | Selection of objective function in genome scale flux balance analysis for process feed development in antibiotic production[END_REF]. [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF] showed that many reactions have fixed flux rate in all optimal solutions. These are determined by flux variability analysis (FVA) [START_REF] Burgard | Minimal reaction sets for escherichia coli metabolism under different growth requirements and uptake environments[END_REF][START_REF] Mahadevan | The effects of alternate optimal solutions in constraintbased genome-scale metabolic models[END_REF]. The remaining variability is due to variability of the fluxes on a number of relatively small subnetworks, which we call flux modules. As an example we use here an artificial network similar to the one presented by Müller and Bockmayr (2013b) in Figure 1.

[Figure 1 about here.] In the example, all stoichiometric coefficients are supposed to be +1, -1 or 0. Assuming an input flux rate 1 of the nutrient, any optimum outputs rate 2 of biomass. The continuous hyperarcs represent the reactions that carry fixed flux in any optimal solution. The various dashed sub-hypernetworks indicate the variability present in various optimal solutions. However, for any of the dashed, dash-dotted or dash-dot-dotted hypernetworks, we notice that the net influx and the net outflux is the same in every optimal solution: e.g. in every optimal solution, the dash-dotted subnetwork consumes 1 unit of metabolite m 4 and produces 1 unit of metabolite m 6 , but there is flexibility in which route this unit flow goes through the dash-dotted network. This allows to see a module as one sort of aggregated reaction, see Figure 2.

[Figure 2 about here.]

The fixed input and output compounds of a subnetwork characterizes the notion of flux-module (Müller and Bockmayr, 2013b) in a mathematically rigorous way. Müller and Bockmayr (2013b) showed that every optimal yield elementary flux mode (EFM) [START_REF] Schuster | On elementary flux modes in biochemical systems at steady state[END_REF]) is

INTRODUCTION

a concatenation of reactions with fixed flux and an elementary mode of each of the flux modules. This is illustrated in Figure 3. There are two ways to go through the dashed module, three ways to go through the dash-dotted module, and two ways to go through the dash-dot-dotted module. Hence 7 sub-paths suffice to define the 12 elementary optimal flux modes. Clearly in large networks this combinatorial explosion can be much more dramatic.

[Figure 3 about here.] While the method by [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF] required the enumeration of exponentially many vertices of a flux polyhedron (which are related to the optimal yield EFMs), Müller and Bockmayr (2013b) showed a way to find the modules without needing to compute all extreme solutions. Their method however relied on many runs of FVA. Although faster than EFM enumeration, the method is very sensitive to numerical instabilities and analyses of genome-scale networks could still take several hours.

The most important result in this paper is an extremely simple method allowing to compute the flux-modules in a few seconds for genome-scale metabolic networks. The method, described in Section 2, is based on the observation that the modules correspond to the separators of the linear matroid defined by the columns of the stoichiometric matrix that belong to reactions with variable optimal flux. We will explain all these technical concepts in Section 2.1. The efficiency of our method is demonstrated in Section 3 by application to several genome-scale metabolic networks.

Flux modularity highly depends on the growth conditions. In particular, interesting flux modules can usually only be found in the optimal flux space. Hence, it is of high importance to understand how the decomposition of modules changes under different growth conditions and objective functions. Since with our new method, module computation has become so fast, we can simply compute and compare modules under many different growth conditions and compare the results. Essential for this is a visualization method that shows the interplay of modules in the context of the whole network. In Sec. 2.3 we present a method that automatically generates such a visualization using a clever compression based on flux modules. Results of that method applied to a set of genome-scale metabolic networks can be found in Sec. 3.2.

Methods

Definitions and Preliminaries

We use M to denote the set of metabolites, R to denote the set of reactions. We abuse the notation for sets also for their size. S ∈ R M×R denotes the stoichiometric matrix. By appropriate remodeling the polyhedron P introduced before can be rewritten as P ⊆ {v ∈ R R : Sv = b}. We observe that b = 0 leads to the standard steady-state assumption. Here, we also allow b = 0 to simplify notation in the context of modules. Furthermore, the space of optimalyield fluxes is again a polyhedron and can be written in this form, too (Müller and Bockmayr, 2013b). We will show that we can reduce the analysis of P to the analysis of flux spaces defined by the kernel of S: ker(S) := {v ∈ R R : Sv = 0}. We use v r to denote flux through reaction r.

The support of flux-vector v is denoted by supp(v) := {r ∈ R : v r = 0}.

We will conduct FVA for each reaction r of the network by solving the following two linear optimization problems, yielding, respectively, the minimal and maximal possible flux rate: max / min{v r : Sv = 0, ℓ ≤ v ≤ u} We will also be interested in the flux through a subset of reactions A ⊆ R. Hence, we write v A to denote the components of v corresponding to the reactions in A and we use S A to denote the stoichiometric matrix that only contains the columns corresponding to the reactions in A. We define the projection pr A (P ) := {v A : v ∈ P }.

Definition 1 [Flux Module, Müller and Bockmayr (2013b)] A ⊆ R is a P -module if there exists a d ∈ R M s.t. S A v A = d for all v ∈ P . We call d the interface flux of the module.

✷

In contrast to the definition in Müller and Bockmayr (2013b), we also allow A = ∅ to be a module, which together with R we call the trivial modules. We present here some useful properties of modules proven in Müller and Bockmayr (2013b). They may also help the reader to get some intuition on the concept of module.

Proposition 1 Properties of Modules.

(i) If disjoint sets A and B are P -modules then A ∪ B is a P -module;

(ii) If A and B are P -modules and B ⊂ A then A \ B is a P -module.

The rest of this section is devoted to an introduction to the relevant concepts from Matroid Theory [START_REF] Oxley | Matroid Theory[END_REF], which is a generalization of graph theory and linear algebra. A matroid is defined by a universe of elements and subsets of them that have some independence structure. Definition 2 Given a universe U and a family A of independent subsets of U. Then {U, A} is a matroid if it satisfies the following conditions.

• ∅ ∈ A; • If A ∈ A and A ′ ⊂ A, then A ′ ∈ A;
• If A, A ′ ∈ A and A ′ contains more elements than A, then there exists an element e ∈ A ′ \A, such that A ∪ {e} ∈ A.

✷ As a very relevant example, a set of vectors in R R , together with their linearly independent subsets form a matroid; a so-called linear matroid. Matroid theory has already been used in the past to describe metabolic networks [START_REF] Oliveira | An algebraiccombinatorial model for the identification and mapping of biochemical pathways[END_REF][START_REF] Beard | Thermodynamic constraints for biochemical networks[END_REF]. Indeed, many concepts from metabolic networks also exist in matroid theory. For example, if all reactions are reversible (P = ker(S)) then flux modes in metabolic networks correspond to cycles in matroid theory; i.e., dependent sets of a matroid. Elementary flux modes correspond to circuits; i.e., minimal dependent sets. Notice that in matroid theory we only talk about the support. I.e., A ⊆ R is a cycle if and only if there exists a flux mode v ∈ ker(S) with A = supp(v). Similarly, a circuit C ⊆ R is a cycle with minimal support.

Matroid theory inherits many powerful concepts from linear algebra like duality and rank (which is important for the proofs in the supplementary material). Also graph theory introduces some further useful concepts into matroid theory. Important for us is the notion of a connected component of a matroid: two elements of a matroid are in the same component if there exists a circuit that contains both. We notice that in graph theory this property characterizes a 2connected component. A separator of the matroid is now any union of connected components, i.e., any of the two sides of a partition of the matroid into two parts A, B such that there exists no circuit intersecting A and B. In Sec. 2.2 we show how the flux modules of a metabolic network correspond one-to-one to the separators of the corresponding matroid. We then use matroid theory to derive a very fast and simple algorithm for finding modules. It is based on a result by [START_REF] Krogdahl | The dependence graph for bases in matroids[END_REF]. The runtime results on a set of genome-scale metabolic networks are presented in Sec. 3.1.

Finding Modules Efficiently

We first show that it is sufficient to analyze modularity as a local property of one point in the inside of the flux space, implying that we can ignore reaction reversibilities and simply analyze a subvector-space (Thm. 1). This allows to describe modularity in terms of matroid separators (Thm. 2), which we then exploit in designing an efficient algorithm to compute modules.

To make the first step, consider a point x inside the flux space and a neighbourhood of it (Fig. 4).

[Figure 4 about here.] This neighbourhood captures all the characteristics needed to analyse modularity of the whole flux space. We only have to deal with the term "inside". Since P ⊆ {v ∈ R R : Sv = b}, it follows that P is of lower dimension in R R . Hence, we will only consider the interior relative to ker(S). However, if we have reactions with fixed flux rate, P will also have lower dimension than ker(S). Therefore, we will restrict to reactions with variable flux rate, which we define by:

V := {r ∈ R : v max r = v min r }, where (1) v max r := sup{v r : v ∈ P } v min r := inf{v r : v ∈ P }
This restriction does not destroy the module property:

Observation 1 It holds for all A ⊆ V that A is P -module ⇔ A is pr V (P )-module.
To guarantee that we can find a x inside the flux space after we restricted to reactions with variable flux rate, we require that P is convex. In a future work we will consider the case of non-convex flux spaces.

Theorem 1 If P ⊆ {v ∈ R R : Sv = b}, is convex, it holds for all A ⊆ R A is P -module ⇔ A ∩ V is ker(S V )-module.
The proof can be found in the supplementary material.

By Thm. 1 we can restrict our attention to the analysis of linear vector spaces. Hence, in the following we will only analyse polyhedra of the form P = ker(S). We will relate modules of ker(S) to separators of the matroid defined by the columns of S. Remember the explanation of a separator in a graph in terms of the non-existence of a flow circulation in Section 2.1 and observe, that every module in ker(S) also has interface flux 0 since 0 ∈ ker(S).

Formally, we obtain the following theorem, the proof of which is deferred to the supplementary material.

Theorem 2 A ⊆ R is a ker(S)-module if and only if A is a separator in the matroid represented by S.

✷

The characterization of modules as separators of matroids allows to compute the flux-modules of a metabolic network efficiently. Since separators and modules are closed under disjoint union, it suffices to describe the set of minimal nontrivial separators (modules).

Definition 3 (Minimal Module) A P -module ∅ = A ⊆ R is called minimal if there exists no P -module B = ∅ with B ⊂ A.

✷

To understand the algorithm for finding the modules, we observe that the minimal non-trivial separators are the connected components of the matroid. Formulated in matroid-terminology we recall from Section 2.1 the following characterization of connected component: For any 2 elements (columns of S in the linear matroid, edges in the graph) in the same connected component there exists a circuit that contains them both. For pairs of elements of different connected components this is not true.

We can now build a graph G = (V, E), where V is the set of reactions defined in (1) and there is an edge between two reactions (columns of S V ) if and only if there exists a circuit that contains both. The connected components (in the graph-theoretic sense) of G will be the minimal separators. However, as the number of circuits explodes exponentially, it is not efficient to enumerate all circuits in order to compute the connected components of the graph G. Indeed, this is also not necessary and it suffices to look at a special set of circuits, so called fundamental circuits [START_REF] Truemper | Partial matroid representations[END_REF].

A set of fundamental circuits is obtained as follows: We start by finding a maximal independent set (also called basis) X of the matroid, which we compute by Alg. 1. Notice that, starting from the empty set, the algorithm grows X by adding elements only if this keeps X independent. Since we try to add all elements to X, it follows that at the end of the algorithm, X will be a basis of the linear matroid represented by S V .

Let Y := V \ X. Clearly, for every r ∈ Y , adding r to X will create a cycle C r ⊆ X ∪ {r}. It is easy to see that C r is actually a circuit, which is called fundamental circuit. In Alg. 1 the fundamental circuits are constructed simultaneously with constructing X. This gives us a so-called partial representation.

We now build, by Alg. 2, the graph G ′ = (V, E ′ ), where two reactions are connected by an edge if there exists a fundamental circuit that contains both. Krogdahl and Cunningham showed that the connected components of G ′ , found by Alg. 2, are precisely the minimal separators of the matroid [START_REF] Cunningham | A combinatorial decomposition theory[END_REF][START_REF] Krogdahl | The dependence graph for bases in matroids[END_REF].

To each circuit C there exists a flux vector v that is unique up to scaling with C = supp(v), Sv = 0. If we enter for every fundamental circuit the corresponding flux vector as a column into a matrix, we obtain a null-space matrix of S. Hence, this approach can be understood as computing a block-diagonalization of the null-space matrix. Approaches like this in the context of stoichiometric matrices have already been studied in [START_REF] Schuster | Detecting strictly detailed balanced subnetworks in open chemical reaction networks[END_REF]. However, [START_REF] Schuster | Detecting strictly detailed balanced subnetworks in open chemical reaction networks[END_REF] does not use matroid theory and it is unclear whether their method will always compute the finest block-diagonalization.

Algorithm 1 Computes a basis X and its set of fundamental circuits of a matroid represented by S function ComputePartialRepresentation(S)

C = ∅ X = ∅ for r ∈ V do check feasibility of S X v = -S r if feasible then C := supp(v) ∪ {r} C := C ∪ {C} else X := X ∪ {r} end if end for return C Algorithm 2 Computes the modules of {v : S V v = 0} function ComputeModules() C = ComputePartialRepresentation(S V ) Build Graph G = {V, E} with (x, y) ∈ E iff there exists C ∈ C with x, y ∈ C. A = find connected components of G (e.g.

using depth-first search). return A

Here we recapitulate all the steps for finding the modules of the optimal flux space of a metabolic network.

1. Determine the optimal value by LP;

2. Set the objective function equal to the optimum value and add it as a constraint;

3. For each reaction r maximize and minimize the flux through r in the optimal flux space; 4. Determine the set V of reactions for which the maximum and the minimum are not equal; 5. Select the set of columns S V corresponding to V of the stoichiometric matrix S and neglect the non-negativity constraints; i.e., irreversibilities, directions of the reactions;

6. Apply Alg. 2 to compute the minimal modules A of {v ∈ R V : S V v = 0}.

7.

A is the set of minimal modules that contain reactions in V . The reactions with fixed flux are all minimal modules by themselves.

We notice that steps 3 (and therefore 4) of the algorithm can be parallelized in a trivial way, reducing the computation times even further.

Visualization

We develop a visualization tool to help us understand how the decomposition of modules changes under different growth conditions and objective functions. By the definition of module, the reactions inside a module have together a fixed function (the interface flux). Hence, we can represent the module by a single reaction with a fixed flux in the genome-scale network. The stoichiometry of the representing reaction is precisely the interface flux of the module.

This way we can create a compressed network that contains all the reactions with fixed flux rates and artificial reactions that represent the modules. This compressed network has the following advantages:

• The number of reactions carrying flux is compressed (a module with many reactions, is represented by a single reaction).

• All the reactions in the compressed network have a fixed flux rate.

Unfortunately, the number of fixed reactions is still very large. This prevents automatic visualization of the network and the role of the modules containing variable reactions is obfuscated. However, reactions that have a fixed flux rate can also be grouped together into modules by Prop. 1.

Theoretically, we could group all reactions with a fixed flux rate into 1 module. This would result in a compressed metabolic network consisting of k + 1 reactions, where k is the number of minimal modules containing reactions with variable flux rates. In particular, the module containing all fixed reactions will likely also contain the biomass-and nutrient-uptake reactions.

If we want to understand the role of the modules for biomass production or nutrient uptake, this is not very useful. Moreover, modules of variable reactions may disconnect reactions with fixed flux rates from each other. Such disconnected reactions are important for the mediation between modules and should also be displayed separately. Hence, we decided to build a compressed network as follows:

1. Given: A collection Mod of interesting modules (selected by the user). Mod has to cover all reactions with variable flux rates. Typically Mod contains all minimal modules of variable reactions, a module containing the biomass reaction and modules containing the nutrient uptake reactions.

2. We compute the set R Mod := {r ∈ R : r ∈ M ∃ M ∈ Mod} of reactions in interesting modules.

3. We compute the set R B := {r ∈ R \ R Mod : v r = 0 ∀v ∈ P } of blocked reactions.

4. We compute the set M Mod := {m ∈ M : ∃ r ∈ R Mod such that m ∈ supp(S r )} of metabolites involved in the interesting modules.

5. We consider the metabolic network, where R Mod , R B and M Mod are removed. It is represented by the stoichiometric matrix S ′ := S M\MMod,R\(RMod∪R B ) .

6. We compute the connected components Mod F of S ′ . We do so by defining the incidence matrix of a bipartite graph, the nodes of which on one side of the bipartition correspond to the rows of S ′ , and the ones on the other side to the columns of S ′ , and there is an edge between row-node i and column-node j if and only if S ′ ij = 0. The column-nodes represent the reactions in R\(R Mod ∪B), and the corresponding reactions of the connected components of this bipartite graph, whence Mod F , forms a partition of R \ (R Mod ∪ B). Clearly, every A ∈ Mod F is a module, since Mod F only contains fixed reactions.

7. We represent each module in Mod, Mod F by a single reaction with the corresponding interface flux. Let M 0 be the set of metabolites that have a net interface flux of 0 in all these modules. We suppress M 0 , since they would just show up as isolated metabolites. We obtain a metabolic network with metabolites M ′ := M \ M 0 and reactions

R ′ := Mod ∪ Mod F .
8. We remove reactions disconnected from the network that contain the target reaction, e.g. because of modules that form thermodynamically infeasible cycles or otherwise have no role in the metabolism.

In practice, this results in medium-scale networks that can automatically be visualized with graph-drawing software like GraphViz [START_REF] Gansner | An open graph visualization system and its applications to software engineering[END_REF].

3 Results

Runtime of Module Finding

With the new method we can compute all flux modules for the optimal flux space of genome scale networks in about the same time as is needed for conventional flux variability analysis. In Table 1 we see that the new method using matroid theory outperforms the previous methods in orders of magnitude. We used the metaopt toolbox (Müller and Bockmayr, 2013a) to solve the flux variability subproblems. Unfortunately, we did not have access to all the runtime data of [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF] which is why some of the data is missing and the reported runtimes may be only from some steps in the pipeline. The computations for the matroid approach were obtained by computations on a 4-core desktop computer.

[Table 1 about here.]

In particular notice that large networks like Human recon 2 can now also be analyzed. In addition, the new method is numerically much more stable. In the method introduced by Müller and Bockmayr (2013b) it often happens that error tolerances are chosen too small or too large, which causes that linear programs that should be feasible are detected as infeasible etc. This then usually caused the algorithm to abort and the tolerance sometimes needed to be adjusted according to the problem instance.

We experienced that the new matroid based method is much more robust in this respect. Our initial tolerances of 10 -20 for the optimization step, 10 -8 for the flux variability and 10 -9 for the final module computation worked in all cases.

Note, that the other two methods are solving slightly different problems. In Müller and Bockmayr (2013b) we were actually looking for modules in the thermodynamically constrained flux space and in [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF], rays and linealities are eliminated prior to module computation.

A comparison between the results of Müller and Bockmayr (2013b) and the new method on E. coli iAF1260 revealed that 7 of the modules coincide, 2 modules from the new method contain additional reactions (which have fixed flux under thermodynamic constraints). The remaining modules are computed by the new method, but not by Müller and Bockmayr (2013b) since they again only contain reactions that have fixed flux by thermodynamic constraints (usually those modules are formed by a splitted pair of forward and backward reactions). The differences seem to be small, but a detailed analysis will be subject to future work.

Visualization

We used the visualization method presented in Section 2.3 to create visualizations of the above mentioned genome scale networks. The results can be found on the supplementary website. In Tab. 2, we compare the original size of the networks with the size of the compressed networks that are used to visualize the interplay of the flux modules with variable flux rates. Each reaction of the compressed network is a flux module. Every minimal flux module containing reactions with variable flux rates is represented by exactly one reaction. Reactions with fixed flux rate are grouped together. It is interesting to see that although the networks have quite different sizes originally, the compressed sizes do not vary very much.

Discussion

Enumeration of Optimal-Yield Pathways

We showed that flux modules [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF]Müller and Bockmayr, 2013b) of genome-scale metabolic networks can be computed efficiently using matroids. We confirmed the previous results that the optimal flux space of most genome-scale metabolic networks decomposes into modules. If we want to compute the set of all optimal yield elementary modes, we theoretically can do this by simply computing the optimal yield elementary modes for each module. Then, we can use the decomposition theorem of Müller and Bockmayr (2013b) and obtain all optimal yield elementary modes of the whole network. A small numerical barrier in practice is that EFM enumeration for each module appears to be numerically very unstable. Hence it is likely that EFMs are missed if not everything is computed using precise rational arithmetic.

We noted that the previous methods [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF]Müller and Bockmayr, 2013b) were computing flux modules on slightly different flux spaces (in [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF] rays and linealities were removed, in Müller and Bockmayr (2013b) we worked on the thermodynamically feasible flux space). These differences seem to be small but could be of significant biological importance. For example, it could be that due to thermodynamic constraints a reaction is blocked and hence, we can refine the modules. In a follow up work we will (mathematically and empirically) analyse the impact of these differences. Also, we want to point out here that, for computing modules, the method by [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF] has to enumerate all the extreme points of the flux polyhedron of optimal fluxes (after some preprocessing), a much harder task. As a result more information than modules is obtained.

The full flux space is usually not decomposable into modules. In a follow up paper we will generalize the notion of module. This will allow us to find interesting modules also for the full flux space. Furthermore, this will have the potential to derive similar decomposition theorems as in Müller and Bockmayr (2013b) that then will work on the full flux space as well. We think this will be a major step towards EFM enumeration of genome-scale networks.

Modularity under different Growth-Conditions

It has been observed that the decomposition into modules depends on the growth condition [START_REF] Kelk | Optimal flux spaces of genome-scale stoichiometric models are determined by a few subnetworks[END_REF]Müller and Bockmayr, 2013b). If we want to understand how the optimal flux space changes if the growth condition is modified, we have to recompute the decomposition into modules. Previously, this was a tedious task. Now it is very simple and fast and it can be done even for very small changes.

We presented a visualization method that shows the interplay of the modules and how they contribute to optimal biomass production. We think that this visualization will be very helpful to detect when a change in a growth condition significantly changes the structure of the optimal flux space.

For the visualization we use the definition of module to lump reactions together. This way we compute a compressed metabolic network that shows the optimal flux distribution with only a small number of reactions. These networks were small enough to be visualized using automated graph drawing tools. Currently, we have only little control on how these networks are drawn, causing the visualization to seem to be very sensitive to changes. In particular it would be interesting if we could get more robust drawing results for small changes in the network.

Conclusion

In this paper we demonstrated the power of matroid theory for metabolic network analysis and used it to present a new method that allows us to compute flux modules very efficiently. This allows us to compute flux modules of many metabolic networks under a large set of different conditions to compare flux modules with existing classical metabolic subsystems like Glycolysis.

Compared to classical metabolic subsystems that, at worst, are arbitrary functional groupings of metabolic reactions/species, flux modules are mathematically well defined. They are structural features only depending on a defined set of conditions (inputs, optimality). This qualifies them as a performance and quality metric for genome-scale metabolic networks. Furthermore, it allows us to investigate the modularity, and simplify genome metabolic networks without the risk of a bias from conventional biological interpretation. unused for optimal yield r glc r bio r byproduct Figure 1: Toy example network. All stoichiometric coefficients are +1, -1 or 0.We assume a fixed nutrient uptake rate of 1 through r glc . For optimal biomass production (flux through r bio ) this implies that no sideproduct is produced (flux through r byproduct is fixed to zero) and a optimal biomass production of 2 is achieved. The continuous hyperarcs represent reactions carrying fixed flux in all optimal solutions. The modules of the network are marked with different dash-styles. Since the reaction r byproduct carries no flux in any optimal solution, it is greyed out. Figure 3: Visualization of all 12 optimal-yield EFMs of the toy network (Fig. 2). By taking one EFM through each module together with the reactions of fixed non-zero flux we obtain an optimal yield EFM of the original network. Furthermore, all optimal yield EFMs of the original network can be obtained this way. For each EFM of a module, the used reactions are marked in black, while the unused are marked in grey. For each of the genome-scale networks a compressed network representing the optimal-yield flux space was computed by compressing flux-modules and sets of reactions with fixed flux to single reactions. All reactions in the compressed network have a unique flux and the metabolites display the interactions between the flux-modules.

  FIGURES

  FIGURES

Figure 4 :

 4 FIGURES

TABLES TABLES Table 1 :

 TABLES1 Comparison of runtimes for computing modules in the optimal flux space of genome scale networks. Network Kelk et al. (2012) Müller and Bockmayr (2013b) using matroids

	TABLES			TABLES
	E. coli iAF1260	133495sec	755sec	6.4sec
	E. coli iJR904	1906sec	162sec	1.9sec
	E. coli iJO1366			8.4sec
	H. pylori iIT341		55.5sec	0.8sec
	H. sapiens recon. 1			153.3sec
	H. sapiens recon. 2			1131sec
	M. barkeri iAF692	1088sec	941sec	1.4sec
	M. tuberculosis iNJ661	9317sec	1623sec	4.3sec
	S. aureus iSB619		127.8sec	1.2sec
	S. cerevisiae iND750			3.0sec

Table 2 :

 2 Size of the compressed networks. No. Metabolites No. Reactions No. Metabolites No.

	Reactions

Acknowledgments

We thank Timo Maarleveld for providing us with runtime data of the method by Kelk et al. (2012). This work was supported by the Berlin Mathematical School in form of a PhD-Stipend and by the Tinbergen Institute.

Author Disclosure Statement

No competing interests exist.

Authors Contributions

AR and LS found the connection to matroid theory and took care of mathematical correctness. AR provided all detailed proofs. LS supervised and organised the collaboration process. AR, FB and BO developed the visualization method. AR, LS and FB worked on the manuscript. All authors read and approved the manuscript.

List of Figures