Solving $AXB + X\odot C = D$ 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.
 A: I am going to reduce your problem to another form. In a truncated version, this modified problem has been discussed here and it seems there is no apparent efficient solution à la Bartels–Stewart. However, it still might be valuable to have an alternative approach at hand, especially because such equations are being actively discussed (see T. Damm, Direct methods and ADI-preconditioned Krylov subspace methods for generalized Lyapunov equations, Numer. Linear Algebra Appl. 15, no. 9 (2008): 853-871) and because iterative approaches based on the standard Lyapunov solvers do exist.
Let us apply the vec-operation to the both sides of your equation using the commutativity of Hadamard product
$$
\mathrm{vec}(AXB)+\mathrm{vec}(C\odot X)=\mathrm{vec}(D).
$$
Now we use the following properties
$$
\text{vec}(ABC)=(C^\mathrm{T}\otimes A)\text{vec}(B)
$$
and 
$$
\text{vec}(A\odot B)=\text{vec}(A)\odot\text{vec}(B).
$$
Thus
$$
(B^\mathrm{T}\otimes A)\text{vec}(X)+\mathrm{diag}\!\left[\mathrm{vec} (C)\right]\mathrm{vec} (X)=\mathrm{vec}(D).
$$
Let us introduce a new matrix $U$:
$$
U=(B^\mathrm{T}\otimes A)+\mathrm{diag}\!\left[\mathrm{vec} (C)\right].
$$
Our goal is now to write it in the form of a Kronecker product.
$$
U=\sum_i \sigma_i V_i^\mathrm{T}\otimes W_i.
$$ 
This is known as the nearest Kronecker product problem. Using SVD decomposition on a permuted version of $U$ as Van Loan (J. Comp. Appl. Math, 123 (2000) 85) proposed, the factors $V_i$ and $W_i$ and the singular values $\sigma_i$ can be determined, and the original equation can be written as
$$
\sum_i \sigma_i\left(V_i^\mathrm{T}\otimes W_i\right)\mathrm{vec} (X)=\mathrm{vec}(D),
$$
which is equivalent to
$$
\sum_i \sigma_i W_i X V_i=D.
$$
Truncating the sum to just 2 terms, a standard Lyapunov equation is obtained. 
