Why do I end up with a highly correlated matrix when I multiply two strictly positive random matrices?

Question

I've been working with some large positive matrices for a machine learning problem I'm trying to solve. The problem involves multiplying a matrix $A$ that is a set of users and items, with a matrix $B$ that is a set of items and attributes. When we look at the columns of $B$, we see minimal correlations. However, when we create a new matrix, $C=A \bullet B$, we see that the columns of $C$ are highly correlated. This is clearly a problem for machine learning.

One thing that we've tried is normalization. It seems like centering and scaling the the features of $B$ to a standard normal distribution fixes the problem of high correlations. Also, in order to maintain the positive constraints, if we $L1$ norm the rows of $A$ and the columns of $B$, then the correlations do not appear in the final product.

The weirdness comes when we try to use $L2$ normalizations, which to me shouldn't be that different than $L1$. However, if we $L2$ norm the rows of $A$ and the columns of $B$, then the correlations still appear (although they're slightly weaker) ... This is intuitively very strange to me and I can't seem to figure out why this is.

Does anyone have any insight into what's going on here? I've included R code below that recreates the issue with random matrices.

    library(corrplot)

    # randomly generate strictly positive matrices
    item_attributes <- matrix(runif(100000, min=0, max=1),1000,10)
    users_items <- matrix(sample(0:1, 10000*1000, rep=T, prob=c(0.95, 0.05)),10000,1000)

    # L1 Norm
    item_attributes_l1 <- t(t(item_attributes)/abs(apply(item_attributes,2,sum)))
    users_items_l1 <- users_items / abs(apply(users_items,1,sum))

    # L2 Norm
    item_attributes_l2 <- t(t(item_attributes)/sqrt(apply(item_attributes^2,2,sum)))
    users_items_l2 <- users_items/sqrt(apply(users_items^2,1,sum))

    # multiply the matrices
    mat_prod <- users_items %*% item_attributes
    mat_prod_l1 <- users_items_l1 %*% item_attributes_l1
    mat_prod_l2 <- users_items_l2 %*% item_attributes_l2

    # check for correlations in the final products
    corrplot(cor(mat_prod), main='\ndot product')
    corrplot(cor(mat_prod_l1), main='\nl1 norm')
    corrplot(cor(mat_prod_l2), main='\nl2 norm')

Answer 1

It helps to think about the underlying meaning of the scores that result from each approach. Denote the binary user-item scores as $X_{hi}$ for user $h=1,...,H$ on item $i=1,...,I$ and the item-attribute scores as $T_{ia}$ for item $i=1,...,I$ and attribute $a=1,...,A$.

First consider using the dot product (without normalization). The score on attribute $a$ for user $i$ is $\sum_{i=1}^I X_{hi} T_{ia}$, which is the sum of the scores on attribute $a$ for all items purchased by user $h$. The scores on the $A$ attributes will be correlated because they are driven in part by the number of items purchased by user $h$.

library(corrplot)

# randomly generate strictly positive matrices
item_attributes <- matrix(runif(100000, min=0, max=1),1000,10)
users_items <- matrix(sample(0:1, 10000*1000, rep=T, prob=c(0.95, 0.05)),10000,1000)

# dotproduct
mat_prod <- users_items %*% item_attributes
corrplot(cor(mat_prod), main='\ndot product')

Next consider taking the L1 norm of the user-item scores (but leaving the item-attribute scores as is):

item_attributes_l1 <- t(t(item_attributes)/abs(apply(item_attributes,2,sum)))
users_items_l1 <- users_items / abs(apply(users_items,1,sum))
mat_prod_l1 <- users_items_l1 %*% item_attributes_l1
# corrplot(cor(mat_prod_l1), main='\nl1 norm')

mat_prod_l1A <- users_items_l1 %*% item_attributes
corrplot(cor(mat_prod_l1A), main='\nl1 norm - items only')

The score on attribute $a$ for user $i$ would then be mean of the scores on attribute $a$ across all items purchased by user $h$: $\frac{1}{n_a} \sum_{i=1}^I X_{hi} T_{ia}$, where $n_a = \sum_{i=1}^I X_{hi}$ is the number of items purchased. Taking the mean rather than the sum removes the correlation driven by the number of items purchased.

Next consider taking the L2 norm of the user-item scores (again, leaving the item-attribute scores as is):

item_attributes_l2 <- t(t(item_attributes)/sqrt(apply(item_attributes^2,2,sum)))
users_items_l2 <- users_items/sqrt(apply(users_items^2,1,sum))
mat_prod_l2 <- users_items_l2 %*% item_attributes_l2
# corrplot(cor(mat_prod_l2), main='\nl2 norm')

mat_prod_l2A <- users_items_l2 %*% item_attributes
corrplot(cor(mat_prod_l2A), main='\nl2 norm - items only')

The L2 norm of a set of $i$ binary items is $\sqrt{\sum_{i=1}^I X_{hi}^2} = \sqrt{\sum_{i=1}^I X_{hi}} = \sqrt{n_a}$, and so the score on attribute $a$ for user $i$ would then be $\sqrt{n_a}$ times the mean of the scores on attribute $a$ across all items puchased by user $h$: $\frac{1}{\sqrt{n_a}} \sum_{i=1}^I X_{hi} T_{ia}$. The scores on the $A$ attributes will still be correlated because they are all driven in part by $\sqrt{n_a}$.

Of these three options, it seems to me that taking the L1 norm on items (i.e., averaging the item-attribute scores instead of summing them) produces the most interpretable results. If total number of items purchased is of interest, then you can always calculate this and include it as a further predictor.

How best to norm the item-attribute scores seems to me to be an entirely separate question. Note that taking the L1 norm of the item-attribute scores amounts to assigning each item a "portion" of the total amount of an attribute, which might or might not make sense depending on context. Rather than taking the L2 norm, it might make more intuitive sense to re-scale the item-attribute scores to have mean 0 and SD 1, i.e.,

item_attributes_Z <- apply(item_attributes, 2, scale)
apply(item_attributes_Z, 2, mean)
apply(item_attributes_Z, 2, sd)