# Solving Kronecker-structured linear equations

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]