enabled complex eigenvectors for lapack

[geoeo]

As mentioned in #1105
complex eigenvector computation is actually implemented but the complex values are simply discarded.

I changed the signature of complex_eigenvalues() to maybe return the eigenvectors.
Important! : I do not return the complex type because it would require an n^2 iteration according to: http://www.netlib.org/lapack/explore-html/d9/d8e/group__double_g_eeigen_ga66e19253344358f5dee1e60502b9e96f.html
If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + iVL(:,j+1) and u(j+1) = VL(:,j) - iVL(:,j+1).

I simply return the eigenvector matrix as is. Any user interested in the complex eigenvectors would have to do this themselves. Usually one would be interested in the real vectors (at least for me). This should probably be documented and is somewhat at odds with the current interface.

Furthermore I replaced the macro calls if lapack_info!with lapacl_panic! due to typing issues

[geoeo]

@Andlon Is there any update on this issue? Im stuck compiling against my local nalgebra version

[geoeo]

@metric-space Maybe you are the better ping since you implemented the current function

[Andlon]

@geoeo: Sorry, I did not mean to ignore you for so long. I'm not too familiar with the nalgebra-lapack code, and certainly not with the complex aspect.

I'm also confused - the original issue in #1105 states that the eigenvalue computation crashes. However, this PR completely changes the signature of the method. Wouldn't it be better to separately fix the computation of the eigenvalues and have a separate routine for the eigenvectors?

[geoeo]

@Andlon Yes. Maybe my attempt to explain this was very convoluted, since I wasnt sure what was going on at the start. I fixed two issues I had with the codebase:

1. The main issue is that complex eigenvectors can not be returned by nalgebra-lapack
2. The crash by eigenvalue

I am completely open to putting this a new function. I just wanted to get some feedback on my code and have this being tracked.

[geoeo]

I was just reviewing the code and I dont think its a good idea to split up the method into two. The problem is that as stated above:
I do not return the complex type because it would require an n^2 iteration according to: http://www.netlib.org/lapack/explore-html/d9/d8e/group__double_g_eeigen_ga66e19253344358f5dee1e60502b9e96f.html

So in order to check if an eigenvector is complex I query the corresponding eigenvalue. Its similar to the lapack method which has the same signature. There is also only one method for that. I can rename it though since the name is no longer descriptive

[geoeo]

@Andlon could you give feedback if this would be mergable or what is still neccessary?

[metric-space]

@geoeo apologies for the radio silence to your ping. Haven't had the time to look at this (still don't ). That said please note I am just a contributor and am not a maintainer like @Andlon or @sebcrozet

[Andlon]

@geoeo: I'm sorry but I don't really know the code base of nalgebra-lapack and don't feel particularly well-suited to review. IIRC someone other than @sebcrozet initially contributed the crate, but has since probably left. I think nalgebra-lapack is technically in dire need of one or more maintainers.

I'd also say the current API is very strange, as I don't understand why complex_eigenvalues is a method on Eigen, which explicitly says that it represents an eigendecomposition of a real square matrix with real eigenvalues.

[sebcrozet]

@geoeo, @Andlon I will take care of this next week!

[geoeo]

@Andlon i think a real square matrix can have complex eigenvales in conjugate pairs

[Andlon]

@geoeo: indeed it can, but the Eigen struct on which the method is defined explicitly says that the eigenvalues are real. But maybe just the docs for Eigen need to be updated, I don't know the full context.

[geoeo]

@sebcrozet hey. Any update on this?

[sebcrozet]

@geoeo Thanks for this PR! I pushed some change to:

• Combine both ::new and ::complex_eigen_decomposition into a single ::new method.
• Get rid of the match to avoid unnecessary near-duplicate code.
• Add the imaginary part of the eigenvalues to the Eigen struct.
• Add a method to check if all the eigenvalues stored in Eigen are real.

Please, let me know if you have any issue with these changes.

[sebcrozet]

Important! : I do not return the complex type because it would require an n^2 iteration according to: http://www.netlib.org/lapack/explore-html/d9/d8e/group__double_g_eeigen_ga66e19253344358f5dee1e60502b9e96f.html
If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + iVL(:,j+1) and u(j+1) = VL(:,j) - iVL(:,j+1).

Now that the complex eigenvalues are stored into Eigen, perhaps we could add a method that computes these complex eigenvectors by running this n² iteration?

[geoeo]

As far as I can tell everything is working fine. I can have a look at the computation of the explicit complex eigenvectors.

[geoeo]

