为什么涉及Inf时，稀疏矩阵乘法与密集矩阵乘法不同？

Question

我注意到，在Julia和Python中，如果涉及无穷大，在稀疏数组和密集数组之间进行矩阵乘法的结果是不同的。请查看以下示例代码：

julia> using SparseArrays

julia> using LinearAlgebra

julia> A = spdiagm(0 => [0, 1])
2×2 SparseMatrixCSC{Int64,Int64} with 2 stored entries:
  [1, 1]  =  0
  [2, 2]  =  1

julia> B = [1 Inf; 1 2]
2×2 Array{Float64,2}:
 1.0  Inf
 1.0   2.0

julia> A * B
2×2 Array{Float64,2}:
 0.0  NaN
 1.0    2.0

julia> Array(A) * B
2×2 Array{Float64,2}:
 0.0  NaN
 1.0  NaN

julia> dropzeros(A) * B
2×2 Array{Float64,2}:
 0.0  0.0
 1.0  2.0

在Python中相同

from scipy.sparse import diags
import numpy as np

A = diags([0, 1])
B = np.array([[1, np.inf], [1, 2]])
print(f"A=\n{A}")
print(f"B=\n{B}")
print(f"sparse mul:\n{A @ B}")
print(f"dense mul:\n{A.toarray() @ B}")

输出

A=
  (1, 1)    1.0
B=
[[ 1. inf]
 [ 1.  2.]]
sparse mul:
[[0. 0.]
 [1. 2.]]
/home/.../TestSparseInf.py:9: RuntimeWarning: invalid value encountered in matmul
  print(f"dense mul:\n{A.toarray() @ B}")
dense mul:
[[ 0. nan]
 [ 1. nan]]

也许这是由于相同的子程序所致，到目前为止还没有检查其他语言。 看起来具有一些未存储条目的乘积始终设置为零，因此不会产生NaN，就像在稠密数组中出现的0 * inf一样。

我没有找到任何提到这种行为的文档。有人知道这是否常见或是否达成共识吗？ 特别是在Julia中，从数学角度来看，我希望dropzeros不会改变结果，但在这种情况下并非如此。 另一方面，scipy会自动删除零，因此我找不到重现第一个Julia乘法（A * B）结果的方法。

- Hugou

Answer 1

简而言之，稀疏矩阵具有巨大的性能优势，因为您不必检查0*x是什么。如果您的条目中有99.9％是零（这通常是情况），那么在另一个矩阵中检查inf值会做很多额外的工作。

- Oscar Smith

Answer 2

使用Python的sparse库：

将两个数组均转换为稀疏矩阵：

In [641]: A = sparse.diags([0,1]).tocsr()                                                            
In [642]: B = sparse.csr_matrix(np.array([[1, np.nan],[1,2]]))                                       
In [643]: A                                                                                          
Out[643]: 
<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 1 stored elements in Compressed Sparse Row format>
In [644]: B                                                                                          
Out[644]: 
<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in Compressed Sparse Row format>

在这种情况下，结果的稀疏性也很高。

In [645]: A@B                                                                                        
Out[645]: 
<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements in Compressed Sparse Row format>
In [646]: _.A                                                                                        
Out[646]: 
array([[0., 0.],
       [1., 2.]])

稀疏结果更明显地表明计算仅使用了A的非零元素，跳过与0有关的任何内容。因此，它避免了0*np.nan（或0*np.inf）问题。

- hpaulj

---

我假设了绕过。我很烦恼，因为这种行为甚至没有在文档中提到。 - Hugou

稀疏矩阵乘法代码中有一条注释引用了一篇数学论文。该算法多年前为有限差分和元素问题而开发。 - hpaulj

此代码的早期探索，请参见 https://stackoverflow.com/a/37540149/901925。 - hpaulj

Answer 3

每个 Julia 库都提供了自己的方法来实现相应功能，详见：

methods(*)

要检查有多少个*方法（或者只看SparseArrays库提供的方法，可以使用methods(*, SparseArrays)）。

由于Julia中的运算符只是一个简单的函数，实际上，对于给定的一组参数，您可以查看其实现。在您的代码中尝试：

@less *(A, B)
# or to actaully open the code in an editor
@edit *(A, B)

你很快就会发现，SparseArrays 的实现在执行操作时只是简单地省略了零。由于在 Julia 中 0 * Inf == NaN，因此您观察到的行为可能应该被视为错误 - 库在执行乘法时可能应该检查这种情况。另一方面，这个角落可能会影响整个库的性能。

- Przemyslaw Szufel

---

谢谢您的意见。如果需要考虑（可选）检查，我将前往Julia社区。 - Hugou