Of monodromy, parallelization, and numerical hell

In a previous blog post, we described a software package that exploits the monodromy action for solving polynomial systems in a parametric family. We’re happy to report that the original paper accompanying this package is now published in the IMA Journal of Numerical Analysis. In this paper, we laid out the basic framework for monodromy solvers and highlight potential strengths of this approach, both by analyzing its complexity (conditional on strong randomness assumptions!) and testing our package against other solvers.

In this post, we will describe the sequel to this paper, presented at ISSAC 2018.

Many of the celebrated approaches to numerically solving polynomial systems employ a methodology in which a start system and a target system are both fixed. Numerical solutions to the target system are realized by numerical continuation along a fixed number of solution paths. This paradigm is ripe for speeeding up via parallel computing.


In the cartoon above, H(t) denotes a homotopy connecting the solutions between two fictitious polynomial systems. For two solution paths to be tracked, zero communication is needed from the two path trackers. Thus, four parallel processors result in a perfect speedup. This phenomenon is known as embarrassing parallelism.

With monodromy solvers, one structural property of systems we wish to exploit is lack of sharpness in the number of paths needed to find all solutions. For instance, the homotopy M(t) in the image above has two diverging solution paths. In some sense, an optimal homotopy for the system¬†M(1) would track only two paths. The goal of our monodromy approach is essentially to produce an optimal start system for a given family of systems. However, such a start system may not be easily accessible by deterministic calculation. Our approach is fundamentally randomized—for a family with generically d solutions, we allow for

O \big(|G| \cdot d \big)

path-tracking steps, where |G| measures the size of an underlying data structure we call the homotopy graph. Our desired output is a generic system with d solutions. In our parlance, we may visualize the homotopy graph connecting some systems in the parameter space and a solution graph of size |G| \cdot d lying above it.


In reality, our goal is to track as few paths along the homotopy graph as possible. With this goal in mind, we may observe that devising a scheduling algorithm for a parallel incarnation of monodromy is not a trivial task. One of our main contributions is a description of such an algorithm and numerous simulations indicating its potential to provide substantial speedups. If you’re curious, check out the paper or the accompanying github repository.

Another factor complicating the aforementioned scheduling algorithm: in reality, path-tracking can fail. In theory, paths in the parameter space which are suitably randomized will avoid the locus of ill-conditioned systems with probability 1. In reality, the core numerics may be unable (or simply unwilling) to distinguish nearly-singular systems from those that truly are. We may imagine, in the worlds inhabited by our graphs, several scattered regions of numerical hell, as depicted below.


In the ISSAC paper, we consider a simplified model of path-tracking failures, and describe its implications for the scheduling algorithm as well as robustness of our approach to a few failures—namely, if our graph is a clique as in two images above, we may tolerate a certain number of failures and still expect to produce the correct number of solutions, under suitable randomness assumptions.

One final word on randomness assumptions—despite being demonstrably false, they turn out to be useful in practice. For a different, but related, perspective, check out this nice lecture on sampling in numerical algebraic geometry.

Null Symbols and CRN Parsing

Having modified the parser to fix issues addressed in the previous post, I thought it might help to make explicit the precise types of inputs accepted — in other words, the language of all valid input strings for the reactionNetwork constructor. This is a regular language, but somewhat easier to describe in Extended Backus-Naur Normal form, a notational device for context-free languages. Basically, this allows you describe your language in terms of primitive strings called “terminals” and valid production rules for making larger strings.


For simplicity, we’ll stipulate that a reaction network in string format should contain a single reaction; additionally, we allow for the option that this reaction is separated by a comma from another reaction network. Symbolically, this is encoded by the production rule below.

reactionNetwork = reaction ,
{ "," reactionNetwork };

Macaulay2 isn’t going to recognize ‘reaction’ as a nonterminal, so we have to go further in specifying the grammar. Here it is.

# important definitions

reaction = complex , delimiter , complex ;

complex = summand , { "+" , complex } ;

species = letter , { letter } , [ "_" ] ,
[ number ] , [ "'" ] ;

summand = nullSymbol | [coefficient] , species ;

coefficient =  digit , { coefficient } ;

# more or less obvious definitions

posDigit = "1" | "2" | "3" | "4" | "5" | "6"
| "7" | "8" | "9" ;

digit = "0" | posDigit ;

posNum = posDigit , { digit } ;

numbers = posNum | "0" ;

capLetter = "a" | "b" | "c" | "d" | "e" | "f" |
"g" | "h" | "i" | "j" | "k" | "l" | "m" | "n" |
"o" | "p" | "q" | "r" | "s" | "t" | "u" | "v" |
"w" | "x" | "y" | "z" ;

lowLetter = "A" | "B" | "C" | "D" | "E" | "F" |
"G" | "H" | "I" | "J" | "K" | "L" | "M" | "N" |
"O" | "P" | "Q" | "R" | "S" | "T" | "U" | "V" |
"W" | "X" | "Y" | "Z" ;

letter = capLetter | lowLetter ;

delimiter = "-->" | "<--" | "<-->" ;

nullSymbol = "" ;

This seems to roughly correspond to the format used by the Control program. The main difference seems to be that Control accepts a variety of strings for the null Symbol (eg. “”, “0”.) Currently, our package only allows one null symbol at a time, but it may be user specified.


In the package, objects of type ReactionNetwork now have attributes “NullSymbol” and “NullIndex”, with respective defaults of “” and -1. The function ‘laplacian’ is a prototype for constructing steady state equations. It could use some improvements, but at least generates results consistent with previous worked examples. This extension also introduced a weird error for the ‘merge’ methods in which the default method is assumed, so I renamed the method as “glue.” I’m not sure what’s wrong here.