Least squares problem with left and right unknowns For $i=1,...,n$, let $b_i$ be a scalar and $A_i$ be an $k\times l$ matrix. Is there a closed form solution for the following problem assuming $n>k+l$?
$$\min_{x\in \mathbb{R}^k ,y\in \mathbb{R}^l} \sum_{i=1}^n (b_i-x' A_i y)^2$$
 A: I don't think a closed-form solution will be forthcoming (since it's a nonlinear least squares problem). As an alternative, you could try a numerical solution.
For example, this Mathematica code gives a quick answer (I tried it on 100 random matrices of size 30x60).
k=30;l=60;n=100;
a=Table[RandomReal[],{i,1,k},{j,1,l},{p,1,n}];b=Table[RandomReal[],{i,1,n}];
xx=Array[x,k];yy=Array[y,l];zz=Flatten[{xx,yy}];
f[zz_]:=Sum[(b[[p]]-Sum[a[[i,j,p]]*xx[[i]]*yy[[j]],{i,1,k},{j,1,l}])^2,{p,1,n}]
NMinimize[f[zz],zz]

Prompted by Fedja's comment I experimented a bit. Mathematica offers a variety of numerical nonlinear global optimization methods. For a large number of variables (as in the example above), I found the Nelder-Mead method and the simulated annealing method to be most efficient. Both use random numbers to search for the global minimum, so one strategy to see how well they are doing is to run each multiple times with different random number seeds, and compare the minima.
My limited experience is that the same minimum does show up repeatedly, giving some confidence that it is a global minimum.
Do[Print[NMinimize[f[zz],zz,Method->{"SimulatedAnnealing","PerturbationScale"->7,"RandomSeed"->i}]], {i,5}]        
Do[Print[NMinimize[f[zz],zz,Method->{"NelderMead","ShrinkRatio"->0.85,"ContractRatio"->0.85,"ReflectRatio"->2,"RandomSeed"->i}]], {i,5}]

