# 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.

• Can you clarify what $\odot$ means? – Steven Landsburg Sep 28 '19 at 3:26
• I think what to do depends on the data $A, B, C$. If $A,B$ are invertible and if $\|A^{-1}\| \|B^{-1}\|\|C\|<1$ (say for the Frobenius norms) then $X\mapsto X + A^{-1}(X\odot C)B^{-1}$ is invertible by a Neumann series, with the usual estimate for the remainder. On the other hand, if all coefficients of $C$ are non-zero, and $\|A \| \|B \|\|C^{-\odot}\|<1$ (here $C^{-\odot}$ denotes the Hadamard inverse) one can invert instead $X\mapsto X + (AXB)\odot C^{-\odot}$. In both cases the above conditions on the norms are of course stronger than needed for the convergence. – Pietro Majer Sep 28 '19 at 7:45
• So I think a relevant point is: are your $A,B,C,D$ totally generic $d\times d$ matrices, or do they have any special feature that may address to a convenient method? – Pietro Majer Sep 28 '19 at 7:56
• I would try some iterative method like conjugated gradient possible with a preconditioner but the choice depends on $A,B,C$ – user35593 Sep 28 '19 at 10:58
• – Federico Poloni Sep 30 '19 at 11:25

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.

• This would just be a preconditioner though, right? – Nick Alger Oct 1 '19 at 3:34
• @NickAlger Why do you say so? – yarchik Oct 1 '19 at 7:48
• It seems to me that only $U \approx V^T\otimes W$, and they cannot be made equal in general (unless I am misunderstanding something). For example imagine $A=B=I$, and $C$ is a matrix filled with iid gaussian random entries. Then $U$ is diagonal, so $V$ and $W$ must be diagonal. But then equality is not possible, by a degree of freedom counting argument ($2d$ vs $d^2$). – Nick Alger Oct 1 '19 at 8:14
• @NickAlger Yes, you are right. – yarchik Oct 1 '19 at 8:46