Some queries from a new users

[TorkelE]

I was happily surprised to discovered that homotopy continuation was implemented in Julia. I have been trying the package out, however there are a few non obvious things which I do not understand (this might partially be due to my lack of experience in the area). I am unsure whenever this is the right place to ask my questions, but it was the best I found (if not I apologise).

First, it seems that when there is a singularity present an additional 1.0 is added to the solution array, e.g.

@polyvar x y z w 
f = [x*y+z^2, y-3x*z+z^2+2, x^3+y^2+z-1];
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
   println(real(sols[i].solution))
end

outputs

[3.62126e-46, 1.50418e-25, 1.0, 8.62497e-7]
[-1.00432e-24, -7.17904e-17, 1.0, -7.61496e-7]
[-1.40608, 2.36735, -1.82447]
[-76.8682, 674.107, -227.634]

where both the first and second solution contains a singularity. However, I did not find anything on what the actual meaning of the 1.0 and its position is.

Next I have a few things which all might be regarding my lack of knowledge in the area, I am no expert. However, when I run the part

sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
   println(real(sols[i].solution))
end

several times I do not all time get the previously mentioned solution, sometimes the solution [-76.8682, 674.107, -227.634] is not found (roughly half of the time). In addition this solution does not actually seem to be a good solution, indeed

sol = [-76.8682, 674.107, -227.634]
for i = 1:length(f)
   println(f[i](x=> sol[1] , y=>sol[2], z=>sol[3]),"   ") 
end

returns

-0.15374140000494663   
-0.10256040000604116   
-1.070417910619625   

is the solving method simply inaccurate for certain cases, not running long enough, or am I doing something wrong? The documentations says

The aim of this project is twofold: establishing a fast numerical polynomial solver in Julia and at the same time providing a highly customizable algorithmic environment for researchers for designing and performing individual experiments.

Then the documentations starts with a simple example of solving a polynomial system (what I want to do and also what I am trying here), and then continue with advanced stuff of defining your own homotopies and pathtracking. I am uncertain where the "fast numerical polynomial solver" (what I am interested in) ends and the "highly customizable algorithmic environment for researchers for designing and performing individual experiments" starts.

If the issue lies in my ability to use the methods and that they require some experience, do you maybe have a suggested reading for understanding more?

A further example where I struggle. For 4 variable systems I have been unable to find proper solutions for even simple linear systems like

@polyvar x y z w 
f = [w+x-4y, z+3-2x, 3w+x, 10-w+y];
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
   println(real(sols[i].solution))
end

This gives

[6.66667, -20.0, -3.33333, -43.0]

However,

sol = [6.66667, -20.0, -3.33333, -43.0]
for i = 1:length(f)
   println(f[i](x=> sol[1] , y=>sol[2], z=>sol[3], w=>sol[4]),"   ") 
end

outputs

43.666669999999996   
-13.66667   
-122.33333   
33.0

which is distinctively non zero.

[saschatimme]

Hi

I am unsure whenever this is the right place to ask my questions, but it was the best I found (if not I apologise).

Absolutely, but sorry for the late response.

first, it seems that when there is a singularity present an additional 1.0 is added to the solution array, [...] However, I did not find anything on what the actual meaning of the 1.0 and its position is.

This is definitely not obvious and should probably be handled better. If you look closely you see that the solutions with the additional 1.0 have 4 entries while the other ones have only 3. The reason for this is that all calculations are executed in projective space (which is beneficial for multiple reason). Let's make this concrete. In your example we have the 3 variables x, y, z, so any solution of this system in in C^3. We then can introduce a new variable a to make the system f homogenous:

f_hom := [x*y+z^2, a*y-3x*z+z^2+2a^2, x^3+y^2*a+z*a^2-a^3];

For any zero (x, y, z) of f the value (1, x, y, z) is a zero of f_hom. And for any zero (a, x, y, z) of f_hom with a != 0 the value (x / a, y / a, z / a) is a zero of f. Therefore, the homogenous system always has at least the same number of solutions as the original system. But it can have more solutions, namely those with a=0. Those solutions, we call solutions at infinity. You filter these away by using

solutions(solve(f), only_real=true, at_infinity=false)

(Currently the default is at_infinity=true which should be changed.)

several times I do not all time get the previously mentioned solution, sometimes the solution [-76.8682, 674.107, -227.634]

This solutions seems to be numerically a little bit tricky, but if I make the path precision more restrictive (solve(f, path_precision=1e-8)) I always get both solutions. Note that the zero is very sensible against small perturbations, i.e. you need to use the exact value and not the abbreviation printed.

