Using rank() and nullspace() with sparse matrix

[Alex Dowling]

When I use rank() on a sparse matrix, I get the following error:

ERROR: LoadError: MethodError: `svdvals!` has no method matching svdvals!(::SparseMatrixCSC{Float64,Int64})

 in rank at linalg/generic.jl:300

Small example:

A = sprand(5,5,0.6)

rank(A)

Is there an alternate way to calculate the rank of a sparse matrix? Same question for calculating the nullspace basis.

Thanks,

Alex

[vav...@uwaterloo.ca]

You can compute the rank of a sparse matrix using the 'squeezed QR factorization' available in SuiteSparse, which has been partly wrapped in Julia.  This method for rank determination is less robust than either the SVD or conventional rank-revealing QR factorizations, but it appears to be better at preserving sparsity, at least according to the authors of SuiteSparse.

Here is a routine I wrote or in some cases lifted from Base (Julia 0.4.6) to extract the R and colperm factors from the sparse squeezed QR factorization of a matrix A using SuiteSparse.  The number of rows in R upon termination is the computed rank of the matrix.  (Feel free to use this code if you wish.  I'm posting it here for whoever wants it under the same license as Julia, that is, the MIT license.  At some point I will write some test routines for it and submit it as PR for the SuiteSparse.jl package.)  I have run this code only on Windows 10 64-bit.

-- Steve Vavasis


module MySparseExtensions


import Base.sparse
import Base.SparseMatrix.CHOLMOD.FactorComponent
import Base.SparseMatrix.CHOLMOD.Factor
import Base.SparseMatrix.CHOLMOD.CHOLMODException
import Base.SparseMatrix.CHOLMOD.common
import Base.SparseMatrix.CHOLMOD.C_Sparse
import Base.SparseMatrix.CHOLMOD.Sparse
import Base.SparseMatrix.CHOLMOD.free_sparse!
import Base.SparseMatrix.CHOLMOD.increment
import Base.SparseMatrix.CHOLMOD.SuiteSparse_long
import Base.SparseMatrix.CHOLMOD.defaults
import Base.SparseMatrix.CHOLMOD.fact_
import Base.SparseMatrix.CHOLMOD.set_print_level
import Base.SparseMatrix.CHOLMOD.common_final_ll


function qrfact_get_r_colperm(A::SparseMatrixCSC{Float64, Int},
                              tol::Float64,
                              ordering = Base.SparseMatrix.SPQR.ORDERING_DEFAULT)
    Aw = Sparse(A,0)
    s = unsafe_load(Aw.p)
    if s.stype != 0
        throw(ArgumentError("stype must be zero"))
    end
    ptr_R = Ref{Ptr{C_Sparse{Float64}}}()
    ptr_E = Ref{Ptr{SuiteSparse_long}}()
    cc = common()
    rk = ccall((:SuiteSparseQR_C, :libspqr), Clong,
               (Cint, #ordering
                Cdouble, #tol
                Clong, #econ
                Cint, #getCTX
                Ptr{C_Sparse{Float64}},  # A
                Ptr{Void}, #Bsparse
                Ptr{Void}, #Bdense
                Ptr{Void}, #Zsparse
                Ptr{Void}, #Zdense
                Ptr{Void}, #R
                Ptr{Void}, #E
                Ptr{Void}, #H
                Ptr{Void}, #HPInv
                Ptr{Void}, #HTau
                Ptr{Void}), #cc
               ordering, #ordering
               tol, #tol
               1000000000, #econ
               0, #getCTX
               get(Aw.p),  # A
               C_NULL, #Bsparse
               C_NULL, #Bdense
               C_NULL, #Zsparse
               C_NULL, #Zdense
               ptr_R, #R
               ptr_E, #E
               C_NULL, #H
               C_NULL, #HPInv
               C_NULL, #HTau
               cc) #cc
    R = ptr_R[]
    E = ptr_E[]
    if E != C_NULL
        colprm = pointer_to_array(E, size(A,2), false) + 1
    else
        colprm = collect(1:size(A,2))
    end
    R1 = unsafe_load(R)
    if R1.stype != 0
        throw(ArgumentError("matrix has stype != 0. Convert to matrix with stype == 0 before converting to SparseMatrixCSC"))
    end
    maxrownum = 0
    for entryind = 1 : R1.nzmax
        maxrownum = max(maxrownum, unsafe_load(R1.i, entryind) + 1)
    end
    R_cvt = SparseMatrixCSC(maxrownum,
                            R1.ncol,
                            increment(pointer_to_array(R1.p, (R1.ncol + 1,), false)),
                            increment(pointer_to_array(R1.i, (R1.nzmax,), false)),
                            copy(pointer_to_array(R1.x, (R1.nzmax,), false)))
    free_sparse!(R)
    ccall((:cholmod_l_free, :libcholmod), Ptr{Void}, (Csize_t, Csize_t, Ptr{Void}, Ptr{Void}),
          size(A,2), sizeof(SuiteSparse_long), E, cc)
    R_cvt, colprm
end


end