Efficiently solve this system of quadratic equations?

Questions
1

Hello,

the problem i need to solve is to find ALL solutions of a system of quadratic equations that has the relatively simple form aijkxjxk - xi = 0, using the Einstein summation convention and with each index running from 1 to 24. The values of the coefficients aijk are known, most of them are zero, all others (about 15 per equation on average) are positive.

I have typed out the quadratic equations with all values already plugged in and left out all summands that were zero. Then I have tried using the “Solve”-command with a Working Precision of 10 and Cubics, Quartics and VerifySolutions set to False. But after about 5 days of calculation it ran out of RAM and terminated. The solutions should be real, but I did not specify that, because I feared that all the “ConditionalExpressions” generated for when a root exists would need even more memory. Is there a better way to solve my problem that uses less memory?

2

Posting the explicit example would help. Could provide it as a sparse array to save on space. Absent that one can only guess, and I will venture that maybe it really does need huge memory, and maybe NSolve will do better.

There might also be a more clever way using linear algebra.

3

Thanks for the reply.
This is how the expression shows up in my notebook:

{-x[1] + 6 x[1]^2 + 16 x[8]^2 + 2 x[14]^2 + 4 x[15]^2 + 20 x[16]^2 + 
  10 x[17]^2 + 10 x[17] x[24] + 5 x[24]^2 + x[28]^2 + 2 x[14] x[29] + 
  x[29]^2,
     -x[3] + 8 x[3]^2 + 4 x[10]^2 + 4 x[15]^2 + 
  20 x[21]^2 + 10 x[22]^2 + 10 x[21] x[26] + 5 x[26]^2,
     -x[4] + 24 x[4]^2 + 36 x[4] x[6] + 3 x[6]^2 + 4 x[11]^2 + 4 x[16]^2 + 
  4 x[21]^2 + 10 x[23]^2 + 2 x[21] x[26] + 4 x[23] x[27],
     -x[5] + 14 x[5]^2 + 28 x[5] x[7] + 3 x[7]^2 + 4 x[12]^2 + 8 x[12] x[13] + 
  4 x[17]^2 + 4 x[22]^2 + 20 x[23]^2 + 4 x[17] x[24] + 8 x[23] x[27],
     -x[6] + 6 x[4] x[6] + 9 x[6]^2 + x[26]^2 + 2 x[27]^2,
     -x[7] + 6 x[5] x[7] + 7 x[7]^2 + 2 x[13]^2 + 2 x[24]^2 + 4 x[27]^2,
     -x[8] + 6 x[1] x[8] + 7 x[8] x[14] + 4 x[10] x[15] + 20 x[11] x[16] + 
  10 x[12] x[17] + 2 x[13] x[17] + 2 x[12] x[24] + 
  x[13] x[24] + 6 x[8] x[28] + 7 x[8] x[29],
     -x[10] + 5 x[3] x[10] + x[10] x[14] + 6 x[8] x[15] + 
  4 x[10] x[15] + 20 x[11] x[21] + 10 x[12] x[22] + 2 x[13] x[22] + 
  4 x[11] x[26] + 3 x[10] x[28] + 2 x[10] x[29] + 5 x[30]^2,
     -x[11] + 21 x[4] x[11] + 9 x[6] x[11] + x[11] x[14] + 6 x[8] x[16] + 
  4 x[11] x[16] + 4 x[10] x[21] + 10 x[12] x[23] + 2 x[13] x[23] + 
  4 x[10] x[26] + 5 x[12] x[27] + x[13] x[27] + 3 x[11] x[28] + 
  2 x[11] x[29] + x[30]^2,
     -x[12] + 11 x[5] x[12] + 7 x[7] x[12] + 
  5 x[5] x[13] + x[7] x[13] + x[12] x[14] + 6 x[8] x[17] + 
  4 x[12] x[17] + 4 x[10] x[22] + 20 x[11] x[23] + 4 x[8] x[24] + 
  2 x[13] x[24] + 5 x[11] x[27] + 3 x[12] x[28] + 2 x[13] x[28] + 
  2 x[12] x[29] + x[13] x[29],
     -x[13] + x[5] x[13] + 5 x[7] x[13] + 
  x[13] x[14] + 4 x[13] x[17] + 2 x[8] x[24] + 4 x[12] x[24] + 
  4 x[13] x[24] + x[13] x[28] + 4 x[30]^2, 
     12 x[8]^2 - x[14] + 6 x[1] x[14] + 2 x[14]^2 + 4 x[15]^2 + 
  20 x[16]^2 + 10 x[17]^2 + 20 x[17] x[24] + x[28]^2 + 2 x[1] x[29] + 
  x[29]^2,
     6 x[8] x[10] + 2 x[10]^2 - x[15] + 3 x[1] x[15] + 
  5 x[3] x[15] + 2 x[14] x[15] + 2 x[15]^2 + 20 x[16] x[21] + 
  10 x[17] x[22] + 10 x[22] x[24] + 10 x[16] x[26] + x[15] x[29] + 
  5 x[30]^2, 
     6 x[8] x[11] + 2 x[11]^2 - x[16] + 3 x[1] x[16] + 21 x[4] x[16] + 
  9 x[6] x[16] + 2 x[14] x[16] + 2 x[16]^2 + 4 x[15] x[21] + 
  10 x[17] x[22] + 10 x[23] x[24] + 2 x[15] x[26] + 2 x[17] x[27] + 
  x[24] x[27] + x[16] x[29] + x[30]^2, 
     6 x[8] x[12] + 2 x[12]^2 + 4 x[8] x[13] + x[13]^2 - x[17] + 
  3 x[1] x[17] + 11 x[5] x[17] + 7 x[7] x[17] + 2 x[14] x[17] + 
  2 x[17]^2 + 4 x[15] x[22] + 20 x[16] x[23] + x[1] x[24] + 
  5 x[5] x[24] + x[7] x[24] + 2 x[14] x[24] + x[24]^2 + 
  4 x[16] x[27] + x[17] x[29],
     4 x[10] x[11] + 4 x[15] x[16] - x[21] + 5 x[3] x[21] + 
  21 x[4] x[21] + 9 x[6] x[21] + 2 x[21]^2 + 10 x[22] x[23] + 
  x[3] x[26] + 5 x[4] x[26] + x[6] x[26] + x[26]^2 + 2 x[22] x[27], 
     4 x[10] x[12] + 4 x[10] x[13] + 4 x[15] x[17] - x[22] + 
  5 x[3] x[22] + 11 x[5] x[22] + 7 x[7] x[22] + 2 x[22]^2 + 
  20 x[21] x[23] + 4 x[15] x[24] + 10 x[23] x[26] + 4 x[21] x[27] + 
  x[26] x[27] + 2 x[30]^2, 
     4 x[11] x[12] + 4 x[11] x[13] + 4 x[16] x[17] + 4 x[21] x[22] - 
  x[23] + 21 x[4] x[23] + 11 x[5] x[23] + 9 x[6] x[23] + 
  7 x[7] x[23] + 2 x[23]^2 + 4 x[16] x[24] + 2 x[22] x[26] + 
  4 x[4] x[27] + 2 x[5] x[27] + x[6] x[27] + x[7] x[27] + x[27]^2, 
     2 x[8] x[13] + 4 x[12] x[13] + 2 x[13]^2 - x[24] + x[1] x[24] + 
  x[5] x[24] + 5 x[7] x[24] + 4 x[17] x[24] + 2 x[24]^2 + 
  x[24] x[29] + 4 x[30]^2,
     -x[26] + x[3] x[26] + x[4] x[26] + 5 x[6] x[26] + 4 x[21] x[26] + 
  4 x[26]^2 + 4 x[30]^2,
     -x[27] + x[4] x[27] + x[5] x[27] + 4 x[6] x[27] + 2 x[7] x[27] + 
  4 x[23] x[27] + 5 x[27]^2 + 2 x[30]^2, 
     16 x[8]^2 + 4 x[10]^2 + 20 x[11]^2 + 10 x[12]^2 + 4 x[12] x[13] + 
  5 x[13]^2 - x[28] + 2 x[1] x[28] + 6 x[14] x[28] + 6 x[28] x[29], 
     16 x[8]^2 + 4 x[10]^2 + 20 x[11]^2 + 10 x[12]^2 + 4 x[12] x[13] + 
  5 x[13]^2 + 5 x[24]^2 + 4 x[28]^2 - x[29] + 2 x[1] x[29] + 
  6 x[14] x[29] + 2 x[29]^2,
     -x[30] + x[10] x[30] + x[11] x[30] + x[12] x[30] + 2 x[13] x[30] +
  2 x[15] x[30] + x[16] x[30] + 2 x[17] x[30] + 2 x[21] x[30] + 
  x[22] x[30] + 2 x[23] x[30] + 2 x[24] x[30] + 4 x[26] x[30] + 
  5 x[27] x[30]}

