It's tricky, but sprank is definitely not competitive with finer
implementations. As you certainly know, rank calculations are very
unstable, so if your matrix has integer entries you should prefer an exact
method, or a calculation in a finite field. (Sadly, Julia does not contain,
yet, sparse matrices with non-real/complex entries). See
http://dx.doi.org.sci-hub.io/10.1145/2513109.2513116 for a recent
discussion of sparse matrix rank.
The following routines compute null spaces. To obtain the rank, you
subtract the nullspace dimension from the row space dimension (i.e. number
of columns). If you want pure Julia, you could try
function julia_null_space{T}(A::SparseMatrixCSC{T})
SVD = svdfact(full(A), thin = false)
any(1.e-8 .< SVD.S .< 1.e-4) && error("Values are dangerously small")
indstart = sum(SVD.S .> 1.e-8) + 1
nullspace = SVD.Vt[indstart:end,:]'
return sparse(nullspace .* (abs(nullspace) .> 1.e-8))
end
If you have MATLAB, the following is faster:
using MATLAB
global m_session = MSession()
function matlab_null_space{T}(A::SparseMatrixCSC{T})
m, n = size(A)
(m == 0 || n == 0) && return speye(T, n)
put_variable(m_session,:A,convert(SparseMatrixCSC{Float64},A))
eval_string(m_session,"N = spqr_null(A,struct('tol',1.e-8));")
eval_string(m_session,"ns = spqr_null_mult(N,speye(size(N.X,2)),1);")
eval_string(m_session,"ns = sparseclean(ns,1.e-8);")
ns = get_mvariable(m_session,:ns)
return jvariable(ns)
end
Finally, if your matrices are over Z, Q or a finite field, do give a look
to linbox (http://www.linalg.org/linbox-html/index.html).
HTH, Laurent
On Tue, 12 Apr 2016 at 20:49 Alex Dowling <[email protected]> wrote:
> Hello,
>
> I am also looking to compute the (approximate) rank of a sparse matrix.
> For my applications I'm consider dimensions of 10k - 100k by 10k - 100k,
> not millions by millions. I've previously done this in MATLAB using the sprank
> command <http://www.mathworks.com/help/matlab/ref/sprank.html>. Does
> anyone have recommendations for similar functions in Julia?
>
> Thanks,
> Alex
>
>
> On Tuesday, November 17, 2015 at 3:08:29 AM UTC-6, [email protected]
> wrote:
>>
>> hello,
>>
>> getting the rank of a huge sparse matrix is mathematically difficult
>>
>>
>> http://math.stackexchange.com/questions/554777/rank-computation-of-large-matrices
>>
>> bests,
>> M.
>>
>> Le lundi 16 novembre 2015 22:12:04 UTC+1, Laurent Bartholdi a écrit :
>>>
>>> Hello world,
>>> I'm new at julia, and trying it out as a replacement for matlab and
>>> other computer algebra systems. I'm a bit stuck with sparse matrices: I
>>> have largeish matrices (10^6 x 10^6) that are very sparse, and
>>> want to know their rank and nullspace. (the matrices decompose by
>>> blocks, so the nullspace should be expressible as a sparse matrix).
>>>
>>> I tried with the latest (github) julia 0.5:
>>>
>>> julia> nullspace(sparse([1],[1],[1]))
>>> ERROR: MethodError: `nullspace` has no method matching
>>> nullspace(::SparseMatrixCSC{Int64,Int64})
>>>
>>> julia> nullspace(full(sparse([1],[1],[1])))
>>> 1x0 Array{Float64,2} # I'm a bit unhappy here, I was hoping to get a
>>> rational answer.
>>>
>>> julia> nullspace([1//1])
>>> 1x0 Array{Float32,2} # yikes! I'm down to 32 bits floats now.
>>>
>>> julia> rank(sparse([1],[1],[1.0]))
>>> ERROR: MethodError: `svdvals!` has no method matching
>>> svdvals!(::SparseMatrixCSC{Float64,Int64})
>>>
>>> julia> rank([1.0])
>>> ERROR: MethodError: `rank` has no method matching
>>> rank(::Array{Float64,1})
>>>
>>> julia> rank([1 0;0 1])
>>> 2 # finally something that works...
>>>
>>> Many thanks in advance! Laurent
>>>
>>> --
Laurent Bartholdi
DMA, École Normale Supérieure, 45 rue d'Ulm, 75005 Paris. +33 14432 2060
Mathematisches Institut, Universität Göttingen, Bunsenstrasse 3-5, D-37073
Göttingen. +49 551 39 7826