@sebcrozet Im having trouble with the generics for a complex number.

    /// Returns a tuple of vectors. The elements of the tuple are the complex eigenvalues, complex left eigenvectors and complex right eigenvectors respectively. 
    /// The elements appear as conjugate pairs within each vector, with the positive of the pair always being first.
    pub fn get_complex_elements(&self) -> (Option<Vec<Complex<T>>>, Option<Vec<OVector<Complex<T>, D>>>, Option<Vec<OVector<Complex<T>, D>>>) where DefaultAllocator: Allocator<Complex<T>, D> {
        match !self.eigenvalues_are_real() {
            true => (None, None, None),
            false => {
                let number_of_elements = self.eigenvalues_re.nrows();
                let number_of_complex_entries = self.eigenvalues_im.iter().fold(0, |acc, e| if !e.is_zero() {acc + 1} else {acc});
                let mut eigenvalues = Vec::<Complex<T>>::with_capacity(2*number_of_complex_entries);
                let mut eigenvectors = match self.eigenvectors.is_some() {
                    true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)),
                    false => None
                };
                let mut left_eigenvectors = match self.left_eigenvectors.is_some() {
                    true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)),
                    false => None
                };

                let eigenvectors_raw = self.eigenvectors;
                let left_eigenvectors_raw = self.left_eigenvectors;

                for mut i in 0..number_of_elements {
                    if self.eigenvalues_im[i] != T::zero() {
                        //Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first.
                        eigenvalues.push(Complex::<T>::new(self.eigenvalues_re[i].clone(), self.eigenvalues_im[i].clone()));
                        eigenvalues.push(Complex::<T>::new(self.eigenvalues_re[i].clone(), -self.eigenvalues_im[i].clone()));

                        if eigenvectors.is_some() {
                            let mut r1_vec = OVector::<Complex<T>, D>::zeros(number_of_elements);
                            let mut r1_vec_conj = OVector::<Complex<T>, D>::zeros(number_of_elements);

                            for j in 0..number_of_elements {
                                r1_vec[j] = Complex::<T>::new(self.eigenvectors.unwrap()[(i,j)].clone(),self.eigenvectors.unwrap()[(i,j+1)].clone());
                                r1_vec_conj[j] = Complex::<T>::new(self.eigenvectors.unwrap()[(i,j)].clone(),-self.eigenvectors.unwrap()[(i,j+1)].clone());
                            }
    
                            eigenvectors.unwrap().push(r1_vec);
                            eigenvectors.unwrap().push(r1_vec_conj);
                        }


                        if left_eigenvectors.is_some() {
                            //TODO: Do the same for left
                        }


                        i += 1;
                    }
                    
                }
                (Some(eigenvalues), left_eigenvectors, eigenvectors)
            }
        }

OVector::<Complex<T>, D>::zeros(number_of_elements) wont compile since zero does not exist for complex

[metric-space]

@geoeo

The error I'm getting is

error[E0599]: the function or associated item `zeros` exists for struct `Matrix<Complex<T>, D, Const<1>, <DefaultAllocator as na::allocator::Allocator<Complex<T>, D>>::Buffer>`, but its trait bounds were not satisfied
   --> nalgebra-lapack/src/eigen.rs:201:78
    |
201 | ...                   let mut r1_vec_conj = OVector::<Complex<T>, D>::zeros();
    |                                                                       ^^^^^ function or associated item cannot be called on `Matrix<Complex<T>, D, Const<1>, <DefaultAllocator as na::allocator::Allocator<Complex<T>, D>>::Buffer>` due to unsatisfied trait bounds
    |
    = note: the following trait bounds were not satisfied:
            `D: DimName

I believe the error has more to do with your dimension constraints than the non-existence of the method

from the below
https://github.com/dimforge/nalgebra/blob/dev/src/base/dimension.rs#L221

nalgebra/src/base/construction.rs

Line 647 in 74c4aa9

impl<T: Scalar, R: DimName, C: DimName> OMatrix<T, R, C>

and Eigen has type constraint of Dim the base type of DimName

One way out is to replace Dim with DimName

It's been a while so I could be wrong

[geoeo]

Im careful with chaning the constraints. That whole system is a but opaque to me. But if its the only way then so be it.

[sebcrozet]

@geoeo Switching to DimName won’t work because it will prevent the use of the eigendecomposition with dynamically-sized matrices. You will need to use ::zeros_generic (which is compatible with just Dim) instead of ::zeros. Here are some lines you can change:

187: let (number_of_elements, _) = self.eigenvalues_re.shape_generic();
202: for mut i in 0..number_of_elements.value() {
209: let mut r1_vec = OVector::<Complex<T>, D>::zeros_generic(number_of_elements, Const::<1>);
210: let mut r1_vec_conj = OVector::<Complex<T>, D>::zeros_generic(number_of_elements, Const::<1>);

[geoeo]

@sebcrozet @metric-space Thanks a lot. I pushed an initial version that compiles. But I would want to write a test for it. Do you know of a "standard" matrix which produces complex eigenvalues that I can check against?

[geoeo]

I found this test matrix from https://textbooks.math.gatech.edu/ila/1553/complex-eigenvalues.html

$$ \left(\begin{array}{cc} \frac{4}{5} & -\frac{3}{5} & 0.0\\ \frac{3}{5} & \frac{4}{5} & 0.0 \\ 1 & 2 & 2 \end{array}\right) $$

Ill park this here

[geoeo]

@sebcrozet Feel free to review. I wrote an example/test for eigenvalues and (right) eigenvectors. I had to change the Cargo.toml to make it build.

[sebcrozet]

Thank you for your perseverance @geoeo. I think we are good to merge this version.