how to solve this multivariate quadratic equation? Any hope to acquire an analytic solution to such equations:
Solve:
$$\sum_{j=1}^n a_{ij} x_i x_j = b_i$$
for $i=1,\ldots,n$, where $a_{ij}$'s and $b_i$'s are known constants and $x_i$'s are unknowns to be solved.
Thanks a lot!
P.S. Thanks for alex.jordan's interesting comment! Let's consider this problem in a positive setting, i.e. let's require all coefficients($a_{ij}$'s and $b_i$'s) to be positive so that the solution seems to exist. Also $a$ is symmetric, i.e. $a_{ij}=a_{ji}$. If there could be any fast numerical solution it's also useful.
P.S. Just for the information of all readers: the original form of this problem is formulated as:
$$ x_i\cdot\left(\sum_{j=1}^n a_{ij}x_j\right) = \sum_{j=1}^n a_{ij}q_{ij} $$ where $q_{ij}\in(0, 1)$. I simplified the RHS into $b_{ij}$ but now it seems better not to do so(sorry for that!) to at least give the existence of real solutions a better chance.
 A: I couldn't fit this into a comment - it's not an answer. Each equation has a different quadratic form on its left. With $n=3$, here is what the system can be made to look like, just for example:
$$
\begin{align}
\vec{x}^T\begin{bmatrix}a_{11}&a_{12}/2&a_{13}/2\\a_{12}/2&0&0\\a_{13}/2&0&0\end{bmatrix}\vec{x}&=b_1\\
\vec{x}^T\begin{bmatrix}0&a_{21}/2&0\\a_{21}/2&a_{22}&a_{23}/2\\0&a_{23}/2&0\end{bmatrix}\vec{x}&=b_2\\
\vec{x}^T\begin{bmatrix}0&0&a_{31}/2\\0&0&a_{32}/2\\a_{31}/2&a_{32}/2&a_{33}\end{bmatrix}\vec{x}&=b_3\\
\end{align}
$$
In each case because the matrices are symmetric, the matrix is (orthogonally) diagonalizable. And generally they have rank $2$ with a null-space of dimension $n-2$. Maybe this can be exploited. It looks like in general, there will be $2^n$ solutions counting multiplicity, but they may be complex solutions.
A: As "everything in sight is positive" write your system in the form
$$\sum_{j=1}^n a_{ij} x_j={b_i\over x_i}\qquad(1\leq i\leq n)$$
and use Newton's method: Given an approximate solution $x=(x_1,\ldots, x_n)$ find a better approximation $x+\Delta$ by solving
$$\sum_{j=1}^n a_{ij} (x_j+\Delta_j)={b_i\over x_i+\Delta_i}\doteq{b_i\over x_i}\left(1-{\Delta_i\over x_i}\right)\qquad(1\leq i\leq n) $$
for $\Delta$. This means solving the linear system
$$\sum_{j=1}^n \left(a_{ij}+\delta_{ij}{b_i\over x_i^2}\right)\>\Delta_j={b_i\over x_i} -\sum_{j=1}^n a_{ij}x_j\qquad(1\leq i\leq n)$$
for the $\Delta_j$.
A: Let $A=\left( a_{ij} \right)_{n \times n} $, $B=\mathrm{diag}\{b_1, b_2, \cdots, b_n\}=\begin{pmatrix} 
b_1 \\
& \ddots \\
& & b_n
\end{pmatrix}$, $C=A^{-1}B$.
There is a possible numerical solution when $A$ is invertible and $b_i\neq 0$ ($i$ = 1, 2, $\cdots$, $n$ ).
From the conditions, it will result
$$ 
\begin{pmatrix} 
x_1 a_{11} & \cdots & x_1 a_{1n} \\
\vdots & \ddots & \vdots \\
x_n a_{n1} & \cdots & x_n a_{nn}
\end{pmatrix}  
\begin{pmatrix} 
x_1 \\
\vdots \\
x_n
\end{pmatrix}  =  
\begin{pmatrix} 
b_1 \\
\vdots \\
b_n
\end{pmatrix}. 
 $$
