🚀 Solving at scale: many cores, many machines

Homotopy continuation is pleasantly parallel. Every path is independent – once the start points are fixed, tracking solution #3 needs nothing from solution #57. So the wall-clock cost of a solve is, in principle, the cost of the slowest single path, no matter how many paths there are, if you have enough workers to go around.

Bertini turns that “in principle” into a one-line change. This tutorial shows the two knobs:

• MPI ranks – a manager-worker pool, possibly spread across many machines. Pass a communicator to solve() and rank 0 hands paths out to the rest.

• threads per rank – each worker can track several paths at once. Set OMP_NUM_THREADS.

The two scale multiplicatively: R worker ranks of T threads each track R × T paths at a time. We demonstrate on two driving problems – a cyclic-n system (hundreds of cheap paths) and an eigenvalue solve (a few dozen heavy, uneven paths) – because they stress the scheduler in opposite ways, and the same manager-worker pool handles both.

Note

The timing tables below were measured on a 16-core Apple M3 Max (12 performance + 4 efficiency cores). Your numbers will differ; the shape (speedup grows with workers, then flattens once the single slowest path dominates the wall) is the point. Run tools/update_scaling_timings.py to repopulate every table for your own hardware – it pins a fixed seed so every layout solves the identical problem.

The manager-worker model

A bertini solve already has a definite shape – track every path before the endgame, run a midpath check, then run the endgame on each path – and the parallel solve keeps exactly that shape. There is no separate “parallel solver”: you call the same solve(), only with a communicator.

• Rank 0 is the manager. It owns the bookkeeping – the list of path indices still to do, and the results as they come back – but it does not track paths itself. It hands out work on demand: when a worker reports a finished path, the manager sends it the next index. That demand-driven dispatch is what gives automatic load balancing – a worker that draws a cheap path simply comes back for more sooner.

• Ranks 1..N-1 are workers. Each asks the manager for a path index, tracks that path through the whole pre-endgame / midpath / endgame pipeline, ships the result back, and asks again, until the manager says “done”.

Because the manager does not track, a distributed run needs at least two ranks (one manager + one worker). The serial baseline is just solve() with no communicator. The minimal program is:

from mpi4py import MPI
import bertini as pb

comm = MPI.COMM_WORLD
solver = pb.ZeroDimSolver(my_system, mptype='double')

if comm.Get_size() > 1:
    solver.solve(communicator=comm)      # rank 0 dispatches; the others track
else:
    solver.solve()                       # serial: a lone manager would have no workers

if pb.parallel.is_manager():             # only rank 0 has the collected solutions
    print(len(solver.all_solutions()), "solutions")

The bertini.parallel helpers – rank(), size(), is_manager(), is_worker() – let each rank know its role. Only the manager ends up holding the solution list, so guard your reporting and your correctness checks with is_manager().

Launch it with mpirun:

python            my_solve.py     # serial baseline (1 process)
mpirun -n 2 python my_solve.py    # 1 manager + 1 worker  (the smallest parallel run)
mpirun -n 9 python my_solve.py    # 1 manager + 8 workers

Warning

mpirun gives stdin only to rank 0, so a heredoc (mpirun -n 4 python << EOF) deadlocks – the workers wait forever for a script they never receive. Always run a file.

Note

To put one worker on each of your cores you want -n (cores + 1) – the extra rank is the near-idle manager. On a machine with exactly that many cores mpirun will refuse to start (“not enough slots”); add --map-by :OVERSUBSCRIBE to let the manager share a core (it barely uses one). For example, to run 8 workers on an 8-core box:

mpirun -n 9 --map-by :OVERSUBSCRIBE python my_solve.py

Scaling a cyclic system

The cyclic-\(n\) system is a classic benchmark: \(n\) polynomials whose total-degree homotopy tracks \(n!\) paths, most of which diverge, leaving a known number of finite solutions. It is the “many cheap paths” regime – cyclic-6 is 720 paths – so it shows off fine-grained load balancing.

The complete, runnable script is solve_cyclic.py:

examples/solve_cyclic.py

#!/usr/bin/env python
"""Solve the cyclic-n system, distributed across MPI ranks (with optional threads per rank).

The cyclic-n polynomials are a classic benchmark for polynomial system solving.  Total-degree
homotopy tracks one path per Bezout count (the product of the degrees), and the work is
pleasantly parallel: each path is independent.  Bertini hands the paths out from rank 0 (the
manager) to the worker ranks, and -- if OMP_NUM_THREADS > 1 -- each worker runs several tracking
threads.

Run it::

    python solve_cyclic.py --n 6                              # serial baseline
    mpirun -n 5 python solve_cyclic.py --n 6                  # 1 manager + 4 workers
    OMP_NUM_THREADS=4 mpirun -n 3 --bind-to none python solve_cyclic.py --n 6   # 2 workers x 4 threads

Only rank 0 prints; it checks the distinct finite-solution count -- taken from the solver's own
solution metadata -- against the known cyclic-n value.
"""

