Below, $M$ and $A$ are $n \times n$ real symmetric matrices. We can define a column vector as
\begin{align} x := \begin{bmatrix} M_{11} \\ M_{12} \\ \vdots \\ M_{1n}\\ M_{21}\\ M_{22}\\ \vdots \\ M_{nn} \end{bmatrix}\in \mathbb{R}^{n^2}\qquad (1) \end{align} So now we can write the matrix in terms of this column vector by $$ M_{ij} =: x_{n(i-1)+j} \ \forall i,j \in [1,n].\qquad (2) $$ I understand this representation as it can trivially be checked for choosing values of $i,j$. But my problem is from here on.
The Kronecker Delta function can then be used to represent a matrix $A$ by writing
$$ \delta_{il}\delta_{jk} =: A_{n(i-1)+j, n(k-1)+l} \ \forall i,j,k,l \in [1,n] \qquad (3) $$ which then leads to $$ x^\top A x = \sum_{i,j,k,l}^n M_{ij}M_{kl}\delta_{il}\delta_{jk}. \qquad(4) $$
How are equations (3), and (4) justified? Is this obvious? I am not able to see why. I am completely clear on equations (1) and (2) though.
I am looking for a thorough explanation or proof as to why (3), (4) are correct. Thanks.
In terms of the Kronecker product and vectorization operator, we could express what's going on as follows:
Let $e_1,\dots,e_n$ denote the canonical basis vectors of $\Bbb R^n$ (so that the $j$th entry of $e_i$ is $\delta_{ij}$). Let $E_{ij}$ denote the square matrix $e_i e_j^T$ (that is, the $p,q$ entry of $E_{ij}$ is $\delta_{ip}\delta_{jq}$).
The Kronecker product is particularly useful since we can now write $$ e_{n(i-1) + j} = e_i \otimes e_j $$
Your matrix $A$ can therefore be written as $$ A = \sum_{i,j=1}^n (e_i \otimes e_j)(e_j \otimes e_i)^T = \sum_{i,j = 1}^n E_{ij} \otimes E_{ji} $$ With that, we now have $$ \begin{align*} x^T A x &= \operatorname{vec}(M^T)^T \left( \sum_{i,j = 1}^n E_{ij} \otimes E_{ji} \right) \operatorname{vec}(M^T) = \sum_{i,j = 1}^n \operatorname{vec}(M^T)^T(E_{ij} \otimes E_{ji})\operatorname{vec}(M^T) \\ & = \sum_{i,j = 1}^n \operatorname{vec}(M^T)^T \operatorname{vec}(E_{ji}M^TE_{ij}) = \sum_{i,j = 1}^n \operatorname{vec}(M^T)^T \operatorname{vec}(E_{ji}M^TE_{ij}) \\ & = \sum_{i,j = 1}^n \operatorname{vec}(M^T)^T \operatorname{vec}(e_j[e_i^TM^Te_i]e_j^T) = \sum_{i,j = 1}^n M_{ii}\operatorname{vec}(M^T)^T \operatorname{vec}(E_{jj}) \\ & = \sum_{i,j = 1}^n M_{ii}\operatorname{vec}\left(\sum_{p,q=1}^n M_{pq} E_{qp}\right)^T \operatorname{vec}(E_{jj}) \\ &= \sum_{i,j = 1}^n \sum_{p,q=1}^n M_{ii} M_{pq}\operatorname{vec}(E_{qp})^T \operatorname{vec}(E_{jj}) = \sum_{i,j = 1}^n \sum_{p,q = 1}^n M_{ii}M_{pq} \delta_{qj}\delta_{pj} \\ & = \sum_{i,j = 1}^n M_{ii} M_{jj} = \operatorname{tr}(M)^2 \end{align*} $$
Here's the computation assuming that you meant to type $$ \delta_{ik}\delta_{jl} =: A_{n(i-1)+j, n(k-1)+l} \ \forall i,j,k,l \in [1,n] $$ as I explain in my comment. In this case, we have $A = I \otimes I$ (which is to say that $A$ is the $n^2 \times n^2$ identity matrix) so that $$ \begin{align*} x^TAx &= \operatorname{vec}(M^T)^T (I \otimes I) \operatorname{vec}(M^T) = \operatorname{vec}(M^T)^T\operatorname{vec}(I\,M^T\,I) \\ & = \operatorname{tr}(MM^T) = \sum_{i,j = 1}^n M_{ij}^2 \end{align*} $$