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.



Efficient algorithm for matrix equation $AXB + BXA = F$
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:


1. 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.
2. 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.
3. 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.


