问题描述
我注意到，在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）的结果。
解决方法
TLDR认为，稀疏矩阵是一项巨大的性能优势，因为您不必检查0*x是什么。如果99.9％的条目为零（通常是这种情况），那么检查另一个矩阵中的inf值会做很多额外的工作。
使用Python稀疏：
使两个数组都变得稀疏：
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）问题。
每个Julia库都为功能提供了自己的一组方法，请参见：
methods(*)
检查您拥有多少个*方法（甚至查看methods(*,SparseArrays)仅查看SparseArrays库带来的方法）。
由于Julia中的运算符只是一个简单的函数，实际上，对于给定的一组参数，您可以看到其实现。在您的代码中尝试：
@less *(A,B)
# or to actaully open the code in an editor
@edit *(A,B)
您会很快发现，SparseArrays的实现在执行操作时只是省略了零。由于在Julia 0 * Inf == NaN中，您观察到的行为可能应该被视为错误-执行乘法时，库可能应该检查这种极端情况。另一方面，这个角落可能会影响整个库的性能。