# Solve efficiently large system of $N$ quadratic equations

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:

• I think your $\lambda$-values are called $Z$-eigenvalues. Mathematica does not have built-in functionality to find them, as far as I know. Maybe your question would be better off at the math StackExchange. Mar 19 at 18:49
• By the way, even for fully symmetric $T$ there are complex eigenvalues $\lambda$. They appear along with real eigenvalue pairs $\pm \lambda$ and possess a similar symmetry $\pm \lambda, \pm \lambda^*$. With my method below, I can reliably get the lowest real eigenvalue. This approach can be generalized to all real $\lambda$. But I do not know how to get the remaining complex ones. Do you know, by the way, how many eigenvalues are expected for a given $N$? Mar 20 at 20:46
• The problem of finding all eigenvalues $\lambda$ for $N=250$ is likely to be unsolvable. There are $1.8\times10^{75}$ eigenvalue pairs. No chance even to store them all :) Mar 21 at 22:45
• Nice plot, your question has quite deep connections. Mar 25 at 16:16

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).

• Not so fast. All you've shown is that if we have found such a vector $\vec{x}$, then the way to compute the corresponding $\lambda$ is what you write. But the question of how to find the allowed vectors $\vec{x}$ remains open. Mar 19 at 18:17
• @Roman see my update. Thanks for being careful and providing a stimulating question ;) Mar 19 at 19:26
• Nice! The OP is looking for the complete spectrum though, not just for the "ground state". Mar 19 at 19:31
• @Roman yes, nice observation again! I will indicate in my answer that I only look for the ground state, the generalization to other states should be possible, but there could be some numerical issues with changing the eigenvalues order during the iteration. This is harder to fix. There is no universal solution. Mar 19 at 19:36
• Hi, this is a very interesting answer, thank you. The issue I have with your algorithm is that how can I be sure that I will reach the actual smallest real $\lambda$? When I compare my naive implementation with your algorithm for $n=5$, the converged $\lambda$ is equal to $\lambda_\text{min}$ only $\sim 75\%$ of the time.
– Matt
Mar 25 at 9:56