# Solving $AXB + X\odot C = D$ matrix equation

Can anyone see a way to solve this equation efficiently?

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

I tried a straightforward solution that involved vectorizing $$X$$ but that turned out too expensive for my application -- my matrices are all $$d\times d$$ with $$d=1000$$ so I'm seeing if there's a way to solve it at a cost comparable to $$d\times d$$ SVD.

$$A$$,$$B$$ are positive-semi-definite and $$C$$ may have a few zeros.

• Have you tried solving the vectorized equation with an iterative method, such as GMRES? You can use the non-vectorized form to compute the matrix-vector product in \$O(d^3), so your total cost would be comparable to SVD if the number of iterations is small. You might need a preconditioner though. Sep 28, 2019 at 14:25
• @YaroslavBulatov I added "matrix-equations", correct me if I am wrong. Sep 29, 2019 at 19:37
• Sep 30, 2019 at 11:24
• Where it appears to have an answer. Please post on one site at a time only. Sep 30, 2019 at 23:46
• I'm starting to think that this problem is not solvable in O(d^3). It comes from approximating T=E[xi,xj,xk,xl] in terms of lower order moments, then trying to solve TX=G. I found that even computing ||T||^2 needs more than O(d^3) operations Oct 2, 2019 at 17:52

$$\def\BR#1{\Big(#1\Big)} \def\LR#1{\left(#1\right)} \def\KR#1{\left[#1\right]} \def\op#1{\operatorname{#1}} \def\trace#1{\op{Tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\c#1{\color{red}{#1}} \def\CLR#1{\c{\LR{#1}}} \def\fracLR#1#2{\KR{\frac{#1}{#2}}}$$Write the problem as a matrix function \eqalign{ F &= \BR{AXB + C\odot X - D} \;\doteq\; f(X) \\ } then use $$F$$ to create a scalar sum-of-squares objective function \eqalign{ \tfrac 12\,\|F\|^2 \;=\; \tfrac 12\,F:F \;\doteq\; \phi(X) \\ } where $$(:)$$ denotes the Frobenius product, which is a concise notation for \eqalign{ A:B &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij} \\ A:A &= \|A\|^2_F \qquad\big({\rm Frobenius\:norm}\big) \\ } Calculate the gradient of $$\phi$$ \eqalign{ d\phi &= F:dF \\ &= F:\LR{A\,dX\,B + C\odot dX} \\ &= \LR{A^TFB^T + C\odot F}:dX \\ \grad{\phi}{X} &= \LR{A^TFB^T + C\odot F} \;\doteq\; g(X) \\ } Use this gradient expression in your favorite gradient descent algorithm. I prefer Barzilai-Borwein for its speed and simplicity.
Initialize with the starting guess \eqalign{ X_0 &= random \\ G_0 &= g(X_0) \\ \phi_0 &= \tfrac 12\,\|f(X_0)\|^2 \\ X_1 &= X_0 - \fracLR{0.05\cdot\phi_0}{G_0:G_0}G_0 \qquad\qquad\qquad \\ k &= 1 \\ } and iterate until the convergence is satisfactory \eqalign{ G_k &= g(X_k) \\ X_{k+1} &= X_k - \fracLR{\LR{X_k-X_{k-1}}:\LR{G_k-G_{k-1}}} {\LR{G_k-G_{k-1}}:\LR{G_k-G_{k-1}}} G_k \\ k &= k+1 \\ }