使用 scipy 稀疏数组上的坐标有效地执行逐元素操作

Question

我想弄清楚如何在 scipy 稀疏数组（csc 格式）上有效地执行以下操作：

逐元素伪代码：

try:
    r = M[i,j] / (V[i] + V[j])
    if r.isfinite():
        M[i,j] = r
    else:
        # Leave at old value
        pass # e.g. 1e200/1e-200 -> inf
except ZeroDivisionError:
    # Leave at old value
    pass # e.g. 1e200/0 -> ZeroDivisionError

请注意，这会将零元素保留为 0，因此不会增加输入矩阵的密度

M

.

对于密集矩阵有一个简单的方法：

# V is a ndarray of shape (N), all values finite and >= 0.0
# M is a ndarray of shape (N,N), all values finite

rows_plus_cols = V[:, np.newaxis] + V[np.newaxis, :]
M = np.divide(M, rows_plus_cols, out=M, where=(rows_plus_cols != 0))

（请注意，此示例无法在所有情况下正确处理下溢。）

不幸的是，这种方法似乎不适用于稀疏矩阵，因为

rows_plus_cols

的中间计算会导致大小为

(N,N)

.

的密集矩阵

典型的输入 M 可能是 100k x 100k，具有约 1M 的非零条目。显然，在不建议具体化一个完整的密集矩阵的情况下。出于测试目的，以下内容与典型输入“足够接近”：

n = 100000
d = 10
# a.k.a. "mostly nonzero, but not always"
V = np.maximum((np.random.normal(size=n)+2), 0)

M = scipy.sparse.csc_array((np.random.normal(size=n*d), np.random.randint(n, size=(2, n*d))), shape=(n,n))
M.sum_duplicates()
M.prune()
M.sort_indices()

小问题：scipy 稀疏数组似乎不允许使用

where

子句进行元素划分。可以解决，但这不是主要问题。我知道您可以使用

M.data

访问 scipy csc 数组的底层数据点；我不知道如何在不具体化密集矩阵的情况下进行元素计算，这是坐标的函数。

建议？

Answer 1

这是您可以使用

coo_array

执行此操作的一种方法。

在此示例中，

M

是一个 8x8

coo_array

，

V

是一个长度为 8 的一维 NumPy 数组。如果

M

可能有重复的条目，您应该在执行以下操作之前调用

M.sum_duplicates()

。

In [106]: M
Out[106]: 
<8x8 sparse array of type '<class 'numpy.float64'>'
    with 16 stored elements in COOrdinate format>

In [107]: M.toarray()
Out[107]: 
array([[0.  , 0.  , 0.44, 0.  , 0.93, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.54, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.2 , 0.21, 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.93, 0.76, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.32, 0.88, 0.76],
       [0.91, 0.72, 0.  , 0.24, 0.87, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.74, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.13, 0.  , 0.  , 0.  , 0.  ]])

In [108]: M.data
Out[108]: array([0.44, 0.93, 0.54, 0.2 , 0.21, 0.93, 0.76, 0.32, 0.88, 0.76, 0.91, 0.72, 0.24, 0.87, 0.74, 0.13])

In [109]: M.row
Out[109]: array([0, 0, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 7], dtype=int32)

In [110]: M.col
Out[110]: array([2, 4, 6, 4, 5, 5, 6, 5, 6, 7, 0, 1, 3, 4, 4, 3], dtype=int32)

In [111]: V
Out[111]: array([0.5, 1. , 2. , 0. , 0. , 0.5, 1. , 0. ])

首先将分母形成为

V[M.row] + V[M.col]

，并将任何出现的0替换为1，因此不除以0，这些位置的值保持不变：

In [112]: denom = V[M.row] + V[M.col]

In [113]: denom[denom == 0] = 1

In [114]: denom
Out[114]: array([2.5, 0.5, 2. , 2. , 2.5, 0.5, 1. , 0.5, 1. , 1. , 1. , 1.5, 0.5, 0.5, 1. , 1. ])

我们称新数组为

Q

。

Q_data

将是

data

的

Q

数组：

In [115]: Q_data = M.data / denom

In [116]: Q_data
Out[116]: array([0.176, 1.86 , 0.27 , 0.1  , 0.084, 1.86 , 0.76 , 0.64 , 0.88 , 0.76 , 0.91 , 0.48 , 0.48 , 1.74 , 0.74 , 0.13 ])

用

coo_array

创建一个新的

Q_data

，使用与

M

相同的形状和相同的行和列索引：

In [117]: Q = sparse.coo_array((Q_data, (M.row, M.col)), shape=M.shape)

In [118]: Q.toarray()
Out[118]: 
array([[0.   , 0.   , 0.176, 0.   , 1.86 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.27 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.1  , 0.084, 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 1.86 , 0.76 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.64 , 0.88 , 0.76 ],
       [0.91 , 0.48 , 0.   , 0.48 , 1.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.13 , 0.   , 0.   , 0.   , 0.   ]])

执行密集版本的计算以验证它是否按预期工作：

In [119]: VV = np.add.outer(V, V)

In [120]: VV[VV == 0] = 1

In [121]: M.toarray() / VV
Out[121]: 
array([[0.   , 0.   , 0.176, 0.   , 1.86 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.27 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.1  , 0.084, 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 1.86 , 0.76 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.64 , 0.88 , 0.76 ],
       [0.91 , 0.48 , 0.   , 0.48 , 1.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.13 , 0.   , 0.   , 0.   , 0.   ]])

您可以获取该 ipython 会话的第 112、113、115 和 117 行并将它们包装为一个函数。