Madupite: Distributed High-Performance MDP Solver

• Madupite is an open-source, distributed solver for exact dynamic programming of finite-state, discounted MDPs using C++20/PETSc/MPI with a Python API.
• It employs inexact policy iteration with configurable Krylov methods to efficiently tackle near-undiscounted regimes and mitigate the curse of dimensionality.
• The solver leverages distributed memory and sparse computation to handle models with millions of states, targeting applications in finance, epidemiology, and traffic control.

Searching arXiv for madupite and directly related iPI papers to ground the article in current sources. arXiv search query: "madupite" Madupite is an open-source, distributed high-performance solver for large-scale discounted infinite-horizon Markov decision processes with finite state and action spaces, implemented as a C++20/PETSc/MPI core with an equivalent Python API via nanobind and packaged for pip installation (Gargiani et al., 20 Feb 2025, Gargiani et al., 30 Jul 2025). It is intended for exact dynamic-programming solution of explicit finite MDP models that are too large for desktop-oriented tools, including problems with more than a million states and near-undiscounted regimes in which classical value-iteration-like methods slow down sharply (Gargiani et al., 20 Feb 2025, Gargiani et al., 30 Jul 2025).

1. Scope and motivation

Madupite was developed to address a specific mismatch between the modeling utility of MDPs and the computational limits of conventional software. The motivating domains named in the papers include finance, epidemiology, and traffic control. In these settings, state cardinality can grow exponentially with problem parameters, producing the familiar curse of dimensionality. The initial presentation emphasizes that approximation methods are often adopted not because they are intrinsically preferable, but because exact solution at scale is otherwise computationally out of reach; in the authors’ wording, this comes “at the price of introducing sub-optimality with respect to the original solution” (Gargiani et al., 20 Feb 2025).

The problem class is deliberately narrow. Madupite targets discounted, infinite-horizon, finite-state, finite-action MDPs, and it is designed for exact solution of the explicit finite model rather than for function approximation, neural approximation, or model-free reinforcement learning. A central practical claim is that distributed-memory execution lets the solver handle MDPs whose transition and cost representations exceed the memory of a single machine, including instances running on multiple nodes of a high-performance computing cluster (Gargiani et al., 20 Feb 2025).

The later technical-design paper sharpens this positioning by adding another axis: ill-conditioning induced by discount factors close to $1$. In that regime, the authors argue that solver design should be understood not only as a software-engineering problem, but also as a numerical linear algebra problem. This motivates madupite’s combination of exact dynamic programming, inexact policy iteration, and distributed sparse computation (Gargiani et al., 30 Jul 2025).

2. Mathematical formulation

The formal model used in the technical-design paper is the finite discounted MDP tuple

$$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$

where $\mathcal{S}=\{1,\dots,n\}$ is the finite state space, $\mathcal{A}=\{1,\dots,m\}$ is the finite action space, $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$ is the transition probability function, $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$ is the stage-cost function, and $\gamma\in(0,1)$ is the discount factor (Gargiani et al., 30 Jul 2025). The formulation is explicitly cost-minimizing rather than reward-maximizing, with

$$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$

For initial state $s$, the objective is to minimize

$$\lim_{T\rightarrow\infty}\mathbb{E}\left[\sum_{t=0}^{T-1}\gamma^t g(s_t,a_t)\,\middle|\, s_0=s\right].$$

A policy is a stationary deterministic mapping

$$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$0

and the number of such policies is finite, with

$$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$1

This finiteness explains why brute-force enumeration is only possible for very small problems.

The matrix-vector representation is central to madupite. For a policy $$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$2, define $$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$3 by $$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$4 and $$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$5 by $$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$6. Policy evaluation is then the linear system

$$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$7

The optimal value function satisfies

$$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$8

With Bellman operator $$\left\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\},$$9 defined by

$\mathcal{S}=\{1,\dots,n\}$0

