Solving system of bilinear equations

Question

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.

Answer 1

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