Best way to remove column from CSC Matrix

Viewed 76

Lest say I have the following Sparse Matrix:

A = sparse([1 4 7; 2 5 8; 3 6 9])

What is the most efficient way to get:

A = sparse([1 0 7; 2 0 8; 3 0 9])

The actual matrix I'm using is actually a lot bigger...

2 Answers

To set column two to zeros you can use:

function set2zero!(A, col)
    @assert 1 <= col <= size(A, 2)
    A.nzval[A.colptr[col]:A.colptr[col+1]-1] .= zero(eltype(A))
    return A
end

which will be fast (especially if A is sparse).

However, note that then 0 entries are stored in the matrix (so there is a tradeoff of speed of update vs later speed of working with such array).

Bogumil's answer most likely will do the trick. But, just for reference, another method to remove a column is using the following function:

function deletecol!(A, col)
    @assert 1 <= col <= size(A, 2)
    V = A.colptr
    deleteat!(A.nzval, V[col]:V[col+1]-1)
    deleteat!(A.rowval, V[col]:V[col+1]-1)
    @views V[col+1:end] .-= V[col+1]-V[col]
    return A
end

This function not only zeroes the column but also updates the sparsity pattern to reflect an empty column. In some usage cases this may have performance benefits down the line. The cost is immediate, as the operation takes longer. In the toy case presented in OP, this is not a huge difference:

julia> bench_1() = ( A = sparse([1 4 7; 2 5 8; 3 6 9]) ; set2zero(A,2) )
bench_1 (generic function with 1 method)

julia> bench_2() = ( A = sparse([1 4 7; 2 5 8; 3 6 9]) ; deletecol!(A,2) )
bench_2 (generic function with 1 method)

julia> using BenchmarkTools

julia> @btime bench_1();
  404.020 ns (5 allocations: 560 bytes)

julia> @btime bench_2();
  421.774 ns (5 allocations: 560 bytes)

The output A is compressed as can be seen by:

julia> A
3×3 SparseMatrixCSC{Int64, Int64} with 6 stored entries:
 1  ⋅  7
 2  ⋅  8
 3  ⋅  9

Ideally if there are many column deletion operations, after using Bogumil's function, a sparse pattern compression can be called (without going through non-sparse matrix). But I couldn't find the function name just now. Perhaps a commenter can (A = sparse(collect(A)) goes through full matrix first).

UPDATE: The function to compress a sparse matrix representation is dropzeros!. So we have:

julia> A
3×3 SparseMatrixCSC{Int64, Int64} with 9 stored entries:
 1  0  7
 2  0  8
 3  0  9

julia> dropzeros!(A)
3×3 SparseMatrixCSC{Int64, Int64} with 6 stored entries:
 1  ⋅  7
 2  ⋅  8
 3  ⋅  9

(Thanks Bogumil!)

The implementation of dropzeros! also points to a more compact version of deletecol! defined as follows:

deletecol!(A,cind) = SparseArrays.fkeep!(A,(i,j,v) -> j != cind)

giving:

julia> A = sparse([1 4 7; 2 5 8; 3 6 9])
3×3 SparseMatrixCSC{Int64, Int64} with 9 stored entries:
 1  4  7
 2  5  8
 3  6  9

julia> deletecol!(A,3)
3×3 SparseMatrixCSC{Int64, Int64} with 6 stored entries:
 1  4  ⋅
 2  5  ⋅
 3  6  ⋅

This definition is more Julian because of higher composability and it does fit on one-line :-). fkeep! can also be used to define other useful manipulations (deleting multiple columns/rows etc.)

Related