It's a usual trick in deep learning. Matrix multiplications can be done in parallel:
$$
(A.B)_{ij} = \sum_{k} A_{ik} B_{kj}
$$
since you have to do this sum for each component, then you can divide your computing time by the number of processes you use - i.e. compute these components several at a time rather than one by one. Assuming you use a naive matrix computation then you have a $O(\dfrac{n^3}{K})$ with K number of processes. That's why GPUs are so much better for deep learning (usually have hundreds if not thousands of cores) which allows you to divide by a factor of several 100s at worst the time to compute the matrix multiplication.
In this particular case, you can use block matrix multiplications (we expend the rules of matrix multiplications to matrices):
$$
\mathbf{𝑊}_{∗(𝑥𝑟)}[𝑥_𝑡,𝑦_{𝑡−1}]= \begin{bmatrix}
\mathbf{𝑊}_{∗𝑥} & \mathbf{0} \\
\mathbf{0} & \mathbf{𝑊}_{∗y}
\end{bmatrix}[𝑥_𝑡,𝑦_{𝑡−1}]
$$
if you apply the traditional matrix multiplications (keeping in mind matrix sizes), you get:
$$
\mathbf{𝑊}_{∗(𝑥𝑟)}[𝑥_𝑡,𝑦_{𝑡−1}]=\mathbf{𝑊}_{∗𝑥}𝑥_𝑡 + \mathbf{0} * 𝑦_{𝑡−1} + \mathbf{0} * 𝑥_𝑡 + \mathbf{𝑊}_{∗y} * 𝑦_{𝑡−1} =\mathbf{𝑊}_{∗𝑥}𝑥_𝑡 + \mathbf{𝑊}_{∗y} * 𝑦_{𝑡−1}
$$