切片稀疏（scipy）矩陣

Question

我將不勝感激任何幫助，以瞭解從scipy.sparse包中切片lil_matrix（A）時的以下行爲。

其實，我想提取一個基於行和列的任意索引列表的子矩陣。

當我使用的代碼，這兩條線：

x1 = A[list 1,:] 
x2 = x1[:,list 2] 

一切都很好，我可以提取合適的子矩陣。

當我試圖做到這一點的一條線，它失敗（返回的矩陣是空的）

x=A[list 1,list 2] 

爲什麼會這樣呢？總的來說，我在matlab中使用了一個類似的命令，它在那裏工作。 那麼，爲什麼不使用第一個，因爲它的工作原理？這似乎相當耗時。由於我必須經歷大量的輸入，我想用一個命令來加速它。也許我使用錯誤的稀疏矩陣類型...任何想法？

Answer 1

切片情況與此語法：

a[1:4] 

用於=陣列（[1,2,3,4,5,6,7,8,9]），結果是

array([2, 3, 4]) 

元組的第一個參數指示要保留的第一個值，第二個參數指示不保留的第一個值。

如果您在兩側都使用列表，這意味着您的數組具有與列表長度一樣多的維度。

所以，用你的語法，你可能會需要像這樣：

x = A[list1,:,list2] 

根據A的形狀

希望它確實幫助你。

Answer 2

您已經使用的方法，

A[list1, :][:, list2] 

似乎是選擇從備件矩陣所需的值的最快方法。請參閱下面的基準。

然而，要回答你如何與一個單一的指標， 你就需要使用選擇的A從任意行值和列問題，所謂的"advanced indexing"：

A[np.array(list1)[:,np.newaxis], np.array(list2)] 

憑藉先進的索引如果arr1和arr2是NDarrays，的A[arr1, arr2]的(i,j)分量等於

A[arr1[i,j], arr2[i,j]] 

因此，您希望arr1[i,j]等於list1[i]所有j和 arr2[i,j]等於list2[j]所有i。

可以在broadcasting（見下文）的幫助下設置 arr1 = np.array(list1)[:,np.newaxis]和arr2 = np.array(list2)。

的arr1形狀(len(list1), 1)而arr2形狀是 (len(list2),)因爲需要時新的軸被添加上自動左 ，其廣播到(1, len(list2))。

每個陣列可以進一步廣播以形成(len(list1),len(list2))。 這正是我們想要的 A[arr1[i,j],arr2[i,j]]是有意義的，因爲我們希望(i,j)運行於形狀爲(len(list1),len(list2))的結果數組的所有可能的索引。

---

下面是這表明A[list1, :][:, list2]是最快的選項微基準爲一個測試用例：

In [32]: %timeit orig(A, list1, list2) 
10 loops, best of 3: 110 ms per loop 

In [34]: %timeit using_listener(A, list1, list2) 
1 loop, best of 3: 1.29 s per loop 

In [33]: %timeit using_advanced_indexing(A, list1, list2) 
1 loop, best of 3: 1.8 s per loop 

---

這裏是我以前爲基準設置：

import numpy as np 
import scipy.sparse as sparse 
import random 
random.seed(1) 

def setup(N): 
    A = sparse.rand(N, N, .1, format='lil') 
    list1 = np.random.choice(N, size=N//10, replace=False).tolist() 
    list2 = np.random.choice(N, size=N//20, replace=False).tolist() 
    return A, list1, list2 

def orig(A, list1, list2): 
    return A[list1, :][:, list2] 

def using_advanced_indexing(A, list1, list2): 
    B = A.tocsc() # or `.tocsr()` 
    B = B[np.array(list1)[:, np.newaxis], np.array(list2)] 
    return B 

def using_listener(A, list1, list2): 
    """https://stackoverflow.com/a/26592783/190597 (listener)""" 
    B = A.tocsr()[list1, :].tocsc()[:, list2] 
    return B 

N = 10000 
A, list1, list2 = setup(N) 
B = orig(A, list1, list2) 
C = using_advanced_indexing(A, list1, list2) 
D = using_listener(A, list1, list2) 
assert np.allclose(B.toarray(), C.toarray()) 
assert np.allclose(B.toarray(), D.toarray()) 

Answer 3

對我來說，unutbu的解決方案效果很好，但速度很慢。

我發現作爲一個快速替代，

A = B.tocsr()[np.array(list1),:].tocsc()[:,np.array(list2)] 

你可以看到行的和COL的獲得單獨切割，但每一個轉換，以最快的稀疏格式，讓指數這段時間。

在我的測試環境中，這個代碼比另一個快1000倍。

我希望我不會說錯或錯誤。

Answer 4

B[arr1, arr2]中的同時索引確實有效，它比我的機器上的listener's solution快。請參閱下面Jupyter示例中的In [5]。要將其與上述答案進行比較，請參閱In [6]。此外，我的解決方案不需要.tocsc()轉換，使其更易讀IMO。

請注意B[arr1, arr2]工作，arr1和arr2必須broadcastable numpy的陣列。

A 更快的解決方案，但是，使用B[list1][:, list2]作爲pointed out by unutbu。請參閱下面的In [7]。

In [1]: from scipy import sparse 
     : import numpy as np 
     : 
     : 

In [2]: B = sparse.rand(1000, 1000, .1, format='lil') 
     : list1=[1,4,6,8] 
     : list2=[2,4] 
     : 
     : 

In [3]: arr1 = np.array(list1)[:, None] # make arr1 a (n x 1)-array 
     : arr1 
     : 
     : 
Out[3]: 
array([[1], 
     [4], 
     [6], 
     [8]]) 

In [4]: arr2 = np.array(list2)[None, :] # make arr2 a (1 x m)-array 
     : arr2 
     : 
     : 
Out[4]: array([[2, 4]]) 

In [5]: %timeit A = B.tocsr()[arr1, arr2] 
100 loops, best of 3: 13.1 ms per loop 

In [6]: %timeit A = B.tocsr()[np.array(list1),:].tocsc()[:,np.array(list2)] 
100 loops, best of 3: 14.6 ms per loop 

In [7]: %timeit B[list1][:, list2] 
1000 loops, best of 3: 205 µs per loop