Merge pull request #455 from JuliaHomotopyContinuation/detect-singula…

…r-jacobian

Detect singular jacobian for invalid start values

docs/src/model_kit.md

@@ -36,6 +36,9 @@ variables(::Expression)
 ## System
 ```@docs
 System
+evaluate(F::System, x, p = nothing)
+jacobian(F::System)
+jacobian(F::System, x, p = nothing)
 degrees(F::System)
 expressions(F::System)
 optimize(::System)

src/endgame_tracker.jl

@@ -81,6 +81,7 @@ The possible states an `EndgameTracker` can be in:
 * `EndgameTrackerCode.terminated_accuracy_limit`
 * `EndgameTrackerCode.terminated_ill_conditioned`
 * `EndgameTrackerCode.terminated_invalid_startvalue`
+* `EndgameTrackerCode.terminated_invalid_startvalue_singular_jacobian`
 * `EndgameTrackerCode.terminated_max_winding_number`
 * `EndgameTrackerCode.terminated_max_steps`
 * `EndgameTrackerCode.terminated_step_size_too_small`
@@ -95,6 +96,7 @@ using ..TrackerCode: TrackerCode
     at_zero
     terminated_accuracy_limit
     terminated_invalid_startvalue
+    terminated_invalid_startvalue_singular_jacobian
     terminated_ill_conditioned
     terminated_max_steps
     terminated_max_extended_steps
@@ -117,6 +119,8 @@ function Base.convert(::Type{codes}, code::TrackerCode.codes)
         return terminated_ill_conditioned
     elseif code == TrackerCode.terminated_invalid_startvalue
         return terminated_invalid_startvalue
+    elseif code == TrackerCode.terminated_invalid_startvalue_singular_jacobian
+        return terminated_invalid_startvalue_singular_jacobian
     elseif code == TrackerCode.terminated_step_size_too_small
         return terminated_step_size_too_small
     elseif code == TrackerCode.terminated_unknown

src/model_kit/symbolic.jl

@@ -943,6 +943,7 @@ struct System
     variables::Vector{Variable}
     parameters::Vector{Variable}
     variable_groups::Union{Nothing,Vector{Vector{Variable}}}
+    _jacobian::Ref{Union{Nothing,Matrix{Expression}}}
 
     function System(
         exprs::Vector{Expression},
@@ -955,7 +956,13 @@ struct System
             vars == reduce(vcat, variable_groups) ||
                 throw(ArgumentError("Variable groups and variables don't match."))
         end
-        new(exprs, vars, params, variable_groups)
+        new(
+            exprs,
+            vars,
+            params,
+            variable_groups,
+            Ref{Union{Nothing,Matrix{Expression}}}(nothing),
+        )
     end
 end
 
@@ -1105,13 +1112,83 @@ function Base.show(io::IO, F::System)
     end
 end
 
+"""
+    evaluate(F::System, x, p = nothing)
+
+Evaluates the system `F` at `(x, p)`.
+
+# Example
+
+```julia
+@var x y
+F = System([x^2 + 3y, (y*x+1)^3])
+evaluate(F, [2, 3])
+```
+```
+2-element Array{Int32,1}:
+  13
+ 343
+```
+"""
+
 evaluate(F::System, x::AbstractVector) = evaluate(F.expressions, F.variables => x)
 function evaluate(F::System, x::AbstractVector, p::AbstractVector)
     evaluate(F.expressions, F.variables => x, F.parameters => p)
 end
 (F::System)(x::AbstractVector, p::Nothing = nothing) = evaluate(F, x)
 (F::System)(x::AbstractVector, p::AbstractVector) = evaluate(F, x, p)
 
+"""
+    jacobian(F::System)
+
+Computes the symbolic Jacobian of the system `F`.
+
+# Example
+
+```julia
+@var x y
+F = System([x^2 + 3y, (y*x+1)^3])
+jacobian(F)
+```
+```
+2×2 Array{Expression,2}:
+             2*x                3
+ 3*y*(1 + x*y)^2  3*x*(1 + x*y)^2
+```
+"""
+function jacobian(F::System)
+    if isnothing(F._jacobian[])
+        F._jacobian[] = differentiate(F.expressions, F.variables)
+    end
+    F._jacobian[]::Matrix{Expression}
+end
+
+"""
+    jacobian(F::System, x, p = nothing)
+
+Evaluates the Jacobian of the system `F` at `(x, p)`.
+
+# Example
+
+```julia
+@var x y
+F = System([x^2 + 3y, (y*x+1)^3])
+jacobian(F, [2, 3])
+```
+```
+2×2 Array{Int32,2}:
+   4    3
+ 441  294
+```
+"""
+function jacobian(F::System, x, p = nothing)
+    if p isa Nothing
+        evaluate(jacobian(F), F.variables => x)
+    else
+        evaluate(jacobian(F), F.variables => x, F.parameters => p)
+    end
+end
+
 function Base.:(==)(F::System, G::System)
     F.expressions == G.expressions &&
         F.variables == G.variables &&

src/tracker.jl

@@ -170,6 +170,7 @@ module TrackerCode
     terminated_accuracy_limit
     terminated_ill_conditioned
     terminated_invalid_startvalue
+    terminated_invalid_startvalue_singular_jacobian
     terminated_step_size_too_small
     terminated_unknown
 end
@@ -194,9 +195,12 @@ is_terminated(S::TrackerCode.codes) = S ≠ TrackerCode.tracking && S ≠ Tracke
     is_invalid_startvalue(code::TrackerCode.codes)
 
 Returns `true` if `code` indicates that the path tracking got terminated since the start
-value was not a zero.
+value was not a regular zero.
 """
-is_invalid_startvalue(S::TrackerCode.codes) = S == TrackerCode.terminated_invalid_startvalue
+function is_invalid_startvalue(S::TrackerCode.codes)
+    S == TrackerCode.terminated_invalid_startvalue ||
+        S == TrackerCode.terminated_invalid_startvalue_singular_jacobian
+end
 
 """
     is_tracking(code::TrackerCode.codes)
@@ -256,8 +260,17 @@ is_success(R::TrackerResult) = R.return_code == :success
     is_invalid_startvalue(result::TrackerResult)
 
 Returns `true` if the path tracking failed since the start value was invalid.
+You can inspect `result.return_code` to get the exact return code. Possible values
+if `is_invalid_startvalue(result) == true` are
+* `:terminated_invalid_startvalue_singular_jacobian` indicates that the Jacobian of the homotopy at
+  the provided start value is singular, i.e., if it has not full-column rank.
+* `:terminated_invalid_startvalue` indicates that the the provided start value is not sufficiently
+  close to a solution of the homotopy.
 """
-is_invalid_startvalue(R::TrackerResult) = R.return_code == :invalid_startvalue
+function is_invalid_startvalue(R::TrackerResult)
+    R.return_code == :terminated_invalid_startvalue ||
+        R.return_code == :terminated_invalid_startvalue_singular_jacobian
+end
 
 """
     solution(result::TrackerResult)
@@ -696,7 +709,15 @@ function init!(
         state.accuracy = μ
         state.μ = max(μ, eps())
     else
-        state.code = TrackerCode.terminated_invalid_startvalue
+        evaluate_and_jacobian!(corrector.r, workspace(jacobian), homotopy, x, t)
+        J = matrix(workspace(jacobian))
+        corank = size(J, 2) - LA.rank(J, rtol = 1e-14)
+        if corank > 0
+            state.code = TrackerCode.terminated_invalid_startvalue_singular_jacobian
+        else
+            state.code = TrackerCode.terminated_invalid_startvalue
+        end
+
         return false
     end
 

test/tracker_test.jl

@@ -98,6 +98,8 @@
 
         @test is_success(track(tracker, S[1], 1, 0))
         @test is_success(track(tracker, S[2], 1, 0))
+
+        @test is_invalid_startvalue(track(tracker, [100, -100], 1, 0))
     end
 
     @testset "Compiled/InterpretedHomotopy" begin
@@ -130,6 +132,90 @@
         end
     end
 
+    @testset "invalid_startvalue_singular_jacobian" begin
+        # https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/issues/454
+        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 = 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 = evaluate(
+            a * xvarz_moving_frame,
+            [Variable(:x, i, k) => p0[i, k] for i in freevertices for k = 1:2]...,
+        )[1]
+        # 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 = 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 = 1:2])
+            eqn += -sum([(p0[i, k] - p0[j, k])^2 for k = 1:2])
+        end
+        fs = [edge_equation(E[m, 1], E[m, 2]) for m = 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 = 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 = 1:2]
+        tracker = Tracker(ParameterHomotopy(G, [b0], [b1]; compile = false))
+        result = track(tracker, startsolutions0, 1, 0)
+        @test is_invalid_startvalue(result)
+        @test result.return_code == :terminated_invalid_startvalue_singular_jacobian
+    end
+
     include("test_cases/steiner_higher_prec.jl")
     include("test_cases/four_bar.jl")
 end