Explicitly solvable subset of equations

[jarlebring]

I have a system where parts of the equations can be solved explicitly.

@var  a[1:10]
Random.seed!(0);
A=randn(10,10);
expr=A*a-randn(10) + a.^2;
expr .+= 3.0*a[10].^4
expr[end]=1-a[10];
expr[end-1]=expr[end-1]+ a[9]^12*(1-a[10]^8);

HomotopyContinuation can solve this system.

julia> expr[end]
1 - a₁₀
julia> @time solve(sys)
Tracking 3072 paths... 100%|████████████████████████████████████████| Time: 0:00:33
  # paths tracked:                  3072
  # non-singular solutions (real):  460 (0)
  # singular endpoints (real):      52 (0)
  # total solutions (real):         512 (0)
 33.789564 seconds (164.15 k allocations: 9.301 MiB, 1.15% gc time)
Result with 512 solutions
=========================
• 3072 paths tracked
• 460 non-singular solutions (0 real)
• 52 singular solutions (0 real)
• random_seed: 0x371d6fcf
• start_system: :polyhedral
• multiplicity table of singular solutions:
╭───────┬───────┬────────┬────────────╮
│ mult. │ total │ # real │ # non-real │
├───────┼───────┼────────┼────────────┤
│   1   │  52   │   0    │     52     │
╰───────┴───────┴────────┴────────────╯

The last equation is an equation in one variable, and if the solution is substituted into the other equations, the number of equations reduces by one, and (more importantly) the degree of the system reduces substantially. If I reduce the system by repeated substitutions it's all much faster.

julia> expr2=greedy_eliminate(expr)
(j, variables(thisvars)) = (10, Variable[a₁₀])
9-element Vector{Expression}:
     5.41756338269563 - 0.231909069576951*a₁ - 0.269884920377936*a₂ - 1.07457699497494*a₃ + 1.72102099090346*a₄ - 0.516791014396629*a₅ + 1.86106051520025*a₆ - 0.577960019067313*a₇ + 0.272619688706261*a₈ + 0.849341440756862*a₉ + a₁^2
      1.65642535758593 + 0.940390210248594*a₁ + 0.578003602234251*a₂ - 1.3819748518823*a₃ - 0.670523187729683*a₄ - 1.55823128446815*a₅ + 0.230189216512471*a₆ - 0.907194592696011*a₇ + 1.23520646835216*a₈ - 0.210716011210666*a₉ + a₂^2
        2.5004031428378 + 0.596761671133522*a₁ + 1.16034424207562*a₂ + 0.951852764202713*a₃ + 0.257971553840438*a₄ - 0.96394519056479*a₅ - 1.63450222498368*a₆ - 1.36705348269359*a₇ + 1.11072340702794*a₈ - 0.143171444497545*a₉ + a₃^2
     2.26469417268335 + 1.99782409305779*a₁ + 0.287887995585263*a₂ - 0.292291265796998*a₃ + 0.641934935288121*a₄ + 1.25945167398595*a₅ + 0.814542054939574*a₆ + 0.660498208562693*a₇ + 1.85157752048501*a₈ - 0.121578262796389*a₉ + a₄^2
 3.02930063010366 - 0.0515618422032261*a₁ - 0.441192823034224*a₂ + 0.328946043058385*a₃ - 0.280151928850431*a₄ - 1.3732762654745*a₅ + 0.0719542149895912*a₆ + 0.0371086197757989*a₇ - 0.131820764083343*a₈ + 0.473540021560705*a₉ + a₅^2
  4.89360632089657 + 1.17224486270216*a₁ + 0.134575710473665*a₂ - 0.495998923989969*a₃ + 0.0211415509550022*a₄ - 0.26477164348276*a₅ + 0.134469080567582*a₆ + 0.877643992842206*a₇ - 0.406776052697414*a₈ - 0.0923910828394263*a₉ + a₆^2
      3.08731284450791 - 1.69681301928184*a₁ + 0.165542550919225*a₂ + 0.222770311442103*a₃ - 1.03024535034903*a₄ + 1.24569386430263*a₅ - 0.895554435418489*a₆ + 1.79295575445521*a₇ - 0.379297365028804*a₈ - 0.550620750668124*a₉ + a₇^2
      2.84114029353575 - 2.11614596856045*a₁ - 1.07209741544676*a₂ - 1.66125038995819*a₃ + 0.456951281758197*a₄ - 0.255171803080274*a₅ + 0.758910549437524*a₆ - 1.35784918731374*a₇ - 0.499402525822446*a₈ - 0.474731751909035*a₉ + a₈^2
       1.76701464747205 + 0.55836560578351*a₁ - 1.26016429538773*a₂ + 1.92523528739124*a₃ + 0.281846208971383*a₄ + 1.32342824202356*a₅ + 0.860579389761707*a₆ - 1.48325950792007*a₇ - 0.635276483651858*a₈ + 0.517733167592105*a₉ + a₉^2