I then set this equal to 0 in the argument of “Solve”. Initially it was a System of 31 quadratic equations in 31 variables, but I already found a way to simplify it, which I applied using “Replace”, hence the indices don’t run from 1 to 24, but there’s 24 different values for them. Would it be more efficient to declare the functions in terms of a sparse array?

For now I will try using “NSolve”.

4

Quickly run your example with random initial guess and FindRoot, I obtained a set of none trivial solution:

捕获.PNG
捕获.PNG879×568 16.8 KB

findtoot.nb (38.8 KB)

5

Thanks, I already found several nontrivial solutions as well. The problem is that I need to find ALL solutions, and scanning a 24-dimensional parameter space with random starting points seems impractical, given that one could never be 100% sure to not have missed a single root.

More precisely the task is to find out whether there is any zero point where all eigenvalues of the Jacobi matrix have positive real part. But checking the latter condition first without plugging in values for the x’s to restrain the parameter space seems even more hopeless.

6

Hello C W,

apparently there are some inconsistencies in your contribution:

  • you claim that the index is running from 1 to 24, however, in the code you posted I can see indices up to 30

  • the code you posted is not an equation. Does it represent the “aijkxjxk” part of the equation?

Regards,
Michael

7

Hello Michael,

I think what I wrote right below the code sample box answers your questions pretty well.

