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.

# Sparse SVD (smallest SVs)

**rleegates**#1

**rleegates**#4

At least on v0.5.1, svds does not support `which`

. `IterativeSolvers.svdl`

documentation says it computes â€śsomeâ€ť singular values, whatever this may mean.

**Ralph_Smith**#5

`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.

**rleegates**#6

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.

**rleegates**#8

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?

**dpo**#9

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.

**rleegates**#12

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.

**Ralph_Smith**#13

Out of curiosity, did you consider sparse QR? If your nullspace is large, that might be competitive.

**rleegates**#14

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
```

**rleegates**#16

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 https://github.com/JuliaLang/julia/issues/12298

**rleegates**#18

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`

.