I need to approximately solve the following underdetermined system of $ n$ linear equations

$ $ y_i=a_i^T X b_i$ $

Where $ X$ is $ d\times d$ unknown matrix, $ a_i$ and $ b_i$ are given vectors and $ n=d$ . Is there a way to do this much faster than vectorizing both sides of each equation and calling `LinearSolve`

?

LinearSolve approach destroys the structure, it reduces to solving a system with $ d^3$ coefficients instead of original one with $ 2d^2$ coefficients.

Below is an example of this approach, it’s too slow for my application where $ d=1000$ . On that scale, `backward`

is 2000x slower than `forward`

, with majority of the time spent in `LinearSolve`

. I was hoping for 100x slower since that seems like typical `LinearSolve`

overhead for unstructured systems on that scale.

`n = d = 50; {a, b} = Table[RandomReal[{-1, 1}, {n, d}], {2}]; X = RandomReal[{-1, 1}, {d, d}]; forward[X_] := MapThread[#1.X.#2 &, {a, b}]; y = forward[X]; backward[Y_] := With[{mat = MapThread[Flatten@Outer[Times, #1, #2] &, {a, b}]}, x = LinearSolve[mat, y]; ArrayReshape[x, {d, d}] ]; Print["forward time is ", forward[X]; // Timing // First]; {timing, error} = Timing[Norm[forward[backward[y]] - y, Infinity]]; Print["backward time is ", timing]; Print["error is ", error] `