Consider a collection of $m$ matrices $A\_i$ of size $n\times n$, and a vector $b$ of size $m$. I want to solve the bilinear system


$$\left\{ x^T A\_i y = b\_i : i = 1,\dots,m \right\}$$


in variables $x,y$. Is there an *efficient* way of doing this?


This is both a theoretical and a practical question: the matrices $A\_i$ I have in mind are sparse and have size in the thousands. I understand that if $m=1$ then one can consider a vector $z=[x y]$ and run a semidefinite solver on $\|z^T [0\; A\_1;A\_1^T\; 0]z-2b\_1\|$; however the $m$ I have in mind is also in the thousands.



Solving system of bilinear equations
We have a system of $m$ bilinear equations in $\mathrm x, \mathrm y \in \mathbb R^n$


$$\begin{aligned} \mathrm x^\top \mathrm A\_1 \,\mathrm y &= b\_1\\ \mathrm x^\top \mathrm A\_2 \,\mathrm y &= b\_2\\ &\vdots\\ \mathrm x^\top \mathrm A\_m \,\mathrm y &= b\_m\end{aligned}$$


If matrices $\mathrm A\_1, \mathrm A\_2, \dots, \mathrm A\_m$ are very sparse, perhaps it would not be utterly hopeless to use symbolic methods (e.g., [Gröbner bases](https://en.wikipedia.org/wiki/Gr%C3%B6bner_basis)). However, we can use numerical methods. Note that


$$b\_i = \mathrm x^\top \mathrm A\_i \,\mathrm y = \mbox{tr} \left( \mathrm x^\top \mathrm A\_i \,\mathrm y \right) = \mbox{tr} \left( \mathrm y^\top \mathrm A\_i^\top \,\mathrm x \right) = \mbox{tr} \left( \mathrm A\_i^\top \mathrm x \mathrm y^\top \right) =: \langle \mathrm A\_i, \mathrm x \mathrm y^\top\rangle$$


Let $\mathrm Z := \mathrm x \mathrm y^\top$. Hence, we have $m$ linear equality constraints in $\rm Z$ and a constraint on its rank 


$$\begin{aligned} \langle \mathrm A\_1, \mathrm Z \rangle &= b\_1\\ \langle \mathrm A\_2, \mathrm Z\rangle &= b\_2\\ &\vdots\\ \langle \mathrm A\_m, \mathrm Z\rangle &= b\_m\\ \mbox{rank} (\mathrm Z) &= 1\end{aligned}$$


Since the nuclear norm is a convex proxy for the rank, we solve the following convex program in $\rm Z$


$$\begin{array}{ll} \text{minimize} & \| \mathrm Z \|\_\*\\ \text{subject to} & \langle \mathrm A\_1, \mathrm Z \rangle = b\_1\\ & \langle \mathrm A\_2, \mathrm Z \rangle = b\_2\\ & \qquad\vdots\\ & \langle \mathrm A\_m, \mathrm Z \rangle = b\_m\end{array}$$


If the optimal solution is rank-$1$, then we have solved the original system of $m$ bilinear equations.


