Numerical Root Finding for Extinction Probabilities

Why Root-Finding Needs a Careful Bracket

For a Galton–Watson process, the extinction probability is characterized by a fixed-point equation for the offspring generating function. That equation also has the fixed point \(1\), because every probability generating function satisfies \(G(1)=1\). In a supercritical process, the extinction probability is smaller than \(1\), so a numerical method that simply searches for any root can return the wrong answer. The key is to use the fixed-point geometry to bracket the smallest root and then preserve that bracket throughout the computation.

Work on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). Let \((Z_n)_{n\ge0}\) be a Galton–Watson process in discrete time, with state space \(S=\{0,1,2,\ldots\}\) and initial population \(Z_0=1\). Each individual has an independent number \(\xi\) of offspring, with \(\mathbb{P}(\xi=k)=p_k\). Its offspring generating function is \(G(s)=\mathbb{E}[s^\xi]=\sum_{k\ge0}p_ks^k\), for \(0\le s\le1\). The extinction probability is \(q=\mathbb{P}(Z_n=0\text{ for some }n\ge0)\). The fixed-point characterization established in the earlier material on extinction says that \(q\) is the smallest solution in \([0,1]\) of \(G(s)=s\).

Put \(f(s)=G(s)-s\). A numerical root finder needs more than a formula for \(f\): it needs a bracket whose endpoints lie on opposite sides of the desired root. In the supercritical, nondegenerate case with \(p_0>0\) and finite mean \(m=\mathbb{E}\xi>1\), fixed-point geometry gives \(q\in(0,1)\), \(f(0)=p_0>0\), and \(f(s)<0\) for \(q<s<1\). Thus a valid bracket has left endpoint \(0\) and right endpoint \(b<1\) with \(f(b)<0\). The endpoint must be strictly below \(1\).

Definition: A sign-changing bracket for the interior extinction root is an interval \([a,b]\) with \(0\le a<b<1\), \(f(a)>0\), \(f(b)<0\), and exactly one zero of \(f\) in the interval. Under the supercritical, nondegenerate fixed-point geometry, that zero is \(q\). A maintained bracket is an interval with these sign conditions after every iteration.

There are important cases in which an interior root finder is unnecessary. If \(m\le1\), the extinction probability is \(1\), except for the deterministic one-child law \(\mathbb{P}(\xi=1)=1\), for which it is \(0\). If \(p_0=0\), every individual has at least one child, so the population never reaches zero and \(q=0\). These classifications, established in the earlier material on extinction, should be made before numerical root finding. In particular, do not try to estimate \(q=1\) by locating the root at the right endpoint of a bracket for an interior root.

Bisection with an Error Certificate

Bisection uses only evaluations of \(f\) and its sign. Given a valid bracket \([a,b]\), evaluate \(f(c)\) at \(c=(a+b)/2\). If \(f(c)=0\), the root has been found. If \(f(c)>0\), replace \(a\) by \(c\); if \(f(c)<0\), replace \(b\) by \(c\). The sign conditions and the root remain bracketed at every step.

Theorem (Bisection guarantee): Suppose \(f\) is continuous and has exactly one zero \(q\) in a bracket \([a_0,b_0]\), where \(f(a_0)>0\) and \(f(b_0)<0\). After \(n\) bisection steps, the maintained interval \([a_n,b_n]\) contains \(q\) and has width \((b_0-a_0)/2^n\). Its midpoint \(\widehat q_n=(a_n+b_n)/2\) satisfies \(|\widehat q_n-q|\le(b_0-a_0)/2^{n+1}\).

Proof. Initially the signs at the endpoints are opposite, and continuity gives a zero between them; uniqueness makes that zero \(q\). At each step, continuity and the sign at the midpoint show that the half-interval selected by the sign rule still has opposite endpoint signs and contains a zero. By uniqueness, it contains \(q\). Its width is half the previous width, so induction gives \(b_n-a_n=(b_0-a_0)/2^n\). Since \(q\in[a_n,b_n]\), its distance from the midpoint is at most half the interval width. This gives the stated error bound. \(\square\)

For an absolute error tolerance \(\varepsilon>0\), stop when \(b_n-a_n\le2\varepsilon\); then the midpoint has error at most \(\varepsilon\). This is a certificate based on the bracket, not a guess based on the size of \(f(\widehat q_n)\). A small residual alone need not imply a small root error when \(f\) is flat near its zero.

Worked Example: Bisection for a Three-Point Offspring Law

Let \(\mathbb{P}(\xi=0)=\tfrac14\), \(\mathbb{P}(\xi=2)=\tfrac12\), and \(\mathbb{P}(\xi=3)=\tfrac14\). The probabilities sum to \(\tfrac14+\tfrac12+\tfrac14=1\), and the mean is \(0\cdot\tfrac14+2\cdot\tfrac12+3\cdot\tfrac14=\tfrac74>1\). The generating function and root function are

G(s)=\tfrac14+\tfrac12s^2+\tfrac14s^3,\qquad f(s)=\tfrac14+\tfrac12s^2+\tfrac14s^3-s.

