Documentation for taylor series and custom homotopies

[AayushSabharwal]

The documentation is sparse on both the requirements of taylor! and how to write custom homotopies.

For the former, it seems my implementation of taylor! for a custom AbstractSystem is called with x::HomotopyContinuation.ModelKit.TaylorVector which is an unfamiliar type. How am I meant to interact with this type? What is its public API? The documentation mentions:

If the system should be used in context of a parameter homotopy it is also necessary to implement

taylor!(u, ::Val{1}, F::AbstractSystem, x, p::TaylorVector{2})

Yet I am unable to find what p::TaylorVector{2} means. Is there something I'm missing?

As for the latter, what are the functions required to implement a custom AbstractHomotopy? I can dive into the code to find these, but it seems apt to have documentation as well.

[PBrdng]

You are right. The documentation for writing custom homotopies is not good. I should write a guide on the website. Can you share what homotopy you have in mind and the code you have so far? We can fix it together.

[AayushSabharwal]

I'm trying to write a SciML interface to HomotopyContinuation.jl, through SciML/NonlinearSolve.jl#523. The custom system and homotopy implementations are in interface_types.jl. The system is designed to be a wrapper over a user provided polynomial function f(u, p), where u is traceable via ModelKit but p isn't. The homotopy is H(x, t) = t * f(x) + (1 - t) * (f(x) - f(x_0)) where x_0 is the user-provided initial guess for x. x = x_0 is the only path this homotopy is designed to trace. It can be reformulated as H(x, t) = f(x) + (1 + t) * f(x_0)

[PBrdng]

For this specific case, you can simply implement a straight line homotopy. Or what am I missing?

[AayushSabharwal]

Using StraightLineHomotopy with the custom system will result in two f-calls since it'll use the H(x, t) = t * f(x) + (1 - t) * (f(x) - f(x_0)) formulation. It also still requires implementing taylor! which I'm having trouble with.

[AayushSabharwal]

The custom system is necessary since it's possible to do a lot more than HC.ModelKit allows and still have the system be polynomial.

[PBrdng]

I understand. I will look into it and code something next week.

[AayushSabharwal]

Great! Thanks for the prompt attention to this issue.

[PBrdng]

Did you take a look at the test folder? There are some examples of custom homotopies.

[AayushSabharwal]

That seems to test existing homotopies in HC. I'm interested in implementing my own AbstractHomotopy.

[AayushSabharwal]

Also worth noting is that using an AbstractSystem in a StraightLineHomotopy requires implementing taylor! for various orders and where x::TaylorVector{N} for various N. TaylorVector isn't particularly clear about how it's meant to be used.

[saschatimme]

A TaylorVector{K} of length n is a vector where each entry is a Tuple representing the first K summands of a Taylor series .
So for K=4 you have in each entry f(t), f’(t),f’’(t)/2,f’’’(t)/6 with f(t)=h(x_current,t)
These derivatives should be computable using AD.

I think you can assume that the taylor! function only needs to be able to compute an output of order K for an input of order K-1.

[AayushSabharwal]

Thanks for the clarification!

I think you can assume that the taylor! function only needs to be able to compute an output of order K for an input of order K-1.

For the signature HC.ModelKit.taylor!(u::TaylorVector{A}, ::Val{B}, sys::AbstractSystem, x::TaylorVector{C}, p), my custom system always seems to gets called with B == C. Granted, the solve fails because the taylor value is wrong, but it doesn't seem to be the case that C == B-1 always.

I have a few followup questions if you don't mind:

1. If my system has m variables, will the argument x be an m-length vector?
2. Can I assume A >= B is always true?
3. Can I assume B >= C is always true?
4. For a system with m variables, the kth term in the multivariate taylor series has >=m terms. Where do the extra values go in u?
5. Is it expected that I only need to fill in the B-order terms and can assume that u contains correct values for terms up to order B-1, or do I need to fill in all the terms of all orders?

[saschatimme]

, my custom system always seems to gets called with B == C. Granted, the solve fails because the taylor value is wrong, but it doesn't seem to be the case that C == B-1 always.

I am not 100% sure but I think in all my wisdom, I treated B as the mathematical order, wher as C is the length of the tuple which is always highest order + 1.
So the first call will be taylor!(u, Val{1}, xt::TaylorVector{1}). Where xt is basically just [(x[1],),....(x[n],))] and the return value should be f'(t) from my previous comment.

Trying to anwser your questions out of my head:

1. If my system has m variables, will the argument x be an m-length vector?

That should be the case

Can I assume A >= B is always true?
Can I assume B >= C is always true?

