Linear Algebra – Solving Linear Matrix Equations of the Form $X = AXA^T + C$

linear algebramatricesproblem solvingsolution-verification

I'm trying to solve this matrix equation:

$$X = AXA^T + C$$

In particular,

$$
X = \begin{bmatrix} 1.5 & 1 \\ -0.7 & 0 \end{bmatrix}
X
\begin{bmatrix} 1.5 & -0.7 \\ 1 & 0 \end{bmatrix}
+ \begin{bmatrix} 1 & 0.5 \\ 0.5 & 0.25 \end{bmatrix}.
$$

I started by writing $X$ in open form.

$$X = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$$

The equation becomes

$$
\begin{bmatrix} a & b \\ c & d \end{bmatrix}
= \begin{bmatrix} 1.5 & 1 \\ -0.7 & 0 \end{bmatrix}
\begin{bmatrix} a & b \\ c & d \end{bmatrix}
\begin{bmatrix} 1.5 & -0.7 \\ 1 & 0 \end{bmatrix}
+ \begin{bmatrix} 1 & 0.5 \\ 0.5 & 0.25 \end{bmatrix} \\
\begin{bmatrix} a & b \\ c & d \end{bmatrix}
= \begin{bmatrix} 1.5a+c & 1.5b+d \\ -0.7a & -0.7b \end{bmatrix}
\begin{bmatrix} 1.5 & -0.7 \\ 1 & 0 \end{bmatrix}
+ \begin{bmatrix} 1 & 0.5 \\ 0.5 & 0.25 \end{bmatrix} \\
\begin{bmatrix} a & b \\ c & d \end{bmatrix}
= \begin{bmatrix} 2.25a+1.5b+1.5c+d & -1.05a-0.7c \\ -1.05a-0.7b & 0.49a \end{bmatrix}
+ \begin{bmatrix} 1 & 0.5 \\ 0.5 & 0.25 \end{bmatrix} \\
\begin{bmatrix} a-1 & b-0.5 \\ c-0.5 & d-0.25 \end{bmatrix}
= \begin{bmatrix} 2.25a+1.5b+1.5c+d & -1.05a-0.7c \\ -1.05a-0.7b & 0.49a \end{bmatrix} \\
\begin{bmatrix} -1 & -0.5 \\ -0.5 & -0.25 \end{bmatrix}
= \begin{bmatrix} 1.25a+1.5b+1.5c+d & -1.05a-b-0.7c \\ -1.05a-0.7b-c & 0.49a-d \end{bmatrix} \\
$$

Equating cells of RHS and LHS matrices I wrote the linear equation set

$$
\begin{bmatrix}
1.25 & 1.5 & 1.5 & 1 \\
-1.05 & -1 & -0.7 & 0 \\
-1.05 & -0.7 & -1 & 0 \\
0.49 & 0 & 0 & -1 \\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
c \\
d
\end{bmatrix} =
\begin{bmatrix}
-1 \\
-0.5 \\
-0.5 \\
-0.25
\end{bmatrix}
$$

I solved this set of equations using two different software tool to find exact same $a$, $b$, $c$ and $d$ values:

$$X = \begin{bmatrix} a & b \\ c & d \end{bmatrix} =
\begin{bmatrix} 18.8802 & -11.3672 \\ -11.3672 & 9.5013 \end{bmatrix}$$

Now remember the equation above:

$$
\begin{bmatrix} a & b \\ c & d \end{bmatrix}
= \begin{bmatrix} 1.5 & 1 \\ -0.7 & 0 \end{bmatrix}
\begin{bmatrix} a & b \\ c & d \end{bmatrix}
\begin{bmatrix} 1.5 & -0.7 \\ 1 & 0 \end{bmatrix}
+ \begin{bmatrix} 1 & 0.5 \\ 0.5 & 0.25 \end{bmatrix}
$$

When I put the $a$, $b$, $c$ and $d$ values in the RHS, I don't find the $X$ matrix on the LHS.

$$
\begin{bmatrix} 1.5 & 1 \\ -0.7 & 0 \end{bmatrix}
\begin{bmatrix} 18.8802 & -11.3672 \\ -11.3672 & 9.5013 \end{bmatrix}
\begin{bmatrix} 1.5 & -0.7 \\ 1 & 0 \end{bmatrix}
+ \begin{bmatrix} 1 & 0.5 \\ 0.5 & 0.25 \end{bmatrix}
= \begin{bmatrix} 163.4421 & -332.7002 \\ -332.7002 & 857.9770 \end{bmatrix}
$$

What am I doing wrong in my calculations?

Best Answer

Your solution for $X$ is correct.

Here is another solution (it is cheating, though, because I used Octave). On he plus side, it shows another approach to solving the equation.

$A$ is diagonalizable, and so has a basis of eigenvectors, $u_1,u_2$. Define the linear operator $L(X) = X-A X A^*$. You want to solve $L(X) = C$. Note that $u_i u_j^*$ forms a basis for the space of $2 \times 2 $ matrices, and $L(u_i u_j^*) = (1-\lambda_i \overline{\lambda_j}) u_i u_j^*$, hence $L$ is diagonalizable.

To express $C$ in terms of $u_i u_j^*$, let $C = \sum_{i,j} [\Gamma]_{ij} u_i u_j^*$, and solve for $\Gamma$. Let $U$ be the matrix whose columns are $u_1 ,u_2$. Note that $\sum_i [\Gamma]_{ij} u_i = U \Gamma e_j$, so we have $C = \sum_{j} U \Gamma e_j u_j^* = \sum_{j} U \Gamma e_j (U e_j)^* = \sum_{j} U \Gamma e_j e_j^T U^* = U \Gamma (\sum_{j} e_j e_j^T) U^* = U \Gamma U^*$. It follows that $\Gamma = U^{-1} C (U^{-1})^*$, and if we solve $L(\sum_{i,j} [\tilde{X}]_{ij} u_i u_j^*) = \sum_{i,j} [\Gamma]_{ij} u_i u_j^*$, we get $[\tilde{X}]_{ij} = \frac{1}{1-\lambda_i \overline{\lambda_j}} [\Gamma]_{ij}$.

Reversing the procedure to obtain $X$ in the standard basis, we have $X = U \tilde{X} U^*$.

Using Octave we obtain the same result as above.

a = [1.5 1 ; -0.7 0 ]
c = [1 0.5 ; 0.5 0.25 ]

[u,d] = eig(a)

l = [d(1,1) ; d(2,2) ]

phi = ones(2,2)-l*l'

gamma = inv(u)*c*inv(u')

x_tilde = (1./phi).*gamma
x = u*x_tilde*u'

# check result...
a*x*a'+c-x

Here is the output:

octave:44> a*x*a'+c-x
ans =

   0.00000 - 0.00000i   0.00000 - 0.00000i
   0.00000 - 0.00000i   0.00000 + 0.00000i

octave:45> x
x =

   18.8802 +  0.0000i  -11.3672 -  0.0000i
  -11.3672 -  0.0000i    9.5013 +  0.0000i
Related Question