At \(a_0=0\), \(f(0)=\tfrac14>0\). At \(b_0=\tfrac12\), \(G(\tfrac12)=\tfrac14+\tfrac18+\tfrac1{32}=\tfrac{13}{32}\), so \(f(\tfrac12)=\tfrac{13}{32}-\tfrac{16}{32}=-\tfrac3{32}<0\). This is a valid bracket strictly below \(1\).

The first midpoint is \(c_1=\tfrac14\). Direct substitution gives \(G(\tfrac14)=\tfrac14+\tfrac1{32}+\tfrac1{256}=\tfrac{73}{256}\), so \(f(\tfrac14)=\tfrac{73}{256}-\tfrac{64}{256}=\tfrac9{256}>0\). The new bracket is \([\tfrac14,\tfrac12]\).

The next midpoint is \(c_2=\tfrac38\). Here \(G(\tfrac38)=\tfrac14+\tfrac9{128}+\tfrac{27}{2048}=\tfrac{683}{2048}\), while \(\tfrac38=\tfrac{768}{2048}\). Thus \(f(\tfrac38)=-\tfrac{85}{2048}<0\), and the new bracket is \([\tfrac14,\tfrac38]\).

The third midpoint is \(c_3=\tfrac5{16}\). Substitution gives \(G(\tfrac5{16})=\tfrac14+\tfrac{25}{512}+\tfrac{125}{16384}=\tfrac{5021}{16384}\), whereas \(\tfrac5{16}=\tfrac{5120}{16384}\). Hence \(f(\tfrac5{16})=-\tfrac{99}{16384}<0\), giving the bracket \([\tfrac14,\tfrac5{16}]\). After three steps its width is \(\tfrac5{16}-\tfrac4{16}=\tfrac1{16}\). The midpoint estimate is \(\tfrac9{32}\), with error at most \(\tfrac1{32}\). The signs verify that the computation has retained the smaller root, not the root \(1\).

Safeguarded Newton Iteration