Both yes. But if I understand correctly you actually don't need to bother about the complicated
signature HC.ModelKit.taylor!(u::TaylorVector{A}, ::Val{B}, sys::AbstractSystem, x::TaylorVector{C}, p).
Just implementing
HC.ModelKit.taylor!(u::AbstractVector, ::Val{B}, sys::AbstractSystem, x::TaylorVector{C}, p)
is enough.
Only this signature is used in the predictor

HomotopyContinuation.jl/src/predictor.jl

Lines 209 to 210 in 925f6ae

	taylor!(u, Val(1), H, x, t)
	u .= .-u

You need the more complicated signature only if you want to enable composition of homotopies or similar which I don't think is the case for you.

So most likely you can only implement B == C where you write into U the Bth summand of the taylor series.

For a system with m variables, the kth term in the multivariate taylor series has >=m terms. Where do the extra values go in u?

This is a univariate Taylor series in t.
Here is a symbolic-numeric method to test your implementation

HomotopyContinuation.jl/test/homotopies_test.jl

Line 26 in 925f6ae

function test_homotopy_taylor(homotopy, symbolic_homotopy)

Is it expected that I only need to fill in the B-order terms and can assume that u contains correct values for terms up to order B-1, or do I need to fill in all the terms of all orders?

Should be covered above

[AayushSabharwal]

Thanks for all the clarification, it is incredibly helpful!

This is a univariate Taylor series in t.

Ah right, I confused myself trying to get this to work

[AayushSabharwal]

I've implemented all the method and tested them with the utility functions. Somehow, the evaluate_and_jacobian! function gets called with x = [NaN + NaN * im, NaN + NaN * im], t = NaN + NaN * im at the beginning of the solve. The last non-NaN time is 1.0. I'd really appreciate help with this if anyone has the time.

The stacktrace:

ERROR: InexactError: Float64(NaN + NaN*im)
Stacktrace:
  [1] Real
    @ ./complex.jl:44 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./genericmemory.jl:243 [inlined]
  [4] unsafe_copyto!(dest::Memory{Float64}, doffs::Int64, src::Memory{ComplexF64}, soffs::Int64, n::Int64)
    @ Base ./genericmemory.jl:153
  [5] unsafe_copyto!
    @ ./genericmemory.jl:133 [inlined]
  [6] _copyto_impl!
    @ ./array.jl:308 [inlined]
  [7] copyto!
    @ ./array.jl:294 [inlined]
  [8] copyto!
    @ ./array.jl:319 [inlined]
  [9] evaluate_and_jacobian!(u::Vector{…}, U::Matrix{…}, sys::NonlinearSolveHomotopyContinuation.HomotopySystemWrapper{…}, x::Vector{…}, p::Nothing)
    @ NonlinearSolveHomotopyContinuation ~/Julia/SciML/NL-hc/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl:45
 [10] evaluate_and_jacobian!(u::Vector{…}, U::Matrix{…}, h::NonlinearSolveHomotopyContinuation.GuessHomotopy{…}, x::Vector{…}, t::ComplexF64, p::Nothing)
    @ NonlinearSolveHomotopyContinuation ~/Julia/SciML/NL-hc/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl:161
 [11] evaluate_and_jacobian!
    @ ~/Julia/SciML/NL-hc/lib/NonlinearSolveHomotopyContinuation/src/interface_types.jl:161 [inlined]
 [12] newton!(x̄::Vector{…}, NC::HomotopyContinuation.NewtonCorrector, H::NonlinearSolveHomotopyContinuation.GuessHomotopy{…}, x₀::Vector{…}, t::ComplexF64, J::HomotopyContinuation.Jacobian{…}, norm::HomotopyContinuation.WeightedNorm{…}; μ::Float64, ω::Float64, extended_precision::Bool, accurate_μ::Bool, first_correction::Bool)
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/newton_corrector.jl:102
 [13] newton!
    @ ~/.julia/packages/HomotopyContinuation/HdN5z/src/newton_corrector.jl:80 [inlined]
 [14] step!(tracker::HomotopyContinuation.Tracker{…}, debug::Bool)
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/tracker.jl:858
 [15] step!(endgame_tracker::HomotopyContinuation.EndgameTracker{…}, debug::Bool)
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/endgame_tracker.jl:360
 [16] #track!#151
    @ ~/.julia/packages/HomotopyContinuation/HdN5z/src/endgame_tracker.jl:308 [inlined]
 [17] track!
    @ ~/.julia/packages/HomotopyContinuation/HdN5z/src/endgame_tracker.jl:296 [inlined]
 [18] track(endgame_tracker::HomotopyContinuation.EndgameTracker{…}, x::Vector{…}, t₁::Float64; path_number::Int64, kwargs::@Kwargs{})
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/endgame_tracker.jl:908
 [19] track (repeats 2 times)
    @ ~/.julia/packages/HomotopyContinuation/HdN5z/src/endgame_tracker.jl:901 [inlined]
 [20] serial_solve(solver::HomotopyContinuation.Solver{…}, starts::Vector{…}, progress::ProgressMeter.Progress, stop_early_cb::typeof(HomotopyContinuation.always_false); catch_interrupt::Bool)
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/solve.jl:570
 [21] solve(S::HomotopyContinuation.Solver{…}, starts::Vector{…}; stop_early_cb::Function, show_progress::Bool, threading::Bool, catch_interrupt::Bool)
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/solve.jl:524
 [22] solve(::NonlinearSolveHomotopyContinuation.GuessHomotopy{…}, ::Vararg{…}; show_progress::Bool, threading::Bool, catch_interrupt::Bool, target_parameters::Nothing, stop_early_cb::Function, transform_result::Nothing, transform_parameters::typeof(identity), flatten::Nothing, target_subspaces::Nothing, kwargs::@Kwargs{})
    @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/HdN5z/src/solve.jl:489
 [23] solve(prob::NonlinearProblem{…}, alg::HomotopyContinuationJL{…}; denominator_abstol::Float64, kwargs::@Kwargs{})
    @ NonlinearSolveHomotopyContinuation ~/Julia/SciML/NL-hc/lib/NonlinearSolveHomotopyContinuation/src/solve.jl:120
 [24] solve(prob::NonlinearProblem{…}, alg::HomotopyContinuationJL{…})
    @ NonlinearSolveHomotopyContinuation ~/Julia/SciML/NL-hc/lib/NonlinearSolveHomotopyContinuation/src/solve.jl:106
 [25] top-level scope
    @ REPL[169]:1
 [26] top-level scope
    @ none:1
