Given a probability space $\Omega,$ the space of square-integrable measurable functions $\Omega \to \mathbb{R}^n$ ("random vectors") can be made a vector space over $\mathbb{R}$ in a natural way. Call this space $V.$ In probability theory, we proceed to define several operators on this space, like the expectation operator $E : V \to \mathbb{R}^n$ given by $(X_1,X_2...,X_n) \mapsto (E(X_1),E(X_2)...,E(X_n))$.
However, going just a bit deeper into the theory, we start to see some properties of $E$ nicer than linearity over $\mathbb{R}$ would alone suggest. For example, for any $k \times n$ matrix $A$, we find that $E(AX) = AE(X).$ Similar occurrences occur with the bilinear covariance operator $\mathrm{Cov} : V \to \mathbb{R}^{n \times n}$. For example, for any $k \times n$ matrices $A$ and $B,$ we find $\mathrm{Cov}(AX,BY) = A\mathrm{Cov}(X,Y)B^T,$ where $B^T$ denotes the transpose of $B.$
On one level, one can just view this as matrix algebra (and this may be all there is to it). But I've always been inclined to look for deeper algebraic structure than just matrix algebra when I see matrices, so I'm wondering if there's a deeper algebraic reason to this. For example, we could have viewed $V$ as a module over $n \times n$ matrices, but this approach doesn't seem to explain the transposes and the generalization to $k \times n$ matrices with $k \neq n.$ So, I'm wondering if there's some algebraic structure to $V$ in which the "matrix linearity" of the form seen in $E$ and $\mathrm{Cov}$ become natural (and hence easy to remember!).