import argparse
import math
import os
import time

import numpy as np
from mpi4py import MPI

import bertini as pb

# distinct finite (nonsingular) solution counts of the cyclic-n system
KNOWN_FINITE = {4: 16, 5: 70, 6: 156, 7: 924, 8: 1152}


def cyclic_system(n):
    """The cyclic-n system: the n-1 cyclic sums of products, plus (prod - 1)."""
    x = [pb.Variable('x{}'.format(i)) for i in range(n)]
    w = x + x  # doubled, to take products that wrap around
    sys = pb.System()
    for length in range(1, n):                       # degree-`length` cyclic sum
        sys.add_function(np.sum([np.prod(w[start:start + length]) for start in range(n)]))
    sys.add_function(np.prod(x) - 1)                 # the product, normalized
    sys.add_variable_group(pb.VariableGroup(x))
    return sys


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--n', type=int, default=6, help='which cyclic-n (default 6)')
    parser.add_argument('--seed', type=int, default=None,
                        help='optional RNG seed for a reproducible run (default: random each run); '
                             'the same seed on every rank fixes the homotopy, so serial and '
                             'distributed runs are directly comparable')
    args = parser.parse_args()

    comm = MPI.COMM_WORLD
    if args.seed is not None:
        pb.random.set_random_seed(args.seed)
    system = cyclic_system(args.n)
    solver = pb.ZeroDimSolver(system, endgame='cauchy', mptype='adaptive', startsystem='binomial')

    # Adaptive precision is what makes this reliable.  cyclic-n has a few near-singular paths that
    # sit right at the edge of double-precision tracking tolerance; in double precision one of them
    # occasionally fails the endgame (MinStepSizeReached) and a genuine root is silently lost -- so
    # the distinct count comes back 1 or 2 short of the known value.  Adaptive precision raises the
    # working precision on exactly those paths, so every run finds all the finite solutions.

    start = time.time()
    if comm.Get_size() > 1:
        solver.solve(communicator=comm)           # manager (rank 0) hands paths to the workers
    else:
        solver.solve()                            # serial: a lone manager would have no workers
    elapsed = time.time() - start

    if not pb.parallel.is_manager():
        return

    # Correctness, not just bookkeeping.  The solver's own report classifies every path; we trust it
    # rather than rolling our own cutoff.  We print it whenever anything is flagged, so problems are
    # surfaced -- but we gate correctness on the ground truth, distinguishing two very different flags.
    report = solver.report()

    print('cyclic-{}:  ranks={}  threads/rank={}  paths tracked={}  finite solutions={}  wall={:.1f}s'.format(
        args.n, comm.Get_size(), os.environ.get('OMP_NUM_THREADS', '1'),
        report.num_paths_tracked, report.num_finite_solutions, elapsed))

    if not report.all_paths_resolved:        # a failed path OR an unresolved crossing -- show what and why
        print(report)

    assert report.num_paths_tracked == math.factorial(args.n), \
        'expected {} paths, got {}'.format(math.factorial(args.n), report.num_paths_tracked)
    if args.n in KNOWN_FINITE:
        # A FAILED path (the tracker gave up) means a genuinely lost root -- a hard error.  An
        # UNRESOLVED CROSSING is a conservative "may be wrong" warning the midpath check couldn't
        # clear; the count can still be exactly right (and here it is), so it is reported above, not
        # asserted.  Gate on no failed paths + the known finite count.
        assert report.num_failed == 0, 'some paths failed to track -- see the report above'
        assert report.num_finite_solutions == KNOWN_FINITE[args.n], \
            'expected {} finite solutions, got {}'.format(KNOWN_FINITE[args.n], report.num_finite_solutions)
        print('  OK: {} paths, {} finite solutions -- matches the known cyclic-{} count'.format(
            report.num_paths_tracked, report.num_finite_solutions, args.n))


if __name__ == '__main__':
    main()