Some type information was truncated. Use `show(err)` to see complete types.

I also logged the inputs and outputs to all function calls:

("JAC", x, t, u, U) = ("JAC", ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], 1.0 + 0.0im, ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], ComplexF64[9.0 + 0.0im 2.0 + 0.0im; -4.5 + 0.0im 4.5 + 0.0im])
("JAC", x, t, u, U) = ("JAC", ComplexF64[1.4901161193847656e-14 + 0.0im, 1.4901161193847656e-14 + 0.0im], 1.0 + 0.0im, ComplexF64[1.6431300764452317e-13 + 0.0im, 0.0 + 0.0im], ComplexF64[8.999999999999865 + 0.0im 1.9999999999999554 + 0.0im; -4.4999999999999325 + 0.0im 4.4999999999999325 + 0.0im])
("EVAL", x, t, u) = ("EVAL", ComplexF64[-3.6384955654695295e-17 + 0.0im, -3.6384955654695295e-17 + 0.0im], 1.0 + 0.0im, ComplexF64[0.0 + 0.0im, 0.0 - 0.0im])
("JAC", x, t, u, U) = ("JAC", ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], 1.0 + 0.0im, ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], ComplexF64[9.0 + 0.0im 2.0 + 0.0im; -4.5 + 0.0im 4.5 + 0.0im])
("TAY", x, t, 1, u) = ("TAY", ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], 1.0 + 0.0im, 1, ComplexF64[0.0 + 0.0im, 0.0 - 0.0im])
(x, t) = (HomotopyContinuation.ModelKit.TruncatedTaylorSeries{2, ComplexF64}[(0.0 + 0.0im, -0.0 + 0.0im), (0.0 + 0.0im, 0.0 + 0.0im)], 1.0 + 0.0im)
("TAY", x, t, N, u) = ("TAY", HomotopyContinuation.ModelKit.TruncatedTaylorSeries{2, ComplexF64}[(0.0 + 0.0im, -0.0 + 0.0im), (0.0 + 0.0im, 0.0 + 0.0im)], 1.0 + 0.0im, 2, ComplexF64[0.0 + 0.0im, 0.0 + 0.0im])
(x, t) = (HomotopyContinuation.ModelKit.TruncatedTaylorSeries{3, ComplexF64}[(0.0 + 0.0im, -0.0 + 0.0im, 0.0 + 0.0im), (0.0 + 0.0im, 0.0 + 0.0im, -0.0 - 0.0im)], 1.0 + 0.0im)
("TAY", x, t, N, u) = ("TAY", HomotopyContinuation.ModelKit.TruncatedTaylorSeries{3, ComplexF64}[(0.0 + 0.0im, -0.0 + 0.0im, 0.0 + 0.0im), (0.0 + 0.0im, 0.0 + 0.0im, -0.0 - 0.0im)], 1.0 + 0.0im, 3, ComplexF64[0.0 + 0.0im, 0.0 + 0.0im])

"JAC" is a call to evaluate_and_jacobian!, "EVAL" to evaluate! and "TAY" to taylor!