# An Iterative Scheme to Compute Size Probabilities in Random Graphs and Branching Processes

Authors: Paolo Serafini
Journal: Scientific Programming
 (2018)

---

## Abstract

We deal with a functional equation that plays an important role in random graphs and in branching processes. In branching processes, the functional equation relates offspring probabilities to population size probabilities, while in random graph it relates degree probabilities to small component size probabilities. We present an iterative scheme that allows computing the size probabilities numerically. It is also theoretically possible to invert the iteration, although this inverse iteration is numerically unstable.


---

## Body

## 1. Introduction

LetG(x) and H(x) be two probability generating functions that are linked through the functional equation(1)Hx=xGHx.Functions of this type occur in branching processes and in random graphs [1–6]. In branching processes, G(x) represents the probabilities of new offspring from a member of the population and H(x) represents the population size probabilities. In the configuration model [4] of random graphs, G(x) represents the excess degree probabilities of a vertex in small components and H(x) represents the small component size probabilities. Note that, in both cases, H(x) can be a defective generating function; that is, H(1)<1.UsuallyG(x) is given and H(x) has to be computed. Only in rare cases is it possible to find an explicit analytic expression of H(x). However, a numerical iteration to compute the coefficients of H(x) is possible. To the best of our knowledge, such a question has not been investigated and it seems that the iteration we propose in this paper is new.Interestingly enough, this iteration can be inverted; that is, from the size distribution, we can infer the degree probabilities. We present also this inverse iteration, although it has to be remarked that the inverse iteration is numerically unstable.The paper is organized as follows. In Section2, we provide the mathematical background by referring to the case of random graphs. Then, in Section 3, we present the main result, that is, the iteration to compute the size probabilities of the small components of the graph. The possibility of inverting this computation is presented in Section 4. Then, in Section 5, we point out how the same iteration can be used for a branching process. Some conclusions are presented in Section 6.


## 2. Mathematical Background

We first present our result by explicitly referring to random graphs in the configuration model for which the picture is more complex. In a later section, we show how to relate the iteration to branching processes. Hence, all definitions in this section and in Sections3 and 4 are related to random graphs.A random graph has assigned degree probabilitiesph, h=0,1,…, and ph is the probability that a randomly selected vertex has degree h. We recall that the degree of a vertex is the number of vertices adjacent to it. The study of random graphs through generating functions is asymptotic; that is, it assumes an infinite number of vertices. Let G0(x) be the probability generating function of the degree distribution; that is, (2)G0x=∑h≥0phxh.Let d=G0′(1) be the average degree and(3)G1x=G0′xd=∑h≥0qhxh,where clearly (4)qh=h+1ph+1d.The qh values are known asexcess degree probabilities. Let H0(x) and H1(x) be two generating functions that can be expressed as power series as (5)H0x=∑k≥1skxk,H1x=∑k≥1rkxk,and they are defined by the equations(6)H0x=xG0H1x,H1x=xG1H1x.Our aim is to compute the coefficients sk and rk.The motivation for the generating functionsH0(x) and H1(x) derives from the analysis of the asymptotic properties of the random graph in the configuration model. If the graph is sufficiently dense, it exhibits the so-calledgiant component, that is, a connected component whose size asymptotically goes to infinity. The giant component, if present, is unique. The rest of the graph consists of an infinite number of finite trees, the so-calledsmall components (see [3], among many possible references).It can be shown that ifph is the probability that a randomly chosen vertex (in the whole graph) has degree h, then sk is the probability that a randomly chosen vertex belongs to a small component of size k, and rk is the probability that, after choosing a random vertex i of degree at least one and then a random vertex j adjacent to i, the vertex j belongs to a small component of size k after removing the edge {i,j}.If the giant component is present, the conditional probabilityp^h of choosing in the small components a vertex of degree h is different from ph and similarly for the excess degree probability q^h. It can be shown that(7)p^h=uhvph,q^h=uh-1qh,where u is the solution of u=G1(u) and v=G0(u) is the fraction of vertices in the small components. We can briefly justify (7) by using Bayes’ formula:(8)p^h=Pr⁡S∣DhphPrS,with S being the random event of choosing a vertex in a small component and Dh being the random event of choosing a vertex of degree h. Clearly Pr{S∣D0}=1 and consequently p^0=p0/v. If h>0, Pr{S∣Dh} is the probability thatall adjacent h vertices belong to a small component once we have removed the corresponding edges, and so its value is uh. This explains the left expression in (7). To justify the right expression, we need to compute the average degree in the small components by taking the derivative of G0(ux)/v (by using (7)) and computing it for x=1; that is, uG0′(u)/v=udG1(u)/v=u2d/v. From this, we immediately get the expression at the right.It turns out that usingp^h instead of ph in the definition of G0 and G1 has the only effect of scaling the sk values by the constant factor v and the rk values by the constant factor u, which correspond to the conditional probability of choosing within the small components. In particular, we have H0(1)=v and H1(1)=u if we use ph and qh in the definition of G0 and G1, respectively, whereas we have H0(1)=1 and H1(1)=1 if we use p^h and q^h.We also define the probabilitytk that a randomly selected small component has size k. Of course(9)skv=ktkt¯,with t¯ being the average size of a small component. Here we have to discount sk because the choice of small components necessarily conditions the choice within the small components.


