initial svds support based on eigs

[jaak-s]

Added support for sparse truncated SVD based on eigs. Keyword inputs are the same as for eigs. Returns left singular vectors, singular values, right singular vectors (and also outputs from eigs). Added a test case.

[StefanKarpinski]

My sense is that truncated SVD should be an option to the svd function and the fact that the svds sparse SVD function was used for this is a weird historical design artifact in Matlab that we should not copy.

[StefanKarpinski]

In fact, I'm not even sure a separate svds function should exist.

[jaak-s]

Sure, I've no preference on the structuring/naming.

Just was currently missing the ability to run Arnoldi (ARPACK) methods to extract first singular vectors of large matrices that have sparse or special structure.

[ViralBShah]

We do need a different name though since svd on sparse matrices will do the full svd when we have it. We could have a named argument, but the interface adds too many options for my liking to combine it all in svd.

[ViralBShah]

@jaak-s could you add the documentation as well?

[ViralBShah]

Cc @Jutho

[tkelman]

Also remove the trailing whitespace from the comment in the tests

[ViralBShah]

Should the trailing whitespace be discussed in CONTRIBUTING.md along with configuration options for emacs/vim and in git?

For example, I now have this in my .emacs thanks to some help in one of the issues:

;; Delete trailing whitespace
(add-hook 'before-save-hook 'delete-trailing-whitespace)

[tkelman]

We should have this conversation in a separate issue, but yes I think we should move it to several places, and not do the check in Travis. Code style nitpicking should not prevent running the build and tests. I think we should pick an arbitrary throwaway CI competitor to Travis, like Drone or CircleCI or Shippable, and do the whitespace check there as a separate status indicator. Or put together something custom on Heroku or whatever, but using an existing service seems easier.

[ViralBShah]

BTW, the bikeshedding should not block this PR. Let's get the functionality in along with tests and docs. The naming issue applies to eigs also, and is a longstanding issue.

[Jutho]

Looks good to me. The tests are only ok for input matrices with real numbers (which is the case for the input matrix A), otherwise it should be sign(...) \ ... or conj(sign(...)) * ...

[andreasnoack]

I think we'd like to parametrize this type on the input element and matrix type (and use capital letters), i.e. something like

type SVDOperator{T,S}
data::S
end

[andreasnoack]

I addition, we'll have to support more than Float64 and the type should reflect that. The type can probably not be subtype of AbstractMatrix either, as the input might not be a subtype of AbstractMatrix.

[ViralBShah]

@andreasnoack Allocation-free version of eigs is discussed in #7907.

[ViralBShah]

Octave takes the approach of using [0 A; A' 0], and there is then some cleanup afterwards, IIRC. I am happy to look it up and describe the method here so that someone else can do a cleanroom implementation - unless you already know what needs to be done.

[andreasnoack]

@ViralBShah I'm not sure about how to handle the cleanup in general. The [0 A; A' 0] operator is singular when m!=n so some zeros would have to be removed when looking for the small singular values.

[jaak-s]

Pushed few commits: manual for svds and improved handling of ritzvec parameter.

The comments from @andreasnoack are spot on.

Regarding type information I don't have a good solution yet because the SVDOperator is only an internal type and not visible to users. It seems some additional tricks are needed. Currently, svds(x) uses duck typing, i.e. x only needs to implement functions size and *. What about adding extra parameter to svds specifying whether Float64 or Complex64 is used for eigs?

