Multiplying `I` to a sparse matrix drops explicit zeros

I just noticed that multiplying I to a sparse matrix drops the explicit zeros stored in the sparse matrix. For example, A is a 2×2 zero matrix, but it has an explicit zero entry at (1,1):

julia> VERSION
v"0.7.1-pre.0"

julia> using LinearAlgebra, SparseArrays

julia> A = spzeros(2,2); A[1,1] = 1; A[1,1] = 0;

julia> A
2×2 SparseMatrixCSC{Float64,Int64} with 1 stored entry:
  [1, 1]  =  0.0

Now, if we multiply I to A, this explicit zero is dropped:

julia> B = I*A;

julia> B
2×2 SparseMatrixCSC{Float64,Int64} with 0 stored entries

I am not sure if this behavior is intended. Even if it is intended, I am not sure if such behavior is desired, provided that we can always drop zeros at will using dropzeros!.

Currently I have a function where a sparse matrix is created in the form of A = C * D * transpose(C), where C is sparse and D is diagonal. When D is the identity matrix, D = I is used. The matrix A is always symmetric for all D. However, for D = I, the C*D part of A drops the explicit zeros of C, making the sparsity pattern of C*D different from the transpose of the sparsity pattern of transpose(C). As a result, the sparsity pattern of A = C * D * transpose(C) is no longer symmetric for D = I.

It seems that UMFPACK uses different algorithms for LU factorization depending on whether the sparsity pattern of A is symmetric or not. As a result, even though A is mathematically symmetric for D = I, UMFPACK does not consider it symmetric and ends up using the slower algorithm.

1 Like

Almost, but it is Julia that checks issymmetric.

FWIW, always wrap in Symmetric or Hermitian before calling the solver,
numerical round off typically makes things that are mathematically symmetric, non-symmetric anyway.

1 Like

Not sure if this is correct. As you said, a mathematically symmetric matrix could be numerically nonsymmetric due to rounding errors, but still I experience a huge performance degradation depending on whether D = I or D = sparse(Diagonal(ones(m)), where A is m×m. So, I thought what matters is not symmetry of A itself, but symmetry of the sparsity pattern of A. I thought this makes sense because UMFPACK performs symbolic factorization based on the sparsity pattern.

Thanks for the suggestion. I think this tip will be useful when numerical symmetry needs to be enforced. But again, even if the matrix is numerically nonsymmetric, I experience a huge performance gain in LU factorization by UMFPACK if the matrix has a symmetric sparsity pattern. This means there are cases where one can perform LU factorization faster by putting zeros explicitly in a sparse matrix to achieve a symmetric sparsity pattern, even though the matrix is nonsymmetric. So, automatically dropping those explicit zeros when multiplying I may not be a good idea, because it breaks symmetry in the sparsity pattern…?

https://github.com/JuliaLang/julia/blob/e20d1ee68683be888af5eb1cc5d3c9bdc18bd1c0/stdlib/SparseArrays/src/linalg.jl#L1254 is what I meant but sure, structurally symmetric matrices are common options for sparse solvers (https://github.com/JuliaSparse/Pardiso.jl#setting-the-matrix-type).

I agree that dropping sparse stored zeros is typically unwanted. Earlier in Julia, stored zeros were aggressively dropped (similar to Matlab) but as things evolved, they have been kept in more and more operations.