This guide explains how to implement a homotopy of the form
$$H(x,t).$$
In HomotopyContinuation.jl
a homotopy is a struct
for which one must implement a few (in-place) functions, namely
ModelKit.evaluate!
ModelKit.evaluate_and_jacobian!
ModelKit.taylor!
As the names suggest, these functions
- evaluate $H(x,t)$ at $(x,t)$,
- evaluate $H(x,t)$ and compute the jacobian at $(x,t)$, and
- compute the taylor expansion of $H(x,t)$ with respect to $t$ at $(x,t)$.
In fact, $H(x,t)$ does not even need to be a polynomial function, neither in $x$ nor $t$. The algorithm only depends the implementation of these three functions and tracks $t$ from $t=1$ to $t=0$.
Implementing an example homotopy
Let us illustrate how this works in an example. We want to implement a homotopy that, given systems of polynomials $F(x)$ and $G(x)$, is given by
$$H(x,t) = \gamma\cdot \sin(t\cdot \pi/2)\cdot G(x) + \cos(t\cdot \pi/2) \cdot F(x), \quad 0\leq t\leq 1,$$
where $\gamma \in\mathbb C$.
In this case, we have $H(x,1)=G(x)$ and $H(x,0)=F(x)$. That is, $G(x)$ is our start system and $F(x)$ our target system.
This is similar to a straight line homotopy $H(x,t) = \gamma \cdot t \cdot G(x) + (1-t)\cdot F(x)$.
We copy the code for a StraightLineHomotopy
, adapt it to our setting and call the result TrigonometricHomotopy
.
In the following, we explain some of the implemented functions, but not all. The complete code is at the end of this guide.
The struct
that defines our TrigonometricHomotopy
could look like this:
using HomotopyContinuation
struct TrigonometricHomotopy{S<:AbstractSystem, T<:AbstractSystem} <: AbstractHomotopy
start::S
target::T
γ::ComplexF64
half_pi::Float64
u::Vector{ComplexF64}
ū::Vector{ComplexF64}
v̄::Vector{ComplexF64}
U::Matrix{ComplexF64}
dv_start::TaylorVector{5,ComplexF64}
dv_target::TaylorVector{5,ComplexF64}
end
Here, start
and target
are start and target system, respectively. We also save $\pi/24 in half_pi
so that we don't need to compute it all the time. The remaining entries are for caching.
The evaluate!
function comes next:
function ModelKit.evaluate!(
u,
H::TrigonometricHomotopy,
x::Vector{ComplexF64},
t,
p = nothing,
)
evaluate!(H.v̄, H.start, x)
evaluate!(H.ū, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ .* sin(t0), cos(t0)
for i = 1:size(H, 1)
@inbounds u[i] = ts * H.v̄[i] + tt * H.ū[i]
end
end
Here, ts
and tt
are the coefficients at $t$ of $G$ and $F$, respectively. Notice that the coefficient of $G$ is indeed $\gamma \cdot \sin(t \cdot \pi/2)$ and the coefficient of $F$ is $\cos(t \cdot \pi/2)$.
The evaluate_and_jacobian!
is similar, but includes the Jacobian:
function ModelKit.evaluate_and_jacobian!(
u,
U,
H::TrigonometricHomotopy,
x::AbstractVector{T},
t,
p = nothing,
) where {T}
evaluate_and_jacobian!(u, U, H.start, x)
evaluate_and_jacobian!(H.u, H.U, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ .* sin(t0), cos(t0)
for i = 1:size(H, 1)
@inbounds u[i] = ts * u[i] + tt * H.u[i]
end
for j = 1:size(H, 2), i = 1:size(H, 1)
@inbounds U[i, j] = ts * U[i, j] + tt * H.U[i, j]
end
nothing
end
Finally, we implement ModelKit.taylor!
:
function ModelKit.taylor!(u, ::Val{1}, H::TrigonometricHomotopy, x, t)
evaluate!(u, H.start, x)
evaluate!(H.u, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ * cos(t0) * H.half_pi, -sin(t0) * H.half_pi
for i = 1:size(H, 1)
@inbounds u[i] = ts * u[i] + tt * H.u[i]
end
u
end
Here, the coefficients are given by differentiating $\gamma \sin(t\cdot \pi/2)$ and $\cos(t\cdot \pi/2)$ at $t$. The input ::Val{1}
indicated that this function returns the first derivative with respect to $t$. We also have to implement higher order derivatives (see the complete code at the end of this guide).
Running our custom homotopy
Let us use our homotopy on the following data:
$$G(x) = \begin{pmatrix} x^2 - 1\\ y^2 - 1\end{pmatrix}, \quad F(x) = \begin{pmatrix} y^2 + xy + 1\\ -x^2 + x + 2y - 2\end{pmatrix}$$
choosing $\gamma\in\mathbb C$ at random.
@var x y
G = System([x^2 - 1; y^2 - 1])
F = System([y^2 + x * y + 1; -x^2 + x + 2y - 2])
γ = exp(rand() * 2 * pi * im)
H = TrigonometricHomotopy(G, F; γ = γ)
Now, H
is a TrigonometricHomotopy
and we can track the zeros of $G$ from $t=1$ to $t=0$ along this homotopy.
julia> starts = [[1,1], [1,-1], [-1,1], [-1,-1]]
julia> S = solve(H, starts)
Result with 4 solutions
=======================
• 4 paths tracked
• 4 non-singular solutions (0 real)
• random_seed: nothing
The computation was successful and all $4$ zeros of $F$ were computed!
Comparing with a straight line homotopy
Let us compare our newly implemented homotopy with a StraightLineHomotopy
.
H_SL = StraightLineHomotopy(G, F; γ = γ)
We wish to plot every single point tracked along the way. For this, we can use the following code:
T = Tracker(H)
T_SL = Tracker(H_SL)
pts = Vector{ComplexF64}()
for (x, t) in iterator(T, starts[1], 1.0, 0.0)
push!(pts, x[1])
end
pts_SL = Vector{ComplexF64}()
for (x, t) in iterator(T_SL, starts[1], 1.0, 0.0)
push!(pts_SL, x[1])
end
Now, pts
contains the list of points tracked along our TrigonometricHomotopy
when starting at starts[1] = [1,1]
. Similarly, pts_SL
contains the points tracked along a StraightLineHomotopy
.
We plot the result:
using Plots
scatter(real(pts), imag(pts), markercolor = :steelblue, markersize = 10, legend = false)
scatter!(real(pts_SL), imag(pts_SL), markercolor = :indianred, markersize = 10)
for i in 2:length(pts)
v = pts[i] - pts[i-1]
a = pts[i-1:i] + [0.25 .* v; -0.25 .* v]
plot!(real(a) , imag(a) , linecolor = :black, linewidth = 2, arrow=true, opacity = 0.75)
end
plot!()
The blue points have been tracked by our TrigonometricHomotopy
while the red points have been tracked by the StraightLineHomotopy
.
As we can see, both homotopies start and arrive at the same point, but the intermediate points are different.
The complete code for our example
struct TrigonometricHomotopy{S<:AbstractSystem, T<:AbstractSystem} <: AbstractHomotopy
start::S
target::T
γ::ComplexF64
half_pi::Float64
u::Vector{ComplexF64}
ū::Vector{ComplexF64}
v̄::Vector{ComplexF64}
U::Matrix{ComplexF64}
dv_start::TaylorVector{5,ComplexF64}
dv_target::TaylorVector{5,ComplexF64}
end
function TrigonometricHomotopy(
start::System,
target::System;
kwargs...,
)
TrigonometricHomotopy(
fixed(start),
fixed(target);
kwargs...,
)
end
function TrigonometricHomotopy(
start::AbstractSystem,
target::AbstractSystem;
γ = 1.0,
gamma = γ,
)
size(start) == size(target) || throw(
ArgumentError(
"Start and target do not have the same size, got $(size(start)) and $(size(target))",
),
)
m, n = size(start)
u = zeros(ComplexF64, m)
ū = zeros(ComplexF64, m)
v̄ = zeros(ComplexF64, m)
U = zeros(ComplexF64, m, n)
half_pi = π / 2
dv_start = TaylorVector{5}(ComplexF64, m)
dv_target = TaylorVector{5}(ComplexF64, m)
TrigonometricHomotopy(
start,
target,
ComplexF64(gamma),
half_pi,
u,
ū,
v̄,
U,
dv_start,
dv_target,
)
end
Base.size(H::TrigonometricHomotopy) = size(H.start)
function Base.show(io::IO, mime::MIME"text/plain", H::TrigonometricHomotopy)
println(io, typeof(H), ":")
println(io, "γ: ", H.γ)
println(io, "\nG: ")
show(io, mime, H.start)
println(io, "\n\nF:")
show(io, mime, H.target)
end
function ModelKit.evaluate!(
u,
H::TrigonometricHomotopy,
x::Vector{ComplexF64},
t,
p = nothing,
)
evaluate!(H.v̄, H.start, x)
evaluate!(H.ū, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ .* sin(t0), cos(t0)
for i = 1:size(H, 1)
@inbounds u[i] = ts * H.v̄[i] + tt * H.ū[i]
end
end
function ModelKit.evaluate!(u, H::TrigonometricHomotopy, x, t, p = nothing)
evaluate!(u, H.start, x)
evaluate!(H.u, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ .* sin(t0), cos(t0)
for i = 1:size(H, 1)
@inbounds u[i] = ts * u[i] + tt * H.u[i]
end
u
end
function ModelKit.evaluate_and_jacobian!(
u,
U,
H::TrigonometricHomotopy,
x::AbstractVector{T},
t,
p = nothing,
) where {T}
evaluate_and_jacobian!(u, U, H.start, x)
evaluate_and_jacobian!(H.u, H.U, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ .* sin(t0), cos(t0)
for i = 1:size(H, 1)
@inbounds u[i] = ts * u[i] + tt * H.u[i]
end
for j = 1:size(H, 2), i = 1:size(H, 1)
@inbounds U[i, j] = ts * U[i, j] + tt * H.U[i, j]
end
nothing
end
function ModelKit.taylor!(u, ::Val{1}, H::TrigonometricHomotopy, x, t)
evaluate!(u, H.start, x)
evaluate!(H.u, H.target, x)
t0 = H.half_pi * t
ts, tt = H.γ * cos(t0) * H.half_pi, -sin(t0) * H.half_pi
for i = 1:size(H, 1)
@inbounds u[i] = ts * u[i] + tt * H.u[i]
end
u
end
function ModelKit.taylor!(
u,
v::Val{K},
H::TrigonometricHomotopy,
tx::TaylorVector,
t,
) where {K}
taylor!(H.dv_start, v, H.start, tx)
taylor!(H.dv_target, v, H.target, tx)
t0 = H.half_pi * t
for i = 1:size(H, 1)
start = H.γ * (cos(t0) * H.half_pi * H.dv_start[i, K] + sin(t0) * H.dv_start[i, K+1])
target = cos(t0) * H.dv_target[i, K+1] - sin(t0) * H.half_pi * H.dv_target[i, K]
u[i] = start + target
end
u
end