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

-
Just because all of the coefficients are positive does not imply that all solutions will be real. Consider $n=2$, with no constraint on the coefficients. Generally, there will be four complex points of intersection. But with all of the coefficients equal to 1 we have two real and two nonreal solutions. – alex.jordan May 23 '13 at 21:09
Yeah, that's true. – yuanz07 May 24 '13 at 1:05

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.

-

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

-
Thanks a lot for your very helpful comment! Any idea on the stability of this method? (Sorry that I don't know much about computational mathematics and hope it doesn't bother you.) – yuanz07 May 24 '13 at 1:17
I have thought about this for a while, however, it is stable when the solutions are not close to 0. And we also know the solution is not unique, then the final solution of iteration method maybe very unstable w/ resp to initial state. – Yimin May 24 '13 at 2:46

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.

-
The real case is when A is almost non-singular(having a large cond number) and the true solutions are close to 0, both of the schemes are not stable, the error will be surely diverging. And there seems no way to avoid this. – Yimin May 25 '13 at 18:12