Inner Product for Geometric Interpretation of Multivariate Random Vectors I was looking into the geometric interpretation of random variables as random vectors in a vector space. The textbook I'm referring to defined $\operatorname {Cov}(X,Y)$ as the inner product for any two zero-mean RVs and claimed that it is symmetric, positive definite and bilinear. Now, I understand why this is the case for univariate random variables. But this is simply not true for multivariate random variables.
For the latter, $\operatorname {Cov}(X,Y) = \operatorname {Cov}(Y,X)^T$ given $X,$ $Y$ are multivariate RVs.
What would be a legitimate inner product for multivariate RVs then?

By definition, a multivariate random variable $X$ is a measurable map from a probability space $(\Omega,\mathfrak{F}, \mathbb{P})$ into a topological vector space $\mathbb{F}$.  When $\mathbb F$ has a norm $||\ ||$ (which determines its topology), then that norm determines a norm-like function on the random variables defined by
$$||X||^2=\int_\Omega ||X(\omega)||^2\,\mathrm{d}\mathbb{P}(\omega).\tag{*}$$
As is the case even for univariate random variables, where the dimension of $\mathbb F$ is $1,$ this can fail to be a norm because $||X||=0$ implies only that $X$ is almost surely zero -- but not that it is identically zero.  We can create a relevant vector space out of this by declaring two random variables to be equivalent when they are almost surely the same.  The set of equivalence classes of random variables is then a vector space and the quotient of on that vector space indeed is a norm.  This space is called $L^2(\Omega,\mathbb{F}).$
When $\mathbb{F}=\mathbb{R}^n,$ you might use the usual Euclidean norm $||(x_1,x_2,\ldots,x_n)||=\sqrt{x_1^2 + \cdots + x_n^2}.$
We are done, because inner products and norms are equivalent via the polarization identity.  Given a normed vector space $(\mathbb{V}, ||\ ||)$ (not of characteristic $2$) and $x,y\in\mathbb V,$ define their inner product as
$$x\cdot y = \frac{1}{4}\left(||x+y||^2 - ||x-y||^2\right).$$
Conversely, any inner product $\cdot\,,$ determines a norm as
$$||x||^2 = x\cdot x.$$
These formulas invert each other, establishing the claimed one-to-one correspondence.

Notice that when you use the usual Euclidean norm on $\mathbb{R}^n,$ linearity of integration gives a nice formula.  Let $X = (X_1,X_2,\ldots,X_n)$ be an $n$-variate random variable and rewrite $(*)$ as
$$\begin{aligned}
||X||^2 &= \int_\Omega ||X(\omega)||^2\,\mathrm{d}\mathbb{P}(\omega) \\
&= \int_\Omega |X_1(\omega)|^2 + |X_2(\omega)|^2 + \cdots + |X_n(\omega)|^2\,\mathrm{d}\mathbb{P}(\omega)\\
&= \int_\Omega |X_1(\omega)|^2 \,\mathrm{d}\mathbb{P}(\omega) + \int_\Omega |X_2(\omega)|^2 \,\mathrm{d}\mathbb{P}(\omega) + \cdots + \int_\Omega |X_n(\omega)|^2 \,\mathrm{d}\mathbb{P}(\omega)\\
&= ||X_1||^2 + ||X_2||^2 + \cdots + ||X_n||^2.
\end{aligned}$$
In this case, the squared norm of the random variable $X$ is the sum of squares of the norms of the component (univariate) random variables.  A comparable formula for the inner product follows immediately from polarization: the inner product of two such variables is the sum of the inner products of their components qua univariate random variables.
