First of all, $Cov(X,Y)=0$ does not implies independence unless $X$ and $Y$ are jointly Gaussian. So your orthogonality is way weaker than independence.
It looks disappointing but it is only because you are looking at it the wrong way. Independence does have an analogy to vector space orthogonality, but you have to treat $X$ and $Y$ as entire subspaces instead of individual vectors. Also, because the independence business happens in entirely in the $\sigma$-algebra world, which has less structures than an inner product space (e.g., there are no dimensions, basis nor spans in $\sigma$-algebras), it will be more fruitful to look for a weakened generalization of the vector space inner product idea rather than trying to find a vector space inner product that captures the entirely less structured $\sigma$-algebra situation.
Probability Theory View
Because independence is defined as $f(x,y)=f_x(x)f_y(y)$, you can also say that $X$ and $Y$ are independent if and only if for all functions $g(x)$ and $h(y)$,
\begin{equation}
\int\int g(x)h(y)f(x,y)dxdy = \left(\int g(x)f_x(x)dx\right)\left(\int h(x)f_y(y)dy\right).
\end{equation}
Without independence, the left hand side integral simply can't be factored out so nicely. In fact, it is easy to check that the left hand side is bilinear in both $g$ and $h$ and the right hand side is two "simple dot products," so we can write
\begin{equation}
\langle g, h\rangle_f = \langle g, f_x \rangle \langle h, f_y \rangle.
\end{equation}
Remember this factorization and we will see how this has something to do with orthogonality in standard linear algebra.
(Note that the integral equation is same as $E(g(X)h(Y))=E(g(X))E(h(Y))$ but it's not really the point here.)
Linear Algebra Analogy
Now let's make an analogy with $\mathbb{R}^n$ vectors. Let $X$ be a subspace of $\mathbb{R}^n$ and $Y$ another subspace that is $X$'s orthogonal complement. Let $g$ and $h$ be any vectors from $\mathbb{R}^n$, and $F$ be an $n\times n$ rank-$r$ matrix which defines a bilinear form $\langle g, h\rangle_F = g^TFh$. Let $\mathcal{X}=\{x_i\}_i$ and $\mathcal{Y}=\{y_j\}_j$ be some orthogonal unit-length basis of $X$ and $Y$ respectively. Denote $\mathcal{E}=\mathcal{X}\cup\mathcal{Y}$ so that $\mathcal{E}$ is a basis that spans the entire space.
Obviously by basis expansion we can write $g$ and $h$ as so $g=\sum_{e_i\in\mathcal{E}}g_ie_i$ and $h=\sum_{e_j \in\mathcal{E}}h_je_j$.
Any rank-$r$ $F$ can be written as a sum
\begin{equation}
F=\sum_{k=1}^r u_k v_k^T
\end{equation}
by the SVD decomposition. Usually, the bilinear form
\begin{equation}
\langle g, h\rangle_F = \sum_{e_j\in\mathcal{E}}\sum_{e_i\in\mathcal{E}} (g_i e_i)^T F (h_j e_j) = \sum_{e_j\in\mathcal{E}}\sum_{e_i\in\mathcal{E}}\sum_{k=1}^r g_ih_j e_i^T u_kv_k^T e_j
\end{equation}
cannot be factorized into the product of two nice separate single sums. But if $F$ happens to be rank-$1$ with $F=uv^T$ such that $u\in X$ and $v\in Y$ then by basis expansion on $u$ and $v$, we can write $F$ as
\begin{equation}
F=\left(\sum_{x_\ell\in\mathcal{X}}{f_x}_{\ell}x_\ell\right)\left(\sum_{y_m\in\mathcal{Y}}{f_y}_{m}y_m^T\right),
\end{equation}
where ${f_x}_\ell$ and ${f_y}_m$ are scalar coordinates. So we can now factorize the double sum into a product of two single sums as follows:
\begin{align}
\langle g, h\rangle_F
&= \sum_{e_j\in\mathcal{E}}\sum_{e_i\in\mathcal{E}} (g_i e_i)^T \left(\sum_{x_\ell\in\mathcal{X}}{f_x}_{\ell}x_\ell\right)\left(\sum_{y_m\in\mathcal{Y}}{f_y}_{m}y_m^T\right) (h_j e_j)\\
&= \sum_{y_j\in\mathcal{Y}}\sum_{x_i\in\mathcal{X}} (g_i {f_x}_i) (h_j {f_y}_j)\\
&= \left(\sum_{x_i\in\mathcal{X}} g_i {f_x}_i\right)\left(\sum_{y_j\in\mathcal{Y}} h_j {f_y}_j\right)\\
&= \langle g, f_x\rangle \langle h, f_y\rangle,
\end{align}
where the right hand side $\langle \cdot, \cdot\rangle$ are the "simple dot product." We have just proved that the bilinear form can be factored out for any arbitrary $g\in \mathbb{R}^n$ and $h\in \mathbb{R}^n$. The first line simplifies into the second line because $e_i^Tx_\ell$ is zero except when $e_i = x_\ell$; and in such case, $x_\ell^Tx_\ell = 1$ by definition.
Notice that in this argument:
- Orthogonality of $X$ and $Y$ together with $F$ being "compatible" with them is essential, because otherwise the second line of the simplification simply won't cancel out so nicely.
- Be cautious that $F$ is defined using $X$ and $Y$ just for illustration. In the probability version we define the joint density $f$ at the beginning when forming a probability space, and we call $X$ and $Y$ independent if the pre-defined $f$ can be factored out.
- Why can $g$ and $h$ be any vectors in the analogy, but in the probabilistic version $g$ is a function of only $X$ and $h$ only $Y$? Actually this is only because I have taken $\langle \cdot, \cdot \rangle_F$ to be $\mathbb{R}^n\times\mathbb{R}^n\to\mathbb{R}$ for illustration. The more general situation is one where $\langle \cdot, \cdot \rangle_F:V\times W\to\mathbb{R}$ where $V$ and $W$ are different spaces. Now $g$ must fit in the left slot and $h$ the right slot. But if $V$ and $W$ were any general vector spaces then we can't write down the decomposition of $F$ so easily any more, neither can we say a subspace of $V$ is "orthogonal" to a subspace of $W$ because they are from different spaces (so "inner product" equals to zero doesn't make sense any more). Nonetheless, because $F$ is pre-defined, we can still get around this problem by defining that some sub-spaces of $V$ and $W$ are "$F$-independent" if $\langle \cdot, \cdot \rangle_F$ can be factored out. This is somehow what you are seeing in the probabilistic definition.
- In fact, in advanced probability theory, statistical independence is defined on $\sigma$-algebras first, and then we say that two random variables are independent if the $\sigma$-algebras they generate are independent. So intuitively it also makes more sense to think of independence as a property about two sets rather than two vector elements.
