Boltzmann oracle for combinatorial systems.

Boltzmann random generation applies to well-defined systems of recursive combinatorial equations. It relies on oracles giving values of the enumeration generating series inside their disk of convergence. We show that the combinatorial systems translate into numerical iteration schemes that provide such oracles. In particular, we give a fast oracle based on Newton iteration.

Keywords: Random generation, Boltzmann generation, combinatorics, Newton iteration

1 Introduction

The recent discovery of Boltzmann samplers by Duchon, Flajolet, Louchard and Schaeffer [5] brought a considerable progress to the area of random generation of combinatorial structures. For wide families of classes of structures defined recursively, it is possible to construct automatically efficient random generators (samplers). These generators can produce large structures with the property that a structure of size n is drawn uniformly at random among the cn structures of size n in its class. They rely on so-called oracles, that return numerical values of the generating series of the sequence [([c.sub.n]).sub.n [greater than or equal to] 0] and related series inside their disk of convergence.

In this work, we provide such oracles for numerical values, thus making Boltzmann sampling effective, for all structures specified by well-founded systems of combinatorial equations. The precise setting is given in [section] 2. In the same way that the Boltzmann sampler is constructed automatically from a recursive combinatorial specification, we give automatic constructions of numerical iteration schemes for the oracles. The key point of the treatment is that these numerical schemes are based on combinatorics: each value computed during the numerical iteration is the evaluation of a convergent power series whose coefficients enumerate combinatorial structures.

A straightforward translation of combinatorial systems yields a first family of oracles. For instance, the class of rooted labeled trees is described by the combinatorial equation T = Z x SET(T). This equation can be interpreted as a fixed point equation on combinatorial classes for an appropriate notion of distance (the contact). Then the iteration [T.sup.[n+1]] = Z x SET([T.sup.[n]]) can be seen to converge to the class of rooted labeled trees, starting with [T.sup.[0]] the empty set. For any complex [alpha] with [absolute value of [alpha]] < [e.sup.-1], the corresponding numerical iteration [t.sup.[0]] = 0, [t.sup.[n + 1]] = [alpha]exp([t.sup.[n]]) converges to the value T([alpha]) = -W (-[alpha]) of the generating series of rooted labeled trees (W is the Lambert W function). Indeed, each iterate [t.sup.[n]] is equal to the value [T.sup.[n]] ([alpha]) of the generating series of the class [T.sup.[n]] at [alpha]. That these series are convergent at [alpha] also follows from combinatorics.

Thus our method decomposes into three steps. Starting with a fixed point combinatorial specification (a system of combinatorial equations), we first provide an iteration scheme that converges to the combinatorial solution; then propagate this iteration scheme at the level of generating series and finally at the numerical level, showing that convergence at each stage is a consequence of convergence at the previous one. This is described in [section] 3.

In [section] 4, we design faster oracles based on Newton iteration. For instance, in order to solve the equation t = [alpha]exp(t), Newton's iteration is [t.sup.[n + 1]] = [t.sup.[n]] + ([alpha]exp([t.sup.[n]]) - [t.sup.[n]])/(1 - [alpha]exp([t.sup.[n]])). Classically, this iteration converges to one of the roots of the equation, provided the starting point is chosen correctly. We show that with [t.sup.[0]] = 0, the iteration converges to the root T([alpha]). Moreover, the sequence of iterates is positive increasing so that each iterate is closer to the solution than the previous one. Again, this numerical scheme is based on a combinatorial Newton iteration that converges quadratically to the solution. This iteration originates in work by Decoste, Labelle and Leroux [4] (see also [1, ch. 3]), for the computation of species of structures solution to a univariate functional equation y = H(Z,y). For systems, we have to give a combinatorial meaning to the inverse of the Jacobian matrix. This is handled by the combinatorial "bloomings" introduced by Labelle in his study of the combinatorial Lagrange inversion [10]. Then, as above, this combinatorial iteration transfers to generating series, and quadratic convergence is preserved at the level of valuations. Finally, the numerical version of the iteration yields an oracle in the form of an unconditionally convergent iteration, whose convergence is asymptotically quadratic. Thus, instead of a blind resolution of a system of equations, the iterates directly find their way to the solution that is combinatorially meaningful.

This phenomenon is illustrated in the figure below. A combinatorial structure [C.sub.0] is defined by a recursive combinatorial specification. Next, for each value of z in an interval, we plot the solutions (z, [C.sub.0]) of the system of generating function equations. Only the curve marked in red corresponds to the actual generating function. To the right, we zoom on the black rectangle. The crosses indicate the successive values of the Newton iteration starting from [C.sub.0] = 0.

[C.sub.0] = Z[C.sub.1][C.sub.2][C.sub.3]([C.sub.1] + [C.sub.2]) [C.sub.1] = Z + Z SEQ([C.sup.2.sub.1][C.sup.2.sub.3] [C.sub.2] = Z + [Z.sup.2] SEQ(Z[C.sup.2.sub.2]SEQ(Z))SEQ)([C.sub.2]) [C.sub.2] = Z + Z(3Z + [Z.sup.2] + [Z.sup.2][C.sub.1][C.sub.3])SEQ ([C.sup.2.sub.1])

[GRAPHIC OMITTED]

Further experiments are described in [section] 5.

In a forthcoming article [13], several extensions are presented: the unlabeled case [8]; specifications defining nonempty structures of size 0; a faster Newton iteration computing much fewer intermediate structures; and a fast algorithm for enumeration sequences.

