Using HomotopyContinuation.jl from Macaulay2

How do I solve my ideal from Macaulay2 with HomotopyContinuation.jl?

Macaulay2 is a great software for symbolic computations and it has many features which are not yet available in Julia. In this guide we show a possible workflow where the symbolic computations are happening in Macaulay2 and the numerical work in Julia.

Setting things up

Our goal is to have a function writeSystem available which accepts as input an ideal I (or a list of polynomials) and a filename f and generates a file f containing the necessary Julia code to compute the solutions of I using HomotopyContinuation.jl.

For example, the Macaulay2 code

R = QQ[x, y]
I = ideal(x^2 + y^2 + 1, x - y - 4)
writeSystem(I, "example1.jl")

should generate a file example1.jl containing

# This file was auto generated from Macaulay2.
# For more information about HomotopyContinuation.jl visit:
#    https://www.JuliaHomotopyContinuation.org
using HomotopyContinuation

@polyvar x y
f = [x^2+y^2+1, x-y-4];

result = solve(f)

The function writeSystem can be achieved by the following code

You can copy and paste this code in your current Macaulay2 session and then you have the writeSystem function available. But since it is annoying to define this for every new session again we recommend to add the code to your init.m2 file.

Acknowledgement: A first version of the code was written by Tim Duff during a Macaulay2 workshop at the Max-Planck-Institute for Mathematics in the Sciences in Leipzig.

An example from Schubert calculus

Let $M = (m_{ij})$ be a $3 \times 5$ matrix where the entries $m_{ij} \in \mathbb{Q}[x,y,z]$ are generic polynomials of degree 2. Schubert calculus tells us that the ideal $I$ generated by the maximal of minors of $M$ is zero dimensional with degree 80.

The code to generate $I$ is

R = QQ[x,y,z];
M = matrix apply(3,u->apply(5,t->random(2,R)+random(1,R)+random(0,R)))
I = minors(3,M);

Now we can use our new writeSystem method to generate the necessary Julia code

writeSystem(I, "minors35.jl")

Now let's start a Julia session and load the newly generated file

julia> include("minors35.jl")
Result with 80 solutions
==================================
• 80 non-singular solutions (4 real)
• 0 singular solutions (0 real)
• 216 paths tracked
• random seed: 488250

and we see that we found the correct number of solutions.