Let's say I have a sparse array and a dense array that has the same number of columns but fewer rows:
from scipy.sparse import csr_matrix
import numpy as np
sp_arr = csr_matrix(np.array([[1,0,0,0,1],[0,0,1,0,0],[0,1,0,0,1],[0,0,0,1,1],[0,0,0,1,0]]))
arr = np.random.rand(10).reshape(2,5)
print(arr)
[[ 0.47027789 0.82510323 0.01321617 0.66264852 0.3618022 ]
[ 0.80198907 0.36350616 0.10254934 0.65209401 0.094961 ]]
I would like to get an array containing all of the submatrices for the indices that contain values for each row of the sparse array. For example, the indices for data in sp_arr are as follows:
0: [0, 4] 1: [2] 2: [1, 4] 3: [3, 4] 4: [3]
My output should look like this:
array([array([[ 0.47027789, 0.3618022 ],
[ 0.80198907, 0.094961 ]]),
array([[ 0.01321617],
[ 0.10254934]]),
array([[ 0.82510323, 0.3618022 ],
[ 0.36350616, 0.094961 ]]),
array([[ 0.66264852, 0.3618022 ],
[ 0.65209401, 0.094961 ]]),
array([[ 0.66264852],
[ 0.65209401]])], dtype=object)
I can create this with the following code, but as the size of the arrays scale up (greatly in my case) it gets really slow.
output = np.empty(sp_arr.shape[0], dtype=object)
for row in range(sp_arr.shape[0]):
output[row] = arr[:, sp_arr[row].indices]
I've thought about vectorizing the process and applying it along the axis, but np.apply_along_axis doesn't work with sparse matrices, and unfortunately while this example is small enough to make dense and then use apply_along_axis my actual sparse matrix is much too large to do so (>100Gb).
I had thought that perhaps there is a fancy way to index or use something like hsplit to accomplish this with already vectorized methods, but so far I haven't had any luck. Is there someway to achieve this that is just escaping me?
Update
Per the answer from @Divakar, which is great, I found another way to implement the same thing with the ever so slightest, and negligible, improvement.
@Divakars best answer was:
def app2(sp_arr, arr):
r,c = sp_arr.nonzero()
idx = np.flatnonzero(r[1:] > r[:-1])+1
idx0 = np.concatenate(( [0] , idx, [r.size] ))
arr_c = arr[:,c]
return [arr_c[:,i:j] for i,j in zip(idx0[:-1], idx0[1:])]
Which increased my performance by 50 - 60x! But its kind of hard to read.
What I discovered is that given the csr_matrix format you can use the indices and indptr attributes to your advantage here.
def final_app():
idx = sp_arr.indptr
arr_c = arr[:, sp_arr.indices]
out = [arr_c[:, i:j] for i, j in zip(idx[:-1], idx[1:])]
return out
In the end the performance is statistically the same (less than 50ms improvement on a sparse matrix 276538 x 33114), but it feels easier to understand. More Importantly this approach does include the rows for which there are no values at all, whereas the previous approach does not. This may seem unimportant, but for my use case it is pretty critical.
Update 2
In response to @EelcoHoogendoorn. The problem is part of a parallel implementation of the alternative least squares with regularization method that I am trying to implement in python. This comes from the oft cited paper Large-scale Parallel Collaborative Filtering for the Netflix Prize The normal way of doing this is by having distributed copies of the Ratings, User, and Items matrices across processes. I thought it might be interesting to see what would happen if we constructed all of the item submatrices up front and just distributed those to the processes. That way the processes need only return the feature columns for either one user or one item respectively and those could be used to update the User and Item matrices.
The above problem was actually the bottleneck in my current implementation. And per your comment, in this situation I don't believe the transposition is critical as part of the algorithm takes the dot product of each submatrix with its transpose.