How to rewrite a matrix expression involving Kronecker product and trace as a quadratic form? For a vector $\mathbf{x} \in \mathbb{R}^n$ and symmetric matrices $\mathbf{A} \in \mathbb{R}^{n\times n}$, $\mathbf{B} \in \mathbb{R}^{n^2\times n^2}$, I want to find an expression for the matrix $\mathbf{C}$ such that
$$
\mathrm{tr}\big(\mathbf{A}(\mathbf{x}^T\otimes\mathbf{I}_n)\mathbf{B}(\mathbf{x}^T\otimes\mathbf{I}_n)^T \big) = \mathbf{x}^T \mathbf{C} \mathbf{x},
$$
where $\mathbf{I}_n$ is the $n\times n$ unit matrix and $\otimes$ denotes the Kronecker product.
As an example, for the case $n=2$, this is given by
$$
\mathbf{C} = 
\begin{bmatrix}
b_{11}a_{11} + b_{22}a_{22} + 2 b_{12}a_{12} & b_{13}a_{11} + b_{24}a_{22} + (b_{23} + b_{14})a_{12} \\\\ b_{13}a_{11} + b_{24}a_{22} + (b_{23} + b_{14})a_{12} &  b_{33}a_{11} + b_{44}a_{22} + 2 b_{34}a_{12}.
\end{bmatrix}
$$
How can I get a general expression for $\mathbf{C}$ from the above expression?
 A: We first consider the case where $B$ is of the form $B_1 \otimes B_2$, with $B_1,B_2$ of size $n \times n$. We find that
$$
\begin{align}
\operatorname{tr}(A(x^T \otimes I_n)B(x^T \otimes I_n)^T) &= 
\operatorname{tr}(A(x^T \otimes I_n)(B_1 \otimes B_2)(x \otimes I_n))
\\ & = \operatorname{tr}(A(x^TB_1x \otimes B_2)) = (x^TB_1 x)\operatorname{tr}(AB_2). 
\end{align}
$$
Now, note that we can write $B$ in the form
$$
B = \pmatrix{B_{11} & \cdots & B_{1n}\\ \vdots & \ddots & \vdots \\ 
B_{n1} & \cdots & B_{nn}}  \implies B = \sum_{i,j} E_{ij} \otimes B_{ij},
$$
where $E_{ij}$ denotes the $n \times n$ matrix with a $1$ as the $i,j$ entry and zeros elsewhere, and each $B_{ij}$ has size $n \times n$.  We find that
$$
\begin{align}
\operatorname{tr}(A(x^T \otimes I_n)B(x^T \otimes I_n)^T) &= 
\sum_{i,j} (x^T E_{ij} x)\operatorname{tr}(AB_{ij})
\\ & = \sum_{i,j} \operatorname{tr}(AB_{ij}) x_i x_j = x^TCx,
\end{align}
$$
where $C$ is the matrix whose $i,j$ entry is $\operatorname{tr}(AB_{ij})$. In the case that $A = I$, $C$ coincides with the "partial trace" ("over the second space") of the matrix $B$. In general, $C$ can be written as the partial trace of the matrix $(I_n \otimes A)B$.
A: $\def\v{{\rm vec}}\def\o{\otimes}\def\t#1{{\rm Tr}\Big(#1\Big)}\def\p{{\partial}}\def\grad#1#2{\frac{\p #1}{\p #2}}$Calculate
the singular value decomposition
$$\eqalign{
B &= \sum_{k=1}^r\sigma_ku_kv_k^T \qquad r = {\rm rank}(B) \\
}$$
and reshape the singular vectors into matrices such that
$$\eqalign{
u_k &= \v\big(U_k\big) \qquad v_k = \v\big(V_k\big)\\
}$$
For typing convenience, use a colon as a product notation for the trace
$$\eqalign{
F:G &= \sum_{i=1}^m \sum_{j=1}^n F_{ij}G_{ij} \;=\; \t{F^TG} \\
}$$
Write the function using the above notations.
$$\eqalign{
\t{A(x^T\o I)B(x\o I)}
 &= A^T:(x^T\o I)B(x\o I) \\
 &= (x\o I)A^T(x^T\o I):B \\
 &= \sum_{k=1}^r(x\o I)A^T(x^T\o I):\sigma_ku_kv_k^T \\
 &= \sum_{k=1}^r\left(xx^T\o A^T\right)v_k:\sigma_ku_k \\
 &= \sum_{k=1}^r\v\Big(A^TV_kxx^T\Big):\sigma_k\v\Big(U_k\Big) \\
 &= \sum_{k=1}^rA^TV_kxx^T:\sigma_kU_k \\
 &= \sum_{k=1}^rx^T\left(\sigma_kU_k^TA^TV_k\right)x \\
}$$
The desired matrix is therefore
$$C = \sum_{k=1}^r \sigma_kU_k^TA^TV_k$$
