@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.
On Tue, 12 Apr 2016 at 22:29 Stefan Karpinski <[email protected]> wrote:
> 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.
>
> On Tue, Apr 12, 2016 at 3:29 PM, Laurent Bartholdi <
> [email protected]> wrote:
>
>> 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
>>
>
> --
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
