janverschelde / PHCpack

The primary source code repository for PHCpack, a software package to solve polynomial systems with homotopy continuation methods.
http://www.phcpack.org
GNU General Public License v3.0
59 stars 21 forks source link

Diagonal solver solution #66

Closed wilrop closed 1 year ago

wilrop commented 1 year ago

I am trying to solve a system of polynomial equations with a divide-and-conquer approach. Essentially, I want to split a system of polynomials into subsystems, compute their witness sets and then start intersecting witness sets until I arrive at the final intersection which represents the solution to the original system. We had a brief email exchange about this a year ago, where you showed me your paper (http://homepages.math.uic.edu/~jan/par_solver.pdf) which is what I now also aim to do (without the parallelization necessarily). I believe, however, that there is a problem with the witness set that is returned by the diagonal solver in phcpy, or I am misunderstanding something and should do something differently.

Consider the following piece of code. The system is actually generated by my algorithm, so these are the types of problems I need to be able to solve.

from phcpy.solver import solve
from phcpy.sets import embed
from phcpy.diagonal import diagonal_solver as diagsolve

polys1 = [
    '+ (3 * p1_0 * p2_0 + 1 * p1_0 * p2_1 + 9 * p1_0 * p2_2 + 9 * p1_1 * p2_0 + 5 * p1_1 * p2_1 + 3 * p1_1 * p2_2 + 5 * p1_2 * p2_0 + 2 * p1_2 * p2_1 + 1 * p1_2 * p2_2) - (5 * p1_0 * p2_0 + 8 * p1_0 * p2_1 + 1 * p1_0 * p2_2 + 5 * p1_1 * p2_0 + 8 * p1_1 * p2_1 + 7 * p1_1 * p2_2 + 2 * p1_2 * p2_0 + 8 * p1_2 * p2_1 + 9 * p1_2 * p2_2);',
    '+ p0_1 + p0_2 - 1;',
    '+ p0_0;']
polys2 = [
    '+ (3 * p0_0 * p2_0 + 5 * p0_0 * p2_1 + 2 * p0_0 * p2_2 + 6 * p0_1 * p2_0 + 8 * p0_1 * p2_1 + 7 * p0_1 * p2_2 + 4 * p0_2 * p2_0 + 7 * p0_2 * p2_1 + 3 * p0_2 * p2_2) - (2 * p0_0 * p2_0 + 1 * p0_0 * p2_1 + 8 * p0_0 * p2_2 + 9 * p0_1 * p2_0 + 7 * p0_1 * p2_1 + 4 * p0_1 * p2_2 + 8 * p0_2 * p2_0 + 2 * p0_2 * p2_1 + 4 * p0_2 * p2_2);',
    '+ p1_0 + p1_2 - 1;',
    '+ p1_1;']
polys3 = [
    '+ (7 * p0_0 * p1_0 + 7 * p0_0 * p1_1 + 5 * p0_0 * p1_2 + 6 * p0_1 * p1_0 + 9 * p0_1 * p1_1 + 8 * p0_1 * p1_2 + 3 * p0_2 * p1_0 + 8 * p0_2 * p1_1 + 8 * p0_2 * p1_2) - (3 * p0_0 * p1_0 + 8 * p0_0 * p1_1 + 1 * p0_0 * p1_2 + 4 * p0_1 * p1_0 + 2 * p0_1 * p1_1 + 4 * p0_1 * p1_2 + 9 * p0_2 * p1_0 + 6 * p0_2 * p1_1 + 3 * p0_2 * p1_2);',
    '+ p2_1 + p2_2 - 1;',
    '+ p2_0;']

precision = 'd'

combined = polys1 + polys2 + polys3
sols = solve(combined, verbose=False, precision=precision)
for sol in sols:
    print(sol)

print("-------------------")

nvar = 9

dim1 = nvar - len(polys1)
emb1 = embed(nvar, dim1, polys1, precision=precision)

dim2 = nvar - len(polys2)
emb2 = embed(nvar, dim2, polys2, precision=precision)

dim3 = nvar - len(polys3)
emb3 = embed(nvar, dim3, polys3, precision=precision)

emb1[0] = '+ (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) - (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) +' + emb1[0]
emb2[0] = '+ (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) - (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) +' + emb2[0]
emb3[0] = '+ (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) - (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) +' + emb3[0]

sols1 = solve(emb1, verbose=False, precision=precision)
sols2 = solve(emb2, verbose=False, precision=precision)
sols3 = solve(emb3, verbose=False, precision=precision)

print("Intersecting")
inter_polys, inter_sols = diagsolve(nvar, dim1, emb1, sols1, dim2, emb2, sols2, verbose=False, prc=precision)

inter_polys = inter_polys[:-3]  # The last three polynomials are just ';' so we can ignore them.
inter_dim = nvar - len(inter_polys)  # Compute the new dimension.
inter_polys = embed(nvar, inter_dim, inter_polys, precision=precision)  # Make an embedding.
inter_sols = reformat_solutions(inter_sols, inter_dim)  # Reformat the solutions.
inter_sols_correct = solve(inter_polys, verbose=False, precision=precision)  # Solve the system correctly
inter_polys[0] = '+ (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) - (p0_0 + p0_1 + p0_2 + p1_0 + p1_1 + p1_2 + p2_0 + p2_1 + p2_2) +' + inter_polys[0]

polys, sols = diagsolve(nvar, dim3, emb3, sols3, inter_dim, inter_polys, inter_sols, prc=precision, verbose=False)

for sol in sols:
    print(sol)

What this code does is:

  1. Compute the solutions to the overall system first, just for checking purposes.
  2. Solve the three subsystems.
  3. Intersect subsystem 1 and 2.
  4. Reformat the witness set that was output in step 3. Here I use the reformat_solutions function, which I made to add the zero slack variables back to the solutions. If this is wrong and creating the problem, that would also be nice to know. I will add the code for that to the bottom of this issue.
  5. Intersect subsystem 3 with the witness set from the previous intersection.

The output from this code is one of four things. Sometimes it raises the following error:

    NaN****** 
         ^
SyntaxError: invalid syntax

Otherwise, it returns either no solutions, one of the two correct solutions and very rarely both correct solutions. The fact that it does sometimes return (part of) the solution leads me to believe that there is a subtle bug involved here.

It is interesting to note that when not calling my reformat_solutions function, but solving the embedded system again from scratch, the final intersection does work. However, this defeats the purpose of using the diagonal solver in the first place and is therefore not what I want.

My implementation for reformat_solutions.

def reformat_solutions(sols, dim):
    from phcpy.solutions import endmultiplicity, diagnostics, coordinates, make_solution

    new_sols = []

    for sol in sols:
        names, values = coordinates(sol)
        err, rco, res = diagnostics(sol)
        t, m = endmultiplicity(sol)

        for slack_var_idx in range(dim):
            names.append(f'zz{slack_var_idx + 1}')
            values.append(complex(0., 0.))

        updated_sol = make_solution(names, values, err=err, rco=rco, res=res, tval=t, multiplicity=m)
        new_sols.append(updated_sol)
    return new_sols

I have been working on and off on my algorithm for a year now, and this problem is a huge roadblock that I just cannot seem to solve. If you need any additional information, I would be very happy to share it!