# Multiplying \`I\` to a sparse matrix drops explicit zeros

**URL:** <https://discourse.julialang.org/t/multiplying-i-to-a-sparse-matrix-drops-explicit-zeros/16130>\
**Category:** General Usage\
**Created:** [October 10, 2018, 6:46pm UTC](https://discourse.julialang.org/t/multiplying-i-to-a-sparse-matrix-drops-explicit-zeros/16130 "2018-10-10T18:46:57Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![wsshin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wsshin/32/360_2.png) [@wsshin](https://discourse.julialang.org/u/wsshin)\
**Post date:** [October 10, 2018, 6:46pm UTC](https://discourse.julialang.org/t/multiplying-i-to-a-sparse-matrix-drops-explicit-zeros/16130/1 "2018-10-10T18:46:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [October 10, 2018, 6:51pm UTC](https://discourse.julialang.org/t/multiplying-i-to-a-sparse-matrix-drops-explicit-zeros/16130/2 "2018-10-10T18:51:57Z")

</div>

> [@wsshin](#):
>
> It seems that UMFPACK uses different algorithms for LU factorization depending on whether the sparsity pattern of `A` is symmetric or not.

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.

---

<div class="post-metadata">

**Author:** ![wsshin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wsshin/32/360_2.png) [@wsshin](https://discourse.julialang.org/u/wsshin)\
**Post date:** [October 10, 2018, 7:58pm UTC](https://discourse.julialang.org/t/multiplying-i-to-a-sparse-matrix-drops-explicit-zeros/16130/3 "2018-10-10T19:58:34Z")

</div>

> [@kristoffer.carlsson](#):
>
> Almost, but it is Julia that checks `issymmetric` .

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.

> [@kristoffer.carlsson](#):
>
> FWIW, always wrap in `Symmetric` or `Hermitian` before calling the solver,  
> numerical round off typically makes things that are mathematically symmetric, non-symmetric anyway.

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…?

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [October 10, 2018, 8:17pm UTC](https://discourse.julialang.org/t/multiplying-i-to-a-sparse-matrix-drops-explicit-zeros/16130/4 "2018-10-10T20:17:24Z")

</div>

> [@wsshin](#):
>
> Not sure if this is correct

[https://github.com/JuliaLang/julia/blob/e20d1ee68683be888af5eb1cc5d3c9bdc18bd1c0/stdlib/SparseArrays/src/linalg.jl#L1254](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 ([GitHub - JuliaSparse/Pardiso.jl: Calling the PARDISO library from Julia](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.