Efficient algorithm for matrix equation $AXB + BXA = F$ For $n\in\mathbb{N}$, let $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times n}$ be two symmetric positive definite matrices and let $F\in\mathbb{R}^{n\times n}$ be arbitrary.
Is there any efficient algorithm for solving the matrix system $AXB + BXA = F$ (with $X\in\mathbb{R}^{n\times n}$ being the unknown) provided it has a unique solution?
I have found similar questions on the site, but they involve more assumptions on $F$, such as symmetry.
Thank you very much.
 A: For dense problems, the standard algorithm is a generalization of the Bartels--Stewart algorithm: see for instance https://doi.org/10.1016/S0024-3795(87)90314-4 and https://people.cs.umu.se/isak/recsy/ for an implementation.
The basic idea of the algorithm is: there are (QZ decomposition) two orthogonal matrices $Q,Z$ such that $QAZ$ and $QBZ$ are both upper triangular (over $\mathbb{C}$, or block upper triangular with $2\times 2$ blocks on $\mathbb{R}$). Hence with a suitable change of basis one can reduce to the case in which $A$ and $B$ are upper triangular. At this point, the entries of $X$ can be computed with back-substitution. Checking complexities, one sees that all the various parts of this algorithm require only $O(n^3)$ floating-point operations. This compares very favorably with $O(n^6)$ for the "big linear system" approach.
Remarks:

*

*the typical case is the slightly more general problem $AXB+CXD=E$; in your case, since some of the coefficients coincide, you can use only one QZ decomposition instead of two, but you need to change the order in which you compute the entries of $X$ by back-substitution.

*You will need to save some intermediate values to compute the RHS of (13) in that paper in linear time; I did not see it mentioned there at a first glance but it is worth noting.

*This algorithm has good stability properties because it relies on orthogonal transformations and (backward stable) back-substitution, and, unlike the algorithm mentioned in the comments, it does not break down if $A$ is singular or ill-conditioned.

A: For $F$ symmetric (which I note you have already looked at) is a mild generalization of the Lyapunov equation, see, e.g. (5.2) of http://www.dm.unibo.it/~simoncin/matrixeq.pdf (published version is at https://doi.org/10.1137/130912839).
For generic $F$, and indeed more general forms of this equation altogether, see section 7 of the paper above.
A: You can treat this as a system of $n^2$ linear equations for the entries of
$X$.  The coefficient of $x_{ij}$ in the equation whose right side is $f_{kl}$ is $a_{ki} b_{jl} + b_{ki} a_{jl}$.  The only way I can think of to make this "efficient" is to introduce some sparseness by taking a basis in which, say, $A$ is diagonal:  then
the entry is $0$ unless $k=i$ or $j=l$.
