The Blackbox Solver

A blackbox solver runs with default settings and values of the tolerances, executing a selection of algorithms which has shown to work for a large collection of problems. As to what the solver computes, two different types of solvers are available

  1. to approximate all isolated solutions, typically of systems with as many equations as unknows; and

  2. to compute a numerical irreducible decomposition of the entire solution set, which includes the isolated solution points, but also all positive dimensional sets (curves and surfaces), factored into irreducible components.

Because the output of the two types is so vastly different and because the complexity of a numerical irreducible decomposition is much higher than computing only the isolated solutions, the user must decide in advance which solver function to call.

approximating all isolated solutions

The input to the solver is a list of strings, with symbolic representations of the polynomials in the system. The ; at the end signifies the = 0 of the equations.

polynomials = ['x^3 + 2*x*y - x^2;', 'x + y - x^3;']

To call the blackbox solver, we import the solve function from the solver module.

from phcpy.solver import solve
solutions = solve(polynomials)

The output of solve is also a list of strings.

for (idx, sol) in enumerate(solutions):
    print('Solution', idx+1, ':')
    print(sol)

has as output

Solution 1 :
t :  0.00000000000000E+00   0.00000000000000E+00
m : 2
the solution for t :
 x :  0.00000000000000E+00   0.00000000000000E+00
 y :  0.00000000000000E+00   0.00000000000000E+00
== err :  7.124E-17 = rco :  0.000E+00 = res :  0.000E+00 =
Solution 2 :
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 x :  1.00000000000000E+00   0.00000000000000E+00
 y :  0.00000000000000E+00   0.00000000000000E+00
== err :  2.100E-74 = rco :  5.348E-01 = res :  0.000E+00 =
Solution 3 :
t :  1.00000000000000E+00   0.00000000000000E+00
m : 1
the solution for t :
 x : -1.50000000000000E+00   0.00000000000000E+00
 y : -1.87500000000000E+00   3.30053725596684E-238
== err :  6.747E-80 = rco :  1.220E-01 = res : 1.073E-237 =

How did this actually work? We can ask to see more output of the solver, giving a value to the verbose level parameter, the vrblvl as argument to the solver.

solutions = solve(polynomials, vrblvl=1)

shows the following

in solve, tasks : 0, mvfocus : 0, precision : d
the polynomials on input :
x^3 + 2*x*y - x^2;
x + y - x^3;
nbr : 4, roco :
total degree : 9
2-homogeneous Bezout number : 6
  with with partition : { x }{ y }
general linear-product Bezout number : 5
  based on the set structure :
     { x }{ x y }{ x }
     { x y }{ x }{ x }
mixed volume : 2
stable mixed volume : 4

The roco is an abbreviation of root count. A root count is an upper bound on the number of solutions. In the example, the two input polynomials are cubics. Therefore, the largest number of isolated solutions equals nine, the product of the degrees of the polynomials. However, not all monomials of degree three or less that could appear with nonzero coefficient are present. The numbers 6 and 5 are better bounds computed on the degrees. Mixed volumes are generically sharp bounds for the number of isolated solutions with all coordinates different from zero. As in the example, there is one double root at (0, 0), which is counted by the stable mixed volume.

If the coefficients of the polynomials are sufficiently independent from each other, then the number of isolated solutions counted with multiplicity will match one of the computed root counts.

options of the solve function

The blackbox solver runs with default values, which can be changed when calling the solve function.

In addition to the list of polynomials, there are six parameters with default values:

  1. tasks equals the number of the threads and is by default set to zero. The solutions can approximated independently to each other using p threads could in the best case speed up the solver by a factor of p.

  2. mvfocus is a flag to apply only mixed volumes as root counts, if set to one. By default, this flag is set to zero and the solver will compute all bounds based on the degrees which may be too time consuming for sparse polynomial systems.

  3. precision is by default double (d). Other values are dd for double double and qd for quad double. While running in higher precision leads to more accurate results, the computational overhead can be significant. The overhead may be compensated (in part) by multithreading.

  4. checkin is the option to run some basic checks on the input. By default True, setting it to False is okay if the input polynomials are automatically generated in the correct way.

  5. dictionary_output is by default False and a list of strings is returned. If True, then the output is a list of dictionaries, often convenient for processing.

  6. vrblvl is the verbose level parameter to make the blackbox solver less black when set to higher values.

Of course, then there is always the help(solve):

help(solve)

which shows the documentation string of the function.

When changing the default values, consider the following.

  1. The number of threads should never be set to a value higher than the number of available cores on the system. To find out the number of available cores, do the following:

    from phcpy.dimension import get_core_count
get_core_count()

    So, use up to the value returned by phcpy.dimension.get_core_count() as the value to assign to the parameter tasks.

  2. The focus on mixed volumes (in the option mvfocus) is automatically applied when the polynomials have negative exponents.

  3. When computing in higher precision, keep in mind that also the coefficients of the polynomials then must also be evaluated in higher precision. Consider 1/3 in double precision:

    1/3

    which evaluates to

    0.3333333333333333

    which is of course not equal to the rational number 1/3.

  4. One of the checks done by default (checkin=True) is whether the number of polynomials in the list equals the number of unknowns. At this stage, if the syntax of the polynomial is incorrect, an error message will be printed as well.

  5. If dictionary_output is wanted after a run, then it can be computed afterwards, the solve() should not be called again, but can be computed with the strsol2dict() of the solutions module. For example:

    from phcpy.solutions import strsol2dict
strsol2dict(solutions[0])

    shows the output

    {'t': 0j, 'm': 2, 'err': 1.17e-16, 'rco': 0.0, 'res': 0.0, 'x': 0j, 'y': 0j}

  6. Higher values of the verbose level vrblvl are mainly meant for debugging purposes as it should procedures are executed. As the solving of a polynomial system could take a long time, the user can see which procedures are currently running if the solver appears to be stuck.