or
$$ \begin{pmatrix} 
x_1 \\
& \ddots \\
& & x_n
\end{pmatrix}  
\begin{pmatrix} 
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn}
\end{pmatrix}  
\begin{pmatrix} 
x_1 \\
\vdots \\
x_n
\end{pmatrix}  =  
\begin{pmatrix} 
b_1 \\
\vdots \\
b_n
\end{pmatrix}. 
 $$
As $\sum\limits _{j=1} ^n a_{ij} x_i x_j=b_i$, i.e., $x_i \sum\limits _{j=1} ^n a_{ij}  x_j=b_i$, since $b_i \neq 0$, then $x_i \neq 0$.
so 
$$  
\begin{pmatrix} 
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn}
\end{pmatrix}  
\begin{pmatrix} 
x_1 \\
\vdots \\
x_n
\end{pmatrix}  =  
\begin{pmatrix} 
x_1^{-1} \\
& \ddots \\
& & x_n^{-1}
\end{pmatrix} 
\begin{pmatrix} 
b_1 \\
\vdots \\
b_n
\end{pmatrix} =  
\begin{pmatrix} 
b_1 \\
& \ddots \\
& & b_n
\end{pmatrix} 
\begin{pmatrix} 
x_1 ^{-1} \\
\vdots \\
x_n^{-1}
\end{pmatrix}.  
 $$
As $A=\left( a_{ij} \right)_{n \times n} $, $B=\mathrm{diag}\{b_1, b_2, \cdots, b_n\}=\begin{pmatrix} 
b_1 \\
& \ddots \\
& & b_n
\end{pmatrix}$, $C=A^{-1}B$, then
$$ A 
\begin{pmatrix} 
x_1 \\
\vdots \\
x_n
\end{pmatrix}   =  
B
\begin{pmatrix} 
x_1 ^{-1} \\
\vdots \\
x_n^{-1}
\end{pmatrix},  
 $$
so
$$  
\begin{pmatrix} 
x_1 \\
\vdots \\
x_n
\end{pmatrix}   =  
A^{-1} B
\begin{pmatrix} 
x_1 ^{-1} \\
\vdots \\
x_n^{-1}
\end{pmatrix},  
 $$
that is to say, 
$$  
\begin{pmatrix} 
x_1 \\
\vdots \\
x_n
\end{pmatrix}   =  
C
\begin{pmatrix} 
x_1 ^{-1} \\
\vdots \\
x_n^{-1}
\end{pmatrix}.  
 $$
Let $  
\begin{pmatrix} 
x_1^{(0)} \\
\vdots \\
x_n^{(0)}
\end{pmatrix}   =  
\begin{pmatrix} 
1 \\
\vdots \\
1
\end{pmatrix}$,
 by the recurrence 
$$  
\begin{pmatrix} 
x_1^{(k+1)} \\
\vdots \\
x_n^{(k+1)}
\end{pmatrix}   =  
C
\begin{pmatrix} 
{\left(x_1^{(k)}\right)} ^{-1} \\
\vdots \\
{\left(x_n^{(k)}\right)}^{-1},
\end{pmatrix},  
 $$ 
we can get a numerical solution. 
If the result does not converge, we can use another one
$$ 
\begin{pmatrix} 
{\left(x_1^{(k+1)}\right)} ^{-1} \\
\vdots \\
{\left(x_n^{(k+1)}\right)}^{-1},
\end{pmatrix}   =  
C^{-1}
\begin{pmatrix} 
x_1^{(k)} \\
\vdots \\
x_n^{(k)}
\end{pmatrix}.  
 $$ 
If you want the analytic solution, you may have a try by using the Groebner-Shirshov Bases. (The software Maple or Mathematica can do it.) But I am not sure the analytic solution can be find easily.