I have not checked how to do post-processing for [0 A; A' 0]. Any pointers? If it easy I can do it. If it is more complex then we can use the current version as a place holder until the other is ready.

[andreasnoack]

Currently, svds(x) uses duck typing, i.e. x only needs to implement functions size and *. What about adding extra parameter to svds specifying whether Float64 or Complex64 is used for eigs?

I don't think it is a limitation to require that the operator type has a parameter for the element type. The linear operators in Base (the matrices) have and I think it would feel natural to have an element type for user defined linear operators as well. If we assume that then we already have the element type in the operator type. We'd also like to have that in order to have fast dispatch for * and \ because sometimes it is necessary with many iterations.

I have not checked how to do post-processing for [0 A; A' 0]. Any pointers?

Not sure, but maybe it is discussed somewhere in the ARPACK manual, but that is not the most joyful read. Let's browse the literature and see what we can find.

[ViralBShah]

I don't think the ARPACK manual sheds light on that, and even if it does, it will be impossible to find it. IIRC, the only place where I see discussion of svd is in the examples, and they do A'*A there, but it is only a demonstration, not a recommended way.

[ViralBShah]

Here is what octave does. I suggest that @jaak-s not look at the code, since it is GPL, and I am looking at it and describing what it does. For others look for svds.m in the octave source.

Octave calls eigs on [0 A; A' 0]. The comment in the codebase says:

    ## We wish to exclude all eigenvalues that are less than zero as these
    ## are artifacts of the way the matrix passed to eigs is formed. There
    ## is also the possibility that the value of sigma chosen is exactly
    ## a singular value, and in that case we're dead!! So have to rely on
    ## the warning from eigs. We exclude the singular values which are
    ## less than or equal to zero to within some tolerance scaled by the
    ## norm since if we don't we might end up with too many singular
    ## values.

The norm that they use is the largest singular value (if the sigma input variable is a character) or normest (uses a power series analysis to estimate the norm). Basically, they return all singular values where abs(s) <= tol.

Other things they do are scaling the input by the largest value (1-norm) to make things a bit more stable. If sigma is specified, that also needs scaling. If sigma is 0, then they compute twice as many eigenvalues. Quote:

      ## Find the smallest eigenvalues
      ## The eigenvalues returns by eigs for sigma=0 are symmetric about 0.
      ## As we are only interested in the positive eigenvalues, we have to
      ## double k and then throw out the k negative eigenvalues.
      ## Separately, if sigma is non-zero, but smaller than the smallest
      ## singular value, ARPACK may not return k eigenvalues. However, as
      ## computation scales with k we'd like to avoid doubling k for all
      ## scalar values of sigma.

[ViralBShah]

@alanedelman What are your thoughts here on doing A'*A vs. [0 A; A' 0]?

[jaak-s]

I will soon add type information to SVDOperator. Users with x::AbstractArray only need to do svds(x). For custom types (duck typing) users probably have to wrap it themselves by SVDOperator.

@ViralBShah Thank you for information. It is now quite clear how to extract largest singular values from [0 A; A' 0]. For the smallest it is not so clear yet for me. Also is there information how to extract singular vectors (left or right) from that solution?

[ViralBShah]

@jaak-s The singular vectors computed just as @Jutho described.

m, n = size(A)
V, s = eigs(SVDOperator(A), ....)
ind = <list of indices into s, after filtering the singular values>
u  = sqrt(2) * V[1:m, ind]
v = sqrt(2) * V[m+1:end, ind]

[ViralBShah]

Bump @jaak-s

[jaak-s]

Thanks for the bump, got busy during the holidays. I'll try to put patches together in the next couple of days.

I tried out the [0 A; A' 0] approach you and @Jutho outlined and it seems to work well.

[jaak-s]

Now svds uses [0 A; A' 0] approach, saves type information in SVDOperator and works also with complex operators. The doc has not been updated. Plan is to finalize the docs after you all are fine with the patch.

Main changes:

• using approach by @ViralBShah to extract singular vectors
• changed input parameter nev to nsv (number of singular values)
• svds only returns largest singular values. The option for the smallest singular values is currently removed because [0 A; A' 0] creates extra small singular values and post-processing is not so clear. If someone knows the exact steps let me know, if not we can leave it for future ToDo.
• When using duck typed X in svds(X) following methods are required: size(X), eltype(X), X * vector, X' * vector.

[ViralBShah]

I think we can just simplify this to use just eltype(X) and avoid getOperatorType() altogether. The checks for whether eltype() is BlasFloat should ideally be in eigs. It should be possible to make SVDOperator free of types then.

[ViralBShah]

For good performance, we should be using A_mul_B!, but that requires some work on eigs. For now, nothing much can be done here.

[tkelman]

travis doesn't like the trailing whitespace here or on line 151

[ViralBShah]

We also need to filter out the singular values that are zero or close to zero.

[jaak-s]

The reason for introducing getOperatorType was to support (sparse) binary and integer matrices. Their eltype is bool/integer but they can be multiplied by float vectors. Basically changing the operator type to Float64 for integer matrix allows eigs to work on it. Note that svd works with integer matrices and I found it is useful to have out-of-the-box support for them in svds too (also because in my work I use svds on binary matrices).

If you wish I can remove getOperatorType.

White-space should be now fixed.

[jaak-s]

Regarding removing close to zero singular values, I would prefer the current behavior. I think users like to get N values when they set nsv=N, even if last of them are 0. This is also the behavior of svds in Octave and in Scipy (scipy.sparse.linalg.svds).

[ViralBShah]

I agree the integer matrix support is good to have. I suggest that instead of getOperatorType in svds, we should fix this for eigs. To support integer matrices, it is better to convert them to float in eigs, since the matrix-vector products are significantly faster (due to use of BLAS in the dense case). Would be great if you can add it to eigs here, or I can do it later.

[jaak-s]

Sounds good, I will replace getOperatorType with eltype then. Since I am not very familiar with code of eigs I'd prefer to leave it to you.

[ViralBShah]

About the filtering out of singular values that are 0 - octave does do it. But, we can do it once we get this basic piece of code checked in - once you remove getOperatorType.

[jaak-s]

getOperatorType is now removed, also added check for nsv.

I wonder what singular values does Octave remove because they return 0-valued singular values (version 3.8.1):

octave:1> svds( zeros(4,3), 2 )
ans =

   0
   0

[ViralBShah]

This was in the octave source that I quoted above. Of course there are legitimate zero singular values. Reading this carefully, they are doing a filter on negative values, and perhaps something in the case where sigma is specified. I will look deeper.

    ## We wish to exclude all eigenvalues that are less than zero as these
    ## are artifacts of the way the matrix passed to eigs is formed. There
    ## is also the possibility that the value of sigma chosen is exactly
    ## a singular value, and in that case we're dead!! So have to rely on
    ## the warning from eigs. We exclude the singular values which are
    ## less than or equal to zero to within some tolerance scaled by the
    ## norm since if we don't we might end up with too many singular
    ## values.

[ViralBShah]

Travis failure on linux appears unrelated.

[ViralBShah]

I think this is a great start with svds and we will refine it as we go along.

[ViralBShah]

Added to NEWS in 47061a5

[jaak-s]

Thank you.

I forgot to update the docs earlier but now pushed a commit that reflects the latest version of svds. Can you merge it again or should I create a new pull request?

[ViralBShah]

It looks like your last commit here is still 8 days ago. Could you create a new PR?

Also, we should have squashed these commits, but it escaped my mind before merging.

[jaak-s]

Created PR #9723 for docs.