incorrect eigs(A,B) results for indefinite B

[fgerick]

I encountered a strange issue when solving the generalized eigenvalue problem in Julia 0.6 (version on mac os, but also on 0.5 and 0.6 in Juliabox). When trying to solve the generalized eigenvalue problem: Ax = λBx with julia, like so:

A = randn(10,10)
B = randn(10,10);
λ = eigs(A,B)[1]
λ2 = eigs(A,B)[1]
λ≈λ2
false

The outcome is actually changing with every call to eigs(A,B) (also with simpler matrices). In the case of B being a diagonal Matrix (so that the generalized eigenvalue problem reduces to the standard eigenvalue problem Ax = λx) the results of eigs(A,B) are consistent. This problem does not occur, e.g. in scipy.linalg .

Cheers, Felix

[stevengj]

If you compare to eigvals(A,B), you'll see that eigs(A,B)[1] is giving results that are just plain wrong.

@andreasnoack, does eigs support the generalized eigenproblem with indefinite B? The definite case, eigs(A,B'B)[1], appears to work fine.

[antoine-levitt]

From ARPACK: c\Name: dnaupd c c\Description: c Reverse communication interface for the Implicitly Restarted Arnoldi c iteration. This subroutine computes approximations to a few eigenpairs c of a linear operator "OP" with respect to a semi-inner product defined by c a symmetric positive semi-definite real matrix B.

So julia should either check that B is spd and throw if not, or transform to B^-A.

[andreasnoack]

The manual recommends transforming the generalized problem to B\A when B is not psd as @antoine-levitt also says. I introduced this problem in https://github.com/JuliaLang/julia/pull/17873 where I removed the isposdef check of B to allow for singular B as long as a shift was provided. The problem is how to fix this issue while still allowing for the singular case.

[antoine-levitt]

Maybe as a post-process? As far as I can see ARPACK returns B-normalized vectors iff B is positive semidefinite. Maybe ARPACK has more info in INFO or RESID.

[stevengj]

do we need am ispossemidef function?

[StefanKarpinski]

isspd?

[stevengj]

isspd seems like it would mean "is symmetric positive-definite (SPD)", not semidefinite.

[stevengj]

Or just a semidef=true argument to isposdef would work too, I guess.

[antoine-levitt]

I thought that too, but I wonder if that can have a meaningful distinction from isposdef with round-off error.

Looking at the results from ARPACK, in the case of a non-semidefinite B matrix, ARPACK reports success: INFO is set to 0, the number of arnoldi iterations is still smaller than the maximum number, but the residuals are rubbish (order 1). In the use case of https://github.com/JuliaLang/julia/pull/17873 the residuals were order 1e-8. How about a check on the residuals, with a warning when they exceed eps^1/3*norm(A) or similar? Or will that flag false positives?

[andreasnoack]

but I wonder if that can have a meaningful distinction from isposdef with round-off error.

My concern as well. For rank deficient problems the criterion would probably show slight indefiniteness so a tolerance would be required.

How about a check on the residuals, with a warning when they exceed eps^1/3*norm(A) or similar? Or will that flag false positives?

If something like that works reliably, I'd have expected ARPACK to report it in an info code.

[andreasnoack]

Continued from #17866

Is the fix to essentially do what the snippet above says

Yes, and transform back after the calculation. The current solver should probably only be chosen when the B matrix is semidefinite and a positive shift has been specified. As mentioned above, a non-symmetric solver should be used on B\A when B otherwise. So there is some branching to handle here and some work to find proper test cases but only really hard part should be detecting the semidefinite case.

[ViralBShah]

We could perhaps allow a keyword argument to specify so as a stopgap?

[jarlebring]

The difference in the results of two calls to eigs in the issue description is probably due to how ARPACK handles starting vectors. It uses random starting vectors by default, AFAIK. Changing the starting vector can change to which eigenvalues ARPACK converges:

A = randn(10,10)
B = randn(10,10)
λ = eigs(A,B,v0=ones(10))[1]
λ2 = eigs(A,B,v0=ones(10))[1]
λ3 = eigs(A,B,v0=[1;zeros(9)])[1]
λ≈λ2
true
λ3≈λ
false

The bugs reported in the issue comments still seem relevant.

[jarlebring]

For completeness: eigs calls ARPACKs dnaupd.f which calls dnaupd2.f which in turn by default initiates a starting vector by calling dgetv0:

c
c\Name: dgetv0
c
c\Description: 
c  Generate a random initial residual vector for the Arnoldi process.
c  Force the residual vector to be in the range of the operator OP.  

[andreasnoack]

Now tracked in https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/15