Customizing ParameterHomotopy to deal with parameters that are really rotation matrices

[GGooeytoe]

I believe this is in a similar vein to #605 .

I have a very large parameterized system (19 equations and variables, mixed volume 135,008) where some of the parameters are really rotation matrices. I would like to use parameter homotopy to help with this system: 1. polyhedral start is slow with such a large mixed volume 2. for randomly generated parameter values (real or complex) the number of solutions found by polyhedral start is not stable due to large numbers of path tracking failures*. However, I observe that running a parameter homotopy between two random instances tends to mostly fail to track or roots go to infinity (there should be around 736 non-singular solutions, but we will usually only recover 100-200).

When we linearly interpolate the rotation matrix structure is destroyed in the middle, which may be causing difficulties. As a test, if I randomly generate 100 parameter values that differ in everything except the rotation matrix parameters, parameter homotopy is quite successful (not 100%, but we vary between 724 and 736 non-singular solutions found). Is a reasonable next step to attempt to implement a parameter homotopy that maintains the rotation matrix structure?

I will confess the fact that this all needs to happen in complex numbers confuses me as I don't really have a handle on the appropriate generalization of a rotation matrix to the complex numbers. I have been using the definition in Sommese and Wampler that we continue to enforce determinant of +1 and transpose yielding the inverse, but I sometimes wonder if there is a trick with quaternions here.

*Due diligence: I have experimented with switching to the conservative path tracking settings and upping the degree of autodifferentiation but it hasn't had a substantial effect either way.

[PBrdng]

The fact that you are linearly interpolating rotation matrices is most certainly the reason why the path tracking fails. A homotopy inside the space of rotation matrices would be much better. In fact, this is "almost" implemented! The shortest path inside the Stiefel manifold ($n\times k$ matrices with orthonormal columns) connecting to matrices corresponding to two linear spaces is here:

HomotopyContinuation.jl/src/linear_subspaces.jl

Line 656 in 6a64320

"""

It should be straight-forward to adapt this to your setting.

See also this test:

HomotopyContinuation.jl/test/homotopies_test.jl

Line 156 in 6a64320

@testset "IntrinsicSubspaceProjectiveHomotopy" begin

that defines the subspace homotopy manually. I think using this as a blueprint is probably the easiest wayto adapt this to your setting.

PS: I move this to discussions, since it's not really an issue.

[GGooeytoe]

Thanks for the pointers, this is helpful! I also found the discussion in #609 and in the tutorial helpful. However, I am still unclear as to how to implement the taylor! function. I don't understand what, mathematically, taylor! is computing to store in u, and what is provided in x?

The relevant call signature seems to be (from predictor.jl):

taylor!(u, v,H,x,t,incr)

Based on test_homotopy_taylor, the inputs are:
u is a length m vector to store the computed quantity, for a system of m equations
v is a Val(int) representing the degree of the Taylor expansion
H is the Homotopy object
x is a TaylorVector{v} (but what is stored in x?)
t is the current parameter value

The discussion in #609 seems to suggest u is the coefficients of the Taylor polynomial, f^(v-1)(x,t)/v!, which is also what I would expect, but the things passed in to the symbolic_homotopy in test_homotopy_taylor seem to be computing something much more complicated...the coefficients of the Taylor series in t about the sampled t, but assuming that x depends in a particular randomly generated polynomial way on t.

So perhaps the key is what exactly is stored in x and what we are supposed to do with it. In the parameter_homotopy.jl, I see that taylor! is implemented by calling the taylor! function of the parametric system, with the fifth argument containing tuples whose first element are the current parameter value and whose second element are the start parameter minus the target parameter. This second element is the derivative of the current parameter with respect to t; and so the correct derivative of the homotopy with respect to t is the derivative of the parametric system times the second element raised to the power k. Are we supposed to push the trailing elements of this argument through the derivative?

Once this gets cleared up I would be happy to make a PR with updates to the documentation to try to clarify this; #609 and the tutorial that resulted from it helped but I think pulling some more of the info into the docs themselves would be good.

EDIT: replace an ambiguous reference to "x" with the index of the argument