Completing the square in the context of a multivariate gaussian distribution is pretty straightforward (here is a good explanation).
But what if we are talking about a matrix normal distribution?
Consider the following matrices: $Q_1$ and $H_1$ are ($k \times k$), $Q_2$ and $H_2$ are ($l \times l$) and all of them are symmetric and positive definite (covariance matrices). $B_1$, $B_2$ and $B_3$ are ($k \times l$).
When I combine the two matrix normal densities, I end up with the following expression:
$tr[Q_1(B_1-B_2)Q_2(B_1-B_2)' + H_1(B_1-B_3)H_2(B_1-B_3)']$.
Now I am looking for a neat expression of the form
$tr[X(B_1-Y)Z(B_1-Y)' + W]$.
Question: how can I rearrange terms in the original expression so that neither $X$, $Y$, $Z$ nor $W$ contain $B_1$?
Note: I tried to use some properties of the $tr[]$ to vectorize the expression and then complete the square the 'traditional' way. But then I ran into the problem of not being able to recover the matrix form of the distribution because of this.
