numpy:有效地添加矩阵的行

2024-03-20

我有一个矩阵。

mat = array([
   [ 0,  1,  2,  3],
   [ 4,  5,  6,  7],
   [ 8,  9, 10, 11]
   ])

我想获得某些索引处的行的总和:例如。

ixs = np.array([0,2,0,0,0,1,1])

我知道我可以将答案计算为:

mat[ixs].sum(axis=0)
> array([16, 23, 30, 37])

问题是ixs可能很长,我不想使用所有内存来创建中间产品mat[ixs],只是用总和再次减少它。

我还知道我可以简单地计算索引并使用乘法代替。

np.bincount(ixs, minlength=mat.shape[0).dot(mat)
> array([16, 23, 30, 37])

但如果我的 ix 稀疏,那会很昂贵。

我了解 scipy 的稀疏矩阵,我想我可以使用它们,但我更喜欢纯 numpy 解决方案,因为稀疏矩阵在多种方面受到限制(例如仅是二维)

那么,在这种情况下,是否有一种纯粹的 numpy 方法来合并索引和求和减少?

结论:

感谢 Divakar 和 hpaulj 非常详尽的回复。我所说的“稀疏”是指大多数值range(w.shape[0])不在 ix 中。使用这个新定义(并且具有更现实的数据大小,我重新运行了 Divakar 测试,并添加了一些新函数:

rng = np.random.RandomState(1234)
mat = rng.randn(1000, 500)
ixs = rng.choice(rng.randint(mat.shape[0], size=mat.shape[0]/10), size=1000)

# Divakar's solutions
In[42]: %timeit org_indexing_app(mat, ixs)
1000 loops, best of 3: 1.82 ms per loop
In[43]: %timeit org_bincount_app(mat, ixs)
The slowest run took 4.07 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 177 µs per loop
In[44]: %timeit indexing_modified_app(mat, ixs)
1000 loops, best of 3: 1.81 ms per loop
In[45]: %timeit bincount_modified_app(mat, ixs)
1000 loops, best of 3: 258 µs per loop
In[46]: %timeit simply_indexing_app(mat, ixs)
1000 loops, best of 3: 1.86 ms per loop
In[47]: %timeit take_app(mat, ixs)
1000 loops, best of 3: 1.82 ms per loop
In[48]: %timeit unq_mask_einsum_app(mat, ixs)
10 loops, best of 3: 58.2 ms per loop 
# hpaulj's solutions
In[53]: %timeit hpauljs_sparse_solution(mat, ixs)
The slowest run took 9.34 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 524 µs per loop
%timeit hpauljs_second_sparse_solution(mat, ixs)
100 loops, best of 3: 9.91 ms per loop
# Sparse version of original bincount solution (see below):
In[60]: %timeit sparse_bincount(mat, ixs)
10000 loops, best of 3: 71.7 µs per loop

获胜者,冠军在本例中是 bincount 解决方案的稀疏版本。

def sparse_bincount(mat, ixs):
    x = np.bincount(ixs)
    nonzeros, = np.nonzero(x)
    x[nonzeros].dot(mat[nonzeros])

替代方案bincount is add.at:

In [193]: mat
Out[193]: 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
In [194]: ixs
Out[194]: array([0, 2, 0, 0, 0, 1, 1])

In [195]: J = np.zeros(mat.shape[0],int)
In [196]: np.add.at(J, ixs, 1)
In [197]: J
Out[197]: array([4, 2, 1])

In [198]: np.dot(J, mat)
Out[198]: array([16, 23, 30, 37])

我认为,通过稀疏性,你的意思是ixs可能不包括所有行,例如,ixs没有 0:

In [199]: ixs = np.array([2,1,1])
In [200]: J=np.zeros(mat.shape[0],int)
In [201]: np.add.at(J, ixs, 1)
In [202]: J
Out[202]: array([0, 2, 1])
In [203]: np.dot(J, mat)
Out[203]: array([16, 19, 22, 25])

J仍然有mat.shape[0]形状。但是add.at应按长度缩放ixs.

稀疏解决方案看起来像:

制作一个稀疏矩阵ixs看起来像:

In [204]: I
Out[204]: 
array([[1, 0, 1, 1, 1, 0, 0],
       [0, 0, 0, 0, 0, 1, 1],
       [0, 1, 0, 0, 0, 0, 0]])

对行求和;稀疏通过矩阵乘法来实现这一点,例如:

In [205]: np.dot(I, np.ones((7,),int))
Out[205]: array([4, 2, 1])

然后做我们的点:

In [206]: np.dot(np.dot(I, np.ones((7,),int)), mat)
Out[206]: array([16, 23, 30, 37])

或者在稀疏代码中:

In [225]: J = sparse.coo_matrix((np.ones_like(ixs,int),(np.arange(ixs.shape[0]), ixs)))
In [226]: J.A
Out[226]: 
array([[1, 0, 0],
       [0, 0, 1],
       [1, 0, 0],
       [1, 0, 0],
       [1, 0, 0],
       [0, 1, 0],
       [0, 1, 0]])
In [227]: J.sum(axis=0)*mat
Out[227]: matrix([[16, 23, 30, 37]])

sparse,当转换自coo to csr对重复项求和。我可以利用这一点

In [229]: J = sparse.coo_matrix((np.ones_like(ixs,int), (np.zeros_like(ixs,int), ixs)))
In [230]: J
Out[230]: 
<1x3 sparse matrix of type '<class 'numpy.int32'>'
    with 7 stored elements in COOrdinate format>
In [231]: J.A
Out[231]: array([[4, 2, 1]])
In [232]: J*mat
Out[232]: array([[16, 23, 30, 37]], dtype=int32)
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

numpy:有效地添加矩阵的行 的相关文章

随机推荐

  • JTable 行之间的标题

    我想知道如何实现这样的事情 row 1 row 2 HEADLINE比如 2011 年 1 月 27 日 row 3 row 4 HEADLINE比如 2011 年 1 月 28 日 row 5 这些 假 行应该只是信息行而不是真实行 我找
  • 如何在C#中获取当前产品版本?

    如何在 C 中以编程方式获取当前产品版本 My code VersionNumber System Reflection Assembly GetExecutingAssembly GetName Version ToString 我得到
  • 什么是textview.setTextSize()?

    在我的 Android 游戏中 有一个文本视图 我使用以下代码设置文本大小 textview setTextSize 30 30以像素为单位 但它到底是什么 是字符的高度还是宽度 是序数号吗 在android java或kotlin 中设置
  • WPF - 宽度和高度必须为非负数

    我的应用程序遇到了一个奇怪的异常 它一直运行良好 直到我升级了我的开发机器 相同的操作系统 VS2010 并尝试再次调试软件 这是我得到的异常 宽度和高度必须为非负数 堆栈跟踪在这里 at System Windows Size ctor
  • 使用 babel 后,'Symbol' 在 IE 中未定义

    我有一个reactjs使用 ES6 标准编写的应用程序 我使用webpack来建造它 这webpack加载js模块使用babel loader 具体来说 我使用以下版本的包 email protected cdn cgi l email p
  • OpenCV - 如何从 Canny 函数的结果中提取边缘?

    我在 OpenCV 中使用 Canny 函数 如下所示 Mat detected edges GetImage Canny detected edges detected edges 20 20 3 kernel size 我的问题是这个函
  • 在 C++ 中将对象数组设置为 null

    假设我有一个 C 中 Foo 类型的对象数组 Foo array 10 在 Java 中 我可以简单地通过以下方式将此数组中的对象设置为 null array 0 null the first one 我怎样才能在 C 中做到这一点 使用指
  • PHPStorm中Code Sniffer触发的Xdebug

    我在安装 PHPStorm 时 xdebug 和 Code Sniffer 都工作得很好 但真正烦人的部分是 调试器现在似乎将 Code Sniffer 错误视为断点 并中断代码 让我知道样式警告 尝试测试代码 如何防止 Code Snif
  • Firestore存储大小限制如何存储大型数组

    我有一个收藏users userID followers该用户 ID 是 firebase UID 因此长度为 29 个字节 字符串大小的计算方式为 UTF 8 编码字节数 1 在每个用户文档中 我都有一个名为 follower 的数组和另
  • 解析 DateFormat 时的 Java 时区

    我的代码解析日期如下 String ALT DATE TIME FORMAT yyyy MM dd T HH mm ss SSSZ SimpleDateFormat sdf new SimpleDateFormat ALT DATE TIM
  • android AppWidget 未添加到 Lollipop 上的主屏幕

    我开发了一个应用程序 可以在主屏幕小部件上显示新闻源 由于以下情况 在 Lollipop 之前的 Android 设备上一切正常 用户进入启动器的小部件屏幕以选择 添加特定的小部件 用户单击 MyNewsWidget 以添加到其主屏幕 调用
  • 让 div 占据另一个 div 后剩余的所有空间

    我有两个并排的 div 第一个包含一个可能相当长的文本字段 另一个包含一个很短的数字 我需要第一个 div 占据所有可用空间 而无需拉伸父级并在必要时进行剪切 Ant 它应该考虑第二个 div 的宽度 因此 如果文本的长度很短 那么两个 d
  • linuxrc 的用途是什么以及 rootfs 中是否需要它?

    Question 我的问题是 什么是linuxrc做 我需要它吗 rootfs 和使用有什么关系吗systemd vs initd 背景 我目前正在尝试建立一个rootfs适用于使用 Yocto 的 ARM 7 处理器 我对原始 BSP 项
  • docx4j - 删除 wml P 元素

    我正在使用 docx4j 来处理 Microsoft Word 模板 我想知道如何删除或隐藏模板中的 P 元素 我能够遍历代码来获取特定的 P 元素 现在我需要知道如何删除或隐藏该 P 元素 有人可以帮忙吗 我使用以下代码获取所有 P 元素
  • 按 Enter 键提交搜索?

    当有人按下 回车 键时 需要做什么才能提交此表单
  • Visual Studio 弹出窗口:“操作无法完成”

    当我尝试打开一个项目时 无论是本地项目还是在团队基础服务器 https en wikipedia org wiki Team Foundation Server TFS 我得到一个模态窗口告诉我 操作无法完成 未指定的错误 或者相同的消息
  • 内存怎么这么大?

    我有一个 1000x1500 像素位图 我想在 Android 中制作一个可变副本 当我运行以下代码时 int width original getWidth 1000px int height original getHeight 150
  • C++中的restrict关键字是什么意思?

    我总是不确定 restrict关键字在C 中意味着什么 这是否意味着赋予函数的两个或多个指针不重叠 还有什么意思呢 在他的论文中 内存优化 https web archive org web 20160422113037 http www
  • 如何在 sunos 中获取附加到特定端口的进程 ID

    我正在尝试在 SunOS 上使用端口 7085 连接进程 我尝试执行以下命令 netstat ntlp grep 7085没有返回任何东西 netstat anop grep 7085也尝试过这个 此开关在 SunOs 中无效 我得到以下输
  • numpy:有效地添加矩阵的行

    我有一个矩阵 mat array 0 1 2 3 4 5 6 7 8 9 10 11 我想获得某些索引处的行的总和 例如 ixs np array 0 2 0 0 0 1 1 我知道我可以将答案计算为 mat ixs sum axis 0