Solve efficiently large system of $N$ quadratic equations

Question

I am trying to find all $\lambda$ given a tensor $T$ such that $\sum_{jk}^N T_{ijk}x_jx_k=\lambda x_i$ with $1\leq i \leq N$ and $\sum_i x_i^2=1$. $T$ is known but $x$ is unknown.

This boils down to solve $N+1$ quadratic equations for $N+1$ variables.

My $T_{ijk}$ are random and fully symmetric. I tried to do a naive implementation but already for $N=5$ my laptop is giving me a hard time. I was hoping to solve the system for at least $N=250$. How could I improve my solving method to gain some time?

(* Generate my symmetric random tensor T*)
RandomSymmetrizedArray[dims_, sym_, dist_] := 
 Normal@SymmetrizedArray[_ :> RandomVariate[dist], dims, sym]

T[l_] := T[l] = 
  RandomSymmetrizedArray[{l, l, l}, Symmetric[All], 
   NormalDistribution[0, 1/l]]
X[l_] := Table[Subscript[x, i], {i, 1, l}]


l := 3
variables := Join[X[l], {\[Lambda]}]
S = NSolve[T[l].X[l].X[l] == \[Lambda]*X[l] && X[l].X[l] == 1, 
   variables];
S // MatrixForm

I am not interested in the values of $x_i$, only in $\lambda$. I also noticed that for every $\lambda_{sol}$, $-\lambda_{sol}$ will be a solution too.

Any remark, or advice is always appreciated, thanks.

Edit: Looking for the smallest real $\lambda_\text{min}$ is very interesting, and since if $-\lambda$ is a solution then $\lambda$ will also be a solution (associated with $-\vec{x}$). However I would still be interested in improving my method to solve this system of quadratic equations. Although $N=200$ would give me an absurd number of $\lambda$ (Around $2^200\approx 10^60$) I am still hoping to be able to compute this for $N=25$ for example.

If I plot my solutions $\lambda$ on the complex plane for $N=15$ here is what I find when $T$ is random:

Answer 1

Well, there is an analytic solution to your problem :)

$$\left\{\begin{array}{c}\sum_{jk}^N T_{ijk}x_jx_k=\lambda x_i,\\ \sum_{i}^N x_i^2=1,\end{array} \right.\Rightarrow \sum_{ijk}^N T_{ijk}x_i x_jx_k=\lambda.$$

This means that for any given vector $x_i$, there exist $\lambda$ that can be determined by the equation above.

In the comment below Roman is right that the solution is not complete unless we give a prescription on how to determine the vector $x$. Below I give such solution. Let me first mention that problems of this kind are very typical for the electronic structure calculation in physics or quantum chemistry. There it goes under the name self-consistent field (SCF) method.

The idea of SCF method

1. Start from a given vector $\vec x$,

2. Determine a new vector $\vec x'$ and corresponding $\lambda$

3. If $||\vec x'-\vec x||<10^{-\epsilon}$ the solution has been found, otherwise update the vector $x'\rightarrow x$ according to some algorithm and go to step 2.

The only question is how to perform step 2, namely, to determine a new vector and a new value $\lambda$. In the given case it is easy because $T$ is fully symmetric tensor. Therefore $h[\vec x]= T\cdot \vec x$ is a symmetric matrix, which I will call Hamiltonian matrix according to physical interpretation. With this crucial observation $\vec x'$ and $\lambda$ are the eigenvectors and eigenvalues of the Hamiltonian. I will focus on the lowest eigen-pair (the ground state) because in this case the solution is very simple.

Now we can pack these ingredients into a MA algorithm

1. Initialize the vector $\vec x$ with $||\vec x||=1$

n=5;
x = RandomReal[{0, 1}, {n}]
x = x/Sqrt[x.x]

2. Do the refinement

tol = 10^-9;
y[0] = x;
Do[
 h = T[n].y[0];
 {eval, evec} = Eigensystem[h];
 ord = Ordering[eval] // First;
 λ = eval[[ord]];
 y[1] = evec[[ord]];
 y[1] = y[1]/Sqrt[y[1].y[1]];
 err = Norm[y[1] - y[0]];
 If[err < tol, Break[],
        y[0] = y[1];
        ], {cyc, 200}
 ]

3. Show the results

Print["convergence error: ",err];
Print["converged vector: ",y[0]];
Print["converged λ: ", λ]

Optimization

The naive algorithm presented above may not converge fast enough for larger matrices. The solution is to improve the update step. One can add the mixing as follows:

cos = y[1].y[0];
If[cos < 0., y[1] = -y[1]];
y[0] = 0.9 y[1] + 0.1 y[0]; 
y[0] = y[0]/Sqrt[y[0].y[0]]; 

With this modification I can converge $n=200$ problem within 792 cycles in 183 seconds (typical run).