`howmany` argument is not working properly

[GiggleLiu]

The following example should return the lowest eigen value according to the docstring, but returns full spectrum. This behavior also apply for SparseMatrixCSC type input.

julia> eigsolve(randn(4,4), 1, :SR)
(Complex{Float64}[-2.312+0.0im, -0.462406+1.1235im, -0.462406-1.1235im, 2.3231+0.0im], Array{Complex{Float64},1}[[-0.100722+0.0im, -0.827104+0.0im, -0.519301+0.0im, 0.189947+0.0im], [-0.872454-0.0771637im, 0.0691546-0.482002im, 0.415766-0.00552139im, 0.280353-0.280638im], [-0.872454+0.0771637im, 0.0691546+0.482002im, 0.415766+0.00552139im, 0.280353+0.280638im], [-0.12619+0.0im, 0.521181+0.0im, -0.999292+0.0im, -0.192739+0.0im]], ConvergenceInfo: 4 converged values after 1 iterations and 4 applications of the linear map;
norms of residuals are given by (3.766032907369262e-32, 5.7580543111933805e-33, 5.7580543111933805e-33, 3.955301283972102e-33).
)

help?> KrylovKit.eigsolve
  eigsolve(A::AbstractMatrix, [howmany = 1, which = :LM, T = eltype(A)]; kwargs...)
  eigsolve(f, n::Int, [howmany = 1, which = :LM, T = Float64]; kwargs...)
  eigsolve(f, x₀, [howmany = 1, which = :LM]; kwargs...)
  eigsolve(f, x₀, howmany, which, algorithm)

  Compute howmany eigenvalues from the linear map encoded in the matrix A or by the function f. Return eigenvalues, eigenvectors and a ConvergenceInfo structure.

[Jutho]

This is deliberate, if more eigenvalues have converged at the same cost, all of them are returned. I thought this was a useful property, but I am open to changes.

If less have converged, the requested number of eigenvalues is returned, also those that have not converged. You should always check info.converged to see how many have converged.

[GiggleLiu]

I would prefer changing this behavior. When I used this function for the first time, I got the lowest eigenvalue as expected, so I turn the first output into scalar using E[]. But soon, I found the program errors at this operation occasionally due to the unknown number of converged eigenvectors. This can be a dangerous feature.

If there are some practical concerns that when one expects k eigenvalues, obtaining more eigenvalues can be better, you can provide a keyword argument to let user notice what they are doing.

[Jutho]

You can of course also just do E[1], which is how I typically use it. Or even

(e,), (vec,), info = eigsolve(...)

For safety, you would anyhow need to check info.converged to see if the eigenvalue is converged, so I don't think you can every be sure to have exactly k converged eigenvalues.

[GiggleLiu]

It would be appreciated if you can provide a using case that more than expected eigenvalues can be useful.

[Jutho]

I have to admit that I find your tone a bit aggressive. I can equally well ask you for a use case or minimal code where any of the solutions that I proposed don't work, or where the current behaviour is really problematic.

My initial motivation was that, if you have a real matrix, eigenvalues can be complex conjugate pairs and both of them converge equally fast. It is absurd to return one and not the other only because this doesn't match with howmany. Say you request one eigenvalue of :LM, but it happens to be a complex conjugate pair, which one should be returned? So at least it should be allowed to return one more eigenvalue as requested for real problems. But then, if more eigenvalues happen to be converged, why not also return them. You anyway have to check how many have converged.

[GiggleLiu]

Sorry for making you feeling bad, I never meant to offend you. I do think this issue is a proper place to explain your motivation.

So you mean user get not only the required number of eignevectors (or more), but also the convergence information. That's not bad, maybe changing the argument name to 'kmin' would fit your motivation best, meanwhile explicit. How do you like it?

BTW: Let me explain my using case a bit, for a length 1 vector, you can get the item (in numpy, a.item()) using either E[1] or E[], the latter will error when the length of vector is not 1. Normally I would prefer the API E[], it is safer to use. This is why I got an error here.

[Jutho]

I certainly agree that the documentation, and maybe the argument name, can be improved. I'll look into that.

Normally I would prefer the API E[], it is safer to use.

Why is E[] safer than E[1] ?

[GiggleLiu]

Thanks for your consideration, I could make a PR if you want.

Why is E[] safer than E[1] ?

