# Fast conversion from Matrix{Union{Missing, Float64}} to Sparse

**URL:** <https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161>\
**Category:** Performance\
**Tags:** question, sparse\
**Created:** [December 29, 2023, 6:11pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161 "2023-12-29T18:11:24Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [December 29, 2023, 6:11pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/1 "2023-12-29T18:11:24Z")

</div>

I have a matrix that contains lots of missing values. I am toying with the idea of using sparsity to represent those missing values instead. However, one hurdle to easy conversion from Matrix to SparseMatrix is that the non-missing values of the matrix are allowed to be zero (but rarely) and I would like to keep those zeros. Instead, I wrote the following function. I am wondering if anyone else has suggestions for improvement.

```julia
using SparseArrays, Missings

function sparsenonmissing(A::AbstractMatrix{T}) where {T}
	N,K = size(A)
	numvals = sum(!ismissing,A)
	rows = Int64[]
	cols = Int64[]
	vals = nonmissingtype(T)[]

	sizehint!(rows, numvals)
	sizehint!(cols, numvals)
	sizehint!(vals, numvals)

	for ind in eachindex(IndexCartesian(),A)
		x = A[ind]
		if !ismissing(x)
			push!(rows,ind[1])
			push!(cols,ind[2])
			push!(vals,x)
		end
	end
	return sparse(rows,cols,vals,N,K)
end

A = missings(Float64, 500,10000)
A[1:5:end] .= randn(Float64,length(1:5:length(A)))
B = sparsenonmissing(A)

```

The degree of sparsity in the MWE is approximately what I encounter (the patterns actually display much more clustering).

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [December 31, 2023, 1:51am UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/2 "2023-12-31T01:51:34Z")

</div>

I’ve been thinking this one over and, as far as performance goes, I think this is actually a pretty decent version. I can think of a few ways to get it done in fewer lines, but they’re a bit hacky and probably slower.

---

<div class="post-metadata">

**Author:** ![Tortar](https://avatars.discourse-cdn.com/v4/letter/t/6bbea6/32.png) [@Tortar](https://discourse.julialang.org/u/Tortar)\
**Post date:** [December 31, 2023, 3:01am UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/3 "2023-12-31T03:01:43Z")

</div>

This is 2x faster on my pc:

```julia
function sparsenonmissing2(A::AbstractMatrix{T}) where {T}
    N,K = size(A)
    numvals = sum(!ismissing,A)
    rows = Vector{Int64}(undef, numvals)
    cols = Vector{Int64}(undef, numvals)
    vals = Vector{nonmissingtype(T)}(undef, numvals)
    i = 0
    @inbounds for ind in CartesianIndices(A)
        x = A[ind]
        if !ismissing(x)
            i += 1
            rows[i] = ind[1]
            cols[i] = ind[2]
            vals[i] = x
        end
    end
    return sparse(rows,cols,vals,N,K)
end

```

(adding @inbounds to your version instead does nothing)

---

<div class="post-metadata">

**Author:** ![Zentrik](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zentrik/32/35409_2.png) [@Zentrik](https://discourse.julialang.org/u/Zentrik)\
**Post date:** [December 31, 2023, 8:16am UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/4 "2023-12-31T08:16:34Z")

</div>

I think inbounds needs to be on the line where you do the indexing not on a for loop.

---

<div class="post-metadata">

**Author:** ![Tortar](https://avatars.discourse-cdn.com/v4/letter/t/6bbea6/32.png) [@Tortar](https://discourse.julialang.org/u/Tortar)\
**Post date:** [December 31, 2023, 12:29pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/5 "2023-12-31T12:29:56Z")

</div>

Should be the same because if I’m not mistaken it is propagated to each expression inside the for loop this way. To confirm this, I tried to do what you say on both version and nothing changes

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 31, 2023, 12:40pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/6 "2023-12-31T12:40:42Z")

</div>

> [@tbeason](#):
>
> I am toying with the idea of using sparsity to represent those missing values instead.

What operations do you want to support on the array? (Most of the SparseMatrixCSC functionality is oriented towards linear algebra, which is useless if there are missing values.)

> [@Tortar](#):
>
> This is 2x faster on my pc:

If you care about performance you should construct the CSC format directly, rather than going through IJV format first.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [December 31, 2023, 12:56pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/7 "2023-12-31T12:56:44Z")

</div>

This method seems to get a little more juice:

```julia
function sparsenonmissing3(A::AbstractMatrix{T}) where {T}
    B = similar(sparse((!ismissing).(A)), nonmissingtype(eltype(A)))
    cptr = SparseArrays.getcolptr(B)
    c = 1
    cend = cptr[c+1]
    @inbounds for i in 1:nnz(B)
        if i == cend
            c += 1
            cend = cptr[c+1]
        end
        nonzeros(B)[i] = A[rowvals(B)[i],c]
    end
    return B
end

```

and it doesn’t look so bad either.

BTW, the `A` in the OP has density 0.2 which is a little high for sparse advantage (considering the overhead of CSC).

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [January 1, 2024, 3:19pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/8 "2024-01-01T15:19:28Z")

</div>

This isn’t fast, but it’s fast to write

```julia
sparse(coalesce.(A,0))

```

Using a mutated and mutant `coalesce!()` goes a little faster

```julia

function coalesce!(m::Matrix{Union{Missing, T}},val::T) where T<:Real
	for i in eachindex(m)
		ismissing(m[i]) && (m[i]=val)
	end
	m
end

sparsenonmissing3(A)==sparse(coalesce.(A,0))==sparse(coalesce!(A,0.))

```

or better

```julia

function SMCSC(A)
	nzA=(!ismissing).(A)
	r,c=size(A)
	cptr=[0;cumsum(count.(eachcol(nzA)))].+1
	nzval= A[nzA]
	rvals=[first(ci.I) for ci in findall(nzA)]
	SparseMatrixCSC(r,c,cptr,rvals,nzval)
end

```

```julia
function cptr_(nzA)
	cptr=Vector{Int64}(undef,10001)
	cptr[1]=1
	for (col,c) in enumerate(eachcol(nzA))
		cptr[col+1]=cptr[col]+count(c)
	end
  	cptr
end

function SMCSC(A)
	nzA=(!ismissing).(A)
	r,c=size(A)
	#cptr=[0;cumsum(count.(eachcol(nzA)))].+1
	cptr=cptr_(nzA)
	nzval= A[nzA]
	rvals=[first(ci.I) for ci in findall(nzA)]
	SparseMatrixCSC(r,c,cptr,rvals,nzval)
end

```

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 1, 2024, 8:13pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/9 "2024-01-01T20:13:15Z")

</div>

> [@rocco\_sprmnt21](#):
>
> `cptr=Vector{Int64}(undef,10001)`

this looks a bit suspicious?  
Also, the functions need to return a non-Union type.

Another method to consider:

```julia
using MappedArrays

sparsenonmissing4(A::AbstractMatrix) = 
  sparse(mappedarray(x->coalesce(x, zero(nonmissingtype(eltype(A)))), A))

```

As the saying goes: why not make it a one-liner? 😉

This is faster than `sparsenonmissing3`.

Tiny issue: if the original matrix contains both zeroes and missings, then the resulting matrix will make them both structural zeroes in the sparse matrix.

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [January 1, 2024, 9:11pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/10 "2024-01-01T21:11:20Z")

</div>

I’m not sure I understand your doubt. If you are referring to the fact that the size of the cptr vector is hardcoded, I assume this can generally be calculated as `size(A)[2]+1`.  
I don’t get the other aspect about the type of elements of the matrix. I still don’t have a clear understanding of how the type system works in Julia (sooner or later I’ll have to start studying it)

> [@Dan](#):
>
> this looks a bit suspicious?  
> Also, the functions need to return a non-Union type.

Here’s what I measure on my laptop

> **Summary**
>
> ```julia
> julia> using SparseArrays, Missings, Random
> 
> julia> A = missings(Float64, 500,10000)
> 500×10000 Matrix{Union{Missing, Float64}}:
>  
> missing missing missing missing missing missing missing
> ⋮ ⋱
> missing missing missing missing … missing missing missing
>  
> missing missing missing missing missing missing missing
> 
> julia> A[1:5:end] .= randn(Float64,length(1:5:length(A)))
> 1000000-element view(reshape(::Matrix{Union{Missing, Float64}}, 5000000), 1:5:4999996) with eltype Union{Missing, Float64}:
> 0.35955623771445555
> 0.7063141021260603
> 0.8332490679598374
> 1.865779624046154
> ...
> 1.342504874329712
> -1.1325925605137113
> 
> julia> shuffle!(A)
> 500×10000 Matrix{Union{Missing, Float64}}:
> -1.28583 missing 0.471741 … missing missing
> ...
> missing missing missing missing -1.2943
> -2.58484 missing missing -1.07261 missing
> -0.0542834 missing missing missing -1.05134
> 
> julia> using MappedArrays
> 
> julia> sparsenonmissing4(A::AbstractMatrix) =
> sparse(mappedarray(x->coalesce(x, zero(nonmissingtype(eltype(A)))), A))
> sparsenonmissing4 (generic function with 1 method)
> 
> julia> function cptr_(nzA)
> nr,nc=size(nzA)
> cptr=Vector{Int64}(undef,nc+1)
> cptr[1]=1
> for (col,c) in enumerate(eachcol(nzA))
> cptr[col+1]=cptr[col]+count(c)
> end
> cptr
> end
> cptr_ (generic function with 1 method)
> 
> julia> function SMCSC(A)
> nzA=(!ismissing).(A)
> r,c=size(A)
> #cptr=[0;cumsum(count.(eachcol(nzA)))].+1
> cptr=cptr_(nzA)
> nzval= A[nzA]
> rvals=[first(ci.I) for ci in findall(nzA)]
> SparseMatrixCSC(r,c,cptr,rvals,nzval)
> end
> SMCSC (generic function with 1 method)
> 
> julia> SMCSC(A) == sparsenonmissing4(A) #true
> true
> 
> julia> using BenchmarkTools
> 
> julia> @btime sparsenonmissing4(A)
> 29.725 ms (31 allocations: 19.64 MiB)
> 500×10000 SparseMatrixCSC{Float64, Int64} with 1000000 stored entries:
> ⎡⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎤
> ⎣⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎦
> 
> julia> @btime SMCSC(A)
> 14.931 ms (13 allocations: 32.15 MiB)
> 500×10000 SparseMatrixCSC{Union{Missing, Float64}, Int64} with 1000000 stored entries:
> ⎡⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎤
> ⎣⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎦
> 
> ```

Okay. Maybe now I understand what you were referring to (although having a Union{Missing+Float} type didn’t give a different result than yours)  
Here is a modified version that allocates less and runs a little faster

```julia
julia> function SMCSC(A)
        nzA=findall(!ismissing,A)
        r,c=size(A)
        cptr=[1;findall(i->last(nzA[i].I)!=last(nzA[i+1].I), 1:length(nzA)-1).+1;length(nzA)+1]   
        nzval = [A[a] for a in nzA]
        rvals=[first(ci.I) for ci in nzA]
        SparseMatrixCSC(r,c,cptr,rvals,nzval)
       end
SMCSC (generic function with 1 method)

julia> @btime SMCSC(A)
  11.594 ms (29 allocations: 31.47 MiB)
500×10000 SparseMatrixCSC{Float64, Int64} with 1000000 stored entries:
⎡⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎤
⎣⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎦

```

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [January 2, 2024, 1:45pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/11 "2024-01-02T13:45:04Z")

</div>

regarding the belief that a good (or lucky) use of the functions defined in the library can be more effective than a theoretically better for … loop, here it happens that S!MCSC2(A) which would seem in theory to perform the necessary operations in an “optimal” way is less performing than S!MCSC1(A).  
I can’t explain why.

> **Summary**
>
> ```julia
> julia> function S!MCSC1(A)
> nzA=findall(!ismissing,A)
> r,c=size(A)
> cptr=[1;findall(i->last(nzA[i].I)!=last(nzA[i+1].I), 1:length(nzA)-1).+1;length(nzA)+1]
> nzval = [A[a] for a in nzA]
> #rvals=[first(ci.I) for ci in nzA]
> rvals=[first(Tuple(ci)) for ci in nzA]
> SparseMatrixCSC(r,c,cptr,rvals,nzval)
> end
> S!MCSC1 (generic function with 1 method)
> 
> julia> @btime S!MCSC1(A)
> 12.054 ms (29 allocations: 31.47 MiB)
> 500×10000 SparseMatrixCSC{Float64, Int64} with 1000000 stored entries:
> ⎡⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎤
> ⎣⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎦
> 
> julia> function S!MCSC2(A)
> r,c = size(A)
> rvals = Int64[]
> cptr = [1]
> nzvals = Float64[]
> sizehint!(rvals, trunc(Int,r*c*0.2))
> sizehint!(cptr, c+1)
> sizehint!(nzvals, trunc(Int,r*c*0.2))
> for (ci,col) in enumerate(eachcol(A))
> nzA=findall(!ismissing,col)
> push!(cptr,length(nzA)+cptr[end])
> append!(rvals,nzA)
> append!(nzvals,@view A[nzA,ci])
> # nzvals = [A[a] for a in nzA]
> # #rvals=[first(ci.I) for ci in nzA]
> # rvals=[first(Tuple(ci)) for ci in nzA]
> end
> SparseMatrixCSC(r,c,cptr,rvals,nzvals)
> end
> S!MCSC2 (generic function with 1 method)
> 
> julia> sparsenonmissing4
> sparsenonmissing4 (generic function with 1 method)
> 
> julia> S!MCSC1(A)==S!MCSC2(A)== sparsenonmissing4(A)
> true
> 
> julia> @btime S!MCSC2(A)
> 17.470 ms (40007 allocations: 66.32 MiB)
> 500×10000 SparseMatrixCSC{Float64, Int64} with 1000000 stored entries:
> ⎡⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎤
> ⎣⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎦
> 
> ```

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [January 2, 2024, 2:10pm UTC](https://discourse.julialang.org/t/fast-conversion-from-matrix-union-missing-float64-to-sparse/108161/12 "2024-01-02T14:10:55Z")

</div>

Thanks for the replies everyone!

> [@stevengj](#):
>
> What operations do you want to support on the array? (Most of the SparseMatrixCSC functionality is oriented towards linear algebra, which is useless if there are missing values.)

It is linear algebra. Basically, this is an unbalanced panel of data (rows are dates, columns are firm observations). The unbalanced nature led to a lot of joins on dates in my first go at the estimation when I kept things in the DataFrame. So, I am trying out the approach where I store the data in this sparse matrix instead, and I also manually handle the joins via clever indexing. Maybe it will not be as good, I’m not sure.