I then set this equal to 0 in the argument of “Solve”. Initially it was a System of 31 quadratic equations in 31 variables, but I already found a way to simplify it, which I applied using “Replace”, hence the indices don’t run from 1 to 24, but there’s 24 different values for them.

So to get the system of quadratic equations the expression just needs to be set equal to 0.

Regards

8

Since this quite complicated system of equaitons, I suggest using Reduce on a subset of the listed equations in your problem first. For example, you will find some insights about how some variables are related to each other with equation 16 and 18 only:

Reduce[Thread[eq == ConstantArray[0, 24]][[{16, 18}]], varlist]

Then gradully build your solution toward the whole set.

9

Interesting suggestion! I hope I followed it correctly. I tried to apply it to a simpler system of 32 equations where i set all but 12 of the variables equal to either 0 or another variable and added the equation 0==0 to get a nice total count. Fist i generated one condition from each equation like so:

cond1 = Table[Reduce[betay[[i]] == 0, Table[y[j], {j, 12}]], {i, 32}];

Then I went on to iteratively combine two conditions to get a more strict condition:

cond2 = Table[Reduce[cond1[[2 i - 1]] && cond1[[2 i]],Table[y[j], {j, 12}]], {i, 16}];

But the third step already took more than 20 minutes so I aborted it. “Solve” was able to find the solutions within 7.5 seconds. Do I need to think more carefully about which conditions to combine?

10

“NSolve” has been running for a few days now as well. At one point it occupied 90% of my RAM. Now it’s back to 70%, but I doubt the task will be completed before running out of memory at some point. I’ve been looking around and I’ve found that the typical algorithm for finding complete sets of solutions of systems of polynomial equations is doubly exponential in complexity, making it seem entirely impossible to get there in a reasonable amount of time. I don’t know which algorithms exactly are implemented in Mathematica for this problem, but the runtime for my problem is probably not that much better. Can someone confirm this?

I’m thinking I should restrict myself to finding subsets of the whole solution set. If i do find a zero point where all eigenvalues of the Jacobi matrix have positive real part, my task is completed. But if I can’t get there after checking a few thousand of them I should probably give up on it. What would you say?

11

My suggestion would be classify the solutions based on how many zeros in x[i]. The more zeros, the easier to find with FindRoot. Based on the generic pattern of your equations,

a*x^2+x == 0

I feel most x’s would be quite small. The solution with least zeros so far I have found is this one:

{Subscript[x, 1]->0.049298,Subscript[x, 3]->0.115546,Subscript[x, 4]->-0.0294936,Subscript[x, 5]->-0.080845,Subscript[x, 6]->0.0596149,Subscript[x, 7]->0.0875315,Subscript[x, 8]->0.000925203,Subscript[x, 10]->0.00173587,Subscript[x, 11]->-0.0186874,Subscript[x, 12]->-0.0164931,Subscript[x, 13]->0,Subscript[x, 14]->0.049294,Subscript[x, 15]->0.0465952,Subscript[x, 16]->0.00177485,Subscript[x, 17]->0.00182401,Subscript[x, 21]->0.00114009,Subscript[x, 22]->0.00127775,Subscript[x, 23]->0.00149794,Subscript[x, 24]->0,Subscript[x, 26]->0,Subscript[x, 27]->0.138165,Subscript[x, 28]->0.080892,Subscript[x, 29]->0.080892,Subscript[x, 30]->0}
12

It is probably not in the double-exponential category but it is still likely to be intractable.