Returns invalid start value when the start value is actually valid
Closed this issue · 8 comments
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."
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.
Hey Alex,
- You can quote code using triple quotes (
```
). Eg
```julia
a = 3
b = a^2
```
produces
a = 3
b = a^2
- What Paul said (I was too slow 😅)
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
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
This looks like a textbook cusp, nice! :)
Pretty cuspy, indeed.
I just merged #455 which should help to point out the problem faster for the next person who runs into this.