Is there any efficient solution of the matrix equation AXB + (AXB)' + cX = D

Question

I am trying to find the symmetric solution $X\in \mathbb{R}^{p\times p}$ of following matrix equation:

$AXB + (AXB)^T + cX = D$

where $A,B\in \mathbb{R}^{p\times p}$ are symmetric positive definite matrix, and $D\in \mathbb{R}^{p\times p}$ will also be symmetric.

e.g.

1. if $c=0$ and one of $A$ or $B$ equals identity matrix $I$, we call this equation Sylvester Equation, and have closed form solution using SVD.

2. The existence and uniqueness of the solution can be proved under the condition that $\beta_i\alpha_j + \alpha_i\beta_j + c \neq 0$ for all $i,j=1,\cdots ,p$ where $\alpha_i$ and $\beta_i$ are eigenvalue of $A,B$ respectively.

Currently, I can't find closed form solution of this equation when $A,B,D$ are arbitrary matrix. It seems like the only way to solve it is transforming this equation into linear equation using kronecker product

$(B\otimes A + A\otimes B + cI)Vec(X) = Vec(D)$

and solve this problem using conjugate gradient descent naively. but it's computational unafforfable when the size $p$ is relative large, since we need to solve a $p^2 \times p^2$ problem.

What if we take advantage of the special structure of the matrix? e.g. $X$ is a symmetric matrix and $A, B$ are symmetric positive definite matrix.

Can someone figure it out? Thanks!

Answer 1

Remark 1. The eigenvalues of $A\otimes B +B\otimes A +cI$ are not the $\alpha_i\beta_j+\beta_i\alpha_j+c$ except when (for example) $AB=BA$.

Remark 2. The complexity of your problem is not in $O(p^6)$. Indeed, there are special methods to solve equations in the form $AXB+CXD=E$ where $A,C\in M_{m,m},B,D\in M_{n,n},E\in M_{m,n}$ are known and $X\in M_{m,n}$ is unknown. cf.

i) http://www.maths.lth.se/na/courses/NUM115/NUM115-09/sylvester.pdf

uses an algorithm that is an extension of the Bartels–Stewart method and the Hessenberg-Schur method. Originally the Bartels –Stewart algorithm was used to solve the Sylvester equation.

ii) http://www.dm.unibo.it/~simoncin/matrixeq.pdf

That is important, is that the previous algorithms have complexity $O(n^3+m^3)$ that is much smaller than the complexity of the Kronecker product method. For example, solving a Lyapounov equation ($AX+XA^T=B$ with $n\times n$ matrices) -with Bartels–Stewart- has the same complexity ($\approx 40 n^3$) as finding eigenvalues and eigenvectors of a $n\times n$ matrix.

EDIT. You can also read this paper

https://pdfs.semanticscholar.org/5c47/c67aabe615a94174c418ebc5029e5630567d.pdf

Answer 2

As @loupblanc states in his answer, there is ample literature on solving equations of the form $AXB+CXD=E$ in $O(p^3)$. Unfortunately, as far as I know, none of these methods extend to cases where there is a third term, i.e., $AXB+CXD+EXF=G$, apart from special cases in which everything can be triangularized simultaneously.

Equations of your form appear in stochastic control; see for instance the book

Damm, Tobias, Rational matrix equations in stochastic control., Bremen: Univ. Bremen, Fachbereich Mathematik und Informatik (Diss.). viii, 200 p. (2002). ZBL1066.93500,

where this kind of equation is discussed.

One case in which you can approach the problem numerically is if one of the terms (say, the last) is much smaller than the others; in this case, you can make a series expansion of that $p^2\times p^2$ matrix inverse as $(M+N)^{-1}=M^{-1}-M^{-1}NM^{-1} + O(\|N\|^2)$, where $M=B^T\otimes A + D^T\otimes C$ and $N=F^T\otimes E$. You can use this approximate solution as a preconditioner for an iterative method such as CG or GMRES. This approach is (as far as I understand) described in

Damm, Tobias; Mena, Hermann; Stillfjord, Tony, Numerical solution of the finite horizon stochastic linear quadratic control problem., ZBL06770161.

This article focuses on a differential problem whose discretization produces equations of this kind.