Hello there,
I’m trying to solve the generalized eigen problem Ax=\lambda Bx with A and B being linear maps. However, I encountered a problem when using a shift-invert method (A-\sigma B)^{-1}Bx=\frac{1}{\lambda-\sigma}x to solve the system: Even when considering a symmetric matrices (see example with Laplacian) the eigen modes have a non-zero imaginary part, that “takes the shape” of the real part. I don’t see why this is happening and how I can avoid this, especially when the problem is more complicated and I can’t just ignore the imaginary part. Here’s the problem:
using DiffEqOperators, LinearMaps
n = 50
B=complex.(sparse(DerivativeOperator{Float64}(2,2,1.,n,:Dirichlet0, :Dirichlet0)))
A=complex.(speye(n))
a,u=eig(full(A),full(B))
@show isreal(u)
@show isreal(a)
σ = a[1]
Alo=LinearMap{Complex128}((y,x)->A_mul_B!(y,A,x),n)
Blo=LinearMap{Complex128}((y,x)->A_mul_B!(y,B,x),n)
C = A-σ*B
C_lu = lufact(C)
Clo = LinearMap{Complex128}((y,x)->A_ldiv_B!(y,C_lu,Blo*x),n)
ae,ue = eigs(Clo,nev=5)
@show isreal(ae)
@show isreal(ue)
@show maximum(imag.(ue))
Thanks for any help, there is probably an easy explanation or workaround!?
Ideally this might be correctable, but you would have to specify that the LinearMap is Hermitian in the constructor.
However, the “correction” would just be postprocessing the computed result to have the desired form. This is because eigs (i.e. Arpack) has no special algorithms for Hermitian complex problems. Thus there are arbitrary complex factors for the eigenvectors, and the eigenvalues will not be exactly real. There has been discussion of adding such postprocessing to eigs, but it hasn’t been done yet AFAIK.
This is in contrast to the case of real symmetric operators (with the appropriate flag in the constructor), where there is a change of algorithm and purely real results.
Also, note that using σ = a[1] makes the problem (nearly) singular so the results will have less accuracy than is possible.
I don’t think that is exactly the problem I’m having. The eigenvalues/vectors are not just complex up to a certain inaccuracy. Here’s what the first eigenvectors look like in this example. Note: the eigenvalues are the same up to a certain accuracy.
ARPACK returns the correct eigenvector, with an arbitrary phase: ue = exp(im*theta)*u.
2 Likes