## 3. Computing the Coefficients ofH0(x) and H1(x)

From the recursive equation,(10)H1x=xG1H1x=x∑h≥0qhH1xh,and from (11)H1x=∑k≥1rkxk(necessarily r0=0), we derive(12)∑k≥1rkxk=xq0+x∑h≥1qh∑j≥1rjxjh=xq0+∑h≥1qhxh+1∑j≥0rj+1xjh,so that(13)∑k≥1rkxk=xq0+∑h≥2qh-1xh∑j≥0rj+1xjh-1.Let akh be the coefficient of xk in (∑j≥0‍rj+1xj)h. Note that ak1=rk+1 and in particular a01=r1. From (13), (14)r1=q0,rk=∑h=2kqh-1ak-hh-1,k≥2.Hence, the computation of rk requires the coefficients ak-21,ak-32,…,a0k-1. In turn, the computation of akh requires the terms r1,…,rk+1 and so to compute rk we only need knowledge of r1,…,rk-1.The recursion works as follows: initiallyr1=q0, and then, for k=2,3,…, the following block is computed:(15)rk=∑h=2kqh-1ak-hh-1,ak-11=rk,ak-1h=∑j=0k-1ajh-1ak-1-j1=∑j=0k-1ajh-1rk-j,h=2,…,k-1,ahk=∑j=0hajk-1rh+1-jh=0,…,k,Note that a0k=a0k-1r1=a0k-1q0 and therefore a0k=q0k, where k is an upper index for a and a power exponent for q.We also derive from(16)H0x=xG0H1x,H0x=∑k≥1skxkthe expression (17)∑k≥1skxk=xp0+∑h=2nph-1xh∑j≥0rj+1xjh-1,so that(18)s1=p0,sk=∑h=2kph-1ak-hh-1,k≥2.In this case, the computation is straightforward, since it involves all previously computed quantities.Theoretically, the generating functionsH0(x) and H1(x) involve an infinite series, but obviously only a finite number of coefficients can be computed. Hence, the computation has to be stopped after having computed the desired number of terms rk and sk. Since each term is computed only once and it is not the result of subsequent smaller and smaller additions, truncating the computation up to a certain index has no effect on the accuracy of the values we compute. In other words, if we compute just a few terms, they are computed with the same accuracy as we had computed all coefficients.It is clear from the definitions and the previous iteration thatH1(x) implies H0(x); that is, once we know the rk values, the sk values are also known and implied by the rk values. It is not difficult to see that the converse is also true. By differentiating (6) and using (3), we get(19)H0′x=H0xx+xdG1H1xH1′x=H0xx+dH1xH1′x,and by integrating (19), we get(20)d2H12x=H0x-∫0xH0ξξdξ;that is,(21)d2∑k≥1rkxk2=∑k≥1skxk-∫0x∑k≥1skξk-1dξ,which leads to the following identities term by term:(22)d2∑h=1k-1rhrk-h=sk1-1k,k≥2.For k=2, we get in particular(23)r1=s2d,and, for k>2, we have(24)d22r1rk-1+∑h=2k-2rhrk-h=sk1-1k,which allows writing(25)rk-1=2sk/d1-1/k-∑h=2k-2rhrk-h2r1,so that all rk values can be recursively computed once we know d. We note that p1>0 implies that q0=r1>0. Hence, the recursion is well defined if p1>0, which is an almost necessary assumption if we investigate the presence of small components. Hence, knowledge of the sk values implies knowledge of the rk values.If we do not knowd, we may still compute d from the recursion. We first note that all rk depend on d through the factor 1/d. Therefore, we initially guess the value d=1 and compute tentative values r~k. Since ∑krk=u, we find the correct value for d as d=(∑kr~k/u)2 and so we have the correct values: (26)rk=r~k∑kr~ku.


