Is there a package which provides a singular value decomposition of a (square, numerically unsymmetric, structurally symmetric, indefinite) sparse matrix and allows for finding the n smallest singular values (instead of the n largest singular values as svds)? In particular, I am looking to compute a (full, but will settle for truncated) basis for the nullspace.
1 Like
Maybe IterativeSolvers.jl?
1 Like
the built-in svds with option which=:SM?
At least on v0.5.1, svds does not support which. IterativeSolvers.svdl documentation says it computes “some” singular values, whatever this may mean.
svds is fairly straightforward code (in base/linalg/arnoldi.jl) which invokes eigs. You could copy the logic and add the which=:SM option to the eigs call.
Yes, I looked at that but was unsure as to whether this alters something which shouldn’t be altered. I assumed that the original author had something in mind when restricting the function to finding the largest singular values.
Thanks for the suggestion of PROPACK, however it looks as though my matrix is too large: I get an InexactError in wrappers.jl line 248:
lwork = Int32(m + n + 10kmax + 5kmax*kmax + 4 + max(3kmax*kmax + 4kmax + 4, nb*max(m, n)))
…while:
julia> m
16992
julia> n
16992
julia> kmax
17002
julia> m + n + 10kmax + 5kmax*kmax + 4 + max(3kmax*kmax + 4kmax + 4, nb*max(m, n))
2312816052
julia> typemax(Int32)
2147483647
Any way to go bigger without risking non-convergence?
That’s a limitation of the Fortran code using 32 bit ints. The short-term solution is to set kmax = 16383, or to try a small value of kmax and hope for the best. The long-term solution would be to rewrite PROPACK in Julia. Andreas Noack started TSVD.jl but I’m not sure if it allows you to compute the smallest singular values.
The mid-term solution is to try the int64 branch created just now. I admit I haven’t really checked that all will work correctly but since the Fortran code only uses standard integer variables, I made a few changes and the tests pass when compiling with fdefault-integer-8 (I’m assuming GFortran).
So if you’d like to give it a go, Pkg.checkout("PROPACK", "int64"); Pkg.build("PROPACK"); Pkg.test("PROPACK"). Better report issues on the tracker than here though.
And the long-term solution is to rewrite the library in Julia.
Be my guest.
2 Likes
I will try the int64 branch in a few days, thanks for the quick response!
I read some stuff on PRIMME_SVDS, which appears to be more reliable and efficient than PROPACK. In particular it appears as though it handles matrices with N in the order of millions quite well. There already are some wrappers for MATLAB and Python, perhaps it’s worth the effort to wrap this in Julia if there is no similarly-performant alternative available as of now.
https://github.com/andreasnoack/Primme.jl
Out of curiosity, did you consider sparse QR? If your nullspace is large, that might be competitive.
Wow, I must have missed that one when googling! I’m not on v0.6 yet, as I’m waiting for things to stabilize in the package ecosystem, so I’ll have to wait on Primme.jl.
No I haven’t tried sparse QR as it appears to fail on v0.5.1 but perhaps I’m doing something wrong:
julia> qr(speye(20))
ERROR: MethodError: no method matching qrfact!(::SparseMatrixCSC{Float64,Int64}, ::Type{Val{false}})
Closest candidates are:
qrfact!{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(::Union{Base.ReshapedArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N}},L}}, ::Type{Val{false}}) at linalg/qr.jl:88
qrfact!(::Union{Base.ReshapedArray{T,2,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray{T,2},SubArray{T,2,A<:Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N}},L}}, ::Type{Val{false}}) at linalg/qr.jl:93
in qrfact(::SparseMatrixCSC{Float64,Int64}, ::Type{T}) at ./linalg/qr.jl:164
in #_qr#16 at ./linalg/qr.jl:176 [inlined]
in (::Base.LinAlg.#kw##_qr)(::Array{Any,1}, ::Base.LinAlg.#_qr, ::SparseMatrixCSC{Float64,Int64}, ::Type{Val{false}}) at ./<missing>:0
in qr(::SparseMatrixCSC{Float64,Int64}) at ./linalg/qr.jl:173
That seems like a bug. You can try qrfact(speye(20))
Then how do I get Q and R from
julia> qrfact(speye(20))
Base.SparseArrays.SPQR.Factorization{Float64}(20,20,Ptr{Base.SparseArrays.SPQR.C_Factorization{Float64}} @0x00007ff237146230)
…couldn’t find any documentation on how to access the output
There are some breadcrumbs here Documentation for SPQR qrfact() · Issue #12298 · JuliaLang/julia · GitHub
Perhaps we should file an issue for the qr problem, don’t really feel like building a proprietary thing to get Q and R from qrfact.
Yes, an issue would be good here.
Issue #21937
1 Like