# Solving for the trace of a matrix

In my research I have commonly had to deal with the following matrix equation (where $\Theta$ is the unknown): $$A \Theta + \Theta A^\text{T} = 2 D$$ known as the continuous time Lyapunov equation. All matrices are $n\times n$ with real coefficients, with $\Theta$ and $D$ also being symmetric positive definite.

This equation is known not to have a closed form solution and must be solved for each particular problem. However, all I really need is the quantity $$F = \text{tr}(A^\text{T} D^{-1} A \Theta)$$

Since all I need is this trace, I was wondering if it was possible to bypass the solution of the Lyapunov equation or, at least, simplify the process in some way.

Thus far I have not had any progress.

Thank you in advance for your time.

Best regards,

Gabriel

-
You mean, $A$ and $D$ are given, and you want all the possible values of $F$ for those $\Theta$'s that satisfy the first equation? – Berci Nov 27 '12 at 22:31
@Berci the matrix $\Theta$ is uniquely determined from the Lyapunov equation. Hence so is $F$. – Gabriel Landi Nov 27 '12 at 23:15

I have considered two different paths to this problem, one requires the diagonalization of $A$ (so may be not so useful depending on $A$) and the other uses the Kronecker product properties. One other note at the end gives an equivalent skew form of the equation. The Kronecker seems to have the most usefulness for finding the desired trace. The rest is maybe something of interest worth looking at.

First, the Kronecker. Using $\otimes$ to denote the Kronecker product, the vectorization of the formula gives

$$A\Theta + \Theta A^\top = 2D \iff \left( I \otimes A + A \otimes I\right) \operatorname{vec}(\Theta) = \operatorname{vec}(2D)$$

That together with a formula for the trace of a product of four matrices (the first line of the following), we get

\begin{align} \operatorname{Trace}(A^\top D^{-1} A \Theta) &= \left( \operatorname{vec}(D^{-\top})\right)^\top\left(A^\top \otimes A\right)\operatorname{vec}(\Theta) \\ &=\left( \operatorname{vec}(D^{-\top})\right)^\top\left(A^\top \otimes A\right)\left(I\otimes A + A \otimes I\right)^{-1}\operatorname{vec}(2D) \end{align} This may be useful in that it is a formula using only the knowns, though it is of higher dimension than the original.

The next path attempts to keep in terms of the original dimensions, though requires the diagonalization of A. Let $D = EE^\top$ since it is symmetric positive definite. This is found with use of the Cholesky factorization. Then we have: \begin{align} A\Theta + \Theta A^\top &= 2EE^\top \\ \end{align}
Let $A$ be diagonalizable with $A = Q \Lambda Q^{-1}$ for some diagonal $\Lambda$

\begin{align} Q \Lambda Q^{-1}\Theta + \Theta \left(Q \Lambda Q^{-1}\right)^\top &= 2EE^\top \\ Q \Lambda Q^{-1}\Theta + \Theta Q^{-\top} \Lambda^\top Q^\top &= 2EE^\top \\ \end{align}

Multiply the entire equation: $Q^{-1}\left( \dots \right) Q^{-\top}$: \begin{align} \Lambda Q^{-1}\Theta Q^{-\top} + Q^{-1}\Theta Q^{-\top} \Lambda^\top &= 2Q^{-1}EE^\top Q^{-\top} \\ \end{align} Use substitution $X = Q^{-1}\Theta Q^{-\top}$ ( noting that $\Lambda^\top = \Lambda$): \begin{align} \Lambda X + X \Lambda &= 2Q^{-1}EE^\top Q^{-\top} \\ \end{align}

Solving for $X$ here requires an element by element inspection. Let $C = 2Q^{-1}EE^\top Q^{-\top}$. Then element $i,j$ has equation: \begin{align} \left[ \Lambda X + X \Lambda \right]_{ij} &= C_{ij} \\ \Lambda_{ii}X_{ij} + X_{ij}\Lambda{jj} &= C_{ij} \\ \end{align} This may be seen since left multiplication by a diagonal scales the rows and right multiplication scales the columns. So for each element of $X$ we have:

$$X_{ij} = \frac{C_{ij}}{\Lambda_{ii} + \Lambda_{jj}}$$

