sparse rank, sparse nullspace, sparse linear algebra over rationals?

[Laurent Bartholdi]

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

[Seth]

well, nullspace appears to be defined for dense arrays only - and rank is defined for abstract arrays but will fail in svdvals!() when passed a sparse array (this is arguably a bug).

[ming...@irt-systemx.fr]

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.

[Alex Dowling]

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. Does anyone have recommendations for similar functions in Julia?

Thanks,

Alex

[Laurent Bartholdi]

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

--

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

[Alex Dowling]

Hello Laurent,

Thank you for the prompt reply. I noticed your pure Julia code converted to a full matrix. Any luck with the sparse SVD routines?

I thought more about my application, and I really want to know the smallest n < 10 singular values and associated vectors for the system. Tight tolerances are not required. I see that eigs() supports specifying the largest or smallest eigenvalues and the number, but SVD does not have a similar option. The sparse matricies I'm considering an not symmetric.

Alex

[Stefan Karpinski]

Julia does have complex sparse matrices:

julia> C = sparse(rand(1:10, 10), rand(1:10, 10), randn(10) + randn(10)im, 10, 10)

10x10 sparse matrix with 9 Complex{Float64} entries:

  [7 ,  1]  =  1.1711-0.587334im

  [5 ,  3]  =  0.940755+1.00755im

  [1 ,  4]  =  1.51515-1.77335im

  [4 ,  4]  =  -0.815624-0.678038im

  [9 ,  5]  =  -0.598111-0.970266im

  [2 ,  6]  =  0.671339-1.07398im

  [7 ,  6]  =  -1.69705+1.08698im

  [10,  7]  =  -1.32233-1.88083im

  [7 , 10]  =  1.26453-0.399423im

How much you can do with these depends on what kinds of special methods are implemented, but you can call eigs on them, which lets you figure out what the rank is:

julia> map(norm,eigs(C)[1])

6-element Array{Float64,1}:

 1.74612

 1.74612

 1.06065

 4.28017e-6

 4.28016e-6

 4.28016e-6

In this case, for example, the matrix is rank 3.

[Laurent Bartholdi]

@Stefan: sorry, I had meant that Julia doesn't yet have sparse non-(real/complex) matrices. However, I see that that's wrong too: I can create sparse matrices with Int64 entries. It would now make a lot of sense to implement the algorithms in linbox so as to compute rank etc. of Int or BigInt sparse matrices.

@Alex: if you're only interested in the 10 lowest singular values of a sparse matrix, then MATLAB is quite good at it. Here's from the doc:

    S = svds(A,K,SIGMA) computes the K singular values closest to the 

    scalar shift SIGMA.  For example, S = svds(A,K,0) computes the K

    smallest singular values.

[Alex Dowling]

Hello Laurent,

Thanks. I've used svds in MATLAB before with some success. Ideally I'd like to have a pure Julia implementation. Do you know of any appropriate packages?

Alex

[James Fairbanks]

Does this help? https://github.com/JuliaLang/IterativeSolvers.jl/blob/master/src/svdl.jl#L192