the optimal value is the fixed point $\mathcal{S}=\{1,\dots,n\}$1. For a fixed policy $\mathcal{S}=\{1,\dots,n\}$2, the evaluation operator is

$\mathcal{S}=\{1,\dots,n\}$3

and $\mathcal{S}=\{1,\dots,n\}$4.

Residuals are used as approximation diagnostics: $\mathcal{S}=\{1,\dots,n\}$5 A policy is greedy with respect to $\mathcal{S}=\{1,\dots,n\}$6 if

$\mathcal{S}=\{1,\dots,n\}$7

The papers do not develop or use a linear-programming formulation; the solver is built around Bellman equations, policy iteration variants, and iterative linear algebra (Gargiani et al., 30 Jul 2025).

3. Algorithmic basis

Madupite is built around inexact policy iteration (iPI). Classical policy iteration alternates between policy improvement and exact policy evaluation. In madupite’s iPI formulation, the policy evaluation stage is solved only approximately, so scalability is gained without replacing the explicit finite MDP by a function-approximation surrogate (Gargiani et al., 30 Jul 2025).

The main practical consequence is that the expensive kernel becomes a configurable sparse linear solve for

$\mathcal{S}=\{1,\dots,n\}$8

The papers emphasize that different inner solvers can materially change performance, especially when $\mathcal{S}=\{1,\dots,n\}$9 is close to $\mathcal{A}=\{1,\dots,m\}$0. The solver repertoire explicitly mentioned includes GMRES, BiCGStab, TFQMR, Richardson iteration, SOR-related options via preconditioning, and CG in compatible settings. The principal claim is not merely that madupite implements iPI, but that it lets the user tune inexactness and choose the inner linear algebra strategy to exploit problem structure.

The reported user-facing controls include -ksp_type, -pc_type, -atol_pi, -alpha, -max_iter_pi, -max_iter_ksp, and PETSc solver-specific options such as Richardson step size. Main defaults are -max_iter_pi = 1000, -max_iter_ksp = 1000, -atol_pi = 1e-8, -alpha = 1e-4, -ksp_type = gmres, and -pc_type = none. Within this framework, madupite can emulate classical dynamic-programming schemes by parameter choice, including Value Iteration, Policy Iteration, Optimistic Policy Iteration, $\mathcal{A}=\{1,\dots,m\}$1-Value Iteration, Gauss-Seidel VI, and Jacobi VI.

A recurring point in both papers is the meaning of “exact.” Madupite is an exact dynamic-programming solver in the sense that it solves the explicit finite MDP model rather than a function-approximation proxy. Numerically, however, policy evaluation is still performed by iterative methods up to prescribed tolerances. This distinction is central: the exactness is model-level rather than symbolic or closed-form (Gargiani et al., 20 Feb 2025, Gargiani et al., 30 Jul 2025).

4. Software architecture and distributed implementation

Madupite’s architecture combines a fully featured C++20 library with an equivalent Python API exposed through nanobind. The interface is intended to serve both native C++ integration and Python-based prototyping, scripting, and experimentation. Users can create an MDP either by reading compliant binary files or by generating local chunks directly from Python callback functions, which enables distributed model construction without first assembling the entire model on a single machine (Gargiani et al., 30 Jul 2025).

The implementation relies on PETSc for distributed sparse data structures and linear solvers, and on MPI for distributed-memory execution. PETSc Mat stores matrices, Vec stores vectors, and IS stores index sets used for policies. Transition probabilities are stored in sparse CSR-like distributed format via MATMPIAIJ, while stage costs are stored in dense format.

Because PETSc does not natively store 3D tensors, the transition tensor

$\mathcal{A}=\{1,\dots,m\}$2

is flattened to a row-stacked matrix of size

$\mathcal{A}=\{1,\dots,m\}$3