Or in matrix form: $$\left[ X_{ij}\right] = \left[ \frac{C_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right] = X$$ In terms of $\Theta$: $$\Theta = Q \left[ \frac{C_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right]Q^\top = Q \left[ \frac{\left[2Q^{-1}EE^\top Q^{-\top}\right]_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right]Q^\top$$

Let us look now at $A \Theta A^\top$, using $A = Q \Lambda Q^{-1}$ $$A \Theta A^\top = Q \Lambda Q^{-1} \left( Q \left[ \frac{\left[2Q^{-1}EE^\top Q^{-\top}\right]_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right]Q^\top\right)\left(Q \Lambda Q^{-1}\right)^\top$$

$$A \Theta A^\top = Q \Lambda \left[ \frac{\left[2Q^{-1}EE^\top Q^{-\top}\right]_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right] \Lambda Q^\top$$

Now for $\operatorname{Trace}(D^{-1}A \Theta A^\top) = \operatorname{Trace}(E^{-\top}E^{-1}A \Theta A^\top) = \operatorname{Trace}(E^{-1}A \Theta A^\top E^{-\top})$ (using $D=EE^\top$ and the cyclic property of trace):

$$\operatorname{Trace}(E^{-1}A \Theta A^\top E^{-\top}) = \operatorname{Trace}(E^{-1}Q \Lambda \left[ \frac{\left[2Q^{-1}EE^\top Q^{-\top}\right]_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right] \Lambda Q^\top E^{-\top})$$

The $\Lambda$ on the left gives row scaling and on the right gives column scaling, for $$\operatorname{Trace}(D^{-1}A \Theta A^\top) = \operatorname{Trace}(E^{-1}Q \left[ \frac{\Lambda_{ii}\Lambda_{jj}\left[2Q^{-1}EE^\top Q^{-\top}\right]_{ij}}{\Lambda_{ii} + \Lambda_{jj}}\right] Q^\top E^{-\top})$$

From here ...? Note that $D=EE^\top$ is congruence transformed, scaled (element $ij$ by $\frac{2\Lambda_{ii}\Lambda{jj}}{\Lambda_{ii} + \Lambda_{jj}}$)then inverse congruence transformed to get the trace equivalent of $D^{-1}A\Theta A^\top$. Not sure exactly how useful that is though.

One other possibly interesting note, since you know that $\Theta$ is symmetric (your statement or also may be seen from above) then $$\Theta = BB^\top$$ for some $B$. This gives \begin{align} ABB^\top + BB^\top A^\top &= 2D \\ \end{align} Using $Y=ABB^\top$ and $Y^\top = BB^\top A^\top$ we have $$Y+Y^\top = 2D$$ and $$Y^\top+Y = 2D^\top = 2D$$ (since $D$ symmetric is given). $Y$ then must be of the form of a skew symmetric plus $D$: $$Y = S + D$$ with $S^\top = -S$. This gives two other forms for $\Theta$ as $$\Theta = BB^\top = A^{-1}Y = A^{-1}\left( S + D\right)$$ and $$\Theta = BB^\top = Y^\top A^{-\top} = \left( -S + D\right)A^{-\top}$$ Plugged into the original equation checks that it works: \begin{align} A\Theta + \Theta A^\top &= 2D \\ AA^{-1}\left( S + D\right) + \left( -S + D\right)A^{-\top}A^\top &= 2D \\ \left( S + D\right) + \left( -S + D\right) &= 2D \\ 2D &= 2D \end{align} Looking at the two forms $\Theta = \Theta$ \begin{align} A^{-1}\left( S + D\right) = \left( -S + D\right)A^{-\top} \\ \left( S + D\right) = A\left( -S + D\right)A^{-\top} \\ \left( S + D\right)A^\top = A\left( -S + D\right) \\ \end{align} and a formula for the skew matrix $S$ is realized: \begin{align} SA^\top + AS &= AD - DA^\top \\ \end{align}

-
This is incredible. Thank you so very very much. – Gabriel Landi Nov 29 '12 at 10:42
Nice answer, but there's a typo in the Kronecker solution, the term $(A^T\otimes A)$ should actually be $(A\otimes A)$. – greg Apr 12 '15 at 6:00
The Kronecker sum can simplify the solution to: $2\,{\rm vec}(D^{-1})^T(A\otimes A) (A\oplus A)^{-1} {\rm vec}(D)$ – greg Apr 12 '15 at 6:04
I think you are correct on both accounts, except I believe $\operatorname{vec}(D^{-1})$ should be $\operatorname{vec}(D^{-\top})$ in your last comment. – adam W Apr 13 '15 at 5:03
But in the problem he says that $D=D^T$. – greg Apr 13 '15 at 17:44

Perhaps a little simplification, using that $$tr.(AB)=tr.(BA)\,\,\,,\,\,tr.(A+B)=tr.(A)+tr.(B)\,\,,\,tr.(A^t)=tr.(A)\;:$$

$$A^tD^{-1}A\Theta=A^t D^{-1}\left(2D-\Theta A^t\right)=2A^t-A^tD^{-1}\Theta A^t\Longrightarrow$$

$$tr.(A^tD^{-1}\Theta A)=2\,tr.(A)-tr.\left((A^t)^2D^{-1}\Theta\right)$$

-
I have come across that manipulation but, unfortunately, it does not remove $\Theta$ so it doesn't really work much. But, in any case, thanks for your help. – Gabriel Landi Nov 27 '12 at 23:16

In my research I have commonly had to deal with the following matrix equation (where Θ is the unknown): $$A \Theta + \Theta A^\text{T} = 2 D$$ known as the continuous time Lyapunov equation. All matrices are $n\times n$ with real coefficients, with $\Theta$ and $D$ also being symmetric positive definite.

This equation is known not to have a closed form solution and must be solved for each particular problem.

Actually, what you have here is a very well known matrix equation, called the continuous Lyapunov equation. It is usually written in the form $$A^T P +P A+Q=0,$$ where $P$ is the unknown and $Q$ is assumed to be positive definite symmetric. When $A$ is stable (eigenvalues in $\mathbb{C}^-$), there is indeed a closed form solution given by $$P=\int_0^\infty \exp(A^T t)Q\exp(At)\,dt.$$ (You can verify this by direct substitution back into the matrix equation. The stability assumption on $A$ is so that the integral converges.)

So, in your context, if $A$ is stable, $$\Theta=\int_0^\infty \exp(At)(-2D)\exp(A^Tt)\,dt.$$ Perhaps this will assist you in getting the trace result you seek.

-