Function Repository Resource:

NSolveByMonodromy

Source Notebook

Numerically solve for isolated roots of a square system of polynomial equations using monodromy

Contributed by: Aravind Baskar

ResourceFunction["NSolveByMonodromy"][eqn,var]

finds isolated roots of the polynomial equation eqn in the variable var.

ResourceFunction["NSolveByMonodromy"][{eqn₁,...},{var₁,...}]

finds isolated roots of the system of polynomial equations eqns in the variables vars.

Details and Options

ResourceFunction["NSolveByMonodromy"] attempts to numerically approximate many isolated roots of a square system of polynomial equations through an iterative numerical continuation algorithm using monodromy.

ResourceFunction["NSolveByMonodromy"] deals primarily with polynomial equations but can handle some analytic transcendental equations with limited success.

The system eqns can be of the following form:

lhs==rhs

equations

expr

functions implicitly set to zero

A single variable or a list of variables matching the dimension of eqns must be specified.

ResourceFunction["NSolveByMonodromy"] finds an initial nonsingular root using FindRoot via a random guess to kick-start the first monodromy iteration. Alternatively, an user-defined initial point may be provided through the option "InitialPoint".

ResourceFunction["NSolveByMonodromy"] takes the following options:

"MaxRoots"

Infinity

Maximum number of roots to return

"Radius"

100

Radius parameter for the construction of monodromy loops

"InitialPoint"

{}

User-defined starting guess for the first loop

"MaxIterations"

Maximum number of iterations of monodromy loops

"MaxSteps"

1000

Maximum number of steps allowed in a single path

"MinStepSize"

10^-16

Minimum size of a single step used during path tracking

"Parallel"

False

Allow parallel computation of monodromy paths

Examples

Basic Examples (2)

Numerically solve a univariate polynomial equation in x:

In[1]:=

Out[2]=

Numerically solve a multivariate system of polynomial equations:

In[3]:=

$eqns = {45*P + 35*S - 165*B - 36, 35*P + 40*Z + 25*T - 27*S, 15*W + 25*P*S + 30*Z - 18*T - 165*B^2, -9*W + 15*P*T + 20*Z*S, W*P + 2*Z*T - 11*B^3, 99*W - 11*S*B + 3*B^2}; vars = {P, S, B, Z, T, W}; solns = ResourceFunction["NSolveByMonodromy"][eqns, vars]; solns // Length$

Out[6]=

Scope (7)

Create a real-valued polynomial objective function f in three variables x, y, z using random coefficient parameters:

In[7]:=

$f = Join[Array[g, 54], {1}] . {1, z, z^2, z^3, z^4, y, y z, y z^2, y z^3, y z^4, y z^5, y z^6, y^2, y^2 z, y^2 z^2, y^2 z^3, y^2 z^4, y^2 z^5, y^2 z^6, y^2 z^7, y^2 z^8, x, x z, x z^2, x z^3, x z^4, x y, x y z, x y z^2, x y z^3, x y z^4, x y z^5, x y z^6, x^2, x^2 z, x^2 z^2, x^2 z^3, x^2 z^4, x^2 y, x^2 y z, x^2 y z^2, x^2 y z^3, x^2 y z^4, x^2 y z^5, x^2 y z^6, x^3, x^3 z, x^3 z^2, x^3 z^3, x^3 z^4, x^4, x^4 z, x^4 z^2, x^4 z^3, x^4 z^4}; SeedRandom[1111]; subs = Thread[ Array[g, 54] -> RandomReal[{-1, 1}, 54, WorkingPrecision -> 2 $MachinePrecision]];$

Find stationary points of the objective f:

In[8]:=

vars = {x, y, z};
eqns = D[f, {vars}];
solns = ResourceFunction["NSolveByMonodromy"][(eqns /. subs), vars];
solns // Length

Out[11]=

Select the real valued ones:

In[12]:=

Out[13]=

Function values at these critical points:

In[14]:=

Out[14]=

Sorted critical points according to the function value:

In[15]:=

Out[16]=

