我注意到,在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)结果的方法。