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
On Tuesday, April 12, 2016 at 4:41:06 PM UTC-5, Laurent Bartholdi wrote:
>
> @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 <ste...@karpinski.org
> <javascript:>> 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 <laurent....@gmail.com
>> <javascript:>> 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 <alex...@gmail.com
>>> <javascript:>> 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,
>>>> ming...@irt-systemx.fr 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
>
