Solution of a linear matrix equation

Question

i'm trying to rewrite a linear matrix equation $$ X=A_1XA_2+A_3X^TA_4+B $$

where $ X\in\mathbb{R}^{m\times n} $ and $ A_1\text{, }A_2\text{, }A_3\text{, }A_4\text{, }B $ are matrices of appropriate sizes.

So far i got to express the equation in terms of the components of $X$: $$ X_{ij}-\sum_{k=1}^{m}\sum_{l=1}^{n}\left(A_{1,ik}\cdot A_{2,lj}+A_{3,il}\cdot A_{4,kj}\right)\cdot X_{kl}=B_{ij} $$ but is there any neat way to express this linear system in terms of a coefficient matrix $C$ and unfolded vectors $x$ and $b$ $$ Cx=b $$

so that i can feed the problem to Matlab? I'd rather avoid for loops, since the matrices are quite big and the linear problem is framed inside a bigger, iterative computation.

Thanks!

Answer 1

As the comment states, we can use vectorization. In Matlab notation, $\operatorname{vec}(X)$ just is $X(:)$, and we can go from $\operatorname{vec}(X)$ to $X$ using the reshape command.

Vectorizing both sides yields $$ X - A_1 XA_2 - A_3 X^TA_4 = B \implies\\ \operatorname{vec}[X - A_1 XA_2 - A_3 X^TA_4] = \operatorname{vec}(B) \implies\\ \operatorname{vec}(X) - A_2^T \otimes A_1 \operatorname{vec}(X) - A_4^T \otimes A_3\operatorname{vec}(X^T) = \operatorname{vec}(B). $$ Here, $\otimes$ denotes the Kronecker product, implemented in Matlab with kron. This is almost in the form that we want, but we need to deal with $X^T$. For that, we need the fact that $\operatorname{vec}(X^T) = P\operatorname{vec}(X)$, where $$ P = \sum_{i=1}^m \sum_{j=1}^n (e_j \otimes e_i)(e_i \otimes e_j)^T = \sum_{i=1}^n \sum_{j=1}^m \operatorname{vec}(E_{ij})[\operatorname{vec}(E_{ji})]^T. $$ $P$ can be built with a quick for-loop and only needs to be put together once (assuming $m,n$ are constant). I don't know of a nicer approch here, unfortunately.

From there, we can finally write $$ \operatorname{vec}(X) - A_2^T \otimes A_1 \operatorname{vec}(X) - (A_4^T \otimes A_3)P\operatorname{vec}(X) = \operatorname{vec}(B) \implies\\ (I - A_2^T \otimes A_1 - (A_4^T \otimes A_3)P)\operatorname{vec}(X) = \operatorname{vec}(B) \implies\\ (I - A_2^T \otimes A_1 - (A_4^T \otimes A_3)P)x = b $$ where $x = \operatorname{vec}(X),b = \operatorname{vec}(B)$.