## 4. Inferring the Degree Probabilities from the Component Size Probabilities

We may also consider the inverse problem of findingG0(x) and G1(x) from H0(x) and H1(x), that is, computing the degree distribution which gives rise to a particular small component size distribution. This problem presents interesting features. Arbitrary degree distributions of rk and sk may not be feasible; that is, there may be no degree distribution that can lead to those values.Formally, the recursion can be easily inverted; that is, knowing therk values, we can compute the pk values and qk values. Indeed, from (14), we have(27)rk=qk-1a0k-1+∑h=2k-1qh-1ak-hh-1;that is,(28)qk-1=rk-∑h=2k-1qh-1ak-hh-1a0k-1.Computing the ahk values is straightforward once we know the rk values. From the qk values, we easily deduce the pk values, apart from the fact that p0 cannot be derived from the qk values. However, p0=s1, and so it is known a priori. Note also that r1>0 implies that a0k-1>0.There is, however, a subtle point to be settled. Let us assume that a giant component may be present but we do not know theu and v values. Then, it is simpler to work with the conditional probabilities within the small components. Starting from the (conditional) sk probabilities, we compute the rk values as explained in the previous section but by using the normalization ∑k‍r~k=1. This way, we actually compute q^h and p^h. Then, from (7), we get ph and qh. The unknowns u and v are computed by imposing ∑h‍ph=1 and ∑h‍qh=1, which is equivalent to solving G1(u-1)=u-1 and G0(u-1)=v-1 with G0 and G1 defined on p^h and q^h.However, the inverse recursion is numerically unstable, and, unless we use exact data, it can produce absurd outcomes, like probabilities outside the range[0,1]. The reason of the instability is clear from (28), where we have a difference in the numerator and the denominator is getting smaller and smaller with q0k. As a simple exercise, suppose that we wonder which degree distribution gives rise to a size distribution of the small components of exponential type; that is, (29)tk=1-ββk-1,k≥1,with 0<β<1. Hence, we have t¯=1/(1-β) and(30)sk=ktkt¯=k1-β2βk-1.We remark that these sk are conditional probabilities. Now we have to compute the rk values from the sk values. As explained in the previous section, we initially fix d=1 and compute from (25) the tentative values:(31)r~k=2β1-ββk,for which ∑k‍r~k=2β. Hence, d=2β, implying the correct values: (32)rk=1-ββk-1.If we carry out the computation in (28) symbolically, we get(33)q^0=1-β,q^1=β,q^k=0,k>1,from which(34)p^0=1-β2,p^1=21-ββ,p^2=β2,p^k=0,k>2.From (7), we have(35)ph=vuhp^h,qh=1uh-1q^h.By imposing ∑h‍qh=1, we get(36)u1-β+β=1⟹u=1,and also v=1. Hence, there is no giant component in this case.Now assume that we have experimental data from which we infer the values:(37)t1=0.667499,t2=0.221782,t3=0.0739131,t4=0.024639,t5=0.00821253,t6=0.00273827,t7=0.000912502,t8=0.000304167.(these data have been generated by slightly perturbing the previous theoretical values with β=1/3). The previous computation leads to (38)q^0=0.667002,q^1=0.333269,q^2=0.0000610268,q3=-0.0000590538,q^4=0.0000833483,q^5=-0.000133577,q^6=0.00018942,q^7=-0.00542309.Not only are there negative values but also the absolute value of qk is increasing with k showing an amplifying effect of error propagation. Therefore, a lot of care should be exerted in order to carry out computations on experimental data. This can be matter of further investigation and is beyond the scope of this paper.We show a second example for the inverse computation. Assume that(39)r2k+1=3k+124k+2Ck,r2k+2=0,k=0,1,…,where Ck are the Catalan numbers. If we carry out the computation in (28) symbolically, we get(40)q^0=34,q^1=0,q^2=14,q^k=0,k>2,from which (41)p^1=910,p^2=0,p^3=110,p^k=0,k>3,so that(42)q0=34u,q1=0,q2=14u,qk=0,k>2,and(43)p1=9v10u,p2=0,p3=v10u3,pk=0,k>3.The normalization yields u=1/3 and v=5/27, so that p1=p3=1/2. Again, we show how perturbed data can lead to strange outcomes. If we perturb the data as (44)r1=0.746705,r3=0.141068,r5=0.0529938,r7=0.02473,r9=0.0130018,r11=0.00728396,r13=0.00430417,r15=0.00261343,r17=0.00163685,r19=0.00104289,with rk=0 for the other indices, we get(45)q^0=0.746705,q^2=0.253005,q^4=-0.000988595,q^6=-0.000552489,q^8=0.00141413,q^10=-0.00223664,q^12=0.00360276,q^14=-0.00597299,q^16=0.0098459,q18=-0.0155087.We see again the same inconsistencies and the amplifying effect. In any case, we may note that the values with odd index are correctly computed as null values.


