# Sparse SVD (smallest SVs)
**URL:** https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645
**Category:** General Usage
**Created:** [May 11, 2017, 8:57am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645 "2017-05-11T08:57:28Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 11, 2017, 8:57am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/1 "2017-05-11T08:57:28Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 11, 2017, 9:42am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/2 "2017-05-11T09:42:15Z")
Maybe IterativeSolvers.jl?
---
### Author:  [@innerlee](https://discourse.julialang.org/u/innerlee)
#### Post date: [May 11, 2017, 10:00am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/3 "2017-05-11T10:00:58Z")
the built-in `svds` with option `which=:SM`?
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 11, 2017, 10:03am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/4 "2017-05-11T10:03:47Z")
At least on v0.5.1, svds does not support `which`. `IterativeSolvers.svdl` documentation says it computes “some” singular values, whatever this may mean.
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [May 11, 2017, 2:02pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/5 "2017-05-11T14:02:07Z")
`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.
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 11, 2017, 3:38pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/6 "2017-05-11T15:38:59Z")
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.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [May 11, 2017, 5:26pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/7 "2017-05-11T17:26:15Z")
[https://github.com/JuliaSmoothOptimizers/PROPACK.jl](https://github.com/JuliaSmoothOptimizers/PROPACK.jl)
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 11, 2017, 8:47pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/8 "2017-05-11T20:47:16Z")
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:
```julia
lwork = Int32(m + n + 10kmax + 5kmax*kmax + 4 + max(3kmax*kmax + 4kmax + 4, nb*max(m, n)))
```
…while:
```julia
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?
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [May 11, 2017, 10:52pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/9 "2017-05-11T22:52:14Z")
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.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [May 12, 2017, 12:14am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/10 "2017-05-12T00:14:07Z")
And the long-term solution is to rewrite the library in Julia.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [May 12, 2017, 12:29am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/11 "2017-05-12T00:29:01Z")
Be my guest.
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 12, 2017, 1:18pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/12 "2017-05-12T13:18:56Z")
I will try the int64 branch in a few days, thanks for the quick response!
I read some stuff on [PRIMME\_SVDS](https://arxiv.org/abs/1607.01404), 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.
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [May 13, 2017, 7:27pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/13 "2017-05-13T19:27:52Z")
[https://github.com/andreasnoack/Primme.jl](https://github.com/andreasnoack/Primme.jl)
Out of curiosity, did you consider sparse QR? If your nullspace is large, that might be competitive.
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 17, 2017, 1:23pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/14 "2017-05-17T13:23:52Z")
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
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 ./:0
in qr(::SparseMatrixCSC{Float64,Int64}) at ./linalg/qr.jl:173
```
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [May 17, 2017, 1:28pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/15 "2017-05-17T13:28:29Z")
That seems like a bug. You can try `qrfact(speye(20))`
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 18, 2017, 7:22am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/16 "2017-05-18T07:22:57Z")
Then how do I get Q and R from
```julia
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
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [May 18, 2017, 8:19am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/17 "2017-05-18T08:19:52Z")
> [@rleegates](#):
>
> qrfact(speye(20))
There are some breadcrumbs here [Documentation for SPQR qrfact() · Issue #12298 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/12298)
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 18, 2017, 11:59am UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/18 "2017-05-18T11:59:56Z")
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`.
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [May 18, 2017, 1:00pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/19 "2017-05-18T13:00:56Z")
Yes, an issue would be good here.
---
### Author:  [@rleegates](https://discourse.julialang.org/u/rleegates)
#### Post date: [May 18, 2017, 1:13pm UTC](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645/20 "2017-05-18T13:13:25Z")
Issue [#21937](https://github.com/JuliaLang/julia/issues/21937)
[Next page](https://discourse.julialang.org/t/sparse-svd-smallest-svs/3645.md?page=2)