Does this help? https://github.com/JuliaLang/IterativeSolvers.jl/blob/master/src/svdl.jl#L192
On Wednesday, April 13, 2016 at 10:30:32 AM UTC-4, Alex Dowling wrote: > > 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 <[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 >> >