The authors report that this layout was chosen after profiling: a more natural $\mathcal{A}=\{1,\dots,m\}$4 arrangement would be less favorable because PETSc’s sparse row access is much more efficient than column access in CSR format. The chosen layout therefore improves greedy-policy extraction and assembly of policy-specific matrices.

State ownership is balanced across $\mathcal{A}=\{1,\dots,m\}$5 MPI ranks. For rank $\mathcal{A}=\{1,\dots,m\}$6, the local state count is

$\mathcal{A}=\{1,\dots,m\}$7

Rank $\mathcal{A}=\{1,\dots,m\}$8 then stores $\mathcal{A}=\{1,\dots,m\}$9 rows of the flattened transition matrix and $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$0 rows of the stage-cost matrix. All action-conditioned transitions for a given state are kept on the same rank, which makes greedy policy extraction largely local and reduces communication (Gargiani et al., 30 Jul 2025).

The main distributed kernels are local greedy policy extraction, PETSc KSP-based distributed solution of

$$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$1

and global convergence tests. Communication overhead is described as relatively low but not negligible; synchronization and communication account for the nonparallel component in the Amdahl analysis. For Krylov-based policy evaluation, madupite often avoids explicit formation of $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$2 by using a matrix shell. An exception is the SOR preconditioner, which requires explicit coefficient entries and therefore usually reduces efficiency. The implementation also replaces the theoretical inner stopping condition based on the infinity norm of the policy residual with a more stringent but cheaper criterion based on the 2-norm of the linear residual, because PETSc already computes that quantity internally (Gargiani et al., 30 Jul 2025).

5. Position relative to prior software

The papers position madupite against two prior tools: pymdptoolbox and mdpsolver. The comparison is qualitative in the initial paper and benchmarked in the later one. pymdptoolbox is characterized as plain Python, with limited scalability and no support for parallel or distributed computing. mdpsolver is described as a C++ solver with a Python API and parallel-computing support, but the initial paper identifies two limitations: methodological restriction and data-structure design (Gargiani et al., 20 Feb 2025).

On the methodological side, the initial paper states that mdpsolver is limited to modified policy iteration and argues that this performs poorly for an important class of practical problems. On the implementation side, it criticizes storage of matrix values and indices in nested std::vectors regardless of sparsity, which, according to the authors, prevents effective use of optimized sparse linear algebra routines. Madupite’s response is to support a wider class of methods, particularly configurable iPI, while delegating sparse distributed storage and solver infrastructure to PETSc.

A broader interpretive point follows from the two papers taken together. Madupite is presented not simply as a faster implementation of familiar value-iteration-like workflows, but as a reorientation of exact MDP solution toward distributed sparse numerical linear algebra. This suggests that its main novelty lies less in a new Bellman-theoretic formulation than in the systems-level combination of iPI, PETSc/MPI, sparse data structures, and a dual Python/C++ interface.

6. Performance evidence and application benchmarks

The evidentiary base is split across the two 2025 papers. The initial introduction states qualitative scalability aims, including solution of MDPs with “more than a million states,” but does not provide benchmark tables, speedup curves, or runtime measurements. The later technical-design paper adds numerical benchmarks, comparative runtimes, and application case studies on the Euler HPC cluster at ETH Zurich and on an Intel Core i7-10750H laptop CPU (Gargiani et al., 20 Feb 2025, Gargiani et al., 30 Jul 2025).

