eigs does not return complex eigenvalues if matrix is not complex

[KCErb]

require "linalg"

m = LA::GMat[[3, -2],[4, -1]]
vals, vecs = m.eigs

res = vecs.inv * m * vecs
res.detect
res.flags.diagonal?

The last line of the above is false because it turns out that the vecs are all Float64 instead of complex as they should be. It seems this is because LAPACK method detection uses the type of the matrix, so only complex matrices are handed off to the complex solver.

By the way, the flags feature is awesome! I love that diagonal? sugar there :tada:

[konovod]

Yes, LAPACK can return eigenvectors either in complex form or as a real matrix - if you check vecs, you can see that imag part is still returned, but on a next cell. Quote from LAPACK docs:

If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR.
If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
v(j+1) = VR(:,j) - i*VR(:,j+1)

I thought that this form have some deep meaning, so decided to keep this behavior for real matrices and return complex values for Complex matrices. Now when i think about it looks like it is just a hack to save space.

So, there are 3 variants of what to return

1. Return Complex | Float vectors depending on whether complex eigenvalues present. Matlab (and perhaps scipy?) go this way, but in a statically typed language this could be annoying - most of time user will have to check result type for further computations.
2. Always return Complex vectors, they just sometimes have zero imaginary part. User will have to convert them to real if he wants.
3. Return Float vectors packing complex values in two sequential cells. Positive side is that if user knows that all eigenvalues are real, he don't need to do anything.

Maybe 2 should be default with 3 provided in something like eish_real. What do you think?

[KCErb]

I know what you mean. Currently the method returns an Array of complex or an array of float for the eigenvalues, and it has been a little annoying to have to handle both.

In my code I was toying with introducing an alias type "EigenSystem" which was just the tuple return from eigs. So I could imagine structs EigenSystemFloat64 and EigenSystemComplex which inherit from a parent with methods float?, complex?, eigenvalues, vals, eigenvectors, vecs and using something like

eigen_system = mat.eigs
if eigen_system.complex?
  work_with_complex(eigen_system.vals)
else
  work_with_float(eigen_system.vals)
end

Then if a user doesn't want to be bothered maybe then can just say mat.eigs(complex: true) and then type is guaranteed. (And complex: false for the eish / packing strategy).

That all said, option 2 is appealing. I wouldn't have been surprised at all if when I first ran eigs it had returned a bunch of 2.3 + 0.0i type numbers. In that case, I'd be looking for some helper methods like mathematica's chop.

vals, vecs = mat.eigs
# vecs is a GMatComplex
vecs.chop!
# vecs is now a GMat if imaginary parts are small

Note

It seems that Mathematica is the only language with such a function, but it's one of the first ones I implement in new projects in other languages. The basic idea is that you just take Type::EPSILON and compare the absolute value. If it's less, then you send it to zero. And if it's complex, and the imaginary part is sent to zero, you cast to a real number.

Edit -- or was it Type::DIGITS ... I think it depends on the use case. i.e. if num.round(Type::DIGITS).zero? just send it to zero.

[konovod]

After trying to tackle with this.

1. complex: true\false approach would look best, but that doesn't work because result type doesn't depend from argument value. So it will still return Matrix(T) | Matrix(Complex) no matter what user pass. I can do eigs returning union and also eigs_complex and eigs_real methods that return guaranted type. But names aren't very intuitive.

2. as for chop - I like the idea. Not sure about naming, but if Mathematica uses it perhaps no need to invent another name. vecs.chop! won't work though (it is impossible to change type of calling object) so it can be used like

   if v = vecs.chop
   # here v is real
   else
   # here v is nil, so use vecs
   end

   This carries small overhead though, as Complex matrices take 2x more memory (4x for Float32) and they would be allocated only to be converted back to float.

[KCErb]

Cool.

I guess I misunderstood the complex: true\false possibility. I thought that LAPACK could return in either complex or real form regardless of input. So we were basically letting the user choose the return form.

Could eigs_complex and eigs_real be accessed via eigs(complex:true) for example? I don't love the name eigs_real but I don't mind if if I know I can get the same behavior via eigs(complex: false) or eigs(real: true) or what have you.

As for chop, good point about the bang method, that would change the type. Is it anti-semantic to do something in the vein of not_nil!? I feel like sometimes in crystal we use ! as an insistence which is the use-case I imagine for calling chop. Like "look, I know the result has to be real" so the result would be

vecs = vecs.chop! and we get a runtime error if an imaginary part exists.

If you are interested in such a direction though, the name real! feels a little nicer.

The main use-case for chop in my experience has had less to do with turning complex into real and more to do with turning 1.499999999999 into 1.5 for logical comparison. So maybe if we want to keep chop as a feature it would just be:

a = Math.cos(2*Math::PI/3) # => -0.4999999999999998 : Float64
a.chop # => -0.5 : Float64

I've actually thought about proposing this as a language feature. And ... in the process of writing this comment I came up with the following example. Can I get your feedback before submitting it over there?

struct Float64
  def chop
    self.round(DIGITS)
  end
end

neg_half = Math.cos(2*Math::PI/3) # one epsilon greater than -0.5
pos_half = Math.cos(Math::PI/3) # one epsilon greater than 0.5

(neg_half + pos_half).abs <= Float64::EPSILON # => false, two epsilons off unfortunately
(neg_half.chop + pos_half.chop).abs <= Float64::EPSILON # => true, in fact you could chop either one

[konovod]

eigs(complex:true) from a library point of view is eigs(complex: Bool?). So it's return type is Matrix(T) | Matrix(Complex), it won't change depending on what value actually was passed by a user. So you will have to do

vals, vecs = eigs(complex: true)
# here vecs is Matrix(T)  | Matrix(Complex), even though we are sure that it is complex
vecs = vecs.as(GMatComplex)
# here we can use vecs as Matrix(Complex)

This can be solved by making eigs a macros, not usual function, but i think this will create more problems than solve, right now all macroses are hidden in internals of library.

not_nil! still can't change a type of variable - it only cut a cases off the union (due to compiler cleverness). So

a = 1+0.i # a is Complex
a.to_real!
# a become Float from now on

is impossible. Well, until you make to_real! a macros, but i'm not sure it's a good idea.

I've actually committed "chop" just now) For a Matrix, I use tolerance that is (relatively arbitrarily) evaluated as (maximum element value)*nrows*ncolumns*10eps. As for numbers... this isn't easy question. 0.1 can't be represented exactly, so i'm not sure about cornercases. Is it possible that it will be chopped to something like 0.099999... after adding two epsilons. Otherwise looks fine, but if you are going to actually propose it to crystal-lang, be prepared to conservative answers like "it can live in a shard" or "we need real use case" 😄

[KCErb]

1. Oh I see now, I guess I assumed there was some way to restrict return types. Thanks for explaining.
2. Good point, I'll do a little more thinking about corner cases before drafting up a proper proposal.