Two things to notice. First, the whole distributed machinery is the single solver.solve(communicator=comm) line – everything else is building the system and checking the answer. Second, the check is real and uses the solver’s own report rather than a hand-rolled cutoff: report.num_finite_solutions is the number of distinct finite solutions – the report sums \(1/\text{multiplicity}\) over the finite endpoints, so a genuine multiple root (several paths converging to one true solution of higher multiplicity) is counted once – and we confirm it equals the mathematically known cyclic-\(n\) value. Just as important, we assert report.all_paths_resolved: had any path failed to track, a genuine root would be silently missing and the count would come up short, so the report catches the loss rather than letting a quiet 69 pass for 70. (A coincident pair of distinct paths is the one other way two endpoints meet – see the note below.)

Note

Two distinct paths landing on the same point is a different matter entirely. Because the homotopy’s \(\gamma\) is random, distinct paths meeting is a probability-0 event: if it happens it is a path crossing – a sign the tracker under-resolved the paths – not a benign duplicate. The solver checks for exactly this at the endgame boundary and re-tracks the offending paths; with the default predictor and tolerances it essentially never triggers. The 🪢 Crossed paths: when two paths become one tutorial provokes one on purpose and shows what to do about it.

Correctness, not just a path count

A parallel solve is only useful if it gives the same answer as the serial one. This is subtler than it sounds: every rank constructs the homotopy independently, and the start system, patch, and \(\gamma\) are all random. If the ranks disagree on those random choices, “path #7” means a different path on each worker and the manager stitches together nonsense – the solve finishes, reports a plausible count, and is silently wrong. Bertini guards against this by broadcasting one rank’s random seed to all of them at the start of the parallel solve, so every rank forms the identical homotopy. The cyclic check above – distinct finite count equals the known value – is what catches a regression here, which is why one of these two examples must verify results and not merely tally paths.

Run the ladder (serial, then 2, 4, 8 workers):

$ python python/examples/solve_cyclic.py --n 6              # serial baseline
$ mpirun -n 3 python python/examples/solve_cyclic.py --n 6  # 2 workers
$ mpirun -n 5 python python/examples/solve_cyclic.py --n 6  # 4 workers
$ mpirun -n 9 python python/examples/solve_cyclic.py --n 6  # 8 workers
$ mpirun -n 13 python python/examples/solve_cyclic.py --n 6 # 12 workers

cyclic-6 (720 paths, 156 finite solutions), measured on a 16-core Apple M3 Max (12 performance + 4 efficiency cores), 2026-06-30

launch	workers	wall-clock (s)	speedup
serial	–	160.1	1.0x
-n 3	2	89.5	1.8x
-n 5	4	46.0	3.5x
-n 9	8	29.4	5.4x
-n 13	12	23.7	6.8x

The speedup grows steadily with the worker count – here to about 6.8x at 12 workers. cyclic-6 is 720 paths, so there is plenty to spread, and it keeps gaining well past the point where the eigenvalue solve below has flattened. It still falls short of the ideal 12x for two reasons: the slowest single path from the opening (a few near-singular paths grind far longer than the rest, so once the workers outnumber the heavy paths the extra ones idle), and the hardware – only 12 of the 16 cores are performance cores, so beyond a dozen workers the slower efficiency cores set the floor. The largest jumps are early (2->4->8 workers), where spreading the few heavy paths off a shared worker buys the most; the 8->12 step gains less. (Note that -n 13 is 13 processes – 12 workers plus the near-idle manager – which fit the 16 cores without oversubscription; the manager, which only dispatches, costs essentially nothing.)

Scaling the eigenvalue solve

The eigenvalue solve from 🧮 Eigenvalues by homotopy continuation is the opposite regime: a few paths, each expensive and of uneven cost. Writing the eigenvector projectively makes the multihomogeneous Bézout number exactly \(n\) – one path per eigenvalue – so an \(n \times n\) matrix gives only \(n\) paths. To have work worth distributing you scale the matrix, not the path count.

The complete script is solve_eigenvalues.py:

examples/solve_eigenvalues.py

#!/usr/bin/env python
"""Find ALL eigenvalues of a matrix as a distributed multihomogeneous homotopy solve.

The eigenproblem ``(A - lam I) x = 0`` is the textbook case where multihomogeneous structure
crushes total degree: write the eigenvector ``x`` in projective space ``P^{n-1}`` and ``lam`` as
an affine unknown, and the multihomogeneous Bezout number is exactly ``n`` -- one path per
eigenvalue -- versus the total-degree bound ``2^n``.  So an ``n x n`` matrix gives ``n`` paths,
and those ``n`` independent paths spread across MPI workers (make ``n`` large enough that there
is real work to hand out).

Run it::

    python solve_eigenvalues.py --size 24                            # serial baseline
    mpirun -n 5 python solve_eigenvalues.py --size 24                # 1 manager + 4 workers
    OMP_NUM_THREADS=4 mpirun -n 3 --bind-to none python solve_eigenvalues.py --size 24

Only rank 0 prints; it checks the recovered eigenvalues against ``numpy.linalg.eigvals`` --
correctness, not just a path count.
"""

