愚蠢，有效地将一个矩阵插入另一个矩阵

Question

我正试图通过删除零元素、进行外积，然后用零行放大得到的矩阵或插入到零矩阵中，来提高两个向量的外积的效率。(使用scipy稀疏矩阵实际上并不起作用，因为转换成本很高，我正在一遍又一遍地进行。(

import numpy
dim = 100
vec = np.random.rand(1, dim)
mask = np.flatnonzero(vec > 0.8)
vec_sp = vec[:, mask]
mat_sp = vec_sp.T * vec_sp # This is faster than dot product
# Enlarge matrix or insert into zero matrix

由于它是两个向量的外积，我知道原始矩阵中的零行和零列，它们是掩码变量中的索引。要看到这一点，

a = np.array(((1,0,2,0))).reshape(1,-1)
a.T * a
>> array([[1, 0, 2, 0],
[0, 0, 0, 0],
[2, 0, 4, 0],
[0, 0, 0, 0]])

我尝试了两种不同的解决方案：一种是使用numpy的insert方法，另一种是对mat_sp变量附加方法。整个过程变成了一个for循环，而且非常缓慢。

for val in mask:
if val < mat_sp.shape[0]:
mat_sp = np.insert(mat_sp, val, values=0, axis=1)
mat_sp = np.insert(mat_sp, val, values=0, axis=0)
else:
mat_sp = np.append(mat_sp, values=np.zeros((mat_sp.shape[0], 1)), axis=1)
mat_sp = np.append(mat_sp, values=np.zeros((1, mat_sp.shape[1])), axis=0)

另一种方法是创建大小为dim x dim的零矩阵，然后通过两个for循环从掩码创建一个巨大的索引向量。然后使用索引向量将矩阵乘法插入到零矩阵中。然而，这也非常缓慢。

任何能够有效解决问题的想法或见解都是很好的，因为稀疏矩阵乘积占用了非稀疏矩阵的2/3的时间。

---

使用@hjpaul的示例，我们得到以下比较代码

import numpy as np
dims = 400
def test_non_sparse():
vec = np.random.rand(1, dims)
a = vec.T * vec
def test_sparse():  
vec = np.random.rand(1, dims)
idx = np.flatnonzero(vec>0.75)
oprod = vec[:,idx].T * vec[:,idx]
vec_oprod = np.zeros((dims, dims))
vec_oprod[idx[:,None], idx] = oprod

if __name__ == '__main__':
import timeit
print('Non sparse:',timeit.timeit("test_non_sparse()", setup="from __main__ import test_non_sparse", number=10000))
print('Sparse:',timeit.timeit("test_sparse()", setup="from __main__ import test_sparse", number=10000))

当然，该代码根据向量的维度和零的数量进行了改进。超过300个维度和大约70%的零给出了适度的速度改进，该速度随着零元素和维度的数量而增加。如果矩阵和掩码一次又一次地相同，那么肯定有可能获得更大的加速。

(我做逻辑索引的错误是做idx而不是idx[:,None](

Answer 1

将一个矩阵插入另一个矩阵的最快方法是使用索引。

用你的外部产品来说明：

In [94]: a = np.array(((1,0,2,0)))
In [95]: idx = np.where(a>0)[0]
In [96]: idx
Out[96]: array([0, 2])
In [97]: a1 = a[idx]

浓缩阵列的外产物：

In [98]: a2 = a1[:,None]*a1
In [99]: a2
Out[99]: 
array([[1, 2],
[2, 4]])

设置结果数组，并使用块索引来识别a2值的去向：

In [100]: res = np.zeros((4,4),int)
In [101]: res[idx[:,None], idx] = a2
In [102]: res
Out[102]: 
array([[1, 0, 2, 0],
[0, 0, 0, 0],
[2, 0, 4, 0],
[0, 0, 0, 0]])

非密集阵列的直接外部：

In [103]: a[:,None]*a
Out[103]: 
array([[1, 0, 2, 0],
[0, 0, 0, 0],
[2, 0, 4, 0],
[0, 0, 0, 0]])
In [104]: np.outer(a,a)
Out[104]: 
array([[1, 0, 2, 0],
[0, 0, 0, 0],
[2, 0, 4, 0],
[0, 0, 0, 0]])

如果a是2d，(n，1(，那么这个外层可以写成np.dot(a.T,a)。dot涉及一个和，在这种情况下大小为1维。

我认为a必须非常稀疏才能从这种额外的索引工作中受益。对于scipy稀疏矩阵，我发现1%数量级的稀疏性具有任何速度优势，即使矩阵是预制的。

---

In [105]: from scipy import sparse
In [106]: A = sparse.csr_matrix(a)
In [107]: A
Out[107]: 
<1x4 sparse matrix of type '<class 'numpy.int64'>'
with 2 stored elements in Compressed Sparse Row format>
In [108]: A.A
Out[108]: array([[1, 0, 2, 0]], dtype=int64)
In [109]: A.T*A           # sparse matrix product, dot
Out[109]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
with 4 stored elements in Compressed Sparse Column format>
In [110]: _.A
Out[110]: 
array([[1, 0, 2, 0],
[0, 0, 0, 0],
[2, 0, 4, 0],
[0, 0, 0, 0]], dtype=int64)