Detecting solutions with multiplicity with custom homotopies

[oskarhenriksson]

What is the most reliable way to deal with solutions with multiplicity when tracing custom homotopies created with the Homotopy constructor? Suppose, for instance, that we have a simple homotopy, and some solutions for the start system at t=1, like these:

@var x t
H = Homotopy([t*im*(x^2-1)+(1-t)*(x-1)^2],[x],t)
start_solutions = [[-1],[1]]

If we use the EndgameTracker constructor to track from t=1 to t=0 in the following way, we find the solution x=1 twice:

start_solutions = [[-1],[1]]
res = track.(EndgameTracker(H),start_solutions)
target_solutions = solution.(res)

julia> target_solutions

2-element Vector{Vector{ComplexF64}}:
 [1.0 - 2.42861286636753e-17im]
 [1.0 + 0.0im]

julia> multiplicity.(res)

2-element Vector{Int64}:
 1
 1

What is the best way to detect which approximate solutions correspond to the same true solution, and obtain the multiplicity of such repeated solutions?

[saschatimme]

Multiplicities are always a best effort currently in the software. I think here the multiplicities are wrong since one of the initial solutions is also a solution of the final system.
The question you are asking is in all generality a very hard question. The general keyword you are looking for is "deflation".
You can take a look at https://arxiv.org/abs/1501.05083 for example if you only care about isolated singular solutions.
Unfortunately I am not aware of a Julia implementation of deflation.

Note that the certification in HC can verify distinctness of solutions once they are reduced to isolated non-singular solution.

[oskarhenriksson]

Thanks! I'll check out the paper!

If I understand the source code right, it seems like the heuristic for checking multiplicity in HC.jl is to filter out the singular solutions (which is checked through the condition number and accuracy from the path tracking), and among them, check how many ones are approximately the same.

I hadn't realized this before, but one can easily enhance the path results with this information through the Result constructor:

julia> R = Result(res)

Result with 2 solutions
=======================
• 2 paths tracked
• 1 non-singular solution (1 real)
• 1 singular solution (1 real)
• random_seed: nothing
• multiplicity table of singular solutions:
╭───────┬───────┬────────┬────────────╮
│ mult. │ total │ # real │ # non-real │
├───────┼───────┼────────┼────────────┤
│   1   │   1   │   1    │     0      │
╰───────┴───────┴────────┴────────────╯

or simply write R = solve(H, start_solutions) directly.

What seems to go wrong in this case is that the condition number is estimated as quite low for one of the paths, which means that this solution isn't counted in the multiplicity calculation:

julia> display(res)

2-element Vector{PathResult}:
 PathResult:
 • return_code → :success
 • solution → ComplexF64[1.0 + 0.0im]
 • accuracy → 0.0
 • residual → 0.0
 • condition_jacobian → Inf
 • singular → true
 • multiplicity → 1
 • winding_number → 1
 • steps → 7 / 1
 • extended_precision → false
 • path_number → 1

 PathResult:
 • return_code → :success
 • solution → ComplexF64[0.9999999999999998 - 5.204170427930421e-18im]
 • accuracy → 2.0163e-16
 • residual → 4.9331e-32
 **• condition_jacobian → 3.6302**
 • steps → 8 / 0
 • extended_precision → false
 • path_number → 2

[PBrdng]

When I run this code, I get condition_jacobian = 2.0588e16 and Inf, respectively. The multiplicity is correctly stated as 2.

Did you update to the latest version? Two months ago or so, I fixed a bug that caused the software to underestimate condition numbers in some cases (see #570 ).

[oskarhenriksson]

This code was indeed run in an environment with version v2.9.3. When I run it with the latest version, I also get a single solution of multiplicity 2, and the following path results, so everything now behaves as expected! Thanks for the help!

2-element Vector{PathResult}:
 PathResult:
 • return_code → :success
 • solution → ComplexF64[1.0 - 2.42861286636753e-17im]
 • accuracy → 1.1365e-16
 • residual → 5.8982e-34
 • condition_jacobian → 2.0588e16
 • singular → true
 • multiplicity → 1
 • steps → 6 / 0
 • extended_precision → false

 PathResult:
 • return_code → :success
 • solution → ComplexF64[1.0 + 0.0im]
 • accuracy → 0.0
 • residual → 0.0
 • condition_jacobian → Inf
 • singular → true
 • multiplicity → 1
 • winding_number → 1
 • steps → 7 / 1
 • extended_precision → false