For the linear system. I cannot reproduce this here. I get as a solution

[-20.0, -3.33333, -43.0, 6.66667]

which is the zero of the system (although the precision is not as good as it should). I think something got messed up with the variable ordering. (In DynamicPolynomials the ordering of the variables is determined by the time they were created.)

I am note sure about your mathematical background, so if I should go into more details just let me know. And I know that everything is a little bit unpolished at the moment, but soon everything should get much better :)

[saschatimme]

Update: I made a mistake and edited the above answer accordingly.

[TorkelE]

Thank you, I understand the addition 1.0 now.

You are probably right about the linear system as well, the reordering indeed works, and if I restart my Kernel and do everything again I get the right order.

I would probably consider myself a mathematician, although I only work with applications these days. I will look more at this, trying out more complicated systems, might return with more queries. Exactly at what point do you need to start thinking about setting up your own homotopies and pathtracking?

This is good stuff. Cheers.

[TorkelE]

Exactly would you recommend working with path_precision? I tried a somewhat more complicated system: g = [x*w^3+x-4y^2, z*x+3-2x*w*z^2, 3w^4-2x, 10-w*y^2-6*z^3+y];, concentrating on real, finite, positive, solutions. If I do not set any path_precision I get these solutions, with corresponding evaluations of the 4 polynomials:

#Solution 1
[1.7336, 0.957339, 1.18591, 1.03685]
1.9984014443252818e-15   
1.3322676295501878e-15   
0.0   
5.329070518200751e-15   

#Solution 2
[7.32158, 2.80021, 0.575777, 1.48638]
-2.6645352591003757e-15   
-8.881784197001252e-16   
-1.7763568394002505e-15   
5.329070518200751e-15   

#Solution 3
[1.28099e7, 5.82889e6, 847.792, 0.0755278]
-1.3590379079403522e14   
-1.3799325627831165e12   
-2.5619887553624373e7   
-2.569777615087e12   

#Solution 4
[9.86097e7, 7.28712e6, 5006.43, 0.0758873]
-2.1240852455788066e14   
-3.746305747334043e14   
-1.9721942159340465e8   
-4.7826684160198e12   

#Solution 5
[1.11219e8, 4.52783e6, 8557.27, 0.0759869]
-8.200494488982812e13   
-1.236754840151639e15   
-2.224378582484602e8   
-5.317552381002262e12   

The first two are correct. However, I do get another few (very) faulty ones. If I set path_precision=1e-8 I only get these solutions:

#Solution 1
[1.7336, 0.957339, 1.18591, 1.03685]
-2.4424906541753444e-15   
-4.440892098500626e-16   
0.0   
0.0   

#Solution 2
[7.32158, 2.80021, 0.575777, 1.48638]
-2.3447910280083306e-13   
1.4210854715202004e-14   
1.7763568394002505e-14   
0.0 

which are the same as the two correct ones as in the first case. It seems like having a lower path_precision can ensure that I get the correct set of solutions?

However; Once I have a set of solutions it is easy to discard the ones that are wrong. What is important is that my generated set of solutions include all the correct ones, without missing any. It feels easier, rather than trusting that I have chosen the right path_precision, to simply solve it and discard the wrong solutions?

[PBrdng]

Hi,

what happens is the following: As Sascha mentioned above we are computing the solutions in projective space. This is beneficial because algebraic geometry provides a collection of genericity statements for projective space, which is useful to exploit in homotopy continuation.

Long story short: The computations are made with an additional variable u that is used to homogenize your system g. Then, after the solutions have been found it is checked whether they were successful by simply plugging them into the polynomial and checking whether you get zero. And here is where the problem starts: For almost singular solutions (note that your solutions 3-5 are almost singular, solutions(solve(g), singular = false) does not return them), the evaluation of the homogenized g at the projective solution [x,y,z,w,u] and the evaluation of g at the affine solution [x,y,z,w] may differ by an order of 1e30. This is why in the next release of the package the :success of a solution is not decided by the projective solution, but by the affine.

Increasing the path-precision will correctly identify solutions 3-5 as non-successful, which is why you don't see them anymore.

PS: Sorry for the late answer.

[PBrdng]

If you are interested in learning more about numerical algebraic geometry, the following summer school may be of interest to you: https://www.mis.mpg.de/calendar/conferences/2018/nc2018.html.

[TorkelE]

Thank you, that is useful. Checking the summer school out. cheers

[saschatimme]

@Gaussia A lot of things improved and I now consistently get the correct solutions in your examples with the current master branch.

[TorkelE]

Sounds good! Really appreciate the work you have put in implementing this in Julia.