Why does matrix multiplication with sparse matrices differ from dense ones if Inf is involved?

Question

I noticed, that in Julia and Python the result of matrix-matrix multiplication is different for sparse and dense arrays, if infinity is involved, see sample code:

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

same in 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}")

prints out

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]]

maybe this is due to the same underling subroutine, haven't checked this with other languages so far. It looks like the product with some unstored entry is always set to zero and therefore no NaN is produced as in case of 0 * inf, which comes up in dense arrays.

I haven't found any documentation mentioning this behavior. Does anyboy know if this is common or somewhere agreed upon? Especially in Julia I would expect, from a mathematical point of view, that dropzeros does not alter the result, which is not the case here. scipy on the other hand drops zeros automatically, therefore I found no way the reproduce the result of the first Julia multiplication (A * B).

Answer 1

The TLDR is that sparse matrices are a massive performance win because you don't have to check what 0*x is. If 99.9% of your entries are zeroes (this is often the case), then checking for inf values in the other matrix is doing a lot of extra work.

Answer 2

With Python sparse:

Make both arrays 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 this case, the result is also sparse:

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.]])

The sparse result makes it more obvious that the calculation has just used the non-zero elements of A, skipping anything involving the 0s. Thus it bypasses the 0*np.nan (or 0*np.inf) issue.

Answer 3

Each Julia library brings its own set of methods to the functionality, see:

methods(*)

to check what how many * methods you have (or even methods(*, SparseArrays) to see only those brought you by SparseArrays library).

Since operator in Julia is just a simple function, in fact, for a given set of parameters you can see its implementation. In your code try:

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

You will find out quickly that the implementation by SparseArrays are simply omitting zeros when performing the operation. Since in Julia 0 * Inf == NaN the behavior you are observing probably should be considered as a bug - the library perhaps should check this corner case when performing the multiplication. On the other hand this corner could affect the performance the entire library.