Solving $AXB + X\odot C = D$

Question

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.

Answer 1

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.