## 5. Branching Processes

Now we defineqh, h=0,1,…, as the probability that a member of the population generates h offspring. We are interested in computing the probability rk that the population will eventually have k members, starting from a population consisting of one member. If G(x) and H(x) are the probability generating functions of qh and rk, respectively, then the following functional equation holds: (46)Hx=xGHx.Hence, the same relations (13) and (14) hold as well as the recursion (15). This time there are no sk coefficients to be computed and the picture is simplified. We may still view a branching process like a random graph. However, while in random graphs we pick up randomly any vertex within the small components, in branching processes, the small components are rooted trees and we pick up the roots randomly. Hence, the rk values we state here for a branching process can be related to the tk values of random graphs.Iteration (15) can also be carried out in exact arithmetic, thus producing results from which closed formulas can be inferred. As a simple example, suppose that q0=q1=q2=1/3 (qi=0 for i>2). Then, by applying (15), we obtain a sequence whose first terms are(47)13,19,227,481,127,7243,17729,1276561,32319683,83559049,2188177147,….We may guess that the denominator grows as the powers of 3, and so if we multiply the nth term of (47) by 3n, we obtain the new sequence (48)1,1,2,4,9,21,51,127,323,835,2188,….By looking at [7], we discover that these are the Motzkin numbers whose nth term is in closed form:(49)Mn=∑k≥0n/2Ckn2k,with Ck being the kth Catalan number. Hence,(50)rn=3-nMn-1,and we have found another combinatorial interpretation of the Motzkin numbers besides the many listed in [7]. The reason for n-1 as subscript is due to the fact that the first index of the sequence rn is n=1, whereas Motzkin numbers in (48) as defined above start from n=0.The same considerations about inferring the probabilitiesqk from the probabilities tk (=rk) can be applied also to branching processes. The example with tk=(1-β)βk-1 is almost trivial if we have in mind a branching process.


## 6. Conclusions

In this paper, we have presented an iterative scheme to compute the coefficients of a generating function that plays an important role in random graphs and in branching processes. The generating function is related to the population size probabilities for a branching process and to the small component size probabilities for random graphs. We also show that the iteration can be inverted; that is, for a branching process, from the population size probabilities, one can infer the offspring probabilities, but the inverse iteration is numerically unstable.



---

*Source: 3791075-2018-04-01.xml*