Two of these (namely, #1 and #3) are local minima with all positive eigenvalues of the local Hessian, while the rest are saddle points:

In[17]:=

Out[18]=

Local minima and saddle points can be also discerned via the following ContourPlot3D visualization:

In[19]:=

styles = Table[Directive[ColorData[97][i], Opacity[0.35], Specularity[White, 100]], {i, Length[realsSort]}];
Show[ContourPlot3D[
f /. subs, {x, -1.75, 1.75}, {y, -1.25, 2.25}, {z, -1.75, 1.75}, Contours -> (f /. subs /. realsSort), Mesh -> None, ContourStyle -> styles, Lighting -> "Neutral", PlotPoints -> 40, PlotLegends -> Table[i, {i, Length[realsSort]}]], ListPointPlot3D[{#} & /@ Values[realsSort], PlotStyle -> PointSize[0.025]], ViewPoint -> {-1.125, 3.175, 0.350}]

Out[20]=

Options (3)

Default execution may not find all isolated roots:

In[21]:=

$eqns = {x1^2 - x1 + x2 + x3 + x4 + x5 - 10, x2^2 + x1 - x2 + x3 + x4 + x5 - 10, x3^2 + x1 + x2 - x3 + x4 + x5 - 10, x4^2 + x1 + x2 + x3 - x4 + x5 - 10, x5^2 + x1 + x2 + x3 + x4 - x5 - 10}; vars = {x1, x2, x3, x4, x5}; SeedRandom[4567]; solns = ResourceFunction["NSolveByMonodromy"][eqns, vars]; solns // Length$

Out[22]=

All the roots may be obtained by increasing or decreasing "Radius" parameter:

In[23]:=

SeedRandom[4567];
solns = ResourceFunction["NSolveByMonodromy"][eqns, vars, "Radius" -> 10];
solns // Length

Out[25]=

"MaxRoots" can be used to compute a subset of the total roots for unwieldy polynomial systems:

In[26]:=

$rpoly[m_(*Coefficient range*), n_(*Max. degree of each variable*), vars_(*Variable list*)] := Module[{dim = Length[vars], poly}, Table[ poly = RandomReal[{-m, m}, Table[n + 1, {i, dim}]]; Do[poly = poly . (vars[[j]]^Range[0, n]);, {j, dim}]; Expand[poly], {k, dim}] ]$

Generate a random large polynomial system of dimension 5 and total degree 100,000:

In[27]:=

SeedRandom[1234];
vars = {v, w, x, y, z};
eqns = rpoly[3, 2, vars];

Compute 50 roots:

In[28]:=

Out[8]=

In some cases transcendental functions can be minimized using "MaxRoots" termination criterion:

In[29]:=

$f = (Sin[10 \[Pi] x ]/(2 E^x)) + (x - 1)^4; var = x; eqn = D[f, var]; solns = ResourceFunction["NSolveByMonodromy"][eqn, var, "MaxRoots" -> 50];$

Visualize the real valued stationary points:

In[30]:=

$reals = Select[solns, Max[Abs[Im[Values[#]]]] <= 10^-16 &]; Show[{Plot[f, {x, -0.5, 2.5}, PlotStyle -> Blue], ListPlot[{x, f} /. reals, PlotStyle -> Red]}]$

Out[31]=

Applications (3)

Direct kinematics of a Gough-Stewart parallel 6-degree of freedom platform:

In[32]:=

$eqns = {t1^2 + t2^2 + t3^2 - 12 t1 - 68, t4^2 + t5^2 + t6^2 - 12 t5 - 68, t7^2 + t8^2 + t9^2 - 24 t8 - 12 t9 + 100, t1 t4 + t2 t5 + t3 t6 - 6 t1 - 6 t5 - 52, t1 t7 + t2 t8 + t3 t9 - 6 t1 - 12 t8 - 6 t9 + 64, t4 t7 + t5 t8 + t6 t9 - 6 t5 - 12 t8 - 6 t9 + 32, 2 t2 + 2 t3 - t4 - t5 - 2 t6 - t7 - t9 + 18, t1 + t2 + 2 t3 + 2 t4 + 2 t6 - 2 t7 + t8 - t9 - 38, t1 + t3 - 2 t4 + t5 - t6 + 2 t7 - 2 t8 + 8}; vars = {t1, t2, t3, t4, t5, t6, t7, t8, t9}; SeedRandom[3456]; solns = ResourceFunction["NSolveByMonodromy"][eqns, vars]; solns // Length$

Out[36]=

Katsura-5 system from magnetism in physics:

In[37]:=

$eqns = {x1^2 - x1 + 2*x2^2 + 2*x3^2 + 2*x4^2 + 2*x5^2 + 2*x6^2, 2*x1*x2 + 2*x2*x3 + 2*x3*x4 + 2*x4*x5 + 2*x5*x6 - x2, 2*x1*x3 + x2^2 + 2*x2*x4 + 2*x3*x5 + 2*x4*x6 - x3, 2*x1*x4 + 2*x2*x3 + 2*x2*x5 + 2*x3*x6 - x4, 2*x1*x5 + 2*x2*x4 + 2*x2*x6 + x3^2 - x5, x1 + 2*x2 + 2*x3 + 2*x4 + 2*x5 + 2*x6 - 1}; vars = {x1, x2, x3, x4, x5, x6}; SeedRandom[1357]; solns = ResourceFunction["NSolveByMonodromy"][eqns, vars]; solns // Length$

Out[41]=

A stationary chemical kinetics problem:

In[42]:=

$eqns = {9*y2^2 - 5.656854249492381*y2 + z2, x3^2 + y3^2 + z3^2 - 1, x4^2 + y4^2 + z4^2 - 1, y5^2 + z5^2 - 0.888888888888889, x3 - 2.828427124746190*y2*x3 + y2*y3 + z2*z3 - 1/3, x3*x4 + y3*y4 + z3*z4 - 1/3, 1/3*x4 + y4*y5 + z4*z5 - 1/3, 8/3 - 2.828427124746190*y2 + x3 + x4, y2 + y3 + y4 + y5 + 0.8888888888888889, z2 + z3 + z4 + z5}; vars = {x3, x4, y2, y3, y4, y5, z2, z3, z4, z5}; SeedRandom[4444]; solns = ResourceFunction["NSolveByMonodromy"][eqns, vars]; solns // Length$

Out[46]=

Properties and Relations (5)

NSolveByMonodromy can find many or all of the isolated roots with or without an initial point, e.g., a bifurcation problem:

In[47]:=

$eqns = {5 z1^9 - 6 z1^5 z2 + z1 z2^4 + 2 z1 z3, -2 z1^6 z2 + 2 z1^2 z2^3 + 2 z2 z3, z1^2 + z2^2 - 0.265625}; vars = {z1, z2, z3}; SeedRandom[6789]; solns = ResourceFunction["NSolveByMonodromy"][eqns, vars, "InitialPoint" -> {0.5, 0.1, 0.01}]; solns // Length$

Out[51]=

FindRoot attempts to find one root at a time based on an initial point:

In[52]:=

Out[52]=

NSolveByMonodromy compares favorably against NSolve in some large systems, e.g., 6R inverse position kinematics:

In[53]:=

$eqns = {z21^2 + z22^2 - 1.25830472585 + 1.05384271933*I, z31^2 + z32^2 + z33^2 - 1, z41^2 + z42^2 + z43^2 - 1, z51^2 + z52^2 + z53^2 - 1, z21*z31 + z22*z32 + (0.642935654806 + 0.819555356316*I)*z33 - 0.266880023988 - 0.452565255666*I, z31*z41 + z32*z42 + z33*z43 - 0.26425551745 - 0.342483846503*I, z41*z51 + z42*z52 + z43*z53 - 0.126010863922 - 0.864590917688*I, (0.352598136811 + 0.116888144319*I)* z51 + (0.539042485525 + 0.687058436892*I)* z52 + (0.391154215376 + 0.128900893182*I)*z53 - 0.179560356712 - 0.8709166566*I, (0.984138451804 + 0.414967172346*I)* z21 + (0.958341609741 + 0.847442419999*I)*z22* z33 + (-0.496254299764 - 0.546020011741*I)* z22 + (0.353268870498 + 0.389909226888*I)* z31 + (0.964759562277 + 0.71397074519*I)*z32* z43 + (0.0783739840935 - 1.33026494666*I)* z32 + (-0.964759562277 - 0.71397074519*I)*z33* z42 + (0.204379350351 + 0.00374294529684*I)* z41 + (0.706319991205 + 0.120097702053*I)*z42* z53 + (-0.706319991205 - 0.120097702053*I)*z43* z52 + (0.907681632783 + 0.405209293447*I)* z51 + (0.148939301127 + 0.182393186752*I)* z52 + (-0.0486385369088 - 0.496934768083*I)*z53 - 0.437312713588 - 0.914780691357*I, (-0.958341609741 - 0.847442419999*I)*z21* z33 + (0.496254299764 + 0.546020011741*I)* z21 + (0.984138451804 + 0.414967172346*I)* z22 + (-0.964759562277 - 0.71397074519*I)*z31* z43 + (-0.0783739840935 + 1.33026494666*I)* z31 + (0.353268870498 + 0.389909226888*I)* z32 + (0.964759562277 + 0.71397074519*I)*z33* z41 + (-0.706319991205 - 0.120097702053*I)*z41* z53 + (0.204379350351 + 0.00374294529684*I)* z42 + (0.706319991205 + 0.120097702053*I)*z43* z51 + (-0.148939301127 - 0.182393186752*I)* z51 + (0.907681632783 + 0.405209293447*I)* z52 + (0.134045297503 + 0.164748774862*I)*z53 - 0.719086796333 - 0.691791591267*I, (0.958341609741 + 0.847442419999*I)*z21* z32 + (-0.958341609741 - 0.847442419999*I)*z22* z31 + (0.964759562277 + 0.71397074519*I)*z31* z42 + (-0.964759562277 - 0.71397074519*I)*z32* z41 + (0.353268870498 + 0.389909226888*I)* z33 + (0.706319991205 + 0.120097702053*I)*z41* z52 + (-0.706319991205 - 0.120097702053*I)*z42* z51 + (0.204379350351 + 0.00374294529684*I)* z43 + (0.0486385369088 + 0.496934768083*I)* z51 + (-0.134045297503 - 0.164748774862*I)* z52 + (0.907681632783 + 0.405209293447*I)*z53 - 0.64863071126 + 0.983034576618*I}; vars = {z21, z22, z31, z32, z33, z41, z42, z43, z51, z52, z53};$

In[54]:=

SeedRandom[4444];
{time, solns} = AbsoluteTiming[ResourceFunction["NSolveByMonodromy"][eqns, vars]];
{time, solns // Length}

Out[17]=

NSolve default execution is 4 times slower and returns duplicate instances requiring an additional filtration step to remove them:

In[55]:=

Out[56]=

In[57]:=

Out[57]=

NSolve execution in high-precision takes even longer than the default execution:

In[58]:=

Out[58]=

Another example, a reduced 8-dimensional system arising in economics:

In[59]:=

eqns = {x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 + x6*x7 - 1*u8,
x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6 + x5*x7 - 2*u8,
x3 + x1*x4 + x2*x5 + x3*x6 + x4*x7 - 3*u8,
x4 + x1*x5 + x2*x6 + x3*x7 - 4*u8,
x5 + x1*x6 + x2*x7 - 5*u8,
x6 + x1*x7 - 6*u8,
x7 - 7*u8,
x1 + x2 + x3 + x4 + x5 + x6 + x7 + 1};
vars = {x1, x2, x3, x4, x5, x6, x7, u8};
SeedRandom[7890];
{time, solns} = AbsoluteTiming[ResourceFunction["NSolveByMonodromy"][eqns, vars]];
{time, solns // Length}

Out[63]=

NSolve default execution misses two isolated roots and returns some duplicate instances instead:

In[64]:=

Out[65]=

In[66]:=

Out[66]=

NSolve in high-precision returns the correct set of roots in a comparable time with NSolveByMonodromy:

In[67]:=

Out[68]=

In[69]:=

Out[69]=

NSolveByMonodromy does not accept a search range as input, unlike NSolve:

In[70]:=

eqn = x^2 - Sin[x];
var = x;
SeedRandom[5555];
solns = ResourceFunction["NSolveByMonodromy"][eqn, var, "MaxRoots" -> 10];

Select the real valued ones:

In[71]:=

Out[71]=

NSolve solves this system if a search range is provided:

In[72]:=

Out[72]=

In[73]:=

Out[73]=

NSolveByMonodromy can only solve square systems:

In[74]:=

$ResourceFunction[ "NSolveByMonodromy"][{x^2 + 3 y - 2, y + 3 z + 10}, {x, y, z}]$

Out[74]=

NSolve computes some witness points for such under-determined systems:

In[75]:=

$NSolve[{x^2 + 3 y - 2, y + 3 z + 10}, {x, y, z}]$

Out[75]=

Possible Issues (2)

NSolveByMonodromy cannot solve non-analytic systems:

In[76]:=

Out[76]=

NSolve solves this system using inverse functions:

In[77]:=

Out[77]=

NSolveByMonodromy may fail to proceed beyond initial point in systems with special monomial structures and/or parameters:

In[78]:=

eqn = Expand[Product[(x - i), {i, 20}]];
var = x;
SeedRandom[4321];
ResourceFunction["NSolveByMonodromy"][eqn, var]

Out[81]=

NSolve finds the approximation of all the 20 roots for the Wilkinson's polynomial:

In[82]:=

Out[82]=

Neat Examples (6)

Static equilibria of a weight suspended by two massless cables, each of different given lengths, in a plane:

In[83]:=

$A0 = {0, 0}(*Fixed pivot*); B0 = {bx, by}(*Fixed pivot*); C0 = A0 + l1 {Cos[\[Theta]1], Sin[\[Theta]1]}(*Moving pivot*); D0 = B0 + l3 {Cos[\[Theta]3], Sin[\[Theta]3]}(*Moving pivot*); f1 = C0 + l2 {Cos[\[Theta]2], Sin[\[Theta]2]} - D0(*Equation of vector loop-closure*); f2 = C0 + RotationMatrix[\[Theta]2] . {px, py} - {x, y}(*Equation of the locus of center of mass*); f3 = F1 {Cos[\[Theta]1], Sin[\[Theta]1]} + F3 {Cos[\[Theta]3], Sin[\[Theta]3]} + W {0, -1}(*Equation of force balance*); f41 = A0 + t1 {Cos[\[Theta]1], Sin[\[Theta]1]} - {x, y0};(*Equation of moment balance - part 1*) f42 = B0 + t3 {Cos[\[Theta]3], Sin[\[Theta]3]} - {x, y0};(*Equation of moment balance - part 2*) eqns = Numerator[ Together[ Flatten[{f1, f2, f3, f41, f42}] /. {Cos[\[Theta]1] -> (T1 + (1/T1))/2, Sin[\[Theta]1] -> (T1 - (1/T1))/(2 I), Cos[\[Theta]2] -> (T2 + (1/T2))/2, Sin[\[Theta]2] -> (T2 - (1/T2))/(2 I), Cos[\[Theta]3] -> (T3 + (1/T3))/2, Sin[\[Theta]3] -> (T3 - (1/ T3))/(2 I)}]];(*Euler parameterization from \[Theta] to T=\[ExponentialE]^(\[ImaginaryI] \[Theta])*) vars = {x, y, T1, T2, T3, F1, F3, t1, t3, y0}; subs = {W -> 1, bx -> 1, by -> -1/8, l1 -> 3/4, l2 -> 7/5, l3 -> 9/10, px -> 3/5, py -> -1/5 };$

This 10-dimensional system admits 12 nonsingular isolated roots in a generic case, as found by NSolveByMonodromy:

In[84]:=

SeedRandom[1234];
{time, solns} = AbsoluteTiming[
ResourceFunction["NSolveByMonodromy"][eqns /. subs, vars, "Radius" -> 0.1]];
{time, solns // Length}

Out[85]=

NSolve default execution takes a long time:

In[86]:=

Out[86]=

The presence of singular and positive-dimensional sets impacts the results of NSolve:

In[87]:=

Out[87]=

Select the real valued ones in x, y:

In[88]:=

$reals = Select[solns, Max[Abs[Im[{x, y} /. #]]] <= 10^-6 &] // Chop; reals // N // MatrixForm$

Out[89]=

Visualize the real equilibria, in which Configuration 1 is the only stable one with both cables under tension:

In[90]:=

$Table[ \[Theta]sub = {\[Theta]1 -> -I Log[T1], \[Theta]2 -> -I Log[ T2], \[Theta]3 -> -I Log[T3]} /. reals[[k]] // Chop; sc = 1/3(*Force scale*); Show[Graphics[ {Line[{A0, C0} /. subs /. \[Theta]sub], Line[{B0, D0} /. subs /. \[Theta]sub], Opacity[0.5], Triangle[{C0, D0, {x, y}} /. subs /. \[Theta]sub /. reals[[k]]], Opacity[1], Blue, Arrow[{{x, y}, {x, y} + sc W {0, -1}} /. subs /. reals[[k]]], Red, Arrow[{C0, C0 + sc F1 {Cos[\[Theta]1], Sin[\[Theta]1]}} /. subs /. \[Theta]sub /. reals[[k]]], Arrow[{D0, D0 + sc F3 {Cos[\[Theta]3], Sin[\[Theta]3]}} /. subs /. \[Theta]sub /. reals[[k]]], Black, Disk[A0, 0.025], Disk[B0 /. subs, 0.025]} ], PlotLabel -> Style[Framed["Config. " <> ToString[k]], Black, Background -> Lighter[Yellow]], PlotRange -> {{-0.5, 1.75}, {-1.5, 1.25}}, AspectRatio -> 1, Frame -> True, GridLines -> {{0}, {0}}] , {k, Length[reals]}]$

Out[90]=

Publisher

Aravind Baskar

Requirements

Wolfram Language 13.0 (December 2021) or above

Version History

1.0.0 – 30 June 2023

Source Metadata

Citation:
- 1. Hauenstein, Jonathan D., and Samantha N. Sherman. "Using monodromy to statistically estimate the number of solutions." In 2nd IMA Conference on Mathematics of Robotics 2, pp. 37-46. Springer International Publishing, 2022.
- 2. Baskar, Aravind, and Mark Plecnik. "Synthesis of six-bar timed curve generators of Stephenson-type using random monodromy loops." Journal of Mechanisms and Robotics 13, no. 1 (2021).
- 3. Brake, Daniel A., Jonathan D. Hauenstein, and Alan C. Liddell Jr. "Decomposing solution sets of polynomial systems using derivatives." In International Congress on Mathematical Software, pp. 127-135. Cham: Springer International Publishing, 2016.
- 4. Timme, Sascha. "Mixed precision path tracking for polynomial homotopy continuation." Advances in Computational Mathematics 47, no. 5 (2021): 75.

Author Notes

1. Encircling singularities: Selecting an appropriate size for the monodromy loops to ensure the encirclement of singularities can be challenging, as assessing the topology of the solution space is not trivial. 2. Missing solution sets: Monodromy in a parameter space limits the roots to the numerically irreducible set the starting seed belongs to, which may result in the omission of other disconnected sets. 3. High curvature points: Randomly constructed monodromy loops can occasionally miss roots with relatively high curvature, where the local Hessian exhibits a large matrix norm. As a result, the algorithm might converge to a root count while missing a small set. 4. Computational efficiency: This iterative method may not be computationally advantageous compared to other algebraic and numerical methods for smaller problems with a simple structure. 5. Termination criterion: A second-order linear trace test is used for polynomial systems as a sufficiency test, using a tolerance setting, to determine whether all the isolated roots in a numerically irreducible set have been computed. However, it is important to note that in some systems, this test may produce a false negative and return a warning, even when all the roots have been accurately computed. Such failures can occur due to certain regularity assumptions not holding in the given context. In such cases, Lincoln-Petersen statistical estimates are used as a fall-back termination criterion.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License

Wolfram Function Repository

NSolveByMonodromy

Details and Options

Examples

Basic Examples (2)

Scope (7)

Options (3)

Applications (3)

Properties and Relations (5)

Possible Issues (2)

Neat Examples (6)

Publisher

Related Links

Requirements

Version History

Source Metadata

Author Notes

License Information

NSolveByMonodromy

Details and Options

Examples

Basic Examples (2)

Scope (7)

Options (3)

Applications (3)

Properties and Relations (5)

Possible Issues (2)

Neat Examples (6)

Publisher

Related Links

Requirements

Version History

Source Metadata

Related Symbols

Author Notes

License Information