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?


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?

  • $\begingroup$ This looks almost like a Stein equation, except you have "+X" instead of "-X" sciencedirect.com/topics/engineering/discrete-lyapunov-equation $\endgroup$ Commented Jul 12, 2022 at 22:37
  • 1
    $\begingroup$ That's not a problem, I can have $A'=-A$ or $B'=-B$ and $E'=-E$. I will look at it, as long as there is a way to factor the equation without inverting the matrix then it's better than what I have now! $\endgroup$
    – PC1
    Commented Jul 12, 2022 at 22:41

2 Answers 2


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

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

  2. 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$


  • $\begingroup$ I thought that the discrete Lyupanov equation is $AXA^T-X=C$. In my case, I don't have $A^T=B$. $\endgroup$
    – PC1
    Commented Jul 13, 2022 at 3:54
  • 2
    $\begingroup$ Is called Stein equation, equation in Eq 6.1 of dm.unibo.it/~simoncin/matrixeq.pdf , but Mathematica calls is Discrete Lyapunov $\endgroup$ Commented Jul 13, 2022 at 4:09
  • $\begingroup$ Thanks for the link. I found in the meantime a "quasi-geometric" solution to the equation but it requires to compute an infinite sum. Is there a simple way to compute this? I added an edit in the question. $\endgroup$
    – PC1
    Commented Jul 13, 2022 at 4:16
  • $\begingroup$ I gave a simple way to compute it in the answer $\endgroup$ Commented Jul 13, 2022 at 4:28
  • $\begingroup$ and here's an implementation that verifies that both methods work -- wolframcloud.com/obj/yaroslavvb/newton/forum-stein-lyapunov.nb $\endgroup$ Commented Jul 13, 2022 at 4:29

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.

  • $\begingroup$ Are there dedicated solvers which work for underconstrained Lyapunov/Sylvester? $\endgroup$ Commented Jul 13, 2022 at 15:23
  • $\begingroup$ Not that I know, just square systems. $\endgroup$ Commented Jul 13, 2022 at 15:28
  • $\begingroup$ Well I get this for square, but the problem becomes underconstrained when A is numerically singular, which is the case for most large covariance matrices $\endgroup$ Commented Jul 13, 2022 at 17:13
  • $\begingroup$ If the coefficients are covariance matrices then they are symmetric, so essentially all these solvers reduce to the closed formula with eigenvalues that you have posted above. With that you can implement all sorts of regularization tricks for cheap: Tikhonov, truncating small eigenvalues, pseudoinverses with 1/0 = 0. $\endgroup$ Commented Jul 13, 2022 at 17:28
  • 2
    $\begingroup$ @YaroslavBulatov Stability. Unsymmetric eigenvalue decompositions can be arbitrarily unstable. If you try an example with ill-conditioned $T$ and $U$ you'll see that you cannot compute a solution with small residual. $\endgroup$ Commented Jul 14, 2022 at 15:48

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.