Because E[] will throw an error immediately if E is not a scalar or a single element array. Throwing an explicit error at early stage is considered as safe in the case when you know what you want is a scalar.

[Jutho]

But, given that E can contain more values (i.e. suppose the documentation was already improved), but if you are interested in only 1, you are guaranteed that it is on position one, then there is nothing wrong or unsafe about E[1] I suppose?

[Jutho]

Actually, the current documentation reads:

*   `vals`: a `Vector` containing the eigenvalues, of length at least `howmany`, but could
    be longer if more eigenvalues were converged at the same cost. Eigenvalues will be real
    if [`Lanczos`](@ref) was used and complex if [`Arnoldi`](@ref) was used (see below).

So I think this is pretty accurate.

[Jutho]

Would you like to see other changes, or can I close this issue?

[GiggleLiu]

How about changing the name howmany to a more accurate one? Since many people do not read the descriptions carefully.

Also, kwargs should error if a keyword argument is not supported, in case one doesn't read the interface description carefully or input some typo. Like me, I once used the interface like eigsolve(A, howmany=1) without getting any warning.

[Jutho]

Any suggestions for a better name? Not something that involves an arbitrary letter like kmin. I do still think howmany is pretty clear and most people won't suffer from having more, especially if you just need one.

Yes, the keyword arguments should error, although, then again, it is sometimes convenient to pass around a whole bunch of keyword arguments even though functions that they don't apply to. But I guess that's not a clean coding practice.

[GiggleLiu]

Maybe at_least_howmany? Since this is not a keyword argument, long name is preferred.

I think the most valuable coding experience is "always writing correct codes", since debugging is the most time consuming and frustrating part of programming. Thanks for your open mind to these interface changes that help users away from mindless errors.

[nrontsis]

Just to mention that howmany is pretty clear to me as well and ArnoldiMethod.jl can also return more eigenvalues than those requested from the user while having a much less informative argument name (nev).

Personally, I find this very useful for e.g. Trust Region Subproblems where the rightmost eigenvalue of a related eigenproblem corresponds to the global solution of the TRS and the second rightmost can correspond to a local solution of the TRS.

[GiggleLiu]

@nrontsis Thanks for your information. I think there is a convention in other communities, and has become a standard, e.g.

Matlab: https://www.mathworks.com/help/matlab/ref/eigs.html

Scipy: https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.sparse.linalg.eigsh.html

A package should either follow the conventions, or state very clearly about the differences, otherwise, many people will get confused. According to your example, returning more eigenvalues can be useful, so this is a good change, but should be clearly stated. nev is bad, howmany is better but not reflecting the difference to k (according the existing standard, kmin is not too had for people switching from python or matlab). Please help us find a better name.

[Jutho]

Maybe at_least_howmany? Since this is not a keyword argument, long name is preferred.

A long name is preferred? That's the first time I hear that style rule. No underscores though is a rule that is currently attempted in Julia Base. This name still shows up in the code, and the documentation string.

I do think you are taking it a bit far; hoping that one can use a package without reading the documentation. Yes, there is something to be said for mimicking existing conventions/APIs. However, some APIs are just bad/ugly. I've never understood what the c is for in ncv. Did somebody think it is Crylov?

And then again, this package aims to be much more general than eigs in matlab/python, by also doing linear solvers and other methods, as well as accepting general vector-like objects instead of restricting to Vectors. If you don't like my changing of existing names, you will be shocked by the interface of linsolve with GMRES. Instead of restart I also call it what it is, namely, krylovdim. I'd rather have a consistent interface than an ancient interface.

The decision for not making howmany and which keyword arguments was that I wanted a consistent set of keyword arguments for the different methods, whereas the normal (positional) arguments specify the problem itself, and the keyword arguments provide settings for the algorithm. This might not have been the best decision, since I also tend to forget the order.

In latest master, you will get errors for unknown keyword arguments, and I will tag a new release soon. Unless there is another viable name suggestion, I will close this issue.

[GiggleLiu]

One last comment, maybe finding a better name is too hard. What about changing 'Compute howmany eigenvalues' to 'Compute at least howmany eigenvalues' in docstring?

[Jutho]

Absolutely, that is a very good suggestion.

[Jutho]

Master has been updated/fixed. I will tag a new minor release later today.

[Jutho]

New release tagged, with some delay, but with new functionality. I think this can be closed