# Generalized eigen problem using LinearMap?

**URL:** <https://discourse.julialang.org/t/generalized-eigen-problem-using-linearmap/12168>\
**Category:** Numerics\
**Created:** [July 4, 2018, 7:47pm UTC](https://discourse.julialang.org/t/generalized-eigen-problem-using-linearmap/12168 "2018-07-04T19:47:21Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [July 4, 2018, 7:47pm UTC](https://discourse.julialang.org/t/generalized-eigen-problem-using-linearmap/12168/1 "2018-07-04T19:47:21Z")

</div>

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:

```julia
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!?

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [July 5, 2018, 12:26am UTC](https://discourse.julialang.org/t/generalized-eigen-problem-using-linearmap/12168/2 "2018-07-05T00:26:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [July 5, 2018, 5:06am UTC](https://discourse.julialang.org/t/generalized-eigen-problem-using-linearmap/12168/3 "2018-07-05T05:06:28Z")

</div>

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.  
 ![ev](https://global.discourse-cdn.com/julialang/original/3X/a/e/aec674af0f1dc283ab9a0c3e9da5191925e4c33e.png)

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [July 5, 2018, 6:40am UTC](https://discourse.julialang.org/t/generalized-eigen-problem-using-linearmap/12168/4 "2018-07-05T06:40:47Z")

</div>

ARPACK returns the correct eigenvector, with an arbitrary phase: ue = exp(im\*theta)\*u.