Decomposing a matrix into two or more independent matrices I am almost certain the solution to the problem is very obvious but I can't see an easy way of doing it.
I would like to implement a solution based on matrix multiplication to compute the following expression:
$$
H = \sum_{i,k}M^2_{i,k} + \sum_{i,k}\sum_{(j,l) \neq (i,k)}\phi_{i,j}\,\rho_{k,l}\,M_{k,i}\,M_{l,j}
$$
For the sake of simplicity, I have bounded the variables $i$ and $k$ to be $i = 0,1$ and $k = 0,1$. I can now write the above expression as a matrix multiplication problem in the following way:
$$
\left( \begin{array}{cccc}
M_{00} & M_{01} & M_{10} & M_{11}\end{array} \right) \times C \times 
\left( \begin{array}{}
M_{00} \\
M_{01} \\
M_{10} \\
M_{11}
\end{array} \right)
$$
where, if I havent made any mistakes, $C$ would be defined as
$$
C = \left( \begin{array}{cccc}
1 & \phi_{01}\,\rho_{00} & \phi_{00}\,\rho_{01} & \phi_{01}\,\rho_{01}\\
\phi_{10}\,\rho_{00} & 1 & \phi_{10}\,\rho_{01} & \phi_{11}\,\rho_{01}\\
\phi_{00}\,\rho_{10} & \phi_{01}\,\rho_{10} & 1 & \phi_{01}\,\rho_{11}\\
\phi_{10}\,\rho_{10} & \phi_{11}\,\rho_{10} & \phi_{10}\,\rho_{11} & 1\end{array} \right) 
$$
My question is:
How can I express the matrix $C$ as a multiplication of 2 or more matrixes such that the coefficients for $\phi$ and $\rho$ are in different matrixes? Basically I would like to have a simple way of storing this data in different matrixes, and then just simply multiply all my matrixes to solve the combinations that give me $H$.
 A: Ignore the summation constraints for a moment, and write the function as
$$\eqalign{
 H &= M:M + M:\rho M\phi^T \cr
   &= {\rm vec}(M):{\rm vec}(M) + {\rm vec}(M):{\rm vec}(\rho M\phi^T) \cr
   &= {\rm vec}(M)^T{\rm vec}(M) + {\rm vec}(M)^T(\phi\otimes\rho){\rm vec}(M) \cr
 &= m^Tm + m^TGm\cr
}$$
where : represents the inner/Frobenius product, $\otimes$ represents the Kronecker product, and ${\rm vec}()$ is the matrix vectorization operator.
The constraint says something about the intrinsic structure of the $G$ matrix. 
From your sample calculation, it appears the constrained matrix has its diagonal elements removed. This can be expressed in a variety of ways
$$\eqalign{
 G^\prime &= G - I\odot G \cr
  &= G - {\rm Diag}(G) \cr
  &= G - D \cr
}$$
where $\odot$ represents the elementwise/Hadamard product, and $I$ is the identity matrix.
Substituting the constained matrix yields
$$\eqalign{
 H &= m^Tm + m^TG^\prime m \cr
   &= m^TIm + m^TGm - m^TDm \cr
   &= m^T\,(I+G-D)\,m \cr
   &= m^T\,C\,m \cr\cr
}$$
You were mainly interested in an expression for the $C$ matrix 
$$\eqalign{
 C &= I \,+ (J-I)\odot(\phi\otimes\rho) \cr
}$$
where $J$ is the all ones matrix, which is the identity element for the Hadamard product.
