JuliaHomotopyContinuation / HomotopyContinuation.jl

A Julia package for solving systems of polynomials via homotopy continuation.
https://www.JuliaHomotopyContinuation.org
MIT License
184 stars 30 forks source link

Missing Solutions #574

Closed ElviraLupoian closed 4 months ago

ElviraLupoian commented 4 months ago

I have a system which is theoretically known to have 364 solutions, that I'm trying to approximate using Homotopy Continuation. However, in most case this fails, and we usually obtain around 350. This seems to be constant across multiple runs.

using HomotopyContinuation

@var x1, x2, x3, x4, x5, x6, x7, x8, x9, x10

e1 = x2*x6^2 - x7*x10^3
e2 = -x1^2 + 2*x2*x5*x6 - 3*x7*x9*x10^2
e3 = -2*x1 + 2*x2*x4*x6 + x2*x5^2 - 3*x7*x8*x10^2 - 3*x7*x9^2*x10
e4 = 2*x2*x3*x6 + 2*x2*x4*x5 - 6*x7*x8*x9*x10 - x7*x9^3 - 3*x7*x10^2 - 1
e5 = 2*x2*x3*x5 + x2*x4^2 + 2*x2*x6 - 3*x7*x8^2*x10 - 3*x7*x8*x9^2 - 6*x7*x9*x10
e6 = 2*x2*x3*x4 + 2*x2*x5 - 3*x7*x8^2*x9 - 6*x7*x8*x10 - 3*x7*x9^2
e7 = x2*x3^2 + 2*x2*x4 - x7*x8^3 - 6*x7*x8*x9 - 3*x7*x10
e8 = -x1^2 + 2*x2*x3 - 3*x7*x8^2 - 3*x7*x9
e9 = -2*x1 + x2 - 3*x7*x8
e10 = -x7 - 1

es = [e1, e2, e3, e4, e5, e6, e7, e8, e9, e10]

sys = System(es)

sols = solve(sys)

julia> sols = solve(sys)
Tracking 736 paths... 100%|█████████████████████████████| Time: 0:00:04
  # paths tracked:                  736
  # non-singular solutions (real):  348 (22)
  # singular endpoints (real):      324 (30)
  # total solutions (real):         672 (52)
Result with 569 solutions
=========================
• 736 paths tracked
• 348 non-singular solutions (22 real)
• 221 singular solutions (18 real)
• random_seed: 0xbd9f1d30
• start_system: :polyhedral
• multiplicity table of singular solutions:
╭───────┬───────┬────────┬────────────╮
│ mult. │ total │ # real │ # non-real │
├───────┼───────┼────────┼────────────┤
│   1   │  164  │   12   │    152     │
│   2   │  23   │   2    │     21     │
│   3   │  22   │   2    │     20     │
│   4   │  12   │   2    │     10     │
╰───────┴───────┴────────┴────────────╯
PBrdng commented 4 months ago

Hi, I can confirm your computation. Do you know if the 364 theoretical solutions are all regular? Because this system is still small enough that I trust the software to give the correct answer.

ElviraLupoian commented 4 months ago

yes, they should all be regular.

PBrdng commented 4 months ago

348 seems the be a constant output. I would assume that the true number is 348 then. Can you point me to a reference?

ElviraLupoian commented 4 months ago

This system should correspond to the 3-torsion points of the Jacobian variety attached to a curve (genus 3, hyper elliptic) and their inverse.

The system is adapted from here : (https://arxiv.org/pdf/2210.02225) where we were initially trying to work with a system of 728 solutions, and had a similar issue of missing solutions. I initially suspected that the degree was too big and that was causing the issue, but this modified system has the same issue.

Our problem essentially reduces to this: we have a degree 7 polynomial which we use to generate a system of 10 equations in 10 variables. We know that the system has exactly 364 solutions, which we want to approximate.

I should say, that there are examples where this works and we get exactly 364. (for instance when f = x^7 -1, which I've included below)

using HomotopyContinuation

@var x1, x2, x3, x4, x5, x6, x7, x8, x9

e1 = -x1^2 + x2*x6^2 + x9^3
e2 = -2*x1 + 2*x2*x5*x6 + 3*x8*x9^2
e3 = 2*x2*x4*x6 + x2*x5^2 + 3*x7*x9^2 + 3*x8^2*x9 - 1
e4 = 2*x2*x3*x6 + 2*x2*x4*x5 + 6*x7*x8*x9 + x8^3 + 3*x9^2
e5 = 2*x2*x3*x5 + x2*x4^2 + 2*x2*x6 + 3*x7^2*x9 + 3*x7*x8^2 + 6*x8*x9
e6 = 2*x2*x3*x4 + 2*x2*x5 + 3*x7^2*x8 + 6*x7*x9 + 3*x8^2
e7 = x2*x3^2 + 2*x2*x4 + x7^3 + 6*x7*x8 + 3*x9
e8 = -x1^2 + 2*x2*x3 + 3*x7^2 + 3*x8
e9 = -2*x1 + x2 + 3*x7

es = [e1, e2, e3, e4, e5, e6, e7, e8, e9]

sys = System(es)

sols = HomotopyContinuation.solve(sys)
AndrewGibbs commented 4 months ago

I am also a little interested in this. Are there any numerical parameters that could be adjusted to improve the likelihood of finding additional solutions? For example, in endgame_options or tracker_options?

PBrdng commented 4 months ago

That is a good point. The missing 16 solutions seem to be almost singular. One could play with endgame_options to see if one can increase them.

PBrdng commented 4 months ago

With the following parameters I was able to get 364:

sols = solve(sys; endgame_options = EndgameOptions(; sing_accuracy = 0.1, singular_min_accuracy = 1e-5))

Note that I had to drastically reduce sing_accuracy to make it work. It seems that some of your zeros are really close to being singular. I could not certify all of them, only 348.

ElviraLupoian commented 4 months ago

Thank you very much.
I will experiment with this on some of my other systems.