SparseArrays.jl

SparseArrays.jl copied to clipboard

Calling reverse on a sparse array introduces structural nonzeros. Here's a minimal example

julia> a = SparseMatrixCSC(2, 2, [1,3,4], [1,2,1], [1,1,1])
2×2 SparseMatrixCSC{Int64, Int64} with 3 stored entries:
 1  1
 1  ⋅

julia> reverse(a,dims=:)
2×2 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 0  1
 1  1

julia> reverse(a,dims=1)
2×2 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  0
 1  1

julia> reverse(a,dims=2)
2×2 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  1
 0  1

i guess we could overload _reverse! but is there a nice inplace way of doing it for dim=2?

@Wimmerer for reverse(dim=2) can we do better than transpose!, reverse!(dim=1), transpose! ?

I wasn't even aware of a reverse function in Julia for non Vectors. I'll have to think about it, I doubt there's anything better that's not more involved than it's worth

Here's an alternative that is simple and more performant:

using SparseArrays


function reverse_sparse( A:: SparseMatrixCSC{Tv, Ti} ) where { Tv, Ti }

    m, n        =   size(A)
    K           =   length(A.colptr) # = n+1    
    rowval      =   m + 1 .- reverse( A.rowval )
    nzval       =   reverse( A.nzval )

    colptr      =   Array{Ti}(undef, K )
    colptr[1]   =   1

    for p in 2:K
        colptr[p]   =   colptr[p-1] + ( A.colptr[K-p+2] - A.colptr[K-p+1] )
    end

    return SparseMatrixCSC( m, n, colptr, rowval, nzval )
end


function benchmark()

    println()
    println("Large matrix")
    println("--------------")
    println()
    A           =   sprand(1000,1000,0.1)

    println("Current implementation:")
    @time Arev1  =   dropzeros(reverse(A))
    println()

    println("Custom implementation:")    
    @time Arev2  =   reverse_sparse(A)
    println()    

    println("Result is correct?")
    println(Arev1 == Arev2)
    println()    

    println("Small matrix")
    println("--------------")
    println()
    A           =   sparse([0 0 1; 1 0 0])

    println("Current implementation:")
    @time Arev1  =   dropzeros(reverse(A))
    println()

    println("Custom implementation:")    
    @time Arev2  =   reverse_sparse(A)
    println()    

    println("Result is correct?")
    println(Arev1 == Arev2)

end

And the result

Large matrix
--------------

Current implementation:
  1.513817 seconds (14 allocations: 7.960 MiB)

Custom implementation:
  0.000461 seconds (7 allocations: 2.291 MiB)

Result is correct?
true

Small matrix
--------------

Current implementation:
  0.000003 seconds (9 allocations: 656 bytes)

Custom implementation:
  0.000002 seconds (7 allocations: 480 bytes)

Result is correct?
true

The reverse_sparse function can easily be modified to reverse just rows or just columns.

function reverse_sparse2(A::AbstractSparseMatrixCSC)
    rowval = reverse(A.rowval)
    rowval .= size(A, 1) + 1 .- rowval

    nzval = reverse(A.nzval)

    colptr = similar(getcolptr(A))
    colptr[firstindex(colptr)] = 1

    @inline for p in eachindex(colptr)[2:end]
        colptr[p] = colptr[p-1] + (getcolptr(A)[end-p+2] - A.colptr[end-p+1])
    end

    return @if_move_fixed A SparseMatrixCSC(size(A)..., colptr, rowval, nzval)
end

is a bit faster, allocates less and more up to code but you should open a pull request.

Do note that we also need reverse! which would be better to implement. i suggest these but honestly i think given how much time that they take (10x time your code), i'm doing something wrong.

function reverse!_col(x)
    @inline for i in 1:size(x, 2)
        s = getcolptr(x)[i]
        t = getcolptr(x)[i + 1] - 1
        reverse!(view(nonzeros(x), s:t))
        v = view(rowvals(x), s:t)
        v .= (size(x, 1) + 1) .- v
        reverse!(v)
    end
    return x
end
function reverse!_row(x)
    y = transpose!(similar(x), x)
    reverse!_col(y)
    transpose!(x, y)
    return x
end

The reverse col is slow? That's surprising, that's what I was going to propose and I don't see a better way immediately. Rows is harder.

In place*

reverse_sparse2 takes like 10 micro second, _col takes 100 micro seconds

The reverse col is slow? That's surprising, that's what I was going to propose and I don't see a better way immediately. Rows is harder.

Can we do temporaries in the row! one? That might be faster than the double transpose. Allocate similar(x) then just copyto! the vectors in reverse order.

E: not sure why this reposted

a = sprandn(10_000, 10_000, 0.00001)
julia> @btime reverse!_col(copy($a));
  123.924 μs (4 allocations: 94.36 KiB)

julia> @btime reverse!_row(copy($a));
  154.990 μs (8 allocations: 188.72 KiB)
julia> @btime reverse_sparse($a);
  11.662 μs (5 allocations: 102.42 KiB)

julia> @btime reverse_sparse2($a);
  11.191 μs (4 allocations: 94.36 KiB)

julia> @btime reverse_sparse!(copy($a));
  269.779 μs (8 allocations: 188.72 KiB)

Can we do temporaries in the row! one? That might be faster than the double transpose. Allocate similar(x) then just copyto! the vectors in reverse order.

not sure i understand.

I'm a little unfamiliar with macros, and I'm getting an error for code that includes @if_move_fixed. Is this special to SparseArrays?

Yes, it's specific to this package. It converts the last argument of the macro to a matrix with fixed pattern if any of the other arguments have fixed pattern

On Wed, Dec 7, 2022, 20:29 Gregory Henselman-Petrusek < @.***> wrote:

I'm a little unfamiliar with macros, and I'm getting an error for code that includes @if_move_fixed. Is this special to SparseArrays?

— Reply to this email directly, view it on GitHub https://github.com/JuliaSparse/SparseArrays.jl/issues/278#issuecomment-1341845352, or unsubscribe https://github.com/notifications/unsubscribe-auth/AARRZZZBICY3SFMDDLQ6LG3WME2YXANCNFSM6AAAAAARKVADYE . You are receiving this because you commented.Message ID: @.***>

Also, are you using julia 1.9? because it's julia 1.9 thing.

Oh that must be it -- I was running 1.8. I'm new to many aspects of development -- is there a way to PR two different versions of the code, one that'd work with the stable 1.6 release and one that works with 1.9?

you have to make pull request to merge into the backport-1.6 branch and one for the main branch.