Sparse matrix multiplication when results' sparsity is known (in python|scipy|cython)
Question: Suppose we want to compute C=A*B for given sparse matrices A,B but are
interested in a very small subset of entries of C, represented by a list of
index pairs:  
rows=[i1, i2, i3 ... ]  
cols=[j1, j2, j3 ... ]  
Both A and B are quite large (say 50Kx50K), but very sparse (<1% of entries is
non-zero).

How can we compute this subset of the multiplication?

Here's a naive implementation that works really slow:

    
    
    def naive(A, B, rows, cols):
        N = len(rows)
        vals = []
        for n in xrange(N):
            v = A.getrow(rows[n]) * B.getcol(cols[n])
            vals.append(v[0, 0])
    
        R = sps.coo_matrix((np.array(vals), (np.array(rows), np.array(cols))), shape=(A.shape[0], B.shape[1]), dtype=np.float64)
        return R
    

even for small matrices this is quite bad:

    
    
    import scipy.sparse as sps
    import numpy as np
    D = 1000
    
    A = np.random.randn(D, D)
    A[np.abs(A) > 0.1] = 0
    A = sps.csr_matrix(A)
    B = np.random.randn(D, D)
    B[np.abs(B) > 0.1] = 0
    B = sps.csr_matrix(B)
    
    X = np.random.randn(D, D)
    X[np.abs(X) > 0.1] = 0
    X[X != 0] = 1
    X = sps.csr_matrix(X)
    rows, cols = X.nonzero()
    naive(A, B, rows, cols)
    

On my machine, naive() finishes after 1 minute, and most of the effort is
spent on structuring the rows/cols (in getrow(), getcol()).  
Of course, converting this (very small) example to dense matrices, the
computation takes about 100ms:

    
    
    A0 = np.array(A.todense())
    B0 = np.array(B.todense())
    X0 = np.array(X.todense())
    A0.dot(B0) * X0
    

Any thoughts on how to **efficiently** compute such matrix multiplication?

  * Note: This question is almost identical to the following question: [Subset of a matrix multiplication, fast, and sparse](http://stackoverflow.com/questions/18792096/subset-of-a-matrix-multiplication-fast-and-sparse)  
However, there, A and B are **full** matrices, and, one of the dimensions is
very low (say, 10) the proposed solutions seem to benefit from both.


Answer: The format of your sparse matrices is important here. You always need a row
form A and a column from B. So, store `A` as a `csr` and `B` as `csc` to get
rid of the `getrow`/`getcol` overhead. Unfortunately, this is only a small
part of the story.

The best solution depends a lot on the structure of your sparse matrix (a lot
of sparse columns/rows, etc), but you might try one based on dictionaries and
sets. For matrix `A` for each row the following are kept:

  * a set with all non-zero column indices on that row
  * a dictionary with the non-zero indices as keys and the corresponding non-zero values as values

For matrix `B` similar dicts and sets are kept for each column.

To calculate element (M, N) in the multiplication result, row M of `A` is
multiplied with column N of `B`. The multiplication:

  * find the set intersection of the non-zero sets
  * calculate the sum of multiplications of the non-zero elements (i.e. the intersection above)

In most cases this should be very fast, as in a sparse matrix the set
intersection is usually very small.

Some code:

    
    
    class rowarray():
        def __init__(self, arr):
            self.rows = []
            for row in arr:
                nonzeros = np.nonzero(row)[0]
                nzvalues = { i: row[i] for i in nonzeros }
                self.rows.append((set(nonzeros), nzvalues))
    
        def __getitem__(self, key):
            return self.rows[key]
    
        def __len__(self):
            return len(self.rows)
    
    
    class colarray(rowarray):
        def __init__(self, arr):
            rowarray.__init__(self, arr.T)
    
    
    def maybe_less_naive(A, B, rows, cols):
        N = len(rows)
        vals = []
        for n in xrange(N):
            nz1,v1 = A[rows[n]]
            nz2,v2 = B[cols[n]]
            # list of common non-zeros
            nz = nz1.intersection(nz2)
            # sum of non-zeros
            vals.append(sum([ v1[i]*v2[i] for i in nz]))
    
        R = sps.coo_matrix((np.array(vals), (np.array(rows), np.array(cols))), shape=(len(A), len(B)), dtype=np.float64)
        return R
    
    D = 1000
    
    Ap = np.random.randn(D, D)
    Ap[np.abs(Ap) > 0.1] = 0
    A = rowarray(Ap)
    Bp = np.random.randn(D, D)
    Bp[np.abs(Bp) > 0.1] = 0
    B = colarray(Bp)
    
    X = np.random.randn(D, D)
    X[np.abs(X) > 0.1] = 0
    X[X != 0] = 1
    X = sps.csr_matrix(X)
    rows, cols = X.nonzero()
    maybe_less_naive(A, B, rows, cols)
    

This is a bit more efficient, the multiplication takes approximately 2 seconds
for the test (80 000 elements). The results seem to be essentially the same.

* * *

A few comments on the performance.

There are two operations performed for each output element:

  * set intersection
  * multiplication

The complexity of set intersection should be O(min(m,n)) where m and n are the
numbers of non-zeros in each operand. This is invariant of the size of the
matrix, only the average number of non-zeros per row/column is important.

The number of multiplications (and dict lookups) depends on the number of non-
zeros found in the intersection above.

If both matrices have randomly distributed non-zeros with probability
(density) p, and the row/column length is n, then:

  * set intersection: O(np)
  * dictionary lookup, multiplication: O(np^2)

This shows that with really sparse matrices finding the intersections is the
critical point. This can also be verified by profiling; most of the time is
spent calculating the intersections.

When this is reflected to the real world, we seem to spend around 20 us for a
row/column of 80 non-zeros. This is not blindingly fast, and the code can
certainly be made faster. Cython may be one solution, but this may be one of
the problems where Python is not the best possible solution. A simple linear
matching (merge sort -type algorithm) for sorted integers should be at least
an order of magnitude faster when written in C.

One important thing to note is that the algorithm can be done in parallel for
several elements at a time. There is no need to settle for a single thread, as
the calculations are independent as far as one thread handles one output
point.