Scenario	Setup	Reported result
Strong scaling synthetic MDP	$$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$3, $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$4, $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$5, about 1.5 billion nonzeros, 1–48 ranks, GMRES, $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$6	Amdahl estimate $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$7; implied maximum speedup $$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$8
Baseline solver comparison	$$P:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\to[0,1]$$9, $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$0, $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$1, tolerance $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$2, 10 runs	pymdptoolbox $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$3 s; mdpsolver $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$4 s; madupite(1) $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$5 s; madupite(4) $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$6 s
Memory feasibility benchmark	Dense synthetic MDP, $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$7, $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$8, varying $g:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$9	At $\gamma\in(0,1)$0, transition matrix alone requires about 74 GiB
SIS epidemiology model	$\gamma\in(0,1)$1, $\gamma\in(0,1)$2, $\gamma\in(0,1)$3, GMRES, 32 cores	Populations up to 1,000,000; largest instance solved in less than 10 s
Inverted pendulum	$\gamma\in(0,1)$4, $\gamma\in(0,1)$5 from 11 to 3163, $\gamma\in(0,1)$6, $\gamma\in(0,1)$7, $\gamma\in(0,1)$8, TFQMR, 32 cores	Up to $\gamma\in(0,1)$9 states; solved in less than 1 hour
Deterministic 2D maze	$$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$0, TFQMR, $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$1, 32 cores	Approximately $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$2 states solved in under 1.6 minutes

Several technical conclusions are explicit. First, strong scaling is substantial but not ideal, with the nonparallel fraction attributed mainly to communication and synchronization. Second, customization matters most in near-undiscounted settings: Richardson and MINRES-based variants degrade strongly as $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$3 increases, whereas Krylov methods such as GMRES remain much less affected. Third, laptop parallelism is not uniformly beneficial; for small problems, 6 cores can be no faster, or slower, than 1 core due to overhead and shared-resource contention.

The application cases illustrate different structural regimes. The SIS model uses a 20-action product space and a superlinear health-cost term $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$4, with the disease-free state absorbing. The inverted-pendulum example begins as a continuous control problem,

$$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$5

then time-discretizes with $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$6 and spatially discretizes $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$7 and $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$8 into a finite MDP. The maze example emphasizes sparse deterministic transitions with a wall penalty implemented as $$\sum_{s'\in\mathcal S} P(s,s',a)=1 \qquad \forall s\in\mathcal S,\ a\in\mathcal A.$$9 and a goal cost typically $s$0. The pendulum case is singled out as particularly difficult, requiring about 130 iPI outer iterations, compared with a more typical average of fewer than 20 (Gargiani et al., 30 Jul 2025).

7. Limitations and practical positioning

Madupite’s scope is intentionally restricted. It is presented as a solver for discounted infinite-horizon finite-state finite-action MDPs. The papers do not claim direct support for continuous-state MDPs, continuous-action MDPs, POMDPs, average-cost MDPs, risk-sensitive or risk-averse MDPs, or model-free RL. Continuous problems must first be discretized, and in that case madupite solves the discretized model exactly rather than the original continuous system (Gargiani et al., 30 Jul 2025).

The curse of dimensionality remains operative. Distributed memory and better algorithms mitigate it, but they do not remove it; very large action spaces and very fine discretizations still imply very large sparse operators. The later paper also notes that generic PETSc data structures were chosen to preserve broad applicability, even though specialized storage formats might improve performance in particular problem classes. Communication and synchronization impose a measurable scaling ceiling, reflected in the estimated nonparallel fraction of about $s$1. The discussion of preconditioners is comparatively limited, and the initial paper does not discuss OpenMP or GPU support (Gargiani et al., 20 Feb 2025, Gargiani et al., 30 Jul 2025).

Within those boundaries, madupite is best suited to problems for which the MDP model is explicit, discounted, finite, and large enough that exact solution quality matters but single-node memory and value-iteration-like methods are inadequate. The papers repeatedly position it as preferable when one wants quantifiable near-optimality for the explicit finite model, when $s$2 is close to $s$3, when sparse distributed exact methods remain feasible, and when HPC resources are available. A plausible implication is that madupite occupies a middle ground between small desktop-oriented tabular solvers and approximation-heavy reinforcement-learning workflows: it expands the range of problems for which exact finite-MDP solution is practical, but it does not supersede approximation when the explicit model is itself too large to store or unavailable (Gargiani et al., 30 Jul 2025).