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]