Constructing, sampling and counting graphical realizations of restricted degree sequences

With the current burst of network theory (especially in connection with social and biological networks) there is a renewed interest on realizations of given degree sequences. In this paper we propose an essentially new degree sequence problem: we want to find graphical realizations of a given degree sequence on labeled vertices, where certain would-be edges are {\em forbidden}. Then we want to sample uniformly and efficiently all these possible realizations. (This problem can be considered as a special case of Tutte's $f$-factor problem, however it has a favorable sampling speed.) We solve this {\em restricted degree sequence} (or RDS for short) problem completely if the forbidden edges form a bipartite graph, which consist of the union of a (not necessarily maximal) 1-factor and a (possible empty) star. Then we show how one can sample the space of all realizations of these RDSs uniformly and efficiently when the degree sequence describes a {\em half-regular} bipartite graph. Our result contains, as special cases, the well-known result of Kannan, Tetali and Vempala on sampling regular bipartite graphs and a recent result of Greenhill on sampling regular directed graphs (so it also provides new proofs of them). The RDS problem descried above is self-reducible, therefore our {\em fully polynomial almost uniform sampler} (a.k.a. FPAUS) on the space of all realizations also provides a {\em fully polynomial randomized approximation scheme} (a.k.a. FPRAS) for approximate counting of all realizations.


Introduction
In the last fifteen years, network theory has been undergoing an exponential grow. One of its more important problems is to algorithmically construct networks with predefined parameters or uniformly sampling all networks with these given parameters. For general background, the interested reader can turn to the now-classic book of Newman, Barabási and Watts ( [14]) or to the more recent book of Newman ([15]).
In this paper we study networks (or graphs as we like to refer them) with given degree sequences. We propose a new degree sequence problem class called restricted degree sequence problem and solve its first instance: we build procedure to decide quickly whether there exists a graph G with a given degree sequence, where G completely avoids a predefined set of forbidden edges, then develop a fast mixing Markov Chain approach to sample almost uniformly all different realizations of the sequence.
We will focus on the technical details of this problem; we do not intend to survey the history of the degree sequence problems in general or their connections to other developments in network theory. One can find detailed information about this specific background in the recent paper of Greenhill ([5]), or in our previous paper [13]. Therefore we will touch only the directly affected definitions and earlier results.
In the next section we discuss some of the known degree sequence problems and introduce our proposed new problem class. In Section 3 we will study in details one instance of this new problem class, describing and proving a fast Havel-Hakimi type greedy algorithm to construct realizations. Then in Section 4 we discuss some known MCMC approaches to sample different kind of degree sequence realizations and state our result about the rapid uniform sampling of this instance of our proposed model. In Sections 5 and 6 we will describe our approach in details. In Section 7 we show that the studied problem is a self-reducible one, therefore our almost uniform sampling method provides good approximate counting of all realizations. Finally in the Appendix we will discuss a necessary technical detail of the required, but very slight generalization of the "Simplified Sinclair's method" introduced in [13].

Degree sequences and restricted degree sequences
Let's fix a labeled underlying vertex set V of n elements. The degree sequence d(G) of a simple graph G = (V, E) is the sequence of its vertex degrees: d(G) i = d(v i ). (Here "simple" means that there are no loops or multiple edges. In general, in cases where multiple edges and/or loops are not excluded the corresponding results are easier. See for example Ryser, [17].) A non-negative integer sequence d = (d 1 , . . . , d n ) is graphical iff d(G) = d for some simple graph G, and then G is a graphical realization of d. The simplest method (a straightforward greedy algorithm) to decide the graphicality of an integer sequence was discovered by Havel ([7]), and was rediscovered independently later by Hakimi ([6]). The validity of their algorithm is proved by the means of the swap operation. It is defined as follows: Let G be a simple graph and assume that a, b, c and d are different vertices. Furthermore, assume that (a, c), (b, d) ∈ E(G) while (b, c), (a, d) ∈ E(G). Then is another realization of the same degree sequence. We call this operation as swap and will denote by ac, bd ⇒ bc, ad.
The analogous notions for bipartite graphs are the following: if B is a simple bipartite graph then its vertex classes will be denoted by U We can define the swap operation for bipartite realizations similarly to (2.1) but we must take some care: it is not enough to assume that (b, c), (a, d) ∈ E(G) but we have to know that a and b are in the same vertex class, and c and d are in the other.
To make clear whether a vertex pair can form an edge in a realization or not we will call a vertex pair a chord if it can hold an actual edge in a realization. Those pairs which cannot accommodate an edge are the non-chords. (Pairs from the same vertex class of a bipartite graph are non-chords.) For directed graphs we consider the following definitions: Let G denote a directed graph (no parallel edges, no loops, but oppositely directed edges between two vertices are allowed) with vertex set X( G) = {x 1 , x 2 , . . . , x n } and edge set E( G). We use the bi-sequence to denote the degree sequence, where d + i denotes the out-degree of vertex x i while d − i denotes its in-degree. A bi-sequence of non-negative integers is called a directed degree sequence if there exists a directed graph G such that (d + , d − ) = dd( G). In this case we say that G realizes our directed degree sequence.
We will apply the following representation of the directed graph G : let B( G) = (U, W ; E) be a bipartite graph where each class consists of one copy of every vertex. The edges adjacent to a vertex u x in class U represent the out-edges from x, while the edges adjacent to a vertex w x in class W represent the in-edges to x (so a directed edge xy is identified with the edge u x w y ). If a vertex has no in-(out-) degree in the directed version, then we delete the corresponding vertex from B( G). (This representation is an old trick, applied already by Gale [4].) There is no loop in our directed graph, therefore there is no (u x , v x ) type edge in its bipartite realization -these vertex pairs are non-chords.
In this paper we will study the following common generalization of all the previously mentioned degree sequence problems: The restricted degree sequence problem d F consists of a degree sequence d and a set F ⊂ V 2 of forbidden edges. The problem, as it is in the original problem, is to decide wether there is a simple graph G on V, completely avoiding the elements of F, which provides the given degree sequence, furthermore to design a uniform sampler of all realizations.
The bipartite restricted degree sequence problem bd F consists of a bipartite degree sequence bd on (U, W ), and a set F ⊂ [U, W ] of forbidden edges. The problem, as it is in the original problem, is to decide wether there is a bipartite graph G on (U, W ) completely avoiding the elements of F, which provides the given degree sequence.
Clearly a bipartite restricted degree sequence problem bd F on (U, W ) is a restricted degree sequence problem It is easy to see that the restricted degree sequence problem is very closely related to Tutte's f -factor theorem ( [19]). However, while Tutte's approach provides a polynomial time algorithm to decide wether a given degree sequence satisfies a d F problem, it does not provide an easy greedy algorithm to do so, furthermore it does not help in the sampling problem.
As we already mentioned the RDS notion is a common generalization of several "classical" degree sequence problems. For instance, the case of bipartite degree sequences is such an example. Another well known example is the case of directed degree sequences: in its bipartite representation, a forbidden 1-factor excludes the directed loops in the original directed graph. Definition 2.1. Let d F be a restricted degree sequence problem and let G be a realization of it. The sequence of vertices C = (x 1 , x 2 , . . . , x 2i ) is a chord-circuit if: (D1) all pair x 1 x 2 , x 2 x 3 , . . . x 2i−1 x 2i , x 2i x 1 are chords; and (D2) each chord occurs only once. A chord-circuit is elementary if (D3) no vertex occurs more than twice; furthermore (D4) when two copies of the same vertex exist, then their distance along the circuit is odd. The chord-circuit C alternates if the chords along C are edges and nonedges in turn (for example x 2j−1 x 2j are edges for 1 ≤ j ≤ i, while the other chords are not edges in G). Deleting the actual edges along C from G and adding the other chords as edges constructs a new graph G ′ which is again a realization of d F . This operation is known as a circular C 2i -swap and denoted by S C . Finally, two different vertices of the alternating chord-circuit C form a potential vertex pair (or PV-pair for short) if (i) this chord does not belong to the circuit; and (ii) the distance of the vertices distance along the sequence (which is the number of chords between them) is odd and it is not 1. If each PV-pair is a non-chord (that is ∈ F), then this circular swap is called a F-compatible swap or F-swap for short.
The F-swap is one of the central notions of this paper. When i = 2 then the circular C 4 -swap coincides with the classical Havel-Hakimi swap. When i = 3 then we get back the notion of the triangular C 6 -swap, which was introduced in paper [2] in connection with directed degree sequences.
We define the weight of the F-compatible circular C 2i -swap as w(C 2i ) = i − 1. This definition is in accordance with the definitions of the weight of the classical HH-swaps, and the weight of a triangular C 6 -swap in paper [2]. Furthermore it is well known (see for example again [2]) that in case of simple graphs (i − 1) Havel-Hakimi swaps are needed to alternate the edges along C 2i . As we will see next the same applies for any (elementary) circular C 2i -swap: Lemma 2.2. Let G be a realization of d F and let the elementary chordcircuit C of length 2i alternate. Then the circular C-swap operation can be carried out by a sequence of F-swaps of total weight i − 1.
More precisely there exists a sequence G = G 0 , G 1 , . . . , G ℓ of realizations such that for each j = 0, . . . , ℓ − 1 there exists an F-swap operation from G j to G j+1 , the difference between G and G ℓ is exactly the alternating circuit C, finally the total weights of those F-swap operations is i − 1.
Proof. We apply mathematical induction for the length of the chord-circuit: assume this is true for all circuits of length at most 2i − 2. (For i = 3 we already saw this.) Then take an alternating elementary chord-circuit C of length 2i in a realization of d F .
If each PV-pair in C is a non-chord, then the circular C 2i -swap itself is a F-swap of weight i − 1. So we may assume that there is a PV-pair uv in C which is a chord. This chord together with the two "half-circuits" of C form chord-circuits C 1 and C 2 using the chords of the original circuit C and twice the chord uv. One of them, say C 1 , is alternating. The length of C 1 is 2j < 2i therefore there exists a F-compatible swap sequence of total weight j − 1 to process it. After the procedure the status of uv will alter into the other status. With this new status of the chord the circuit C 2 becomes an alternating one, so it can be processed with 2i+2−2j 2 − 1 total weight -and after this procedure the chord uv is switched back to its original status. So we found a swap sequence of total weight i − 1 which finishes the proof.
The space of all realizations: Consider now the set of all possible realizations of a restricted graphical degree sequence d F . Let G and H be two different realizations. The natural question is whether G can be transformed into H using F-swaps?
For classical degree sequences this problem was solved affirmatively already in 1891 by Petersen in the (by now almost completely forgotten) paper [16]. Havel's paper [7] also provided (an implicit) solution. For bipartite graphs (with possible multiple edges but no loops) this was done by Ryser ([17]). For simple bipartite graphs it was folklore. Finally for directed graphs it was done in [11] (and later rediscovered in [3]). Let's recall that a circuit in a simple graph G is a sequence of vertices v 0 , . . . , v 2t , where v 0 = v 2t s.t. the consecutive vertices are adjacent and each edge can be used at most once. Note that there can also be other indices i < j such that v i = v j . A circuit is called a cycle, if it is simple, i.e., for any i < j, v i = v j only if i = 0 and j = 2t. A circuit is alternating for G and H if the edges come in turns from E(G) and E(H). When this is the case then the corresponding chord-circuit in realization G (as well in H) is alternating.
We can find a decomposition that no circuit contains a vertex v twice s.t. their distance δ (the number of edges between them) is even. Indeed, since G is simple, therefore δ is at least four, consequently the vertex v can split the original circuit into two smaller, but still alternating circuits. So for example an alternating circuit decomposition with maximal number of circuits has this property. It also implies that no circuit may contain a vertex three times otherwise at least two copies of the vertex would be of even distance from each other.
The application of Lemma 2.2 proves that each circuit C can be processed with |C|/2 − 1 total weight. This finishes the proof.
In the paper [2] the following formula was developed for the required minimum weight of transforming one realization in to an other one of the same (unrestricted) degree sequence: Consider the realizations G and H and the symmetric difference E∆F of their edge sets. Denote maxC(G, H) the maximum possible number of the alternating circuits in such a decomposition.
Theorem 2.4 (Erdős -Király -Miklós [2]). For any pair G, H of realizations of a given degree sequence the weight of the shortest swap sequence between the two realizations is Now we show that the analogous result holds for any restricted degree sequence.
Proof. We show that The first equality is just Theorem 2.4. To check the second inequality assume that the F-swap sequence S C 1 , . . . , S Cn (2.5) transfers G to H. If the length of C j is 2k j , then, applying Theorem 2.4 again, we have that S C j can be obtained as a composition of k j − 1 many "standard" swaps S C j,1 , . . . , S C j,k j −1 . So the "standard" swap sequence transfers G to H, and the weight of this sequence is n j=1 (k j − 1), which is exactly the weight of the F-swap sequence S C 1 , . . . , S Cn from (2.5). Thus the second inequality holds.
Finally the third inequality holds by Lemma 2.2.

The star+factor restricted degree sequences
We turn our interest now for the following specialized restricted degree sequence problem: d F is called a Star -1-Factor Restricted Degree Sequence problem (or star+factor problem for short), if (Ψ) the set F of the forbidden edges is a bipartite graph where the edges are the union of an 1-factor and a star with center s.
Similarly, if bd is a bipartite degree sequence, and (Ψ) holds for F, then bd F is called a Bipartite Star -1-Factor Restricted Degree Sequence problem (or bipartite star+factor problem for short).
As we already mentioned, Tutte's f -factor theorem can always be utilized to find actual graphical realizations of our star+factor restricted degree sequence. However in this special case we can apply a Havel-Hakimi type greedy algorithm to construct such realizations. Consider our star+factor degree sequence problem d F . For a given vertex x ∈ V denote C(x) the set of those vertices in V which form chords together with x. If (i) for each y ∈ C(x) has at most one vertex, denoted by y F , such that chord yy F belongs to F, furthermore if y F = z F then y = z, then we say that C(x) is normal, and we fix an order ≺ x on C(x) such that Lemma 3.1. Let G be a graphical realization of our star+factor RDS d F , let x ∈ V and assume that C(x) is normal with the order ≺ x . Assume furthermore that y ≺ x z and xz ∈ E while xy ∈ E. Then there exists an alternating chord-circuit C = (y, x, z, v 1 , . . . , v i ) with i = 1 or 3. In the acquired new realization we have This statement is actually almost the same as Lemma 4 in [3] and its proof could be recalled. However here we give a complete proof. On one hand this keeps this paper self-contained, on the other hand paper [3] uses a different language.
Proof. We have xz ∈ E but xy ∈ E. At first assume that there exists a vertex u = x, y, z, such that uy ∈ E, and uz ∈ E but u = z F . When such vertex exists then C = (x, z, u, y) is a suitable alternating chord-circuit as the xz, uy ⇒ xy, uz swap shows.
Definition 3.2. From now on the notation xz, uy ⇒ xy, uz always means that xz, uy are edges, xy, uz are non-edge chords, and we create a new realization with the indicated swap.
We continue the proof: when d(y) > d(z) then this happens automatically since y belongs to at most one forbidden pair. However, if d(y) = d(z) then it can happen that z F y ∈ E and It is important to observe that in this case y F z ∈ E, otherwise some u would not satisfy (3.1) (in order to keep d(y) = d(z)).
So the only case when we do not find automatically an appropriate swap with x, y and z is when d(y) = d(z), yz F is an edge and zy F is a chord but not an edge. In this case, we can find a u = y, z such that is the required alternating chord circle. See the figure below. The three line types denotes the edges, the chords which are non-edges, finally the forbidden non-chords.
Lemma 3.1 provides the following easy Havel-Hakimi type greedy algorithm to construct at least one graphical realization of our restricted degree sequence problem.
HH-algorithm for a star+factor degree sequence problem. (Recall that the center of the forbidden star is denoted with s. Since a star can be empty, we can assume that s is always defined): (H 1 ) take an ordering ≺ s on C(s) (which is normal) and connect s to the first d(s) vertices (with respect to ≺ s ) of C(s). Delete s and update the degrees of the used vertices accordingly.
(H 2 ) take the remaining vertices one by one and repeat the process. The d F star+factor restricted degree sequence is graphical if and only if the previous greedy HH-algorithm provides a realization.
Proof. As usual, it is enough recursively apply Lemma 3.1. Since we start our algorithm with the vertex s therefore C(s) satisfies condition (i) therefore it is normal. When s is deleted, then the remaining forbidden edges are from the original 1-factor, so the lemma applies automatically at each subsequent step.
In the case of a bipartite degree sequence d the situation is very similar: we define the normality of any C(x) formally the same way as before. The definition of the order ≺ x is also the same as before.
Here it is interesting to observe, that if x is in class U than C(x) is subset of class W, the vertices y F and z F belong again to class U , so u ∈ W (and, finally, the forbidden edges define those vertices are elements of F). Furthermore, as in any bipartite degree sequence problem, if we determine all edges adjacent to the vertices in U then we automatically placed all edges adjacent to the vertices in W as well. Proof. Indeed, by definition, the center s of the forbidden star belongs to U . Also by definition, in the proof of the Lemma the vertices we are considering are from the vertex class U only: x ∈ U. Therefore C(x) is a subset of W. Consequently the vertices y F and z F belong to vertex class U therefore the vertex u must belong to vertex class W again. Therefore all forbidden edges considered in the proofs belong to F. So the proofs apply without changes for the bipartite case as well.

Sampling uniformly half-regular, bipartite star+factor restricted degree sequences
In this section we describe the main result of our paper which is an MCMC algorithm for (almost) uniform sampling of the space of all realizations of the Half-regular, Bipartite Star -1-Factor Restricted Degree Sequence problem: Let d be a bipartite degree sequence with a star+ 1-factor type forbidden edge set F, where the center of the star is denoted by s (and belong to U ). Furthermore let (Φ) the degree sequence d is a half-regular bipartite one: it requires a B = (U, W ; E) bipartite graph where each vertex u ∈ U -with the one possible exception s -has the same degree. We will write V for U ∪ W.
We will apply a Markov Chain Monte Carlo method for almost uniform sampling of all possible realizations of our d F . Originally the MCMC method for realizations' sampling was proposed by Kannan, Tetali and Vempala (1999, [9]). They conjectured that their process is rapidly mixing on the realizations of any (unrestricted) degree sequences, i.e., starting from an arbitrary realization of the degree sequence, the process reaches a completely random realization in reasonable (i.e., polynomial) time. They managed to prove the result for bipartite regular graphs. Their conjecture was proved for arbitrary regular graphs by Cooper, Dyer and Greenhill (2007, [1]). (Their result does not automatically generalize the previous result, since their version does not allow forbidden edges.) An analogous theorem was proved by Greenhill on regular directed graphs ( [5]). Miklós, Erdős and Soukup proved in [13] that this Markov process is also rapidly mixing on each bipartite half-regular degree sequence (here there is no exceptional vertex s). In this paper we will prove that the analogous Markov process is rapidly mixing on the half-regular star+factor bipartite restricted degree sequence problem.
The state space of our Markov chain is the graph G = (V (G), E(G)) where V (G) consists of all possible realizations of our problem, while the edges represent the possible swap operations: two realizations (which will be indicated by upper case letters like X or Y ) are connected if there is a valid F-swap operation which transforms one realization into the other one (and the inverse swap transforms the second one into the first one as well).
The transition (probability) matrix P is defined as follows: let the current realization be G. Then (a) with probability 1/2 we stay in the current state (that is, our Markov chain is lazy); (b) with probability 1/4 we choose uniformly two-two vertices u 1 , u 2 ; v 1 , v 2 from classes U and W respectively and perform the swap if it is possible; (c) finally with probability 1/4 choose three -three vertices from U and W and check wether they form three pairs of forbidden chords. If this is the case then we perform a circular C 6 -swap if it is possible.
Here cases (b) and (c) correspond to Lemma 3.1. The swaps moving from G to its image G ′ is unique, therefore the probability of this transformation (the jumping probability from G to G ′ = G) is: and (More precisely these are the probabilities that these vertex sets will be checked against making a swap.) The probability of transforming G to G ′ (or vice versa) is time-independent and symmetric. Therefore P is a symmetric matrix, where the entries in the main diagonal are non-zero, but (probably) different values. As we discussed it earlier (see Theorem 2.3), our Markov chain is irreducible (the state space is connected), and it is clearly aperiodic, since it is lazy. Therefore, as it well known, our Markov process is reversible with the uniform distribution as the globally stable stationary distribution.
Our main result is the following: Theorem 4.1. The Markov process defined above is rapidly mixing on each bipartite half-regular degree sequence with a forbidden star and a forbidden 1 -factor.
This result supersedes the previously mentioned three results ( [9,5,13]) (except that this does not care of the actual mixing time, it just proves that it is polynomial). We want to add, however, that the friendly path method, described in [13], was not intended to handle half-regular bipartite degree sequences only but all bipartite degree sequences. Therefore we think that that method should not be completely neglected.
There Sinclair's multicommodity flow method defines a bunch of paths (consecutive sequences of swaps) for each realizations pair X and Y which transform realization X into Y.
So consider two realizations X ∈ G and Y ∈ G, and consider the symmetric difference ∆ = E(X)∆E(Y ). In the bipartite graph Θ = (U, W ; ∆) for each vertex v the number of adjacent X-edges (= E(X) \ E(Y )) and the number of the adjacent Y -edges are the same. Therefore, due to Euler classical reasoning, it can be decomposed into alternating circuits.
The simplified Sinclair's multicommodity path method consists of two phases: In Phase 1 we decompose the symmetric difference ∆ into alternating circuits on all possible ways. In each cases we get an ordered sequence W 1 , W 2 , . . . , W κ of circuits. (Usually there are a huge number of different decompositions.) Each circuit is endorsed with a fixed cyclic order.
In Phase 2 each circuit W i from the (ordered) decomposition derives one unique alternating cycles decomposition: This decomposition is fully determined by the circuit and its well defined edge order. (Both construction algorithms are fully described in Section 5 of the paper [13], we do not discuss them here.) The ordered circuit decomposition together with the ordered cycle decompositions of all circuits altogether provide a well defined ordered cycle decomposition C 1 , . . . C ℓ of ∆.
This ordered cycle decomposition determines ℓ−1 realizations H 1 , . . . H ℓ−1 with the following property: if we use the notations H 0 = X and H ℓ = Y then for each j = 0, . . . , ℓ − 1 we have E(H j )∆E(H j+1 ) = C j+1 . (It is important to recognize that till this point we did not process even one swap operation! We just identified ℓ − 1 realizations which will be along our canonical path.) We will define a unique canonical path from X to Y determined by this circuit decomposition which uses these realizations H j as milestones along the path. The canonical path will be X = G 0 , . . . , G i , . . . , G m = Y where each G i can be derived from G i−1 with one valid swap operation, where we must have the following property: there are some increasing indices 0 < n 1 < n 2 < · · · < n ℓ such that we have G n i = H i . This, together the definitions of H i means that The canonical path we are looking for has two important further properties: for each i < ℓ the constructed path H i = G ′ 0 , G ′ 1 , . . . , G ′ m ′ = H i+1 between G n i and G n i+1 must satisfies that (Θ) m ′ ≤ c · |C i | for a suitable constant c; (Ω) for each j there is where the notations M G stands for the usual bipartite adjacency matrix of G (this will be defined in details at the beginning of the next section), and d stands for the Hamming distance of two matrices of the same form, finally Ω 2 is a small constant.
The current value of the auxiliary matrix M X + M Y − M G ′ j together with the symmetric difference ∆, furthermore a small (polynomial) size parameter set, finally the vertices in G on which the canonical path under investigation goes through uniquely determine the vertices X, Y and the path itself. Therefore it can be used to control certain features of the canonical path system. If the overall number of these auxiliary matrices are small (their number is smaller than a small polynomial of n multiplied with the number of possible realizations -as it is ensured by (Ω)), then -as it was proved in [13] -our Markov chain is rapidly mixing.
So in Phase 2 we have to build up our swap sequence between H i and H i+1 for all values i taking care for conditions (Θ) and (Ω). This will happen in the next Section.

The construction of swap sequences between consecutive "milestones"
Now we are going to implement our plan described above. At first we introduce some shorthand. Instead of H i and H i+1 we will use the names G and G ′ . These two graphs have almost the same edge set. More precisely Of course E(G)∆E(G ′ ) = C i also holds. We refer for the elements of C i ∩ E(X) as X-edges while the others are the Y -edges. We denote the cycle itself as C, it has 2ℓ edges and its vertices are u 1 , w 1 , u 2 , w 2 , . . . , u ℓ , w ℓ . Since C has at least four vertices, therefore we may assume that u 1 = s (so u 1 is not the center of the forbidden star). Finally w.l.o.g. we may assume that the chord u 1 w 1 is an Y -edge (and, of course, w ℓ u 1 is an X-edge).
We are going to construct one by one the realizations G ′ j . We build our canonical path from G toward G ′ and at any particular point the last constructed realization is denoted by Z. (At the beginning of the process we have Z = G.) We are looking for the next realization, denoted by Z ′ .
Before we continue the discussion of the canonical path system, we have to introduce our control mechanism, mentioned in condition (Ω). This auxiliary structure originally was introduced by Kannan, Tetali and Vempala in [9]: For any particular realization G from G the matrix M G denotes the adjacency matrix of the bipartite realization G where the columns and rows are indexed by the vertices of U and W resp. (Therefore the column sums are the same in each realization, except perhaps at column s.) Our indexing method is a bit unusual: the columns are numbered from left to right while the rows are numbered from bottom to the top. (Like in the Cartesian coordinate system.) This matrix is not necessarily symmetric, and elements M i,i can be different from 0.
For example if we consider the submatrix in M G spanned by u 1 , . . . , u ℓ and w 1 , . . . , w ℓ then we have M G (i, i) = 0 for i = 1, . . . , ℓ, while M G (i, i − 1) = 1 (for i = 2, . . . , ℓ) and M G (1, ℓ) = 1. (So the first value gives the column, the second one gives the row.) The non-chords between vertices in the same vertex class are not considered at all, while non-chords which are forbidden are denoted by . As it is clear from the previous sentence, we will identify each chord or non-chord with the corresponding position in the matrix.
Our auxiliary structure is the matrix By definition, each entry of a bipartite adjacency matrix is 0 or 1 (or ). Therefore only −1, 0, 1, 2 can be the "meaningful" entries of M . An entry is −1 if the edge is missing from both X and Y but it exists in Z. This is 2 if the edge is missing from Z but exists in both X and Y. It is 1 if the edge exists in all three graphs (X, Y, Z) or it is there only in one of X and Y but not in Z. Finally it is 0 if the edge is missing from all three graphs, or the edge exists in exactly one of X and Y and in Z. (Therefore if an edge exists in exactly one of X and Y then the corresponding chord in M is always 0 or 1.) One more important, but easy fact is the following: Next we will determine the swap sequence between G and G ′ through an iterative algorithm. At the first iteration we check, step by step, the positions (u 1 , w 2 ), (u 1 , w 3 ), . . . , (u 1 , w ℓ ) and take the smallest j for which (u 1 , w i ) is an actual edge in G. Since (u 1 , w ℓ ) is an edge, therefore such i always exists. So we may face to the following configuration: chord, non-edge non-chord unknown We will call this (u 1 , w i ) chord as start-chord of our current sub-process and u 1 w 1 is the end-chord. We will sweep the alternating chords along the cycle from the start-edge w i u i (non-edge), u i w i−1 (an edge) toward the end-edge w 1 u 1 (non-edge) -switching their status in twos and fours. We check positions u 1 w i−1 , u 1 w i−2 (all are non-edges) and choose the first chord among them, we will call it the current-chord. (Since u 1 = s therefore we never have to check more than two edges to find the first chord, and we need only one times to check two, since there is at most one non-chord adjacent to u 1 .) Case 1: As we just explained the typical situation is that the currentchord is the "next" one, so when we start this is typically u 1 w i−1 . Assume that this is a chord. Then we can proceed with the swap operation We just produced the first "new" realization in our sequence, this is G ′ 1 . For the next swap operation this will be our new current realization. This operation will be called a single-step.
In a realization Z we will call a chord bad, if its current status (being edge or non-edge) is different from its status in G (or, what is the same, in G ′ , since they differ only on the chords along the cycle C). After the previous swap, we have two bad chords in G ′ 1 , namely u 1 w i−1 and w i u 1 . Consider now the auxiliary matrix M (X + Y − Z) (here Z = G ′ 1 ). As we saw earlier, for each position outside the chords in C the status of that particular position in Z is the same as in X or Y or in both. Accordingly, the corresponding matrix value is 0 or 1. Case 2: If the position "below" the start-chord is a non-chord, then we cannot produce the previous swap. Then, however, the non-edge u 1 w i−2 is the current-chord. For sake of simplicity we assume that i − 2 = 2 so we are in Figure 1. Consider now the alternating C 6 cycle: u 1 , w 2 , u 3 , w 3 , u 4 , w 4 . It has altogether three vertex pairs which may be chords. We know already that u 1 w 3 is a non-chord. If none of the three is chord, then this is an F-compatible circular C 6 -swap -and accordingly to the definitions we can swap it in one step. Again, we found the valid swap w 2 u 3 , w 3 u 4 , w 4 u 1 ⇒ u 1 w 2 , u 3 w 3 , u 4 w 4 . After that we again have 2 bad chords, namely u 1 w 2 and w 4 u 1 , and together we have at most two bad positions in the new M (X + Y − Z) with at most one 2-value and at most one −1-value.
Finally if one position, say w 2 u 4 , is a chord then we can process this C 6 with two swap operations. If this chord is, say, an actual edge, then we swap w 2 u 4 , w 4 u 1 ⇒ u 1 w 2 , u 4 w 4 . After this we can take care for the w 2 , u 3 , w 3 , u 4 cycle. Along this sequence we never create more, than 3 bad chords: the first swap makes chords w 2 u 4 , w 4 u 1 and u 1 w 2 bad ones, and the second "cures" w 2 u 4 but does not touch u 1 w 2 and w 4 u 1 . So along this swap sequence we have 3 bad chords, at the end we have only 2. On the other hand, if the chord w 2 u 4 is not an edge, then we can swap w 2 u 3 , w 3 u 4 ⇒ u 3 w 3 , u 4 w 2 , creating one bad edge, then taking care the four cycle u 1 , w 2 , u 4 , w 4 we "cure" w 2 u 4 but we switch u 1 w 2 and w 4 u 1 into bad chords. We finished our double-step along the cycle.
In a double-step we make at most three bad chords. When the first swap uses three chords along the cycle then we may have at most one bad chord (with M -value 0 or −1) and then the next swap switches back the chord into its original status, and makes two new bad chords (with at most one 2-value and one −1-value). When the first swap uses only one chord from the cycle, then it makes three bad chords (changing two chords into non-edge and one into edge), therefore it may make at most two 2-values and one −1-value. After the second swap there will be only two bad chords, with at most one 2-value, and at most one −1-value.
When only the third position corresponds to a chord in our C 6 then after the first swap we may have two −1-values and one 2-value. However, again after the next swap we will have at most one of both types.
Remark 5.2. When two realizations are one swap apart (so they are adjacent in G) then we say that their auxiliary matrices are at swap-distance one. Since one swap changes four positions of the matrix, therefore the Hamming distance of these matrices is 4. Finishing our single-or double-step the previous current-chord becomes the new start-chord and we look for the new current-chord. Then we repeat our procedure. There is only one important point to be mentioned: along the step, the start-chord switches back into its original status, so it will not be a bad chord anymore. So even if we face a double-step the number of bad chords never will be bigger than three (together with the chord w i u 1 which is still in the wrong status, so it is bad), and we have always at most two 2-values and at most one −1-value in M (X + Y − Z).
When our current-chord becomes to w 1 u 2 then the last step will switch back the last start-chord into its correct status, and the last current-chord cannot be in bad status. So, when we finish our sweep from u 1 w i to w 1 u 1 at the end we will have only one bad chord (with a possible 2-value in M ). This concludes the first iteration of our algorithm.
For the next iteration we seeks a new start-chord between w i u 1 and w ℓ u 1 and chord w i u 1 becomes the new end-chord. We will repeat our sweeping process for this setup, and we will repeat it as long as all chords will be processed, so we fond the entire realization sequence from G to G ′ . If in the first sweep we had a double-step, then it will never occur later, so altogether with the bad (new) end-chord we never have more than three bad chords, with at most two 2-values and at most one −1-value.
However, if the double-step occurs sometimes later, for example in the second sweep, then we face to the following situation: if we perform a circular C 6 -swap, then there cannot be any problem. So we may assume that there is a chord in our C 6 , suitable for a swap. If this chord is a non-edge, then the swap around it produces one bad chord, and at most one bad position in M . The only remaining case when that chord is an edge. After the first swap there will be four bad chords, and there may be at most three 2-values and at most one −1 value. However after the next swap (finishing the double step) we annihilate one of the 2-values, and after that swap there are at most two 2-values and at most −1-value along the entire swap sequence. When we finish our second sweep, then chord w i u 1 will be switched back into its original status, it will not be bad anymore.
We apply iteratively the same algorithm, and after at most ℓ sweep sequence, we will process the entire cycle C. This finishes the construction of the required swap sequence (and the required realization sequence).
Meanwhile we also proved the following important observation: Lemma 5.3. Along our procedure each occurring auxiliary matrix M (X + Y − Z) is at most swap-distance one from a matrix with at most three bad positions: with at most two 2-values and with at most one −1-value in the same column, which does not coincide with the center of the forbidden star.

The analysis of the swap sequences between "milestones"
What remains is to show that the defined swap sequences between H i and H i+1 satisfy the properties (Θ) and (Ω) of the simplified Sinclair's method. The first one is easier to see, since we can process a cycle of length 2ℓ in ℓ − 1 swaps. Therefore the derived constant c in (Θ) is actually 1.
We introduce the switch operation on 0/1 matrices with forbidden positions: we fix the four corners of a submatrix (none of them is forbidden), and we add 1 to two corners in a diagonal, and add −1 to the corners on the other diagonal. This operation clearly does not change the column and row sums of the matrix. For example if we consider the matrix M G of a realization of our d F and make a valid swap operation, than it looks like as a switch in this matrix. The next statement is trivial but very useful: Lemma 6.1. If two matrices have switch-distance 1, then their Hamming distance is 4. Consequently if the switch-distance is c then the Hamming distance is bounded by 4c.
We will prove now that property (Ω) holds for our auxiliary matrices: Proof. Consider now a certain M which is not necessarily a valid adjacency matrix of a realization. We will show pictures about the submatrix in this matrix which describes the current alternating cycle C. We choose a submatrix, where the center s of the forbidden star is in the first column. (We choose this submatrix as an illustration tool, but we still consider the entire matrix to work with.) We know that this matrix contains at most two 2-values and at most one 1-value. All these positions are adjacent to the center u 1 of our sweeping sequence (see Figure 1), so they are in the same column.
For simplicity from now on we will denote the center of the sweep as well the column u. The forbidden positions are denoted with . Any column (except column 1) may contain at most one of them, and any row may contain at most two of them. Finally in our pictures the character ⋄ stands for a character which we are not interested in. That is, it can be 0 or 1 or .
We will distinguish cases, depending on the occurring of values 2 and −1.
Case 1. Column u has one bad position, which can be −1 or 2, or it has two 2-values. Consider at first the subcase when M [uw] = −1. By definition that means that chord uw is an edge in Z but non-edge in both X and Y. So vertex w ∈ W has at least one adjacent edge, therefore the row-sum in its row is at least 1. Therefore there are at least two positions in row w with entries 1. They are in column u 1 and u 2 . At least one of them, say u 1 , differs from s. Since the column sums are constant, therefore there exists at least two rows w 1 such that M     If this is not the case then we consider the following: assume that M [u 1 w 1 ] = 0. The column sums are the same, and we assumed that M [u 1 w 3 ] = 0 or . Therefore the difference between column sums in u and u 1 is 1 due to rows w 1 and w 3 , and the difference increase at least 1 for row w 2 , where against a 2-value in column u there is either 1 or 0 in column u 1 . Therefore there exists at least two further rows, where there is a 1 in column u 1 against a 0 or in column u. Since column u can contain at most one , one of the rows must contain a 0. Let it be denoted by w 4 . Hence M [u 1 w 4 ] = 1 while M [uw 4 ] = 0. Then we can switch off this 2-value without making a new bad position. After that we are back to Case 2. Altogether this requires at most three switches. We finished the proof of Lemma 6.3.
In turn this proves Theorem 6.2, so our Markov chain is rapidly mixing as Theorem 4.1 stated.

Self-reduced counting problem
A decision problem is in NP if a non-deterministic Turing Machine can solve it in polynomial time. An equivalent definition is that there exists a witness proving the yes answer to the question which witness can be verified in polynomial time. A counting problem is in #P if it asks for the number of those witnesses of a problem from NP that can be verified in polynomial time (it might happen that not all witnesses are verifiable in polynomial time).
Two complexity classes, FPRAS and FPAUS, concern the approximability of counting problems. Here we give only narrative descriptions of these complexity classes, the detailed definitions can be found, for example, in [8].
A counting problem from #P is in FPRAS (Fully Polynomial Randomized Approximation Scheme) if the number of solutions can be quickly estimated with a randomized algorithm such that the estimation has a small relative error with very high probability.
A counting problem from #P is in FPAUS (Fully Polynomial Almost Uniform Sampler) if the solutions can be sampled quickly with a randomized algorithm that generates samples following a distribution being very close to the uniform one.
It is easy to see that a counting problem is in FPAUS if there is a rapidly mixing Markov chain for which • a starting state can be generated in polynomial running time; • one step in the Markov chain can be conducted in polynomial running time; and • the relaxation time of the Markov chain grows only polynomially with the size of the problem.
The Markov chain we gave the star+factor problem satisfies all these requirements. Jerrum, Valiant and Vazirani proved that any self-reducible counting problem is in FPRAS iff it is in FPAUS [8]. A counting problem is selfreducible if the solutions for any problem instance can be generated recursively such that after each step in the recursion, the remaining task is another problem instance from the same problem, and the number of possible branches at each recursion step is polynomially bounded by the size of the problem instance.
Clearly, a graph with prescribed degree sequence can be built recursively by telling the neighbors of a node at each step, then removing the node in question and reducing the degrees of the selected neighbors. However, this type of recursion does not satisfy all the requirement for being selfreducible since there might be exponentially many possibilities how to select the neighbors of a given vertex.
On the other hand, the degree sequence problem with a forbidden one factor and star tree is a self-reducible counting problem. Indeed, consider the center of the (possibly empty) star, s ∈ U , and the vertex v ∈ V with the smallest index for which (s, v) is a chord. Any solution for the current problem instance belongs to one of the following two cases: • The chord (s, v) is not present in the solution. In that case, extend the size of the star by adding chord (s, v) to the forbidden set, and do not change the degrees. This is another problem instance from the star+factor problem, whose solutions are the continuations of the original problem belonging to this case.
• The chord (s, v) is present in the solution. In that case, extend the size of the star by adding chord (s, v) to the forbidden set, and decrease both d s and d v by one. The new degree sequence is still a half-regular, bipartite star+factor restricted degree sequence, and the solutions of this new problem extended with the previously decided step provide solutions of the original problem.
Since the star+factor counting problem is a self reducible counting problem, it is in FPRAS as it is in FPAUS. [13]. To do so properly it requires a slight generalization of the original method.
The method takes two realizations X and Y of the same degree sequence. It considers all possible ordered circuit decompositions of the symmetric difference of the edge sets, then it uniquely decomposes each such decomposition into an ordered sequence C = C 1 , . . . , C m of oriented cycles. Based on this latter decomposition the method determines a well defined unique path between X and Y in the Markov chain G.
For that end the method defines first a sequence of "milestones". These are different realizations X = H 0 , H 1 , . . . , H m−1 , H m = Y of the degree sequence where the edge set of any two consecutive realizations H i−1 , H i differ exactly in the edges along the cycle C i . (Until this point no swap operation happened.) In the next phase for any particular i = 0, . . . , m − 1 the method determines a sequence of valid swap operations transforming H i−1 into H idescribing a unique path Z 0 , Z 1 , . . . Z ℓ between H i−1 and H i in the Markov chain G. This sequence of course heavily depends on the available swap operations. In paper [13] these are the usual (bipartite) swap operations. In the current paper these are the restricted swap operations. These operations, while exchanging chords in the realizations along the alternating cycle C i , also use some further chords. Therefore the edge set of any Z i is not completely contained by E(X) ∪ E(Y ), there exist a small number of edges in Z i which are non-edges in X and in Y , or non-edges in Z i but edges in X and Y. If Z i is between the milestones H m and H m+1 , then C j for j = m alternates in Z i , and C i alternates with a "small error": there is a very small number of vertices where the alternation does not hold.
Along the process the simplified Sinclair's method requires (see the paper [13], Section 5, (F)(c) ) that this number must be small. In the original application this number is actually one. Here, as we saw in Section 5, this number is three: that many bad chords may occur after any particular RSO. As we saw all these chords are adjacent to the same vertex u 1 .
These numbers are used by the method to determine the size of a parameter set B. This parameters set must have a polynomial size. When we have one bad chord, then it is determined by its end points -there are at most n 2 possibilities for them. This provides an n 2 multiplicative factor to the size of B. When we have at most three bad chords, then they can be chosen at most n 4 ways: point u 1 is fixed (n different choices), while the other three end points can be chosen at most n 3 independent ways. Altogether it provides an at most n 4 multiplicative factor to the size of B. This remark finishes the proof of the simplified Sinclair's method for the case of these restricted swap operations.