CODEWITHSUNDEEP
pythonnumpyscipysparse-matrix
ironv
ironv

Reputation: 1058

Correlation between array and sparse matrix

I have a sparse matrix (x) and an array (y). I would like to compute the correlation between each column in the matrix and the array. Shown below is a very simple approach which is slow. I was hoping somebody would have a faster/better approach.

import numpy as np
from scipy.sparse import rand as r1
from numpy.random import rand as r2

np.random.seed(1000)

nrow,ncol = 50000,4000
x = r1(nrow, ncol, format='csr', density=.05)
y = (r2(nrow)<=.6).astype(int)

correl = [(n,np.corrcoef(np.asarray(x[:,n].todense()).reshape(-1), y)[0,1]) for n in xrange(ncol)]
print correl[:10]

Upvotes: 4

Views: 1433

Answers (1)

Paul Panzer
Paul Panzer

Reputation: 53089

Using sparsity you can easily gain a speedup of >50x:

import numpy as np
from scipy.sparse import rand as r1
from numpy.random import rand as r2
from time import time

np.random.seed(1000)

nrow,ncol = 5000,4000
x = r1(nrow, ncol, format='csc', density=.05)
y = (r2(nrow)<=.6).astype(int)

t = []
t.append(time())
correl = [np.corrcoef(np.asarray(x[:,n].todense()).reshape(-1), y)[0,1] for n in range(ncol)]

t.append(time())

yy = y - y.mean()
xm = x.mean(axis=0).A.ravel()
ys = yy / np.sqrt(np.dot(yy, yy))
xs = np.sqrt(np.add.reduceat(x.data**2, x.indptr[:-1]) - nrow*xm*xm)

correl2 = np.add.reduceat(x.data * ys[x.indices], x.indptr[:-1]) / xs

t.append(time())

print('results equal --', np.allclose(correl, correl2))
print('run time (sec) -- OP: {}, new: {}'.format(*np.diff(t)))

Sample output:

results equal -- True
run time (sec) -- OP: 1.38134884834, new: 0.0178880691528

Explanation: To be able to take advantage of sparsity we standardise y which is dense anyway. And then compute the raw correlation between x and y. Because y is already zero-mean at this point the mean of x is nixed. It therefore remains to divide by the standard deviation of x. Here too we can avoid going through a dense matrix by calculating the raw 2nd moment and subtracting the squared mean.

Implementation detail: Please note that I have taken the liberty of switching to csc which is more suitable here (if needed, add x=x.tocsc() at the start). We use np.add.reduceat to perform the sums along the 'ragged' columns in a vectorized fashion. indices from the csc representation of the sparse matrix is convenient for selecting the elements of y corresponding to nonzero elements in x.

Upvotes: 7

Related Questions