Explicit transform for eigs - bugfix generalized eigenvalue problems

[jarlebring]

Here is an attempt to resolve the bug in JuliaLang/LinearAlgebra.jl#488 (or rather the bug reported in the comments in JuliaLang/LinearAlgebra.jl#488). As pointed out in the issue discussion, eigs is currently broken for generalized eigenvalue problems unless B is positive or negative definite, e.g., it normally does not work for non-symmetric B.

The code introduces an explicittransform kwarg which provides the possibility to do shift-and-invert in the julia code, or (as before) let ARPACK do it. If you do shift-and-invert in the julia code the problem does not occur. The explicittransform kwarg defaults to :AUTO which tries to determine if it is a good idea to let ARPACK do the shift-and-invert or do it in julia-code. This way the previous behaviour is accessible to any user using explicittransform=false.

I have tried to keep the code such that the meaning of which and sigma is independent if explicittransform is true or false, which created a bit more if-statements than one might expect, but the advantage is that you get consistent meaning like this:

julia> A=sprandn(100,100,0.1); e1=eye(100,1); A=A+100*e1*e1';

julia> C=I+0.0001*sprandn(100,100,0.1);

julia> B=C*C';  # symm. pos. definite 

julia> λ1=eigs(A,B,explicittransform=false,sigma=1.0)[1]
6-element Array{Complex{Float64},1}:
0.8853783560528526 - 0.2251061711212125im
0.8853783560528526 + 0.2251061711212125im
1.4198361576177376 + 0.0im               
0.5786082353083483 - 0.5069352320198903im
0.5786082353083483 + 0.5069352320198903im
1.6872265416932886 + 0.0im               

julia> λ2=eigs(A,B,explicittransform=true,sigma=1.0)[1]
6-element Array{Complex{Float64},1}:
0.8853783560528521 + 0.2251061711212136im
0.8853783560528521 - 0.2251061711212136im
1.4198361576177394 + 0.0im               
0.5786082353083479 + 0.5069352320198891im
0.5786082353083479 - 0.5069352320198891im
1.6872265416932901 + 0.0im               

After hopefully some discussion on this here I can add it to the function description.

[jarlebring]

I forgot explanation: I had to change

matvecA!(y, x) = mul!(y, A, x)

to

matvecA! = (y, x) -> mul!(y, A, x)

because I wanted to change the function at a later point. Maybe someone knows if this change can have implications on performance?

[jarlebring]

It seems it passes all unit tests except this:

srand(127) 
n = 10
k = 3
A = randn(n,n); A = A'A
B = randn(n,k); B = B*B'
@test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k])

It fails when it tries to sort complex numbers. With explicittransformation=true we get eltype Complex128 and with explicittransformation=false we get eltype Float64.

It can be easily "fixed" by changing the test to

@test sort(eigs(A, B, nev = k, sigma = 1.0,explicittransform=false)[1]) ≈ sort(eigvals(A, B)[1:k])

I do realize that returning different types based on kwarg-values is text-book type instability. The only fix I can think of is to take another parameter to eigs is a bit intrusive. Something like this:

eigs(A,B;params)=_eigs(Complex128,A,B;params...);
eigs(T,A,B;params)=_eigs(T,A,B;params...);
function _eigs{T}(::Type{T},A,B;nev::Integer=6, ncv::Integer=max(20,2*nev+1), which=:LM,
            tol=0.0, maxiter::Integer=300, sigma=nothing, v0::Vector=zeros(eltype(A),(0,)),
            ritzvec::Bool=true, explicittransform=:AUTO)
...

Not sure how severe is this problem is. Ideas appreciated.

[haampie]

eigs is already having this type instability due to the ritzvec = true/false parameter, so I don't think it's a real problem.

[jarlebring]

Ah. Yes indeed. Also, the older discussion in JuliaLang/LinearAlgebra.jl#279 suggests that making it type stable is non-trivial.

[StefanKarpinski]

If the values of this argument are true, false and :AUTO it seems like it would perhaps be better to do true, false and nothing where nothing chooses automatically?

[StefanKarpinski]

Slightly better to put a Union{Bool, Nothing} type annotation on this then 😀

[jarlebring]

Sounds better indeed. Should be fixed now.

I'm not so experienced in contributing in this way, want me to add tests? Squash commits?

[andreasnoack]

@jarlebring Thanks for preparing this PR Sorry for the lag of response. It's on my todo but I've had to look into other tasks before getting to this. Hopefully, I'll find the time for a review later this week.

[andreasnoack]

@jarlebring As you might have noticed, we have moved ARPACK bindings to Arpack.jl (which partly has contributed to the delay in reviewing this PR) such that developement of our ARPACK wrappers won't be tied to the developement cycles of the main repo. Would you mind moving this PR to Arpack.jl instead. Then we should be able to merge it shortly.

[jarlebring]

Sounds good. See new PR in JuliaLinearAlgebra. Hope I understood you correctly.