Equation solving with GroebnerBasis

Question

I have a system of quadratic equations defined below

nx = 1;
ny = 2;
ne = 1;
nt = 2;
n = nx + ny;
ns = 2;

    For[ii = 1, ii <= ns, ii++,
      G[ii] = Table[g[ii][i, j], {i, ny}, {j, nx + ne + 1}];
      H[ii] = Table[h[ii][i, j], {i, nx}, {j, nx + ne + 1}];
      ];

    eqns = {-1.1111111111111112` g[1][1, 
          1] + (-0.9703966865736877` + 1.` g[1][1, 1] + 0.9` g[1][2, 1] + 
           0.11111111111111109` g[2][1, 1] + 
           0.09999999999999998` g[2][2, 1]) h[1][1, 1], 
       10.` g[1][1, 1] - 
        161.` g[1][2, 
          1] + (143.451` g[1][2, 1] + 15.938999999999995` g[2][2, 1]) h[
           1][1, 1], 
       0.7763173492589502` + 0.6199999999999999` g[1][2, 1] - 
        0.9703966865736877` h[1][1, 1], -1.1111111111111112` g[2][1, 
          1] + (-0.9703966865736877` + 0.11111111111111109` g[1][1, 1] + 
           0.09999999999999998` g[1][2, 1] + 1.` g[2][1, 1] + 
           0.9` g[2][2, 1]) h[2][1, 1], 
       10.` g[2][1, 1] - 
        161.` g[2][2, 
          1] + (15.938999999999995` g[1][2, 1] + 143.451` g[2][2, 1]) h[
           2][1, 1], 
       0.7763173492589502` + 0.13999999999999996` g[2][2, 1] - 
        0.9703966865736877` h[2][1, 1]};

    polys = eqns /. {0. -> 0};

This system has some structure.

First, I find the variables

vars = Variables[polys]; 

Deconstruct the polynomials (I know it is a quadratic system). There are 2 nonzero constant terms, 12 nonzero linear terms, and 12 nonzero 2nd-order terms

{m0, m1, m2} = CoefficientArrays[polys, vars]

Confirm the reconstructed polynomials are the same as the originals (just checking)

polys - (m0 + m1.vars + m2.vars.vars) // Expand // Chop 

Of the 12 nonzero 2nd-order terms, of which there are 8 distinct cross-products ...

