I need to solve the following equation for $X$ with $d$-by-$d$ matrices $A,B,C,D$ and Hadamard product $\odot$

$$AXB + X\odot C = D$$

Vectorizing all terms gives a solution with $O(d^6)$ complexity, which is intractable since $d\approx 1000$ in my application. Is there something I can do to get an estimate in $O(d^3)$ time?

To add more information about the structure, A,B,C are moment matrices. Specifically, for random variables $X,Y$ they are

$$A_{ij}=E[X_iX_j]$$ $$B_{ij}=E[Y_iY_j]$$ $$C_{ij}=E[X_iY_j]$$ $$E_t[D_{ij}]=C_{ij}$$

To give even more background, this comes up in problem of speeding up neural network training which for a single layer can be viewed as the following problem

$$\text{minimize}_{W} E[(x'Wy)^2]$$

Here $x$ and $y$ are random variables with shape $(d,1)$, the gradient is $E[x_iy_j]$ and the curvature is $E[x_i x_j y_k y_l]$. The goal is to use curvature information to obtain Newton-like correction to a noisy estimate of the gradient. Applying Isserlis theorem we can approximate curvature rank-4 tensor in terms of rank-2 covariance tensors which leads to matrices $A,B,C$ above. Substituting small sample estimate of gradient into $D$ and large sample estimate of curvature into $A,B,C$, then solving for $X$, gives us a preconditioned gradient step.

Incidentally, using the same sample to estimate curvature and gradient, gives slightly simplified problem:

$$AXB + X\odot C = C$$

$A$ and $B$ are known to be ill-conditioned -- $n$-th eigenvalue is approximately $1/n$ and majority of eigenvalues are numerically 0, possibly 90% of all eigenvalues. Coincidentally, this ill-conditioning is what allows neural networks to generalize. This implies that $C$ is also singular.