scipy稀疏矩阵中每行或每列的最大值索引

Question

scipy.sparse.coo_matrix.max 返回每行或列的最大值，给定一个轴。我想知道的不是最大值，而是每行或列最大值的索引。目前还没找到一种有效的方法来实现这个功能，所以我欢迎任何帮助。

- Jimmy C

---

你可以尝试这个：https://dev59.com/gGDVa4cB1Zd3GeqPc2Tl#9337071 它会单独处理每一行。 - Guillem Cucurull

Answer 1

我建议学习以下代码：

moo._min_or_max_axis

其中 moo 是一个 coo_matrix。

mat = mat.tocsc()  # for axis=0
mat.sum_duplicates()

major_index, value = mat._minor_reduce(min_or_max)
not_full = np.diff(mat.indptr)[major_index] < N
value[not_full] = min_or_max(value[not_full], 0)

mask = value != 0
major_index = np.compress(mask, major_index)
value = np.compress(mask, value)
return coo_matrix((value, (np.zeros(len(value)), major_index)),
                      dtype=self.dtype, shape=(1, M))

根据偏好轴，它更倾向于使用csc而不是csr。我还没有时间分析这个问题，但我猜想在计算中包括argmax应该是可能的。

---

这个建议可能行不通。关键是mat._minor_reduce方法，它在某种程度上进行了精细处理：

ufunc.reduceat(mat.data, mat.indptr[:-1])

这段内容涉及到IT技术，它将ufunc应用于data数组的块，使用indptr定义块。 np.sum，np.maxiumum都是可以使用这种方法的ufunc。我不知道有没有相当于argmax ufunc的功能。

通常，如果您想对csr矩阵（或csc的col）按“行”进行操作，则必须迭代行，这相对昂贵，或使用此ufunc.reduceat在平面mat.data向量上执行相同的操作。 在numpy中对分区索引执行group argmax / argmin尝试执行argmax.reduceat。那里的解决方案可能适用于稀疏矩阵。

- hpaulj

---

是的，https://dev59.com/mX3aa4cB1Zd3GeqPYh7v是逻辑等效的;只需应用于CSR的“indices”和“indptr”。 - joeln

Answer 2

从scipy版本0.19开始，csr_matrix和csc_matrix都支持argmax()和argmin()方法。

- S.Mohsen sh

Answer 3

如果A是您的scipy.sparse.coo_matrix，那么您可以按如下方式获取最大值的行和列：

I=A.data.argmax()
maxrow = A.row[I]
maxcol=A.col[I]

要获取每行最大值的索引，请参见下面的编辑：

from scipy.sparse import coo_matrix
import numpy as np
row  = np.array([0, 3, 1, 0])
col  = np.array([0, 2, 3, 2])
data = np.array([-3, 4, 11, -7])
A= coo_matrix((data, (row, col)), shape=(4, 4))
print A.toarray()

nrRows=A.shape[0]
maxrowind=[]
for i in range(nrRows):
    r = A.getrow(i)# r is 1xA.shape[1] matrix
    maxrowind.append( r.indices[r.data.argmax()] if r.nnz else 0)
print maxrowind 

r.nnz 是显式存储值的计数（即非零值）

- xecafe

---

这样做会不会只产生一个单一的值，而不是每行或每列都有一个值？ - Jimmy C

Answer 4

在@hpaulj和@joeln的答案基础上，结合numpy中对分区索引求最大值/最小值的群组提供的代码，此函数将计算CSR列的argmax或CSC行的argmax：

import numpy as np
import scipy.sparse as sp

def csr_csc_argmax(X, axis=None):
    is_csr = isinstance(X, sp.csr_matrix)
    is_csc = isinstance(X, sp.csc_matrix)
    assert( is_csr or is_csc )
    assert( not axis or (is_csr and axis==1) or (is_csc and axis==0) )

    major_size = X.shape[0 if is_csr else 1]
    major_lengths = np.diff(X.indptr) # group_lengths
    major_not_empty = (major_lengths > 0)

    result = -np.ones(shape=(major_size,), dtype=X.indices.dtype)
    split_at = X.indptr[:-1][major_not_empty]
    maxima = np.zeros((major_size,), dtype=X.dtype)
    maxima[major_not_empty] = np.maximum.reduceat(X.data, split_at)
    all_argmax = np.flatnonzero(np.repeat(maxima, major_lengths) == X.data)
    result[major_not_empty] = X.indices[all_argmax[np.searchsorted(all_argmax, split_at)]]
    return result

如果任何行（CSR）或列（CSC）完全稀疏（即，在X.eliminate_zeros()之后完全为零），则它将返回-1作为argmax。

- elhuhdron

Answer 5

最新版本的numpy_indexed软件包（免责声明：我是该软件包的作者）可以以高效而优雅的方式解决这个问题：

import numpy_indexed as npi
col, argmax = group_by(coo.col).argmax(coo.data)
row = coo.row[argmax]

在这里，我们按列进行分组，因此它是列上的argmax；交换行和列将为您提供行上的argmax。

- Eelco Hoogendoorn

Answer 6

正如其他人提到的那样，现在有内置的argmax()用于scipy.sparse矩阵。然而，我发现它对于大矩阵来说速度相当慢，所以我看了一下源代码。逻辑非常聪明，但它包含一个python循环，使事情变慢。例如，将源代码缩减为每行的argmax（牺牲所有通用性、形状检查等等，以简化代码），并使用numba进行装饰可以获得一些良好的速度改进。

这是函数：

import numpy as np
from numba import jit


def argmax_row_numba(X):
    return _argmax_row_numba(X.shape[0], X.indptr, X.data, X.indices)

@jit(nopython=True)
def _argmax_row_numba(shape, indptr, data, indices):
    # prep an array to hold the indices
    ret = np.zeros(shape)
    # figure out which lines actually contain data
    nz_lines, = np.diff(indptr).nonzero()
    # loop through the lines
    for i in nz_lines:
        p, q = indptr[i: i + 2]
        line_data = data[p: q]
        line_indices = indices[p: q]
        am = np.argmax(line_data)
        ret[i] = line_indices[am]

    return ret

生成测试矩阵：

from scipy.sparse import random
size = 10000
m = random(m=size, n=size, density=0.0001, format="csr")
n_vals = m.data.shape[0]
m.data = np.random.random(size=n_vals).astype("float")


# the original scipy implementation reformatted to return a np.array
maxima1 = np.squeeze(np.array(m.argmax(axis=1)))
# calling the numba version
maxima2 = argmax_row_numba(m)

# Check that the results are the same
print(np.allclose(maxima1, maxima2))
# True

计时结果：

%timeit m.argmax(axis=1)
# 30.1 ms ± 246 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit argmax_row_numba(m)
# 211 µs ± 1.04 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

- L_W