Rank of a joint distribution matrix $\newcommand{\P}{\mathbf P}$$\newcommand{\rank}{\operatorname{rank}}$Suppose I have discrete random variables $X$ and $Y$ supported on $\{1,\dots,m\}$ and $\{1,\dots,n\}$ respectively, and WLOG I'll assume $m \geq n$. Let $\P$ be the $m\times n$ matrix with $\P_{ij} = P(X=i, Y=j)$ and $f$ and $g$ are the marginal distributions of $X$ and $Y$ respectively. I know $\rank\P = 1 \iff X\perp Y$. In general how does the rank of $\P$ interact with the kinds of dependencies that are possible? I'm also interested in relating more continuous measures of $\P$'s closeness to being rank $1$, like $\frac{d_1}{\sum_i d_i}$ where the $d_i$ are the singular values of $\P$, to the lack of dependence between $X$ and $Y$. I've explored trying to bound $\P - fg^T$ in terms of these quantities as well but to no avail so far.

I'd also be interested in exploring the simpler case of $m=n$ and $X\stackrel{\text{d}}= Y$ so $P$ is square, the marginals are the same, and also I'll take $P(X=i,Y=j) = P(X=j, Y=i)$ so $P$ is symmetric. Then $P = Q\Lambda Q^T$ by the spectral theorem but maybe this isn't helpful. Perhaps nonnegative matric factorization provides more insight?
 A: Here's an idea:
The higher is the rank, the lower the probability that the two variables are almost independent
In the picture below I plot the empirical cumulative distribution of mutual information between two random variables given the rank of the matrix P, assuming that matrix P has shape 3x4. Obviously, for rank=1 mutual information is always zero. For rank=2 it is almost zero most of the time, but for rank=3 it is almost zero much less frequently.

Here's the code to get this result. It is extremely inefficient, as I have to generate random matrices of a given rank until I get a matrix for which all elements are non-negative, so it can be normalized to a 2D probability distribution. This simulaton could be extended to higher ranks if this inefficiency can be overcome in some way, but so far I have no ideas.
from time import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import ortho_group  # Requires version 0.18 of scipy

def rand_posdiag_matrix(m, n, rk):
    D = np.zeros((m, n))
    diag = np.zeros(np.min([m, n]))
    diag[:rk] = np.random.uniform(-1, 1, rk)
    np.fill_diagonal(D, diag)
    
    U = ortho_group.rvs(dim=m)
    V = ortho_group.rvs(dim=n)
    return U.dot(D.dot(V.T))
    
def rand_pos_matrix(m, n, rk):
    M = rand_posdiag_matrix(m, n, rk)
    while np.min(M) < 0:
        M = rand_posdiag_matrix(m, n, rk)
    return M / np.sum(M)

def H1D(p):
    return -p.dot(np.log(p))

def MI(p2D):
    pX = np.sum(p2D, axis=0)
    pY = np.sum(p2D, axis=1)
    return H1D(pX) + H1D(pY) - H1D(p2D.flatten())
    
def cdf(x):
    return sorted(x), np.linspace(0, 1, len(x))
    
nTrial = 100
nRank = 3
rezMat = np.zeros((nRank, nTrial))

for iRank in range(nRank):
    for iTrial in range(nTrial):
        t = time()
        
        M = rand_pos_matrix(3, 4, iRank+1)
        rezMat[iRank, iTrial] = MI(M)
        
#         print(iRank, iTrial, time() - t)
        
plt.figure()
for iRank in range(nRank):
    plt.plot(*cdf(rezMat[iRank]), label=str(iRank+1))
plt.legend()
plt.show()

A: I am not aware whether there is a relation between the rank of P and any kind of dependency between $X$ and $Y$. But, there is a nice characterization of the dependence between $X$ and $Y$ with the corresponding nonnegative rank.
Let us assume for simplicity that both, $X$ and $Y$ take values in $\{1, \ldots, n\}$ (the case $n \neq m$ is analogous).
Let $Z$ be a random variable taking values in $\{1, \ldots, r\}$. Then $X$, $Y$ are independent conditioned on $Z$, if $P$ can be written as
$$ P(X = i, Y = j) = \sum_{k=1}^{r} P(X=i|Z=k) \cdot P(Y=j|Z=k) \cdot P(Z = k).\quad \quad \quad \quad(1)$$
On the other hand we can define the nonnegative rank of a matrix in the following way: The nonnegative rank is the minimal integer $r$ such that there exists a decomposition
$$ P = \sum_{k=1}^{r} \mathbf{v}_k \cdot \mathbf{w}_k^T $$
where $\mathbf{v}_k, \mathbf{w}_k \in \mathbb{R}_{+}^{n}$. Hence the matrix admits a decomposition into rank-1 matrices generated by nonnegative vectors (instead of general vectors as in the case of the usual rank).
Now the following assertions are equivalent:
(a) $X$ and $Y$ are independent conditioned on $Z$ which takes $r$ values.
(b) $\textrm{nn-rank}(P) = r$.
To show (a) $\Longrightarrow$ (b) just set $(\mathbf{v}_k)_i := P(X=i|Z=k)$ and $(\mathbf{w}_k)_j := P(Y=j|Z=k) \cdot P(Z=k)$ and the independency relation leads to a nn-rank decomposition of rank $r$.
To show (b) $\Longrightarrow$ (a) set
$$P(X=i|Z=k) := \frac{(\mathbf{v}_k)_i}{\Vert \mathbf{v}_k \Vert_1}$$
and
$$P(Y=j|Z=k) := \frac{(\mathbf{w}_k)_j}{\Vert \mathbf{w}_k \Vert_1}$$
Dividing by the the 1-norm guarantees the normalization of the defined probability distributions. Now it only remains to define $P(Z=k)$.
We set
$$P(Z = k) := \Vert \mathbf{v}_k \Vert_1 \cdot \Vert \mathbf{w}_k \Vert_1$$
This definition gives again a probability distribution, since all vectors are nonnegative and
$$ \sum_{k=1}^{r} P(Z=k) = \sum_{k=1}^{r} \Vert \mathbf{v}_k \Vert_1 \cdot \Vert \mathbf{w}_k \Vert_1 = \sum_{k=1}^{r} \sum_{i,j=1}^{n} (\mathbf{v}_k)_i \cdot (\mathbf{w}_k)_j = \sum_{i,j=1}^{n} P(X=i, Y=j) = 1$$
Plugging in these definitions into the rank-decomposition yields to the sum (1) which shows the statement.
A: I doubt anything will come out of this. Especially if you are not able to be more specific than "the kinds of dependencies that are possible". Take the most simple example: two Bernoulli variables or a 2x2 matrix. Given the probabilities of $X$ and $Y$ all dependency is contained in a single parameter. Nevertheless, all these joint distributions - except the independent one - have a rank 2 matrix.
To make further progress you should in my opinion:

*

*Study simple cases like 2x3 or 3x3

*Make up your mind what kind of dependency you are really interested in.

