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
On Friday, August 26, 2016 at 4:46:08 PM UTC-4, Alex Dowling wrote:
>
> 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
>