julia> sys2=System(expr2);

julia> @time solve(sys2)
Tracking 512 paths... 100%|████████████████████████████████████████| Time: 0:00:00
  # paths tracked:                  512
  # non-singular solutions (real):  512 (0)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         512 (0)
  0.641606 seconds (82.14 k allocations: 3.657 MiB)
Result with 512 solutions
=========================
• 512 paths tracked
• 512 non-singular solutions (0 real)
• random_seed: 0xbb2fafd6
• start_system: :polyhedral

I did the reduction by a greedy substition approach:

function greedy_eliminate(expr)
    modified=true;
    while modified
        modified=false;
        for j=reverse(1:size(expr,1))
            thisvars=variables(expr[j])
            if (size(thisvars,1)==1)
                res=solve(System([expr[j]]));
                real_sol=real_solutions(res);
                if (size(real_sol,1)>0)
                    x=real_sol[1]
                else
                    x=solutions(res)[1];
                end
                @show j,variables(thisvars)
                expr=subs(expr,thisvars => x);
                modified=true;
                deleteat!(expr,j);
            end
            if (modified)
                break
            end
        end
    end
    return expr
end

In my application, the system was more complicated, and it was not possible to solve the system without first doing a substitution trick.

Does this functionality exist in HomotopyContinuation.jl or companion packages? Or some other way to solve sys efficiently without the substitution trick?

If the answer is "No. It doesn't fit this package.", I think this issue can be closed directly.

[PBrdng]

This functionality does not exist in HomotopyContinuation.jl, but it is a powerful trick if applicable.

In fact, it is always possible to transform a system in an equivalent system that has the properties above. Namely, that the first equation is in one variable, the second in two and so on, so that one can solve it by substitution. This transformation is obtained by computing a Gröbner bases for the Lex order. See, for instance, Example 3.5 in this book.

However, computing Gröbner bases is costly and usually to slow for bigger systems. A Julia implementation is available through Oscar.

Notice that, next to solve, you can also use regeneration or monodromy. It somewhat an art to decide which method is the best in what situation. I do knot know a satisfying way to automatize this.

Does this answer your request? I'm happy to discuss more.

[PBrdng]

I know that people have combined both packages for their research. But I do now know a reference from the top of my head. @oskarhenriksson could know?

I know that a "friends of Oscar" package is currently being developed, where an interface to HC.jl should be part of it.

[oskarhenriksson]

The datatypes are indeed different, but it's not so hard to export a HC.jl system to an Oscar system or vice versa. For example, we have some very basic code for doing this at the bottom of this file: https://github.com/oskarhenriksson/TropicalHomotopies.jl/blob/main/src/numerical_solving.jl

[oskarhenriksson]

But as was pointed out above, Gröbner bases will typically not play well with systems where the coefficients are floating point numbers.

One way to get around this could be to try to replace the floating point numbers with parameters $p_1,\ldots,p_k$ and view your system as a parametric system

$$F=(f_1,\ldots,f_n)\in\mathbb{Q}(p_1,\ldots,p_k)[x_1,\ldots,x_n]^n$$

with coefficients that are rational functions with rational coefficients in the parameters. Then you can analyze the system with Gröbner bases in Oscar over the field of rational functions $\mathbb{Q}(p_1,\ldots,p_k)$, export it back to a parametric system in HC.jl, and specialize the parameters to whatever floating point numbers you wanted (provided that they are "generic enough" so that no illegal divisions by zero would have been made if they were used in the Gröbner basis calculation).

In the example above, you could view the entries of the matrix A as parameters, and then at the very end specialize these parameters to randn(10,10).

[oskarhenriksson]

Just let me know if you want to discuss this in more detail, or if you want help with more specific conversions between Oscar and HC.jl systems with/without parameters!

[jarlebring]

Great! 🙏 I think I will can get back to you and these questions once the current paper we are working on is ready. For our purposes we "solved" the specific issue with greedy substitution. For that paper it is enough for us to be sure that there is no other simple way to do the same thing. But indeed, it sounds like more complicated situations can be approaches with these tools.