JuliaHomotopyContinuation / HomotopyContinuation.jl

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

HC fails to solve sanity check #543

Open votroto opened 1 year ago

votroto commented 1 year ago

Hello, I do not understand why HC is unable to solve the system x = 0.

Solving it

using HomotopyContinuation

@var x
system = System([x])

res = solve(system)

results in 0 solutions. Changing it to x+eps()=0 fixes the issue. Is this something homotopy methods should not be expected to solve, or is this a bug? Should I randomly add tiny noise to my systems?

Thanks

PBrdng commented 1 year ago

The problem here is that your system of equations is homogeneous. The software detects this and then tries to find zeros in projective space. However, if there is only 1 variable this is not meaningful, and we should catch this case.

If you add another variable, it is well-defined and works:

@var x y
solve(System([x], variables = [x;y]))
saschatimme commented 1 year ago

I guess we can require for the homogenous check that the system has at least two variables

votroto commented 1 year ago

I'll give you the context, in case it helps your decision.

I'm using HC in a randomized iterative algorithm as a fallback method to find Nash equilibria. Edge cases are therefore likely. Presumeably, HC is not the algorithm of choice for homogenous systems, and it would be silly to expect a "good" solution. However, what is important to me is that HC be reliable enough to return something reasonable for well defined problems.

saschatimme commented 1 year ago

Presumeably, HC is not the algorithm of choice for homogenous systems

Actually, HC is very good at solving homogenous systems. If a system is homogenous, then the system has naturally at least a one-dimensional solution space (if a $y$ is a solution then also $\lambda y$ for all $\lambda \in \mathbb{C}$). Since we can only deal with isolated solutions we therefore opt for automating "fixing" the problem by computing solutions in an associated affine solution space.

However, what is important to me is that HC be reliable enough to return something reasonable for well defined problems.

We hope that HC is helpful to you and for your research but I like to remind you that this is free and open source software. If you find ways to improve the software you are most welcome to contribute.

votroto commented 1 year ago

Sure, I'll give it a shot if you want; though I doubt my ability to meaningfully contribute to the project as I have no experience with implementing homotopy methods.

votroto commented 1 year ago

BTW @saschatimme, what did you mean by "at least two variables?" This works res = solve(System([-x+y, -y+x])), but this res = solve(System([x, y])) does not.

I tried, but failed. It does not seem it's a straightforward fix. Checking for more than one variable does not solve the problem, and simply checking whether the system is over-constrained fails other tests. This was my quick test set:

@testset "Edge Cases" begin
    function _random_homogeneous_system(variable_count)
        C = rand(variable_count, variable_count) .+ 1 # avoid zeros
        x, = @var x[1:variable_count]

        C, x
    end

    @testset "HC finds trivial solution of tiny homogeneous system ($start, v=$(v))" for start in [:polyhedral, :total_degree], v in [1, 2]
        C, x = _random_homogeneous_system(v)

        sys = System(C * x)
        result = solve(sys; 
            show_progress = false,
            start_system = start,
            compile = false)

        expected = zeros(size(x))
        actual = real_solutions(result)

        is_expected = isapprox(expected; atol = 1e-8)
        @test !isempty(actual)
        @test any(is_expected, actual)
    end
end
votroto commented 1 year ago

While looking into the issue i found broken interpolation in utils.jl at line 152

Should be:

function Base.showerror(io::IO, E::KeywordArgumentException)
    print(io, "Invalid argument $(E.given) for $(E.key). ", E.msg)
end

I reckon it is pointless opening a PR just for that.