2 Combinatorial structures

This work uses the combinatorial framework described in detail in the recent book by Flajolet and Sedgewick [6], of which we use the notations and definitions. In this section we recall the basic vocabulary on admissible constructions and specifications that let us express those combinatorial systems we are dealing with. Since not all specifications are combinatorially meaningful, we isolate precisely the class of well-founded specifications. These specifications have been considered in the framework of combinatorial classes in [7], and here we show how a combinatorial condition due to Joyal [9] gives an algorithmic characterisation. Finally, we recall what we need from the symbolic method leading to generating series equations.

2.1 Combinatorial specifications

The admissible constructions we consider are given in Table 1 (column 1). They form the vocabulary used to express combinatorial specifications. Columns 3 and 4 show how these constructions act on two important special classes: the empty set ([empty set]) that does not contain any structure and the neutral class [epsilon] consisting of a single structure of size 0.

Definition A combinatorial specification for an m-tuple y = ([y.sub.1], ..., [y.sub.m]) of classes is a system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

each [H.sub.i] denoting a term built from the [y.sub.i]'s and the initial class Z (atomic) using the constructions of Table 1. We refer to this system as y = H (Z,y); boldfaced letters denote vectors. Each non-terminal [y.sub.i] stands for a combinatorial class and the objects derived by this generalized grammar, starting with axiom [y.sub.i], are called [y.sub.i]-structures. The size of a structure is the number of its terminals Z. The number of elements of size n in a vector of classes y is denoted [[absolute vale of y].sub.n].

Tab. 1: Admissible constructions, their derivatives and generating functions. The character "-" stands for "undefined". The translation to generating functions of Cycle and Set applies only in the labeled universe, whereas Union, Product and Sequence can be seen in both labeled and unlabeled universes.

2.2 Well-founded combinatorial systems

2.2.1 Definition

We construct oracles for all combinatorial specifications that make sense in the following way. Definition The combinatorial specification y = H(Z,y) is well founded if and only if, for all n [greater than or equal to] 0, it derives only finitely many structures of size n. This is denoted by [[absolute value of y].sub.n] < [infinity].

2.2.2 Derivative

We now introduce an additional admissible construction, which is classical in species theory [1].

Derivative The derivative of a term H(Z,[y.sub.1], ..., [y.sub.m]) with respect to the non-terminal [Y.sub.j] is denoted by [partial derivative]H/[partial derivative][y.sub.j], or H' when m = 1. It is defined recursively, using the rules given in Table 1. The recursion stops at terminals and non-terminals of the term: if A is either [epsilon], Z or [y.sub.k] with k [not equal to] j, then [partial derivative]A/[partial derivative][y.sub.j] = [empty set], while [partial derivative][y.sub.j]/[partial derivative][y.sub.j] is a new terminal of size 0, denoted by [y.sup.*.sub.j].

The only difference with the pointing operation considered in [6] is at the end of the recursion, where instead of a new terminal of size 0, the pointing operation is equal to the cartesian product of [y.sub.i] with such a terminal of size 0.

Example. The derivative of the term H(Z,y) := Z x SEQ(y) with respect to y is

[partial derivative]H/[partial derivative]y = [partial derivative]Z/[partial derivative]y x SEQ (y) + Z x SEQ(y) x [partial derivative]y/[partial derivative]y x SEQ(y) = Z x SEQ(y) x [y.sup.*] X SEQ(y). (1)

Cartesian product This terminal [y.sup.*.sub.j] is called a bud by Labelle [10]. An induction shows that each non-empty structure of the class [partial derivative]H/[partial derivative][y.sub.j](Z, [y.sub.1], ..., [y.sub.m]) contains exactly one bud. This gives a combinatorial interpretation of the cartesian product with a derivative: for any class C, we have the isomorphism 

[partial derivative]H/[partial derivative][y.sub.i] x C ~ [partial derivative]H/[partial derivative][y.sub.j][|.sub.[y.sup.x.sub.j] = c' (2)

where on the right-hand side, we use the substitution of the one bud by C. From now on, we always use this version of the cartesian product with a derivative of one of the [H.sub.i]'s. The interpretation is that the structures of the cartesian product are obtained by grafting a structure of class C in the place of a former [y.sub.j] in a structure of H, in all possible ways.

Example. The class C can be a derivative itself. Thus, the square of the derivative from (1) is

[[partial derivative]H/[partial derivative]Y).sup.2] = Z x SEQ(y) x (Z x SEQ(y) x [y.sup.*] x SEQ(y)) x SEQ(y).

A typical structure of this class is depicted below: black dots correspond to locations of a Z, blue dots to the y's, the bud [y.sup.*] is marked in red.

[ILLUSTRATION OMITTED]

More generally, a sequence SEQ([partial derivative]H/[partial derivative]Y) consists of trees built up by iterating the previous process.

Example. For any term H, the following inclusion holds:

H(y + c) [contains] H(y) + H' (y) x c. (3)

The interpretation is as follows: the structures on the right-hand side are either in H(y), made up only of y-structures; or in H'(y) x C, structures of H(y) with one of the y-structures replaced by a C-structure. These are clearly disjoint sets, both included in H(y + C). These are the first terms of a complete Taylor formula that has been developed by Labelle [11].

