如果涉及 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) 的结果。

Answer 1

TLDR 是稀疏矩阵是一个巨大的性能胜利，因为您不必检查是什么0*x。如果 99.9% 的条目是零（通常是这种情况），那么检查inf另一个矩阵中的值需要做很多额外的工作。