The trick is that the $y_i$ were carefully chosen so that $y_1^2+y_2^2+y_3^2$, $y_1^2y_2^2 + y_1^2y_3^2 + y_2^2y_3^2$, and $y_1y_2y_3$ are symmetric polynomials in the $x_i$s. This means that the value of these polynomials do not change if you permute the $x_i$. For example: \begin{align*} y_1^2+y_2^2+y_3^2 &= a^2(x_1 + x_2 - x_3 - x_4)^2 + a^2(x_1 - x_2 + x_3 - x_4)^2 + a^2(x_1 - x_2 - x_3 + x_4)^2 \\ \text{permute $x_1$ and $x_2$:}&\mathrel{\phantom=} a^2(x_2 + x_1 - x_3 - x_4)^2 + a^2(x_2 - x_1 + x_3 - x_4)^2 + a^2(x_2 - x_1 - x_3 + x_4)^2 \\ &= a^2(x_1 + x_2 - x_3 - x_4)^2 + a^2(x_1 - x_2 - x_3 + x_4)^2 + a^2(x_1 - x_2 + x_3 - x_4)^2 \\ &= y_1^2+y_2^2+y_3^2 \end{align*}
By the fundamental theorem of symmetric polynomials every symmetric polynomial can be written in terms of the elementary symmetric polynomials, which are in fact exactly the expressions derived for $b$, $c$, $d$, and $e$ (up to a factor of $\pm a$). Therefore by this theorem there is some expression for $y_1^2+y_2^2+y_3^2$ in terms of $a$, $b$, $c$, $d$, and $e$ only. There are various ways of computing what that expression actually is; for example, the three identities I gave were computed with the following Maple code using Gröbner bases:
T := lexdeg([y1,y2,y3], [x1,x2,x3,x4], [b,c,d,e]):
G := Groebner[Basis]([
y1-a*(x1+x2-x3-x4),
y2-a*(x1-x2+x3-x4),
y3-a*(x1-x2-x3+x4),
b-(-a*(x1+x2+x3+x4)),
c-(a*(x1*x2+x1*x3+x1*x4+x2*x3+x2*x4+x3*x4)),
d-(-a*(x1*x2*x3+x1*x2*x4+x1*x3*x4+x2*x3*x4)),
e-(a*(x1*x2*x3*x4))],
T):
Groebner[Reduce](y1^2+y2^2+y3^2, G, T);
Groebner[Reduce](y1^2*y2^2+y1^2*y3^2+y2^2*y3^2, G, T);
Groebner[Reduce](y1*y2*y3, G, T);