Inconsistent results for simple polynomial

[AayushSabharwal]

Consider the simple polynomial system:

julia> @var x
julia> F = System([x*x-4x + 4])
System of length 1
 1 variables: x

 4 - 4*x + x^2

Repeatedly solving this in a loop, we get very different results seemingly randomly:

julia> solve(F; threading=false).path_results
2-element Vector{PathResult}:
 PathResult:
 • return_code → :success
 • solution → ComplexF64[2.0000000102978506 + 6.501381957665607e-9im]
 • accuracy → 2.291e-9
 • residual → 1.339e-16
 • condition_jacobian → 1.0264e7
 • singular → true
 • multiplicity → 1
 • winding_number → 2
 • steps → 32 / 3
 • extended_precision → false
 • path_number → 1

 PathResult:
 • return_code → :terminated_accuracy_limit
 • solution → ComplexF64[1.9704407177083874 - 0.006666101338130582im]
 • accuracy → 0.0038685
 • residual → 0.00091819
 • condition_jacobian → NaN
 • steps → 70 / 200
 • extended_precision → true
 • path_number → 2

julia> solve(F; threading=false).path_results
2-element Vector{PathResult}:
 PathResult:
 • return_code → :success
 • solution → ComplexF64[2.0000000000016565 - 1.32418736690043e-10im]
 • accuracy → 1.651e-10
 • residual → 8.8818e-16
 • condition_jacobian → 7.5512e9
 • singular → true
 • multiplicity → 2
 • winding_number → 2
 • steps → 32 / 2
 • extended_precision → true
 • path_number → 1

 PathResult:
 • return_code → :success
 • solution → ComplexF64[1.999999999855291 + 1.311915897327494e-8im]
 • accuracy → 2.03e-9
 • residual → 3.7969e-18
 • condition_jacobian → 37217.0
 • singular → true
 • multiplicity → 2
 • winding_number → 2
 • steps → 38 / 5
 • extended_precision → true
 • path_number → 2

julia> solve(F; threading=false).path_results
2-element Vector{PathResult}:
 PathResult:
 • return_code → :terminated_step_size_too_small
 • solution → ComplexF64[1.999999989084243 + 3.8097523992359725e-8im]
 • accuracy → 2.1831e-17
 • residual → 2.7343e-23
 • condition_jacobian → NaN
 • steps → 46 / 299
 • extended_precision → true
 • path_number → 1

 PathResult:
 • return_code → :success
 • solution → ComplexF64[2.000000000132282 + 3.431008043173774e-10im]
 • accuracy → 6.8238e-11
 • residual → 8.8818e-16
 • condition_jacobian → 5.4389e9
 • singular → true
 • multiplicity → 1
 • winding_number → 2
 • steps → 28 / 3
 • extended_precision → true
 • path_number → 2

Is there a way to make this more consistent? Why do the results vary so wildly for a simple problem?

[PBrdng]

Notice that your polynomial factors as $(x-2)^2$. So you are trying to compute a singular zero. For this one uses special methods, called endgames, to detect them. In both cases, $x=2$ is detected as a singular zero, so that is consistent. And that is what the software is supposed to do.

However, the multiplicity is correctly determined in only one case. Computing multiplicities is a very complicated numerical problem. To make it more consistent we need better algorithms I think.

[PBrdng]

In fact, higher dimensional components of multiplicity $>2$ are currently giving me a headache in another part of the implementation. This is a highly nontrivial and unsolved problem.

[AayushSabharwal]

It does find at least one root, but the accuracy varies quite a bit. Sometimes checking if a root is approximately 2.0 with atol = 1e-10 works, other times it's only accurate up to 1e-6ish. Is there a way to configure this to solve to a given accuracy?

[PBrdng]

No, because that root is singular.

[AayushSabharwal]

I see, thanks!