[Paper Review] Computing Least Fixed Points of Probabilistic Systems of Polynomials

Javier Esparza, Andreas Gaiser+1 more|arXiv (Cornell University)|Dec 21, 2009

Polynomial and algebraic computation7 references3 citations

TL;DR

This paper presents a strongly polynomial algorithm to decide whether the least fixed point of a probabilistic system of polynomials (PSP) equals the vector of all ones, and introduces the first practical algorithm for computing reliable, linearly converging lower and upper bounds on the least fixed point using a hybrid of exact and inexact arithmetic. The method efficiently handles numerical instability while ensuring correctness, outperforming prior LP-based approaches in experiments on branching process models.

ABSTRACT

We study systems of equations of the form X1 = f1(X1, ..., Xn), ..., Xn = fn(X1, ..., Xn), where each fi is a polynomial with nonnegative coefficients that add up to 1. The least nonnegative solution, say mu, of such equation systems is central to problems from various areas, like physics, biology, computational linguistics and probabilistic program verification. We give a simple and strongly polynomial algorithm to decide whether mu=(1, ..., 1) holds. Furthermore, we present an algorithm that computes reliable sequences of lower and upper bounds on mu, converging linearly to mu. Our algorithm has these features despite using inexact arithmetic for efficiency. We report on experiments that show the performance of our algorithms.

Motivation & Objective

• To develop a strongly polynomial algorithm to decide whether the least fixed point of a probabilistic system of polynomials (PSP) equals (1,…,1), which is critical for determining extinction probabilities in branching processes and termination in probabilistic programs.
• To design a practical algorithm that computes reliable sequences of lower and upper bounds on the least fixed point, converging linearly to the true value, even under inexact floating-point arithmetic.
• To overcome numerical instability in Newton’s method for PSPs by combining exact arithmetic for correctness with floating-point arithmetic for efficiency.
• To enable efficient computation of the critical radius in neutron branching processes and similar models, where the threshold for explosion or non-termination lies at the boundary of consistency.
• To demonstrate the practicality and superiority of the proposed method over existing LP-based approaches through experiments on real-world PSP instances.

Proposed method

• The consistency check reduces the problem to solving a system of linear equations, enabling a strongly polynomial algorithm that avoids the computational overhead of exact linear programming.
• For computing bounds, the algorithm uses a novel iterative procedure based on Newton’s method, but enhances it with exact arithmetic to maintain correctness while using floating-point arithmetic for speed.
• The method maintains guaranteed lower and upper bounds at each iteration by combining interval arithmetic with exact rational arithmetic for critical checks.
• It dynamically adjusts precision (via Maple’s Digits) during computation, increasing only when necessary to maintain correctness, thus balancing efficiency and reliability.
• The algorithm applies a contraction mapping principle to ensure linear convergence of both lower and upper bounds toward the least fixed point.
• Experiments use discretized neutron branching processes and stochastic context-free grammars to evaluate performance on real PSP instances.

Experimental results

Research questions

• RQ1 Can the consistency of a probabilistic system of polynomials (i.e., whether its least fixed point is (1,…,1)) be decided in strongly polynomial time?
• RQ2 Can a practical algorithm be designed to compute reliable, linearly converging lower and upper bounds on the least fixed point of a PSP, even when using inexact arithmetic?
• RQ3 How does the performance of the new consistency algorithm compare to existing LP-based methods in terms of speed and scalability on real-world PSPs?
• RQ4 To what extent can hybrid exact/inexact arithmetic be used to balance numerical efficiency and correctness in fixed-point computation for stochastic models?
• RQ5 What is the critical radius in neutron branching processes, and can it be computed efficiently using the proposed algorithms?

Key findings

• The proposed consistency algorithm runs up to one order of magnitude faster than the LP-based method on neutron branching process instances, with speedups of about 10× in tested cases.
• For PSPs with n=150 variables, the critical radius was computed to lie in the interval [2.981, 2.991], closely matching Harris’s estimate of 2.9.
• The consistency algorithm completed in less than 1 second for n=20, 50, and 100, while the exact LP method took up to 258 seconds for n=100, demonstrating superior scalability.
• The bound computation algorithm achieved linear convergence, with the number of accurate bits increasing linearly with iteration count, and required at most two precision increases (to Digits=20) in all tested cases.
• In all tested instances, the algorithm computed an upper bound strictly less than 1, proving inconsistency, within the first few iterations, confirming correctness and efficiency.
• The method successfully avoided the false positives common in floating-point Newton’s method, such as incorrectly reporting μ ≥ 1 when μ < 1, even for small systems with 7 variables.