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_mend{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^toprangle$$

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 Zrangle &= b_2\ &vdots\ langle mathrm A_m, mathrm Zrangle &= b_m\ mbox{rank} (mathrm Z) &= 1end{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\ & qquadvdots\ & langle mathrm A_m, mathrm Z rangle = b_mend{array}$$

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