Jacobian matrix As in the classical case, the Jacobian matrix of H(Z, y) with respect to y, denoted by [partial derivative]H/[partial derivative]y is the matrix whose entries are [partial derivative][H.sub.i](Z, y)/[partial derivative][y.sub.j]. The product of a matrix by a vector or by a matrix is obtained by the classical formulas in sums of products forms, sums being interpreted as disjoint unions and products as cartesian products, themselves obtained by grafting at a bud following Eq. (2).

Example. Series-parallel graphs are specified by ([y.sub.1], [y.sub.2]) = SP(Z, [y.sub.1], [y.sub.2]), with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [SEQ.sub.[greater than or equal to] k] is the disjoint union of [SEQ.sub.i] for i [greater than or equal to] k, and similarly for [SET.sub.[greater than or equal to] k]. Linearity of derivative implies that

[partial derivative][SEQ.sub.[greater than or equal to] k](y)/[partial derivative]y = [SEQ.sub.[greater than or equal to] k-1] (y) x [y.sup.*] x SEQ(y) + [y.sup.*] x [SEQ.sub.[greater than or equal to] k-1](y), [[partial derivative][SET.sub.[greater than or equal to] k](y)/[partial derivative]y = [SET.sub.[greater than or equal to] k-1](y)x[y.sup.*].

The Jacobian matrix of our example is therefore

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Jacobian matrix [partial derivative]H/[partial derivative]y, or any of its powers, is a function of the variables Z and y. In particular, evaluating it at (Z, y) = ([empty set], [empty set]) can be interpreted combinatorially as the matrix whose entry (i, j) is the class of structures obtained by replacing one of the [[y.sub.j]'s in [H.sub.i] by a [y.sup.*.sub.j], replacing all the other terminals and non-terminals by the empty set and performing the simplifications indicated in Column 3 of Table 1.

Example. The evaluation of the previous matrix at (Z, [y.sub.1], [y.sub.2]) = ([empty set], [empty set], [empty set]) gives ([MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]).

Nilpotence As in classical linear algebra, a matrix of combinatorial classes is nilpotent if one of its powers is empty (all its entries are the empty set). The order of nilpotence (the minimal power such that emptyness is reached) is bounded by the dimension of the matrix.

Example. For series-parallel graphs, the Jacobian matrix is nilpotent of order 1. An extended system with a third equation [y.sub.3] = [y.sub.1] + [y.sub.2] to define graphs that are either series or parallel has a Jacobian at (Z, [y.sub.1], [y.sub.2], [y.sub.3]) = ([empty set], [empty set], [empty set], [empty set]) with third row ([y.sup.*.sub.1] [y.sup.*.sub.2] [empty set]). This is nilpotent of order 3.

2.2.3 Characterization of well-founded systems

We can now state an effective criterion for a system to be well founded.

Proposition 1 A combinatorial specification y = H(Z, y) such that W([empty set], [empty set]) = [empty set] is well founded if and only if the Jocabian matrix [partial derivative]H/[partial derivative]y([empty set], [empty set]) is nilpotent.

Proof: First, we prove by contradiction that nilpotence implies that there is a finite number of y-structures of any size. Let n be the smallest size for which there is an infinite number of y-structures. By definition, an arbitrary y-structure y of size n decomposes as H (Z,[[gamma].sub.1], ..., [[gamma].sub.m]), where each [[gamma].sub.i] is a [y.sub.i]-structure of size at most n. If none of the [[gamma].sub.i] has size equal to n, there is a finite number of such decompositions. Otherwise, only one of the [[gamma].sub.i] has size n and all the other ones have size 0. The condition H([empty set], [empty set]) = [empty set] implies that the only y-structure of size 0 is [empty set], thus [gamma] is a structure of [partial derivative]H/[partial derivative]y([empty set], [empty set]) x [beta], with [beta] a vector of y-structures of size n. If p is the order of nilpotence of [partial derivative]H/[partial derivative]y ([empty set], [empty set]), this reasoning cannot be iterated more than p times. Thus there is only a finite number of y-structures of size n.

Conversely, let [gamma] be a y-structure of size n, with n = min{k, [[absolute value of y].sub.k] [not equal to] 0}. Suppose that the matrix [partial derivative]H/[partial derivative]y ([empty set], [empty set]) is not nilpotent. Then, for all q [member of] N, there exists a nonempty structure [[beta].sub.q] [member of] [partial derivative]H/[partial derivative]y[([empty set], [empty set])).sup.q] of size 0 so that all [[beta].sub.q] x [gamma] are y-structures of size n, a contradiction.

2.3 Generating functions

The symbolic method gives a systematic translation of combinatorial constructions into generating functions [6, Thm. II.2]. The basic dictionary is given in Table 1. If a calligraphic letter denotes a class, we use the corresponding roman letter for its generating function.

Test for well-founded specifications To check whether the system is well founded, it is sufficient to: (i) translate the combinatorial system into a system over generating series; (ii) compute the Jacobian matrix J; (iii) compute J[(0,0).sup.m], with m the dimension. The system is well founded if and only if this matrix is 0. This improves in simplicity over Zimmermann's algorithm [14].

Definition A combinatorial specification H(Z, y) is called analytic when the generating series H(z, Y) is analytic in (z, Y) in the neighborhood of (0,0), with nonnegative coefficients.

This restriction on combinatorial specifications is necessary for our method. We note that all combinatorial specifications based on the constructions of Table 1 are analytic, by virtue of the analyticity at the origin of +, x, exp, 1/(1 - z), log(1/(1 - z)) and their compositions. Positivity of the coefficients at the origin also follows by induction. Examples of specifications that are not analytic in this sense are the constructions of unlabeled multiset, powerset or cycles.

3 Iteration and oracle computation

Given a well-founded specification y = [PHI](Z, y) with the condition [PHI]([empty set], [empty set]) = [empty set], the natural iteration [y.sup.[n + 1] = [PHI](Z,[y.sup.[n]]) converges to the solution y, starting from the vector [y.sup.[0]] with empty coordinates (see Theorem 1). In this section, we show that this iteration translates at the level of generating functions, and then at the numerical level when the specification is analytic.

3.1 Example

We first illustrate the processing chain on the simple example of general trees, with specification y = Z x SEQ(y), so that [PHI](Z, y) = Z x SEQ(y).

Check well-founded specification In view of (1), [PHI]([empty set], [empty set]) = [empty set] x [epsilon] = [empty set] and [partial derivative][PHI]/[partial derivative]y ([empty set], [empty set]) = [empty set].

Combinatorial iteration The following picture shows the structures produced in the first six iterations of [y.sup.[n + 1]] = Z x SEQ([y.sup.[n]]).

[ILLUSTRATION OMITTED]

Rectangles enclose structures when for this size, all structures of the limit y have been produced: for example iteration [y.sup.[4]] contains all structures of the solution up to size 4; and iteration [y.sup.[5]] contains all structures of the solution up to size 5.

Generating series The specification translates into Y(z) = z/(1 - Y(z)), while the iteration translates into [Y.sup.[n+1]] (z) = z/(1 - [Y.sup.[n]](z)). The following series are produced by this iteration. By construction, they are the generating series of the classes [y.sup.[0]], ..., [y.sup.[5]] from above. Boldfaced numbers show the coefficients that coincide with those of the series solution Y(z).

[Y.sup.[0]](z) = 0 [Y.sup.[1]](z) = z

[Y.sup.[2]](z) = z + [z.sup.2] + [z.sup.3] + [z.sup.4] + [z.sup.5] + [z.sup.6] + [z.sup.7] + [z.sup.8] + [z.sup.9] + ...

[Y.sup.[2]](z) = z + [z.sup.2] + [2z.sup.3] + [4z.sup.4] + [8z.sup.5] + [16z.sup.6] + [32z.sup.7] + [64z.sup.8] + [128z.sup.9] + ...

[Y.sup.[4]](z) = z + [z.sup.2] + [2z.sup.3] + [5z.sup.4] + [13z.sup.5] + [34z.sup.6] + [89z.sup.7] + [233z.sup.8] + [610z.sup.9] + ...

[Y.sup.[5]](z) = z + [z.sup.2] + [2z.sup.3] + [5z.sup.4] + [14z.sup.5] + [41z.sup.6] + [122z.sup.7] + [365z.sup.8] + ...

Numerical values The last step is numerical iteration. In the case of general plane trees, the radius of convergence of the series is 1/4, so that an evaluation at [alpha] = 0.1 is possible. The iteration translates into [y.sup.[n + 1]] = 0.1/(1 - [y.sup.[n]]) and leads to the following values [y.sup.[i]], which again by construction, are the values of the convergent series [Y.sup.[i]](z), when evaluated at z = [alpha], and have for limit Y([alpha]) = 0.11270166537925831148207346002176 ....

[y.sup.[0]] = 0 [y.sup.[1]] = 0.1

[y.sup.[2]] = 0.11111111111111111111111111111111 ...

[y.sup.[3]] = 0.11250000000000000000000000000000 ...

[y.sup.[4]] = 0.11267605633802816901408450704225 ...

[y.sup.[5]] = 0.11269841269841269841269841269841 ...

Thus the translation of the combinatorial specification into a numerical iteration provides an oracle for general plane trees. The generating series are not used but their existence and convergence properties play a role in the proof that the numerical iteration converges to the desired value.

3.2 Transfer of Convergence

We now turn to the general case and provide a first oracle based on iterating the specification numerically. We first define the notions of convergence that are needed. Convergences The notion of contact is classical in species theory [1]. We define its counterpart in the framework of combinatorial classes of [6]. Two combinatorial classes F and G have contact of order k, denoted by F [=.sub.k] G, when their structures of size up to k are identical. A sequence of classes [([y.sup.[n]]).sub.n[member of]N] converges to a class y if for all k [greater than or equal to] 0, there exists N [greater than or equal to] 0 such that for all n [greater than or equal to] N, [y.sup.[n]] [=.sub.k] y. This is denoted by [lim.sub.n [right arrow] [infinity]] [y.sup.[n]] = y. Convergence of vectors of combinatorial classes is defined as componentwise convergence.

Recall that the valuation of a power series S(z), denoted by val(S(z)), is the exponent of the first nonzero coefficient of the series. A metric is classically deduced by defining the distance between two power series by d(F(z), G(z)) = [2.sup.-val(F(z)-G(z))]; convergence follows.

Theorem 1 (Transfer of Convergence) Let y = F(Z, y) be well founded and F([empty set], [empty set]) = [empty set].

1. The iteration [y.sup.[n+1]] = F(Z, [y.sup.[n]]), with [y.sup.[0]] = [empty set], converges to the combinatorial class y, solution of y = F(Z, y).

2. The iteration [Y.sup.[n+1]] (z) = F(z, [Y.sup.[n]] (z)), with [Y.sup.[0]](z) = 0, converges to the generating series Y(z) of the class y.

3. If F is an analytic specification, then Y has positive radius of convergence [rho] and for all [alpha] such that [absolute value of [alpha]] < [rho], the iteration [y.sup.[n+1]] = F([alpha], [y.sup.[n]]), with [y.sup.[0]] = 0, converges to Y([alpha]).

Proof. 1. Combinatorics The first statement is a consequence of Joyal's proof of his implicit species theorem [9, Thm. 6].

It is useful to note that the sequence [y.sup.[n]] is monotonic, in the sense that [y.sup.[n]] [subset] [y.sub.[n+1]] for all n. This is proved by induction. For n = 0, this comes from the definition of the class [empty set]. Then the inclusion [y.sup.[n]] [subset] [y.sub.[n+1]] is preserved by F since the structures of F(Z, [y.sup.[n]]) are naturally structures of F(Z, [y.sup.[n+1]]). By definition of the iteration, this means [y.sup.[n+1]] [subset] [y.sup.[n+2]].

2. Power series By the symbolic method, the power series [Y.sup.[n]] (z) are the generating series of the classes [y.sup.[n]]. Convergence of combinatorial classes translates at the level of generating series into val([Y.sup.[n]] (z) - Y (z)) [right arrow] [infinity], which gives the convergence of generating series.

3. Numerical values Since F is analytic and the Jacobian matrix (I - [partial derivative]F/[partial derivative]Y)(0, 0) is invertible, the implicit function theorem asserts that Y is analytic at 0 (see e.g., [2, Ch. IV]).

The point is to show that for all [alpha] such that 0 [less than or equal to] [absolute value of [alpha]] < [rho], [Y.sup.[n]] ([alpha]) converges to Y([alpha]), and [Y.sup.[n]] ([alpha]) = [y.sup.[n]] (the evaluation of [Y.sup.[n]] at [alpha] is equal to the value obtained by numerical iteration).

The monotonicity of the combinatorial sequence implies the inequality [[z.sup.k]][Y.sup.[n]](z) [greater than or equal to] [[z.sup.k]]Y(z) for all n and k. Consequently, the [Y.sup.[n]]'s are analytic for [absolute value of z] < [rho], by absolute convergence, and also the tails are bounded by the tail of Y. Thus [Y.sup.[n]] ([alpha]) converges to Y([alpha]). Let r be such that [absolute value of [alpha]] [less than or equal to] r < [rho]. Assuming F (z,Y) to be analytic in a polydisc [absolute value of (z, Y)] [less than or equal to] (r, Y(r)) with componentwise inequality, the vector F([alpha], [Y.sup.[n]]([alpha])) is well defined, and thus by induction [y.sup.[n+1]] = F ([alpha], [y.sup.[n]]) = F ([alpha], [Y.sup.[n]]([alpha])) = [Y.sup.[n+1]]([alpha]).

We now prove the required analyticity of F(z, Y). Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and Y(z) = [summation] [c.sub.k][z.sup.k]. For each coordinate h [member of] {1, ..., m}, extracting the coefficient of [z.sup.k] (k = 0, ..., N) in the identity F(z, Y(z)) = Y(z) leads to an inequality of the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where first indices denote coordinates. The coefficients being positive, the first sum converges to [Y.sub.h](r) as N [right arrow] [infinity]. This proves the convergence of [F.sub.h](z, Y) for [absolute value of (z, Y)] [greater than or equal to] (r, Y(r)) and therefore that of F(z, Y) which concludes the proof.

4 Newton iteration

As the computation for general plane trees in [section] 3.1 exemplifies, the numerical convergence of the previous iteration is typically linear. It is therefore very tempting to use a faster iteration, provided we can show that it converges to the desired solution. The aim of this section is to show that Newton's iteration is appropriate. The proof is based on a combinatorial lifting of Newton's iteration. We first show how this works on the same example.

4.1 Newton Iteration for General Plane Trees

Our oracle for general plane trees is simply the Newton numerical iteration for the equation y = [alpha]/(1 - y), with initial point 0, i.e.,

[y.sup.[n+1]] = [y.sup.[n]] + [alpha]/(1 - [y.sup.[n]]) - [y.sup.[n]]/1 - [alpha]/[(1 - [y.sup.[n]]).sup.2], [y.sup.[0]] = 0. (4)

For [alpha] = 0.1, the iterates are

[y.sup.[0]] = 0

[y.sup.[1]] = 0.11111111111111111111111111111111 ...

[y.sup.[2]] = 0.11270125223613595706618962432916 ...

[y.sup.[0]] = 0.11270166537923032259476392887392 ...

[y.sup.[0]] = 0.11270166537925831148207345989331 ...

showing very fast (quadratic) convergence. Note that the equation y = 0.1/(1 - y) has another real positive solution, namely y = 0.88729833462074168852.... We show below that starting with [y.sup.[0]] = 0 ensures convergence to the desired solution.

As in the direct iteration, the proof is based on showing that the values [y.sub.i] are the evaluations of generating series at 0.1, these generating series enumerating combinatorial classes that are defined by a combinatorial lifting of the Newton iteration.

For the case of one equation, such a combinatorial Newton iteration has been introduced by Decoste, Labelle and Leroux [4], who showed that the equation y = F(Z, y) is solved by

[y.sub.[n+1]] = [y.sup.[n]] + SEQ ([partial derivative]F/[partial derivative]y ([y.sup.[n]]) x (F(Z, [y.sup.[n]]]) - [y.sup.[n]]), [y.sup.[0]] = [empty set].

In the case of general plane trees, F(Z, y) = Z x SEQ(y), so that the iteration becomes

[y.sub.[n+1]] = [y.sup.[n]] + SEQ (Z x SEQ [([y.sup.[n]]).sup.2]) x (Z x SEQ([y.sup.[n]]) - [y.sup.[n]]),

which should be compared with Eq. (4). The first classes are as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The fast convergence can also be seen at the level of generating series:

[Y.sup.[0]] (z) = [Y.sup.[1]] (z) = z + [z.sup.2] + [z.sup.3] + [z.sup.4] + [z.sup.5] + [z.sup.6] + [z.sup.7] + [z.sup.8] + [z.sup.9] + ...

[Y.sup.[2]] (Z) = z + [z.sup.2] + 2[z.sup.3] + 5[z.sup.4] + 14[z.sup.5] + 42[z.sup.6] + 131[z.sup.7] + 417[z.sup.8] + ...

[Y.sup.[3]] (z) = z + [z.sup.2] + 2[z.sup.3] + 5[z.sup.7] + 14[z.sup.5] + 42[z.sup.6] + ... + 742900[z.sup.14] + ...

[Y.sup.[4]] (z) = z + [z.sup.2] + 2[z.sup.3] + 5[z.sup.3] + 14[z.sup.5]+ 42[z.sup.6] + ... + 1002242216651368[z.sup.30] + ...

4.2 General Case

Our main result is the following.

Theorem 2 (Newton Oracle) Let y = H(Z, y) be a well-founded analytic specification with H([empty set], [empty set]) = [empty set]. Let a be inside the disk of convergence of the generating series Y(z) of y. Then the following iteration converges to Y([alpha]).

[y.sup.[n+1]] = [y.sup.[n]] + [(I - [partial derivative]H/[partial derivative] ([alpha], [y.sup.[n]])).sup.-1] x (H([alpha], [y.sup.[n]]) - [y.sup.[0]] = 0.

This is the classical Newton iteration for systems. The main point here is that given the initial point 0, not only does the iteration converge, but the limit is the desired value. Our proof relies on an extension of the Newton combinatorial iteration [4] to the case of systems.

Proposition 2 Let y = H(Z, y) be a well-founded specification with H([empty set], [empty set]) = [empty set]. Let [N.sub.H] be the operator defined by

[N.sub.H](Z, y) = y + [[(I- [partial derivative]H/[partial derivative]y (Z, y)).sup.-1] x (H(Z, y) - y).

Then the sequence defined by [y.sup.[0]] = [empty set], [y.sub.[n+1]] = [N.sub.H](Z, [y.sup.[n]]) (n [greater than or equal to] 0) converges to y. Moreover this convergence is quadratic. [y.sup.[n]] [=.sub.k] [??] [y.sub.[n+1]] [=.sub.2k+1] y.

Concerning the combinatorial meaning of the iteration, a few words of explanation are in order. Following Labelle [10], the inverse of the Jacobian matrix is given by:

[(I - [partial derivative]H/[partial derivative]y (Z, y)).sup.-1] = [summation over k[greater than or equal to] 0] [([partial derivative]H/[partial derivative]y (Z, y)).sup.k].

The subtraction H(Z, y) - y is defined when y [subset] H(Z, y) as the set difference of the classes.

Example. For series-parallel graphs, starting with [y.sup.[0]] = ([empty set], [empty set]), the first step of the combinatorial Newton iteration produces the combinatorial class

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For arbitrary n, the numerical Newton iteration is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The singularity is at 2 - [square root of 5] + ln ((1 + [square root of 5])/2) [approximately equal to] 0.245. With [alpha] = 0.24, the first few iterates are

[y.sup.[1]] -(0.1230510663209943063722 ..., 0.06462664750711721439535 ...)

[y.sup.[2]] = (0.1627000389319615796926 ..., 0.09201293266034877734970 ...)

[y.sup.[3]] = (0.1724333307003245710686 ..., 0.09798441803578338336038 ...)

[y.sup.[4]] = (0.1730460965507535353574 ..., 0.09836831514307466499845 ...)

[y.sup.[5]] = (0.1730486392973095133433 ..., 0.09836989917963665326450 ...)

[y.sup.[6]] = (0.1730486393408452105149 ..., 0.09836989920678769126015 ...)

Proof of Proposition 2: The proof is based on the same three steps as that of the combinatorial Newton iteration in the univariate case [4].

1. The iteration is well defined The subtraction is possible only if [y.sup.[n]] [subset] H(Z, [y.sup.[n]]). The proof of this inclusion is by induction. For n = 0 this is a consequence of [y.sup.[0]] being the empty set. If the property is satisfied for n, then we use Eq. (3) with A = [y.sup.[n]] and B = [y.sup.[n+1]] - [y.sup.[n]]. This latter subtraction is justified since the first summand of [N.sub.H](Z, [y.sup.[n]]) is [y.sup.[n]]. Thus we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where in the second line we use the induction hypothesis to rewrite H(Z, [y.sup.[n]]) and the definition of the iteration to rewrite [y.sup.[n+1]].

2. The iteration is not ambiguous All the structures of .[N.sub.H](Z, [y.sup.[n]]) are distinct: this comes from an induction using the fact that the final grafting of an element of H (Z, [y.sup.[n]]) -[y.sup.[n]] cannot occur anywhere else in a structure built on [y.sup.[n]]'s only.

3. Convergence and its quadratic behaviour Assume that [y.sup.[n-1]] [=.sub.k] [y.sup.[n]] and consider an element [gamma] of size at most 2k+1 in the class [y.sup.[n+1]] - [y.sup.[n]]. Thus [gamma] belongs to [([partial derivative]H/[partial derivative]y (Z, [y.sup.[n]])).sup.l x (H (Z, [y.sup.[n]]) - [y.sup.[n]]) for some l [member of] N. Exactly one of the elements of [y.sup.[n]] composing [gamma] has size larger than k: two of them would give [gamma] a size larger than 2k + 1 and none of them would make all these elements of [y.sup.[n]] elements of [y.sup.[n-1]] and therefore [gamma] would belong to [y.sup.[n]]. Moreover, this element, say [beta], belongs to the final H(Z, [y.sup.[n]]) since otherwise this final element, say [delta], would belong to H(Z, [y.sup.[n-1]]) [subset] [y.sup.[n]]. By definition, [beta] itself belongs to [([partial derivative]H/[partial derivative]y (Z, [y.sup.[n-1]])).sup.l]' x (H(Z, [y.sup.[n-1]]) - [y.sup.[n-1]]) for some l' [member of] N. Thus [delta] belongs to [([partial derivative]H/[partial derivative]y (Z, [y.sup.[n-1]])).sup.l'+1] (H(Z, [y.sup.[n-1]]) - [y.sup.[n-1]]) which is part of [y.sup.[n]], and that is not possible. Thus there is no [gamma] of size at most 2k + 1 in [y.sup.[n+1]] - [y.sup.[n]], or equivalently [y.sup.[n+1]] [=.sub.2k+1] [y.sup.[n]].

The proof is concluded by showing that the limit is the solution of y = H(Z, y). This follows from rewriting the iteration in the form

H(Z, [y.sup.[n]]) - [y.sup.[n]] + [partial derivative]H/[partial derivative]y (Z, [y.sup.[n]]) x ([y.sup.[n+1]] - [y.sup.[n]]) = [y.sup.[n+1]] - [y.sup.[n]]

and observing that since [y.sup.[n+1]]] - [y.sup.[n]] converges to [empty set], so does H(Z, [y.sup.[n]]) - [y.sup.[n]].

Proof of Theorem 2: Since the limit of [y.sup.[n+1]] = [N.sub.H](Z, [y.sup.[n]]) is the solution of y = H(Z, y), which is well founded, there are only finitely many elements in [y.sup.[n]] of any size in y and this makes the specification y =.[N.sub.H](Z, y) well founded too. It is analytic by the analyticity of H. The proof is completed by application of part (3) of Theorem 1 with F = [N.sub.H]

5 Applications

Figure 1 displays the timings of a straightforward Maple implementation of the Newton oracle for our example of series-parallel graphs. Each point corresponds to a computation of the oracle at a specified a and precision (i). The curves display jumps corresponding to precisions where the number of Newton iterations is increased by one.

We have also implemented a prototype of an optimized Newton iteration [13]. Our prototype first decomposes the grammar into strongly connected components and each Newton iteration is run on these components separately. Indications on the performance are given in Table 2. In these tests, we generated random grammars for different values of the number of equations and the number of constructions per equation. The average size of the largest strongly connected component of the system is given in the third row, since it has a significant impact on the performance. Timings are given in seconds. In applications to Boltzmann generation, the location of the point a at which the oracle is invoked is related to the expected size of the objects produced by the random generator, which increases with the closeness to the dominant singularity. Thus, we give our timings at points on a scale given by the dominant singularity [rho]. In the last row, we give the average value of the expected size of the objects produced by the Boltzmann generators at 0.999999[rho], where the average is taken over 100 grammars in each column.

[FIGURE 1 OMITTED]

Besides these random examples, our prototype has also been used to compute oracles and generate random structures for systems appearing in real applications. Darrasse has applied the method to generate XML documents [3]: he developed a program that takes any RELAX NG grammar and produces huge uniform random trees that are valid XML documents. These documents can be used, for instance, for testing vulnerability of web services. With context-free grammars such as MathML or DocBook, leading to combinatorial systems of hundreds of equations, he could compute the oracles in a few seconds, and produce random documents of 10,000 nodes in less than one minute on a standard PC.

In a preliminary work with Oudinet, concerned with paths in large graphs modeling concurrent systems [12], we could deal with a first example comprising 1,183 equations in a few seconds, at a point which leads to an expected size of approximately 100,000 nodes.

References

[1] F. Bergeron, G. Labelle, and P. Leroux. Combinatorial species and tree-like structures, volume 67 of Encyclopedia of Mathematics and its Applications. Cambridge University Press, Cambridge, 1998. ISBN 0-521-57323-8. Translated from the 1994 French original.

[2] H. Cartan. Elementary theory of analytic functions of one or several complex variables. Dover Publications, 1995. ISBN 0486685438 (pbk.). URL http://www.loc. gov/catdir/description/dover031/95013507. html. Reprint of the 1973 ed, translated from the 1961 French original.

[3] A. Darrasse. Random XML sampling the Boltzmann way. Technical Report 0807.0992, arXiv, 2008. URL http: //arxiv.org/abs/0807.0992.

[4] H. Decoste, G. Labelle, and P. Leroux. Une approche combinatoire pour l'iteration de Newton-Raphson. Advances in Applied Mathematics, 3:407-416, 1982.

[5] P. Duchon, P. Flajolet, G. Louchard, and G. Schaeffer. Boltzmann samplers for the random generation of combinatorial structures. Combinatorias, Probability and Computing, 13(4-5): 577-625, 2004. Special issue on Analysis of Algorithms.

[6] P. Flajolet and R. Sedgewick. Analytic Combinatorias. Cambridge University Press, 2008. To be published, available from P. Flajolet's web page.

[7] P. Flajolet, B. Salvy, and P. Zimmermann. Automatic average-case analysis of algorithms. Theoretical Computer Science Series A, 79(1):37-109, Feb. 1991.

[8] P. Flajolet, E. Fusy, and C. Pivoteau. Boltzmann sampling of unlabelled structures. In D. Applegate, G. S. Brodal, D. Panario, and R. Sedgewick, editors, Proceedings of the Ninth Workshop on Algorithm Engineering and Experiments and the Fourth Workshop on Analytic Algorithmics and Combinatorias, volume 126 of SIAM Proceedings in Applied Mathematics, pages 201-211. SIAM, 2007. Workshops held in New Orleans, LA, January 2007.

[9] A. Joyal. Une theorie combinatoire des series formelles. Advances in Mathematics, 42(1): 1-82, 1981.

[10] G. Labelle. Eclosions combinatoires appliquees a l'inversion multidimensionnelle des series formelles. Journal of Combinatorial Theory. Series A, 39(1):52-82, 1985. ISSN 0097-3165.

[11] G. Labelle. Derivees directionnelles et developpements de Taylor combinatoires. Discrete Mathematics, 79(3):279-297, 1990. ISSN 0012-365X.

[12] J. Oudinet. Uniform random walks in very large models. In M.-C. Gaudel, J. Mayer, and R. Merkel, editors, Proceedings of the Second International Workshop on Random Testing (RT 2007), pages 26-29, Atlanta, GA, USA, November 2007. ACM Press.

[13] C. Pivoteau, B. Salvy, and M. Soria. Newton iteration for combinatorial systems with applications to enumeration and random generation. In preparation, 2008.

[14] P. Zimmermann. Series generatrices et analyse automatique d'algorithmes. PhD thesis, Ecole polytechnique, Palaiseau, France, 1991.

(i) The program was run using Maple 11 on an Intel processor at 3.2 GHz with 2 GB of memory.

Carine Pivoteau (1) and Bruno Salvy (2) ([dagger]) and Michele Soria (l) ([double dagger])

(1) LIP6 - UPMC, 104 avenue du president Kennedy, F-75016 Paris, France,

(2) Algorithms Project, Inria Paris-Rocquencourt, France [e-mail], [e-mail], Michele. [e-mail]

([dagger]) This work is supported in part by the French Research Agency (ANR Gecko).

([double dagger]) This work is supported in part by the French Research Agency (ANR Gamma), no BLAN07-2-195422.

Tab. 1: Admissible constructions, their derivatives and generating
functions. The character "--" stands for "undefined". The translation
to generating functions of Cycle and Set applies only in the labeled
universe, whereas Union, Product and Sequence can be seen in both
labeled and unlabeled universes.

construction             notation                B =
                                              [empty set]

Disjoint union           A + B                     A
Cartesian product        A x b                [empty set]
Sequence                 SEQ (B)               [epsilon]
Sequence of card k > 0   [SEQ.sub.k] (B)      [empty set]

Cycle                    CYC (B)              [empty set]
Cycle of card k > 0      [CYC.sub.k] (B)      [empty set]
Set                      SET (B)               [epsilon]
Set of card k > 0        [SET.sub.k] (B)      [empty set]

construction                   B =          derivative
                           [epsilon]

Disjoint union            A + [epsilon]     A' + B'
Cartesian product               A           A' x B + A x B'
Sequence                       --           SEQ (B) x B' x SEQ (B)
Sequence of card k > 0         --           [[summation].sup.k-1
                                            .sub.i=0] [SEQ.sub.i]
                                            (B) x B' x [SEQ.sub
                                            .k-1-i] (B)

Cycle                          --           SEQ (B) x B'
Cycle of card k > 0            --           [SEQ.sub.k-1] (B) x B'
Set                            --           SET (B) x B'
Set of card k > 0              --           [SET.sub.k-1] (B) x B'

construction             generating function

Disjoint union           A(z) + B(z)
Cartesian product        A(z) x B(z)
Sequence                 1/(1 - B(z))
Sequence of card k > 0   [B.sup.k](z)

Cycle                    log(1/1 - B(z)))
Cycle of card k > 0      1/k [B.sup.k](z)
Set                      exp(B(z))
Set of card k > 0        1/k! [B.sup.k](z)

Tab. 2: Experimental Results

# equations                4           10                  50

# constructions/eqn       10           10           10           50

avg size largest scc     2.47         3.42         7.95         18.62

time (0.99[rho])         0.05         0.11         0.17         0.47

time (0.999999[rho])     0.08         0.16         0.19         0.56
avg expected size      4.1 [10.     1.4 [10.     2.2 [10.     1.0 [10.
                        sup.14]      sup.7]       sup.5]       sup.5]

# equations                     100                 500

# constructions/eqn       10           50           50

avg size largest scc     10.93        67.18       339.10

time (0.99[rho])         0.23         7.29         61.73

time (0.999999[rho])     0.25         8.11         61.86
avg expected size      1.2 [10.     5.0 [10.     3.3 [10.
                        sup.6]       sup.4]       sup.4]