Flatten[Most[ArrayRules[#]] & /@ m2][[All, 1]] // Union // Length

... none of which are purely quadratic

Diagonal /@ m2 // Union

So, it is a quadratic system but nothing purely quadratic

If I apply GroebnerBasis

GBordering = 
 Flatten[{H[1][[1 ;; nx, 1 ;; nx]], H[2][[1 ;; nx, 1 ;; nx]], 
   G[1][[1 ;; ny, 1 ;; nx]], G[2][[1 ;; ny, 1 ;; nx]]}]
GB = GroebnerBasis[polys, GBordering];
GB // MatrixForm
soln = Solve[Thread[GB == 0]];
soln // MatrixForm

I get the right solutions. Anyone interested, I will send the code. But you can copy this into your computer and it will give you 9 solutions. Everything OK up to here. But if I make the system a larger (only a little)

nx = 2;
ny = 3;
ne = 1;
nt = 1;
n = nx + ny;
ns = 2;
For[ii = 1, ii <= ns, ii++,
  G[ii] = Table[g[ii][i, j], {i, ny}, {j, nx + ne + 1}];
  H[ii] = Table[h[ii][i, j], {i, nx}, {j, nx + ne + 1}];
  ];
eqns = {9.169301858525786` - 
    1.` g[1][3, 
      1] + (-19.453328842995205` + 8.169847955946475` g[1][2, 1] + 
       0.9077608839940525` g[2][2, 1]) h[1][1, 
      1] + (0.` + 8.169847955946475` g[1][2, 2] + 
       0.9077608839940525` g[2][2, 2]) h[1][2, 1], 
   0.` - 1.` g[1][3, 
      2] + (-19.453328842995205` + 8.169847955946475` g[1][2, 1] + 
       0.9077608839940525` g[2][2, 1]) h[1][1, 
      2] + (0.` + 8.169847955946475` g[1][2, 2] + 
       0.9077608839940525` g[2][2, 2]) h[1][2, 2], -1.` g[1][3, 
      1] + (0.` - 1.` g[1][1, 1] + 0.9` g[1][3, 1] - 
       0.11111111111111109` g[2][1, 1] + 
       0.09999999999999998` g[2][3, 1]) h[1][1, 
      1] + (1.0999999999999999` - 1.` g[1][1, 2] + 0.9` g[1][3, 2] - 
       0.11111111111111109` g[2][1, 2] + 
       0.09999999999999998` g[2][3, 2]) h[1][2, 1], -1.` g[1][3, 
      2] + (0.` - 1.` g[1][1, 1] + 0.9` g[1][3, 1] - 
       0.11111111111111109` g[2][1, 1] + 
       0.09999999999999998` g[2][3, 1]) h[1][1, 
      2] + (1.0999999999999999` - 1.` g[1][1, 2] + 0.9` g[1][3, 2] - 
       0.11111111111111109` g[2][1, 2] + 
       0.09999999999999998` g[2][3, 2]) h[1][2, 
      2], -178.88888888888889` g[1][1, 1] - 
    9.` g[1][3, 
      1] + (159.39` g[1][1, 1] + 17.709999999999994` g[2][1, 1]) h[1][
      1, 1] + (159.39` g[1][1, 2] + 17.709999999999994` g[2][1, 2]) h[
       1][2, 1], -178.88888888888889` g[1][1, 2] - 
    9.` g[1][3, 
      2] + (159.39` g[1][1, 1] + 17.709999999999994` g[2][1, 1]) h[1][
      1, 2] + (159.39` g[1][1, 2] + 17.709999999999994` g[2][1, 2]) h[
       1][2, 2], 
   0.` - 0.21199999999999997` g[1][1, 1] + 
    0.99` h[1][2, 1], -0.792` - 0.21199999999999997` g[1][1, 2] + 
    0.99` h[1][2, 2], 0.` + 1.` g[1][2, 1] - 1.` h[1][1, 1], 
   0.` + 1.` g[1][2, 2] - 1.` h[1][1, 2], 
   9.169301858525786` - 
    1.` g[2][3, 
      1] + (-19.453328842995205` + 0.9077608839940525` g[1][2, 1] + 
       8.169847955946475` g[2][2, 1]) h[2][1, 
      1] + (0.` + 0.9077608839940525` g[1][2, 2] + 
       8.169847955946475` g[2][2, 2]) h[2][2, 1], 
   0.` - 1.` g[2][3, 
      2] + (-19.453328842995205` + 0.9077608839940525` g[1][2, 1] + 
       8.169847955946475` g[2][2, 1]) h[2][1, 
      2] + (0.` + 0.9077608839940525` g[1][2, 2] + 
       8.169847955946475` g[2][2, 2]) h[2][2, 2], -1.` g[2][3, 
      1] + (0.` - 0.11111111111111109` g[1][1, 1] + 
       0.09999999999999998` g[1][3, 1] - 1.` g[2][1, 1] + 
       0.9` g[2][3, 1]) h[2][1, 
      1] + (1.0999999999999999` - 0.11111111111111109` g[1][1, 2] + 
       0.09999999999999998` g[1][3, 2] - 1.` g[2][1, 2] + 
       0.9` g[2][3, 2]) h[2][2, 1], -1.` g[2][3, 
      2] + (0.` - 0.11111111111111109` g[1][1, 1] + 
       0.09999999999999998` g[1][3, 1] - 1.` g[2][1, 1] + 
       0.9` g[2][3, 1]) h[2][1, 
      2] + (1.0999999999999999` - 0.11111111111111109` g[1][1, 2] + 
       0.09999999999999998` g[1][3, 2] - 1.` g[2][1, 2] + 
       0.9` g[2][3, 2]) h[2][2, 2], -178.88888888888889` g[2][1, 1] - 
    9.` g[2][3, 
      1] + (17.709999999999994` g[1][1, 1] + 159.39` g[2][1, 1]) h[2][
      1, 1] + (17.709999999999994` g[1][1, 2] + 159.39` g[2][1, 2]) h[
       2][2, 1], -178.88888888888889` g[2][1, 2] - 
    9.` g[2][3, 
      2] + (17.709999999999994` g[1][1, 1] + 159.39` g[2][1, 1]) h[2][
      1, 2] + (17.709999999999994` g[1][1, 2] + 159.39` g[2][1, 2]) h[
       2][2, 2], 
   0.` - 0.14799999999999996` g[2][1, 1] + 
    0.99` h[2][2, 1], -0.792` - 0.14799999999999996` g[2][1, 2] + 
    0.99` h[2][2, 2], 0.` + 1.` g[2][2, 1] - 1.` h[2][1, 1], 
   0.` + 1.` g[2][2, 2] - 1.` h[2][1, 2]};
polys = eqns /. {0. -> 0};

but with exactly the same structure

Find the variables

vars = Variables[polys]

Deconstruct the polynomials (I know it is a quadratic system). There are 10 nonzero constant terms, 228 nonzero linear terms, and 324 nonzero 2nd-order terms

{m0, m1, m2} = CoefficientArrays[polys, vars]

Confirm the reconstructed polynomials are the same as the originals (just checking)

polys - (m0 + m1.vars + m2.vars.vars) // Expand // Chop 

Of the 324 nonzero 2nd-order terms, of which there are 180 distinct cross-products ...

Flatten[Most[ArrayRules[#]] & /@ m2][[All, 1]] // Union // Length

... none of which are purely quadratic

Diagonal /@ m2 // Union

Mathematica cannot handle it. I am not sure whether it is a problem of memory, a problem of speed or something like that. I wounder whether I could take advantage of the quadratic structure without purely quadratic components to simplify the problem and make it possible for Mathematica to handle it.

Answer 1

"[O]nly a little"?? You just went from a system of 6 equations in 6 variables to a 20 x 20 system. I have to wonder what a big enlargement might be.

You might try

soln = NSolve[polys]

But that also could well hang.

A better possibility in terms of likelihood of completing might be to use local methods such as FindRoot, provided you have some idea of where the relevant (for your purposes) solutions might live.

Followup:

NSolve gave a result after a couple of hours.

In[13]:= Timing[soln = NSolve[polys];]

NSolve::sfail: Subsystem could not be solved for 
    152533 g[1][1, 1]   14327 g[1][1, 2]   171802 g[1][2, 1]
    ----------------- + ---------------- + ----------------- - 
         122535              17505              122535
     113492 g[1][2, 2]   24775 g[1][3, 1]   1475 g[1][3, 2]
     ----------------- + ---------------- - --------------- + <<11>> + 
          122535              24507              1167
     38732 h[2][1, 2]   22301 h[2][2, 1]   145171 h[2][2, 2]
     ---------------- - ---------------- + ----------------- at value 
          24507              17505              122535
                                        17        -20
    6.6504175107872416688377351747276 10   + 0. 10    I
    . The likely cause is failure to detect zero due to low precision. The
     likely effect is the loss of one or more solutions. Increasing
     WorkingPrecision might prevent some solutions from being lost.

Out[13]= {7651.893734, Null}

In[14]:= soln//Length

Out[14]= 110

I believe some of the solutions are numeric gibberish. Others seem plausible though I've not tested them all. Here is one that appears to be reasonable.

{g[1][2, 1] -> 0.877177885546412, g[1][2, 2] -> -0.731787453222285, 
 g[2][2, 1] -> 1.514048096489015, g[2][2, 2] -> 4.056151212079698, 
 g[1][3, 1] -> -0.2520687271834896, g[2][3, 1] -> 0.19110251804564718, 
 g[1][3, 2] -> 5.450124929865966, g[2][3, 2] -> 5.157856942595859, 
 g[1][1, 1] -> -0.3067458899627449, g[2][1, 1] -> 0.11152252700132087, 
 g[1][1, 2] -> 1.4198327350185431, g[2][1, 2] -> 0.9644105684773703, 
 h[1][2, 1] -> -0.06568699865868881, h[1][2, 2] -> 1.1040449897211426, 
 h[2][2, 1] -> 0.016672054541611605, h[2][2, 2] -> 0.94417450922692, 
 h[1][1, 1] -> 0.877177885546412, h[1][1, 2] -> -0.731787453222285, 
 h[2][1, 1] -> 1.514048096489015, h[2][1, 2] -> 4.056151212079698}

Residuals are on the order of 10^(-14) and smaller.

Answer 2

The Buchberger algorithm (the main approach used in computing a GroebnerBasis) strongly depends on the number of variables. It has double exponential dependence. So, if in the first example you have only six variables, while in the latter one there are 20, in the worst case you could expect up to

  2^(2^20)/2^(2^6) // N

3.654379384547651*10^315633

worse timing.

There is however hope when dealing with different orderings and different options like for example GroebnerWalk.