import argparse
import os
import time

import numpy as np
from mpi4py import MPI

import bertini as pb
from bertini import ZeroDimSolver

OK = int(pb.SuccessCode.Success)


def random_symmetric_integer_matrix(n, seed, spread=9):
    """A reproducible n x n symmetric integer matrix (real, generically distinct spectrum).

    The same ``seed`` builds the same matrix on every rank -- essential, since each rank
    constructs its own copy of the system to solve.
    """
    rng = np.random.default_rng(seed)
    M = rng.integers(-spread, spread + 1, size=(n, n))
    return M + M.T                       # symmetric => real eigenvalues


def eigen_system(A):
    """(A - lam I) x = 0 with x a PROJECTIVE group and lam affine.

    Projective ``x`` lives in ``P^{n-1}`` natively, so the eigenvector scale is already fixed --
    no normalization equation, and the multihomogeneous Bezout number is exactly ``n``.
    """
    n = A.shape[0]
    x = np.array(pb.variables('x', n), dtype=object)
    lam = pb.Variable('lam')
    sys = pb.System()
    sys.add_functions(A @ x - lam * x)                  # the rows of (A - lam I) x
    sys.add_hom_variable_group(pb.VariableGroup(list(x)))    # eigenvector in P^{n-1}
    sys.add_variable_group(pb.VariableGroup([lam]))
    return sys


def recovered_eigenvalues(solver):
    """lam is the last dehomogenized coordinate; collect it from each successful, finite path.

    We trust the solver's own classification: a real eigenvalue is the lam coordinate of a path
    whose endgame succeeded and that the library calls finite (is_finite) -- no hand-rolled cutoff.
    """
    sols = solver.all_solutions()
    md = solver.solution_metadata()
    return [
        complex(sols[i][len(sols[i]) - 1]).real
        for i in range(len(sols))
        if int(md[i].endgame_success_code) == OK and md[i].is_finite and len(sols[i]) > 0
    ]


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--size', type=int, default=24, help='matrix size n (default 24)')
    parser.add_argument('--seed', type=int, default=20260614,
                        help='RNG seed for the matrix AND the homotopy (default fixed, so runs are '
                             'reproducible and serial/distributed are directly comparable)')
    args = parser.parse_args()

    comm = MPI.COMM_WORLD
    pb.random.set_random_seed(args.seed)        # reproducible homotopy (same on every rank)
    A = random_symmetric_integer_matrix(args.size, args.seed)
    solver = ZeroDimSolver(eigen_system(A), endgame='cauchy', mptype='adaptive', startsystem='mhom')

    start = time.time()
    if comm.Get_size() > 1:
        solver.solve(communicator=comm)         # manager (rank 0) hands paths to the workers
    else:
        solver.solve()                          # serial: a lone manager would have no workers
    elapsed = time.time() - start

    if not pb.parallel.is_manager():
        return

    got = recovered_eigenvalues(solver)
    expected = sorted(np.linalg.eigvals(A).real)

    # Distinct recovered eigenvalues: a path can occasionally duplicate one, so dedupe (the
    # eigenvalues are real and well separated) before counting -- the meaningful quantity is the
    # number of distinct eigenvalues, which is deterministic.
    distinct = []
    for g in sorted(got):
        if not distinct or abs(g - distinct[-1]) >= 1e-6:
            distinct.append(g)

    print('eigen-{}:  ranks={}  threads/rank={}  paths={}  recovered={}  wall={:.1f}s'.format(
        args.size, comm.Get_size(), os.environ.get('OMP_NUM_THREADS', '1'),
        args.size, len(distinct), elapsed))

    # Correctness: every numpy eigenvalue is matched by a homotopy-recovered one (closest-match,
    # robust to ordering and tiny imaginary parts), and we recover all n distinct eigenvalues.
    assert len(distinct) == args.size, \
        'expected {} distinct eigenvalues, recovered {}'.format(args.size, len(distinct))
    max_err = max(min(abs(g - e) for g in distinct) for e in expected)
    assert max_err < 1e-6, \
        'recovered eigenvalues disagree with numpy by {:.2e}'.format(max_err)
    print('  OK: all {} eigenvalues match numpy.linalg.eigvals (max error {:.2e})'.format(
        args.size, max_err))


