In my case I assume $A,B,C,X \in \mathbb{C}^{n \times n}$ are $n \times n$ square matrices. The Sylvester equation (Wikipedia) is
$$
A X + B X = C
$$
and given $A,B,C$ we want to solve for $X \in \mathbb{C}^{n \times n}$.
The theoretical and naive approach for solving it is to use Kronecker product
$$
A \otimes B = (a_{i,j} B) \in \mathbb{C}^{n^2 \times n^2}
$$
and write the Sylvester equations as a set of linear equation by casting $X$ as a vector. Denote by $\operatorname{vec}(X)$ the operation of taking the columns of $X$ and stacking them vertically so to create a $n^2$ column vector. Then an equivalent equation to the Sylvester equation is
$$
(I \otimes A + B^T \otimes I) \operatorname{vec}(X) = \operatorname{vec}(C) .
$$
The computational complexity of solving it is $o(m^3) = o\left( (n^2)^3 \right) = o(n^6)$. If one uses Gauss elimination to a upper-triangular form and then backward substitutions, with $m = n^2$, then:
$$
\frac{2 m^3 + 3 m^2 – 5 m}{6} + \frac{m^2 + m}{2} = \frac{2m^3 + 6 m^2 – 2m}{6} \mbox{ MACS}
$$
and $m + (m-1)$ divisions (here $m = n^2$). MACs are multiplication, addition and accumulation action, each takes $o(1)$. This is very expansive in terms of computation time.
I guess there are more efficient algorithms to solve the Sylvester equation. I use the Python's SciPy function of $\texttt{solve_sylvester}$. According to SciPy reference:
Computes a solution to the Sylvester matrix equation via the Bartels- Stewart algorithm. The A and B matrices first undergo Schur decompositions. The resulting matrices are used to construct an alternative Sylvester equation (RY + YS^T = F) where the R and S matrices are in quasi-triangular form (or, when R, S or F are complex, triangular form). The simplified equation is then solved using *TRSYL from LAPACK directly.
The question is what is the computational complexity of solving Sylvester equation used by Python using Schur decomposition and the *TRSYL algorithm? Is it still $o(n^6)$ or is it $o(20n^3)$ (when $m=n$) as written here?
Is there a more exact expression? If possible I would like to separate multiplication and addition operations (MACs) from division and inversion operations in the counting.
Best Answer
Python uses Bartels-Stewart algorithm described in Wikipedia. First it uses real Schur decomposition to transform the equation into a triangular one, and the solves using forward or back substitution.
Using the QR algorithm, the real Schur decompositions in step 1 require approximately $ 10(m^{3}+n^{3}) $ flops, so that the overall computational cost is $$10(m^{3}+n^{3})+2.5(mn^{2}+nm^{2}) $$ and when $m=n$, the cost is $25 n^3$. This is quite expensive but far better than the naive $o((n^2)^3) = o(n^6)$.