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.

  • 3
    $\begingroup$ 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. $\endgroup$ Sep 28, 2019 at 14:25
  • 2
    $\begingroup$ @YaroslavBulatov I added "matrix-equations", correct me if I am wrong. $\endgroup$
    – Anton Menshov
    Sep 29, 2019 at 19:37
  • 1
    $\begingroup$ Cross-posted on MO. $\endgroup$ Sep 30, 2019 at 11:24
  • 1
    $\begingroup$ Where it appears to have an answer. Please post on one site at a time only. $\endgroup$
    – Richard
    Sep 30, 2019 at 23:46
  • 1
    $\begingroup$ 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 $\endgroup$ Oct 2, 2019 at 17:52

1 Answer 1


$ \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 \\ }$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.