Interesting possible failure of monodromy solving

[orebas]

Below is an MWE. I noticed that when solving the very simple equation x^3=8, monodromy solving doesn't seem to find any new solutions beyond the one it's fed. The default solver doesn't seem to have this issue. I will try using some different parameters to see if it helps.

Example output:

julia> include("HC-MWE.jl")

=== Regular solving ===
All solutions from regular solve:
3-element Vector{Vector{ComplexF64}}:
 [2.0 + 0.0im]
 [-1.0 + 1.7320508075688772im]
 [-1.0 - 1.7320508075688774im]

=== Monodromy solving ===
1-element Vector{Vector{ComplexF64}}:
 [-1.0 + 1.7320508075688772im]

Solutions from monodromy:
1-element Vector{Vector{ComplexF64}}:
 [-1.0 + 1.7320508075688772im]

Verification for monodromy solutions (should all be close to 0):
x = -1.0 + 1.7320508075688772im => residual = 1.9860273225978185e-15

julia> include("HC-MWE.jl")

=== Regular solving ===
All solutions from regular solve:
3-element Vector{Vector{ComplexF64}}:
 [2.0 + 0.0im]
 [-1.0 + 1.7320508075688774im]
 [-1.0 - 1.7320508075688772im]

=== Monodromy solving ===
1-element Vector{Vector{ComplexF64}}:
 [2.0 + 0.0im]

Solutions from monodromy:
1-element Vector{Vector{ComplexF64}}:
 [2.0 + 0.0im]

Verification for monodromy solutions (should all be close to 0):
x = 2.0 + 0.0im => residual = 0.0

MWE:

using HomotopyContinuation

# Create a simple system with the same structure as your problem
function MonodromyTest()
	@var x p

	# Define the system: x³ - 8 = 0 (when p = 0)
	F = System([x^3 - 8.0 - p], parameters = [p])

	println("\n=== Regular solving ===")
	# Solve with p = 0
	result = solve(F, target_parameters = [0.0])

	# Display all solutions
	println("All solutions from regular solve:")
	display(solutions(result))

	println("\n=== Monodromy solving ===")
	# Now try with monodromy
	# First find a start pair
	found_start_pair = false
	pair_attempts = 0
	newx = nothing
	param_final = [0.0]

	while (!found_start_pair && pair_attempts < 50)
		testx, testp = find_start_pair(F)
		newx = solve(F, testx, start_parameters = testp, target_parameters = param_final,
			tracker_options = TrackerOptions(automatic_differentiation = 3))
		startpsoln = solutions(newx)
		pair_attempts += 1
		if (!isempty(startpsoln))
			found_start_pair = true
		end
	end

	display(solutions(newx))

	# Now do monodromy solve
	result_mono = monodromy_solve(F, solutions(newx), param_final,
		show_progress = true,
		timeout = 120.0,
		max_loops_no_progress = 20,
		unique_points_rtol = 1e-6,
		unique_points_atol = 1e-6,
		tracker_options = TrackerOptions(automatic_differentiation = 3))

	println("\nSolutions from monodromy:")
	display(solutions(result_mono))

	# Verify these are actually solutions
	println("\nVerification for monodromy solutions (should all be close to 0):")
	for sol in solutions(result_mono)
		residual = abs(sol[1]^3 - 8)
		println("x = ", sol[1], " => residual = ", residual)
	end
end

MonodromyTest()

[orebas]

Based on a long chat with chatgpt, see https://chatgpt.com/share/67831bf4-7f2c-8007-895c-34ae1f2b74bd

I tried replacing the parametrization x^3-8-p=0 with the parametrization q*(x^3-8)-p=0, and indeed all 3 solutions are recovered. (code below.) I guess it would be nice if there were a way for the solver to "notice" when this happens programatically, or maybe even correct for it. I would love to hear other thoughts on this. (Acknowledging there may be nothing wrong in the software, more of a thing for users to watch out for)

using HomotopyContinuation

function MonodromyTest()
	# Create a simple system with two parameters
	@var x p q

	# Define the system: q*(x³ - 8) - p = 0
	F = System([q * (x^3 - 8.0) - p], parameters = [p, q])

	println("\n=== Regular solving ===")
	# Solve with p = 0, q = 1
	result = solve(F, target_parameters = [0.0, 1.0])

	# Display all solutions
	println("All solutions from regular solve:")
	display(solutions(result))

	println("\n=== Monodromy solving ===")
	# Now try with monodromy
	# First find a start pair
	found_start_pair = false
	pair_attempts = 0
	newx = nothing
	param_final = [0.0, 1.0]  # Target parameters: p = 0, q = 1

	while (!found_start_pair && pair_attempts < 50)
		testx, testp = find_start_pair(F)
		newx = solve(F, testx, start_parameters = testp, target_parameters = param_final,
			tracker_options = TrackerOptions(automatic_differentiation = 3))
		startpsoln = solutions(newx)
		pair_attempts += 1
		if (!isempty(startpsoln))
			found_start_pair = true
		end
	end

	# Now do monodromy solve with modified parameters
	result_mono = monodromy_solve(F, solutions(newx), param_final,
		show_progress = true,
		timeout = 300.0,
		max_loops_no_progress = 50,
		target_solutions_count = 3,
		unique_points_rtol = 1e-8,
		unique_points_atol = 1e-8,
		trace_test = true,
		trace_test_tol = 1e-10,
		min_solutions = 3,
		tracker_options = TrackerOptions(automatic_differentiation = 3))

	println("\nSolutions from monodromy:")
	display(solutions(result_mono))
	display(certify(F, result_mono))

	# Verify these are actually solutions
	println("\nVerification for monodromy solutions (should all be close to 0):")
	for sol in solutions(result_mono)
		residual = abs(sol[1]^3 - 8)
		println("x = ", sol[1], " => residual = ", residual)
	end
end

MonodromyTest()

[PBrdng]

I tried replacing the parametrization x^3-8-p=0 with the parametrization q*(x^3-8)-p=0, and indeed all 3 solutions are recovered. (code below.) I guess it would be nice if there were a way for the solver to "notice" when this happens programatically, or maybe even correct for it.

We do not understand the theory of monodromy enough to predict when the loops we generate are not enough to find all solutions. In your example, $x^3 - 8 - p$ is in fact enough to find all solutions if you can set up the appropriate loops. For instance, I generated by hand a loop in $p$ that encircles $p=8$ and it finds a new solutions. In general, the set of all loops that encircle the branch locus generate all solutions (ChatGPT is false here: the fact that the polynomial is linear in $p$ is irrelevant). But we do not know how to generate these loops given the parameterization.

[PBrdng]

For instance, this finds all 3 solutions (because the loops are now big enough):

S = monodromy_solve(F, parameter_sampler = p-> 10 .* randn(ComplexF64, length(p)))