scipy 稀疏矩阵中每行或每列的 Argmax
Argmax of each row or column in scipy sparse matrix
scipy.sparse.coo_matrix.max
returns 给定一个轴,每行或每列的最大值。我想知道的不是值,而是每行或列的最大值的索引。我还没有找到一种有效的方法,所以我很乐意接受任何帮助。
如果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
是显式存储值(即非零值)的计数
我建议研究
的代码
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])
也就是将ufunc
应用于矩阵data
数组的块,使用indptr
来定义块。 np.sum
、np.maxiumum
是 ufunc
这行得通的地方。我不知道等效的 argmax
ufunc.
一般来说,如果你想通过 'row' 为 csr 矩阵(或 csc 的 col)做事,你要么必须遍历行,这是相对昂贵的,要么使用这个 ufunc.reduceat
在平面 mat.data
向量上做同样的事情。
group argmax/argmin over partitioning indices in numpy
尝试执行 argmax.reduceat
。那里的解决方案可能适用于稀疏矩阵。
最新发布的numpy_indexed包(免责声明:我是它的作者)可以高效优雅地解决这个问题:
import numpy_indexed as npi
col, argmax = group_by(coo.col).argmax(coo.data)
row = coo.row[argmax]
这里我们按列分组,所以它是列上的 argmax;交换 row 和 col 将为您提供行上的 argmax。
扩展@hpaulj 和@joeln 的答案并按照建议使用group argmax/argmin over partitioning indices in 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
它returns -1 对于完全稀疏的任何行 (CSR) 或列 (CSC) 的 argmax(即,在 X.eliminate_zeros()
之后完全为零)。
从scipy版本0.19开始,csr_matrix
和csc_matrix
都支持argmax()
和argmin()
方法。
正如其他人提到的那样,scipy.sparse
矩阵现在内置了 argmax()
。但是,我发现它对于大型矩阵来说非常慢,所以我看了一下 the source code。逻辑非常聪明,但它包含一个 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)
scipy.sparse.coo_matrix.max
returns 给定一个轴,每行或每列的最大值。我想知道的不是值,而是每行或列的最大值的索引。我还没有找到一种有效的方法,所以我很乐意接受任何帮助。
如果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
是显式存储值(即非零值)的计数
我建议研究
的代码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])
也就是将ufunc
应用于矩阵data
数组的块,使用indptr
来定义块。 np.sum
、np.maxiumum
是 ufunc
这行得通的地方。我不知道等效的 argmax
ufunc.
一般来说,如果你想通过 'row' 为 csr 矩阵(或 csc 的 col)做事,你要么必须遍历行,这是相对昂贵的,要么使用这个 ufunc.reduceat
在平面 mat.data
向量上做同样的事情。
group argmax/argmin over partitioning indices in numpy
尝试执行 argmax.reduceat
。那里的解决方案可能适用于稀疏矩阵。
最新发布的numpy_indexed包(免责声明:我是它的作者)可以高效优雅地解决这个问题:
import numpy_indexed as npi
col, argmax = group_by(coo.col).argmax(coo.data)
row = coo.row[argmax]
这里我们按列分组,所以它是列上的 argmax;交换 row 和 col 将为您提供行上的 argmax。
扩展@hpaulj 和@joeln 的答案并按照建议使用group argmax/argmin over partitioning indices in 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
它returns -1 对于完全稀疏的任何行 (CSR) 或列 (CSC) 的 argmax(即,在 X.eliminate_zeros()
之后完全为零)。
从scipy版本0.19开始,csr_matrix
和csc_matrix
都支持argmax()
和argmin()
方法。
正如其他人提到的那样,scipy.sparse
矩阵现在内置了 argmax()
。但是,我发现它对于大型矩阵来说非常慢,所以我看了一下 the source code。逻辑非常聪明,但它包含一个 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)