Solution to Sylvester-like equation

Question

I have an equation that is a bit similar to a Sylvester equation.

The equations is $AXB^T+X=E$, where all variables are matrices.

I could try to inverse $B$ and rewrite the equation as $AX+XB^{-T}=EB^{-T}$, which is then the standard form for the Sylvester equation, but I need to invert $B$, which is usually not a good idea.

Is there a more standard approach to work with this type of equation?

Edit:

I can see by doing simple algebra that $X=E-AEB^T+A^2EB^{2T}-A^3EB^{3T}+\dots$ seems to be a solution of the equation $AXB^T+X=E$, if the term $A^nEB^{nT}\to0$ as $n\to\infty$. Would there be a simple way to compute this infinite sum?

Edit 2:

I worked on some potential factorization and I realized that in my problem, $A$ and $B$ are actually symmetrical and positive-definite matrices. The matrix $E$ is not symmetrical (and I don't expect the solution $X$ to be symmetrical either).

So with this information, I rewrite the matrices as $A=P_AD_AP_A^T$ and $B=P_BD_BP_B^T$, which means: $$P_AD_AP_A^TXP_BD_BP_B^T+X=E$$ $$D_AYD_B+Y=F$$ with $Y=P_A^TXP_B$ and $F=P_A^TEP_B$. So the matrix $Y$ can be computed as: $$y_{ij}(1+a_ib_j)=f_{ij}$$ where $a_i$ and $b_j$ are the diagonal entries of $D_A$ and $D_B$ respectively. This equation has stable solutions as $a_i>0$ and $b_j>0$. Then $X=P_AYP_B^T$. So it requires 2 diagonalisations and 4 matrix multiplications, plus one scalar equation per entry.

Is there another matrix factorization approach that could be more efficient here?

Answer 1

Working with a more standard form of a Stein equation, sometimes called (generalized?) Discrete Lyapunov Equation

$$AXB - X = C$$

1. Rewrite it in standard linear equation form, use least squares solver, reshape $X$ back into matrix.

$$(B^T\otimes A - I\otimes I)\text{vec}X=\text{vec}C$$

2. Use dedicated solver for Stein equation. They tend to be robust but slow compared to eigenvector solution.

3. Write solution in terms of eigenvalues/eigenvectors.

$$X=T\left( \frac{T^{-1} C B^{-1} U}{\lambda(A)-\lambda(B^{-1})^T}\right)U^{-1}$$ where

• $\lambda(Y)$ is a column vector of eigenvalues of $Y$
• $T$ is a matrix of eigenvectors of $A$ as columns
• $U$ is a matrix of eigenvectors of $B^{-1}$ as columns
• $X/Y$ is pointwise division
• $-$ is subtraction using numpy broadcasting rules

In the underconstrained case, you will have some zeros in denominator. You can simply skip dividing by those entries. I haven't seen this approach in literature, but in experiments this gave me the same solution as the expensive "least squares" form with kronecker products, but much cheaper.

Also note that there are no standalone inverses. Eigenvalues/eigenvectors of $B^{-1}$ can be computed from $B$. Everything else is of the form $A^{-1}B$ or $AB^{-1}$ which is equivalent to solving $AX=B$ or $XA=B$

notebook

Answer 2

EDIT: this answer is essentially useless after the question was updated. If the coefficients $A,B$ are symmetric, these approaches reduce to the simpler closed formula with the eigenvalue decompositions that has already been suggested elsewhere. I don't think there is a simpler / more efficient solution than that.

---

This answer expands on the "dedicated solvers" mentioned in Yaroslav Bulatov's answer.

For these equations there is an algorithm based on reduction to triangular form via Schur factorizations + back-substitution, analogous to the Bartels-Stewart algorithm for Sylvester equations. It is rarely spelled out explicitly in books and sources, but it follows from the same ideas.

The algorithm is "unstructured backward stable", in that it guarantees a low residual $\|AXB -X - C\|$; this follows, for instance, from results for a much more general equation in this paper of mine, but it was surely known before that.

The back-substitution being a for-loop, it's better left to high-performance Fortran or C code. Two packages that implement it are RECSY and MEPACK, but both seem hard to find. In addition, there is the difficulty of knowing how to link Fortran libraries to whatever language you are using. You can try writing to the authors of MEPACK if you want to go this route. It might be easier to find code for the more general equation $AXB+CXD=E$.

EDIT: also SLICOT, with a direct download link to the code.