# Sparse SVD (smallest SVs)

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