if __name__ == '__main__':
    main()

Again the distributed solve is one line, and again it checks correctness – here against numpy.linalg.eigvals, demanding that every numpy eigenvalue is matched by a homotopy-recovered one. Because the matrix is symmetric with integer entries its spectrum is real and (generically) distinct, so a closest-match comparison is unambiguous.

Run the same ladder on a sizable matrix:

$ python python/examples/solve_eigenvalues.py --size 24              # serial baseline
$ mpirun -n 3 python python/examples/solve_eigenvalues.py --size 24  # 2 workers
$ mpirun -n 5 python python/examples/solve_eigenvalues.py --size 24  # 4 workers
$ mpirun -n 9 python python/examples/solve_eigenvalues.py --size 24  # 8 workers
$ mpirun -n 13 python python/examples/solve_eigenvalues.py --size 24 # 12 workers

eigenvalues of a 24x24 symmetric matrix (24 paths), measured on a 16-core Apple M3 Max (12 performance + 4 efficiency cores), 2026-06-30

launch	workers	wall-clock (s)	speedup
serial	–	12.6	1.0x
-n 3	2	7.2	1.8x
-n 5	4	4.5	2.8x
-n 9	8	3.3	3.8x
-n 13	12	3.4	3.7x

With only 24 paths the dynamic, demand-driven dispatch earns its keep: under adaptive precision the per-path cost varies a lot – a near-collision of eigenvalues makes one path linger at high precision – and a static split would leave most workers idle while one grinds. Demand-driven dispatch keeps the others busy, but it cannot split a single path: once enough workers are on hand, that one slowest eigenvalue sets the floor (here ~|tw-eigen-floor-wall| s, the 3.8x plateau). The plateau is unmistakable in the table – 12 workers is no faster than 8 (3.7x vs 3.8x): past the point where every path has a worker, handing out more workers changes nothing, because there is no 25th path to give them. That is the sharp contrast with cyclic-6 above, which has 720 paths and so keeps gaining out to 12 workers. Same lesson either way – distributing helps right up until you hit the longest path, which is exactly the bound the opening promised.

Threads within ranks: the hybrid model

Each worker rank can itself track several paths concurrently, one per thread, set with the OMP_NUM_THREADS environment variable (default 1). This stacks with MPI: R worker ranks of T threads each give R × T paths in flight.

The launch needs one extra flag, --bind-to none. By default mpirun pins each rank to a single core, which would strangle that rank’s threads onto one core; --bind-to none lets a rank’s threads spread across the machine:

OMP_NUM_THREADS=4 mpirun -n 3 --bind-to none python python/examples/solve_cyclic.py --n 6
#               4 threads  x  (3 ranks = 1 manager + 2 workers)  =  8 paths at a time

Why have two knobs for the same cores? Reach versus footprint. MPI ranks can live on different machines – that is the only way past one box’s core count – but each rank is a full process with its own copy of the system. Threads share one process’s memory and stay on one machine, so they add parallelism without the per-rank memory cost. The usual recipe: one rank per machine (or per NUMA node), threads to fill that machine’s cores.

At a fixed budget of 12 worker-cores, here is the same cyclic-6 solve split between ranks and threads:

cyclic-6 at 12 worker-cores, ranks x threads, measured on a 16-core Apple M3 Max (12 performance + 4 efficiency cores), 2026-06-30

layout	wall-clock (s)	speedup
-n 13, 12 workers x 1 thread	23.7	6.8x
-n 7, 6 workers x 2 threads	21.8	7.3x
-n 5, 4 workers x 3 threads	22.9	7.0x
-n 3, 2 workers x 6 threads	22.7	7.1x

On a single shared-memory machine the four layouts land within a few percent of each other – you are using the same twelve worker-cores either way, and they all bottom out at the same slowest-path floor. The reason to prefer more ranks shows up only across machines (more reach) or under memory pressure (threads share, ranks duplicate); the reason to prefer more threads is to keep the per-machine process count – and memory footprint – low.

Where to go from here

• The scripts take a --n / --size argument: push them up until a single machine is no longer enough, then spread the ranks across a cluster with a hostfile (mpirun --hostfile hosts -n ...).

• Combine with 🔁 Parameter homotopy: solve once, re-solve many times: solve once at a generic parameter, then distribute the re-solves across parameter values – each is an independent, already-cheap track, ideal for the manager-worker pool.

• The same solve(communicator=...) works for every zero-dim solver class (total degree, multihomogeneous, user homotopy), so anything you can solve serially, you can solve distributed.