Newton's method can converge quickly when its derivative-based proposal is good, but an unconstrained Newton step can leave the interval or converge to the root \(1\). A safeguard keeps every update inside a bracket. For a differentiable \(f\), start at the midpoint \(x=(a+b)/2\) and form \(y=x-f(x)/f'(x)\) if \(f'(x)\ne0\). Accept \(y\) only if it lies in the central portion \([a+(b-a)/4,\ b-(b-a)/4]\); otherwise use the midpoint. Evaluate \(f\) at the chosen point and update the endpoint with the same sign as that evaluation. If the evaluation is zero, stop.

Theorem (Safeguarded bracket guarantee): Suppose \(f\) is continuous and has a unique zero \(q\) in a strict sign-changing bracket \([a,b]\). At each iteration, choose a test point in \([a+(b-a)/4,\ b-(b-a)/4]\), using a Newton proposal when it lies there and the midpoint otherwise. After the sign-based endpoint update, the bracket still contains \(q\), and its width is at most \(\tfrac34(b-a)\).

Proof. The chosen test point \(x\) is inside the bracket, so continuity and the opposite endpoint signs imply that replacing the endpoint with the same sign as \(f(x)\) preserves a sign-changing interval containing a zero. Uniqueness identifies that zero as \(q\). Because \(x\) lies at least one quarter of the old width from each endpoint, either possible updated interval has width at most three quarters of the old width. If \(f(x)=0\), the exact root is found and no further update is needed. Repeating the argument preserves the root and reduces the width by at least the stated factor at each nonterminal step. \(\square\)

The safeguard gives a reliable width-based stopping rule even if Newton's proposals are poor. The derivative is used to suggest a point, not to certify the answer. For a generating function with finite offspring support, \(G\) is a polynomial and \(f'(s)=G'(s)-1\) can be evaluated directly.

Worked Example: One Safeguarded Newton Step

Let \(\mathbb{P}(\xi=0)=\tfrac13\) and \(\mathbb{P}(\xi=2)=\tfrac23\). The probabilities sum to \(1\), and \(m=0\cdot\tfrac13+2\cdot\tfrac23=\tfrac43>1\). Here \(G(s)=\tfrac13+\tfrac23s^2\) and \(f(s)=\tfrac13+\tfrac23s^2-s\). Use the bracket \([0,\tfrac34]\): \(f(0)=\tfrac13>0\), and \(f(\tfrac34)=\tfrac13+\tfrac38-\tfrac34=-\tfrac1{24}<0\).

The midpoint is \(x=\tfrac38\). Its function value is \(f(\tfrac38)=\tfrac13+\tfrac3{32}-\tfrac38=\tfrac5{96}\), and \(f'(s)=\tfrac43s-1\) gives \(f'(\tfrac38)=-\tfrac12\). The Newton proposal is therefore \(y=\tfrac38-(\tfrac5{96})/(-\tfrac12)=\tfrac{23}{48}\). The central portion of the bracket is \([\tfrac3{16},\tfrac9{16}]\), and \(\tfrac{23}{48}\) lies in this interval, so the proposal is accepted.

To determine which endpoint to update, substitute into \(f\). The factorization \(f(s)=\tfrac13(2s-1)(s-1)\) gives \(f(\tfrac{23}{48})=\tfrac13(-\tfrac1{24})(-\tfrac{25}{48})=\tfrac{25}{3456}>0\). The updated bracket is \([\tfrac{23}{48},\tfrac34]\). As a check, the factorization shows the roots are \(\tfrac12\) and \(1\); the maintained bracket contains \(\tfrac12\) and excludes \(1\), as required. The exact answer here is \(q=\tfrac12\), so the first safeguarded update has correctly narrowed around the smaller root.

Finding a Valid Bracket and Choosing a Stop Rule

When \(m>1\), \(p_0>0\), and \(m\) is finite, fixed-point geometry implies that \(f(s)<0\) for values below \(1\) sufficiently close to \(1\). One practical search is to test \(b_k=1-2^{-k}\) for increasing positive integers \(k\), stopping at the first certified evaluation with \(f(b_k)<0\). Since \(b_k<1\) and \(f(0)>0\), this supplies the required bracket \([0,b_k]\). The existence of such a point follows from the negative values of \(f\) just below \(1\). The search is not a license to evaluate \(f(1)\): that value is zero and does not bracket the smaller root against a negative value.

In exact arithmetic, the sign tests above are decisive. In floating-point arithmetic, cancellation can make a computed value have the wrong sign, especially when \(b\) is close to \(1\) and \(f(b)\) is small. For a certified computation, evaluate \(f\) with enough precision or use interval arithmetic: accept a positive sign only when the computed interval lies wholly above zero, and a negative sign only when it lies wholly below zero. If the interval includes zero, increase precision or choose another test point. For infinite offspring support, a truncated sum for \(G(s)\) also needs a bound on the omitted tail before its sign can be certified.

Use the bracket width for the primary stopping criterion. If the requested absolute error is \(\varepsilon\), stop when the width is no more than \(2\varepsilon\) and return the midpoint. A relative tolerance can be added when \(q\) is known to be bounded away from zero, but an absolute tolerance is essential when \(q=0\) or could be very small. A tiny residual \(|f(\widehat q)|\) is not, by itself, a bound on \(|\widehat q-q|\).

Worked Example: Classify Before Iterating

Suppose \(\mathbb{P}(\xi=0)=\tfrac12\) and \(\mathbb{P}(\xi=1)=\tfrac12\). The probabilities sum to \(1\), and the mean is \(0\cdot\tfrac12+1\cdot\tfrac12=\tfrac12<1\). The offspring generating function is \(G(s)=\tfrac12+\tfrac12s\), so \(f(s)=\tfrac12-\tfrac12s\). The subcritical classification gives \(q=1\); there is no interior root to locate.

Substitution checks the answer: \(G(1)=\tfrac12+\tfrac12=1\), so \(f(1)=0\). In contrast, when \(\mathbb{P}(\xi=1)=1\), the mean is \(1\) and \(G(s)=s\) for every \(s\in[0,1]\). The process never becomes extinct from one initial individual, so \(q=0\); the fixed-point equation alone has infinitely many solutions and cannot identify \(q\) without the smallest-root rule. These checks illustrate why classification and fixed-point selection precede numerical iteration.

A Common Mistake

Wrong calculation: For the law with probabilities \(\tfrac13\) at zero and \(\tfrac23\) at two, evaluate \(f(1)=G(1)-1=0\) and report \(q=1\). This finds a fixed point but not the smallest fixed point. The correct calculation uses \(f(s)=\tfrac13(2s-1)(s-1)\), whose roots are \(\tfrac12\) and \(1\). Since \(f(0)=\tfrac13>0\) and \(f(\tfrac34)=-\tfrac1{24}<0\), the bracket \([0,\tfrac34]\) contains the smaller root \(\tfrac12\) and excludes the trivial root \(1\).

The general rule is to classify the process first, then construct a bracket strictly below \(1\) whenever the extinction probability is interior. Bisection provides a simple, derivative-free error certificate; safeguarded Newton iteration can improve the proposals while retaining that certificate. Neither method should be allowed to turn the known fixed point \(1\) into the answer by default.

Check Your Understanding

For each problem, identify the relevant root or bracket before deciding what a numerical method should return.

For \(G(s)=\tfrac14+\tfrac12s^2+\tfrac14s^3\), verify the signs at \(0\) and \(\tfrac12\), then state the midpoint error bound after four bisection steps.
Why is \([0,1]\) not a valid strict sign-changing bracket for an interior extinction root?
In safeguarded Newton iteration, what should be done if \(f'(x)=0\) or the Newton proposal lies outside the central portion of the bracket?
For \(\mathbb{P}(\xi=0)=\tfrac13\), \(\mathbb{P}(\xi=2)=\tfrac23\), verify the smaller root and explain why returning \(1\) is incorrect.
Why can a small value of \(|f(\widehat q)|\) fail to certify that \(\widehat q\) is close to \(q\)?