scipy.sparse for numpy.random_multivariate_normal

Question

I wanted to save a bit of memory, and thought I'd create a scipy.sparse identity matrix (dim is in the thousands, not terrible, but also not frugal). Notice its shape passes the assert:

cov = sigma_0 * sparse.identity(dim, dtype=np.float32)
assert (dim, dim) == cov.shape
result = np.random.multivariate_normal(mu, cov)

E   ValueError: cov must be 2 dimensional and square

The following, however, works fine:

cov = sigma_0 * np.identity(dim, dtype=np.float32)
assert (dim, dim) == cov.shape
result = np.random.multivariate_normal(mu, cov)

Did I miss it, somewhere, in the docs to say that sparse covariance matrices are expected fail with a ValueError?

Answer 1

What's happening here is that in np.random.multivariate_normal the input array is cast to an array:

cov = np.array(cov)

which ends up creating a scalar array of dtype object since numpy doesn't know anything about sparse matrices.

In [3]: cov = sparse.identity(100, dtype=np.float32)

In [4]: cov.shape
Out[4]: (100, 100)

In [5]: np.array(cov)
Out[5]: 
array(<100x100 sparse matrix of type '<type 'numpy.float32'>'
        with 100 stored elements (1 diagonals) in DIAgonal format>, dtype=object)