Solve a system of 4-th order polynomial equations (numerically)

Question

I have a function defined on $ S^6 \times S^6 $ (two spheres embedded in $ \mathbb{R}^7 $ individually). Let us call this vector $ \vec{f}(\vec{x},\vec{y}) $, where $ \vec{x} $ and $ \vec{y} $ are points in the two spheres. The target space of the function if $ \mathbb{R}^{13} $, namely $ \vec{f}: S^6 \times S^6 \rightarrow \mathbb{R}^{13} $. The vectors $ \vec{x}, \vec{y} $ live in the spheres, so $ \vec{x} \cdot \vec{x} = 1 $ and $ \vec{y} \cdot \vec{y} = 1 $.

This is the problem I am trying to solve. I want to find at least two points $ {\vec{x}_0,\vec{y}_0} $ such that

$$ \vec{f}(\vec{x}_0, \vec{y}_0) = \vec{0} $$

Each entry $ \vec{f} $ is a polynomial of the 4-th order in the entries of the vectors and can be written formally as

$$ f^i(\vec{x}, \vec{y}) = \sum_{jknm} a^{i}_{jknm} x^j x^k y^n y^m $$ where $ a^{i}_{jknm} $s are real coefficients.

I have saved the vector $ \vec{f} $ in this file .wbx (download it) that you can import with the command

Uncompress @ Import["vector.wbx", "String"]

In this file, the entries of the vectors $ \vec{x} $ and $ \vec{y} $ are labelled with $ x_i $ and $ y_i $ for $ i = 1, 2, ..., 7 $.

Given the explicit values of the coefficients $ a^i_{jknm} $, I have tried to solve the system with Mathematica. Obviously, the system is very difficult and an analytic solution is hard to find. So, I have tried to find the solution numerically, but Mathematica takes too long and does not give an answer.

I have tried with

Solve[{f == Table[0, {i, 1, 13}, x.x == 1, y.y == 1}, {x, y} // Flatten]
NSolve[{f == Table[0, {i, 1, 13}, x.x == 1, y.y == 1}, {x , y} // Flatten]

and also I have looked if Mathematica can find just a solution with

FindInstance[{f == Table[0, {i, 1, 13}, x.x == 1, y.y == 1}, {x, y} // Flatten]

but still it takes too long.

Do you have any suggestion in order to solve this problem with Mathematica?

Answer 1

Perhaps the easiest solution, inferred from a comment by the OP below:

obj = Uncompress@Import["/tmp/vector.wdx", "String"];

obj /. {x3 -> 0, x4 -> 0, x5 -> 0, x6 -> 0, x7 -> 0,
        y3 -> 0, y4 -> 0, y5 -> 0, y6 -> 0, y7 -> 0}
(*  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}  *)

So x1, x2, y1, y2 may take on any values.

---

Original numerical approach

One could minimize the squares of the residuals. Use higher WorkingPrecision for a more accurate result.

obj = Uncompress@Import["/tmp/vector.wdx", "String"];
vars = Variables@obj;
sumOfSq = Flatten@{obj, #.# - 1 & /@ Partition[vars, 7]} // #.# &;

SeedRandom[0] (* for reproducibility *)

{acc, sol1} = FindMinimum[
  sumOfSq, Thread@{vars, Normalize@RandomReal[{-1, 1}, 14]}]
(*
{1.87883*10^-32,
 {x1 -> 0.770225, x2 -> 0.637772, x3 -> -1.07917*10^-10, x4 -> -1.40295*10^-9, 
  x5 -> 1.1682*10^-9, x6 -> 2.92351*10^-8, x7 -> -1.02405*10^-8,
  y1 -> 0.626869, y2 -> -0.779124, y3 -> 1.54472*10^-8, y4 -> -5.33267*10^-12, 
  y5 -> -1.04004*10^-9, y6 -> -5.10446*10^-10, y7 -> -2.39731*10^-10}}
*)

{acc, sol2} = FindMinimum[
  sumOfSq, Thread@{vars, Normalize@RandomReal[{-1, 1}, 14]}, 
  WorkingPrecision -> 32, MaxIterations -> 500]
(*
{7.6898897201977925978778380482460*10^-64,
 {x1 -> -0.90297104074904814666962995570081, 
  x2 -> 0.42970140745473574888183421474225, 
  ..., 
  y1 -> 0.65977018314385699231516512470880, 
  y2 -> 0.75146743471312275002694856178381, 
  ..., 
  y7 -> -2.6423411962015994155963564116132*10^-18}}
*)

Flatten@{obj, #.# - 1 & /@ Partition[vars, 7]} /. sol1 // Abs // Max
Flatten@{obj, #.# - 1 & /@ Partition[vars, 7]} /. sol2 // Abs // Max
(*
1.09623*10^-16
0.*10^-32
*)

Note that sol1 is pretty accurate at machine precision, just as accurate as sol2 when it is converted to machine precision:

Flatten@{obj, (#.# - 1) & /@ Partition[vars, 7]} /. N@sol2 // Abs // Max
(*  1.38778*10^-16  *)