$\newcommand{\diag}{\operatorname{diag}}$Let $X = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$, $Y= \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}$
I have a vector:
$$V = \begin{bmatrix} x_1(y_1 - \sum\limits_{i \in \{1,2\}} x_iy_i)\\ x_2(y_2 - \sum\limits_{i \in \{1,2\}} x_iy_i) \end{bmatrix}$$
I re-arranged the above to put into a matrix equation:
$$V= \diag(x_1,x_2)\Bigg(\begin{bmatrix} y_1 \\ y_2 \end{bmatrix} - \begin{bmatrix} 1 \\ 1 \end{bmatrix} \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}\Bigg)$$
where $\diag(x_1,x_2) = \begin{bmatrix} x_1 & 0 \\ 0 & x_2 \end{bmatrix}$
Which I then write $V = \diag(x_1,x_2)(Y-\mathbf{1}X^TY)$
But unfortunately I think is confusing and not very visually pleasing. Is there some way to rewrite the expression of $V$ so that it is more compact. I am thinking some way to use Kronecker product, or something like that but I am not sure
Your $V$ is $$ V = \begin{pmatrix} -x_1^2 y_1 - x_1 x_2 y_2 + x_1 y_1 \\ -x_2^2 y_2 - x_1 x_2 y_1 + x_2 y_2 \\ \end{pmatrix} $$ thus a cubic polynomial for each component.