Returns invalid start value when the start value is actually valid

[alexheaton2]

I have an example where I know a valid start solution, yet my parameter homotopy gives the return_code :terminated_invalid_startvalue. This seems wrong. What am I doing wrong, or could there be a better return_code message? I am 100% sure the start solution is valid, since the polynomials are constructed from the start solution, basically. Also, I used the evaluate(...) function and checked they give zeros for all equations. What am I doing wrong?

Here is the code...

using HomotopyContinuation
p0 = [0 0; 1 0; 1 1; 2 1; 3/2 1/2;
    3.0 1.0; 4 1; 5 1; 9/2 1/2; 5 0; 6 0];
E = [1 2; 2 3; 2 5; 3 4; 3 5; 4 5; 4 6; 5 9; 6 7; 7 8;
    7 9; 8 9; 8 10; 9 10; 10 11];
pinnedvertices = [1; 6; 11]
freevertices = [2;3;4;5; 7;8;9;10]

@var x[1:size(p0)[1], 1:size(p0)[2]] # x[1:11, 1:2] # create all variables
xvarz_moving_frame = [Variable(:x,i,k) for i in freevertices for k in 1:2 ];

# create random, real-valued, linear equation in the moving frame variables
# created so that it passes through the initial configuration p0
a = randn(1, length(xvarz_moving_frame)) # random coefficients
a = [0 1.0 0 1.0 0 1.0 0 1.0 0 1.0 0 1.0 0 1.0 0 1.0]

b0 = subs(a * xvarz_moving_frame, [Variable(:x,i,k) => p0[i,k] for i in freevertices for k in 1:2 ]...)[1]
b0 = to_number(b0) # Float64(b0)

# the parameters to move linear "slice" later in a parameter homotopy, just move the constant term "b"
bvarz = [Variable(:b)]

# the linear equation with parameters
L = (a*xvarz_moving_frame)[1] .- Variable(:b)

ε = 1e-3 #1e-10
p1 = p0 + ε*randn(size(p0))

b1 = subs(a * xvarz_moving_frame, [Variable(:x,i,k) => p1[i,k] for i in freevertices for k in 1:2]...)[1]
b1 = to_number(b1) # Float64(b1) throws an error

function edge_equation(i,j)
    eqn = sum([(x[i,k] - x[j,k])^2 for k in 1:2])
    eqn += -sum([(p0[i,k] - p0[j,k])^2 for k in 1:2])
end
fs = [edge_equation(E[m,1],E[m,2]) for m in 1:size(E)[1]]

# pin the vertices by substitution
gs = [subs(fij, [Variable(:x,i,k) => p0[i,k] for i in pinnedvertices for k in 1:2]...) for fij in fs];

G = System(vcat(gs,L); variables=xvarz_moving_frame, parameters=bvarz)

startsolutions0 = [p0[i,k] for i in freevertices for k in 1:2]
result1 = solve(G, [startsolutions0]; start_parameters=[b0], target_parameters=[b1])
path_results(result1)

1-element Array{PathResult,1}:
 PathResult:
 • return_code → :terminated_invalid_startvalue
 • solution → Complex{Float64}[1.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 2.0 + 0.0im, 1.0 + 0.0im, 1.5 + 0.0im, 0.5 + 0.0im, 4.0 + 0.0im, 1.0 + 0.0im, 5.0 + 0.0im, 1.0 + 0.0im, 4.5 + 0.0im, 0.5 + 0.0im, 5.0 + 0.0im, 0.0 + 0.0im]
 • accuracy → NaN
 • residual → 0.0
 • condition_jacobian → NaN
 • steps → 0 / 0
 • extended_precision → true
 • path_number → 1

In the docs at https://www.juliahomotopycontinuation.org/HomotopyContinuation.jl/stable/tracker/#HomotopyContinuation.is_invalid_startvalue-Tuple{TrackerResult} it says "Tracking terminated since the provided start value was invalid."

[PBrdng]

Hi Alex,

the problem in your example is that the solution you provide is singular:

julia> using LinearAlgebra
julia> J = differentiate(vcat(gs,L), xvarz_moving_frame)
julia> J₀ = evaluate(J, xvarz_moving_frame => startsolutions0, bvarz => [b0])
julia> minimum(svdvals(J₀))
0.00

The algorithm checks if a solution is a valid start solution by applying Newton's method. If the method converges, then the solution is labelled a valid start solution. In your case, it does not converge, because the jacobian is singular.

[saschatimme]

Hey Alex,

1. You can quote code using triple quotes (```). Eg

   ```julia
   a = 3
   b = a^2

produces

 a = 3
 b = a^2

2. What Paul said (I was too slow 😅)

[saschatimme]

But I think it would make sense to improve the return_code / error message to distinguish between a non-solution and a singular start solution

[alexheaton2]

Okay, I see. Thanks for the help. Is there a standard way to deal with this situation? I can't get rid of it by randomization, can I? Is there something I can do to my system of equations? Probably not... I guess I should just compute a witness set for the 1-dimensional component.

The solution set should be 1-dimensional nearby that point, and I want to see if it is a cusp point or not. I was thinking to sample the variety nearby at many points, and then project all the sampled points into some 2 or 3 dimensional subspace. If the singular point I started with was a cusp, then it should still look like a cusp after projection.

Any thoughts on dealing with any of these issues? Thanks! -alex

[alexheaton2]

Okay, I computed a witness set, repeatedly sampled, and projected the cusp from 16 dimensional space into a random 2-dimensional subspace. Here it is...

[saschatimme]

This looks like a textbook cusp, nice! :)

[PBrdng]

Pretty cuspy, indeed.

[saschatimme]

I just merged #455 which should help to point out the problem faster for the next person who runs into this.