nalgebra icon indicating copy to clipboard operation
nalgebra copied to clipboard

Eigendecomposition for square matrix

Open demurgos opened this issue 2 years ago • 7 comments

Hello, I have a non-symmetric square metric and would like to get an eigendecomposition. As far as I see, only symmetric matrices are supported (with symmetric_eigen).

Is eigendecomposition for regular square matrices supported? Am I missing something?

If you need an example, I'm trying to get a decomposition for the following stochastic matrix:

[0.47  0.53  0   ]
[0.52  0     0.48]
[0     0     1   ]

demurgos avatar Jun 07 '23 16:06 demurgos

What you want is the Schur decomposition yes?

https://docs.rs/nalgebra/latest/nalgebra/base/struct.Matrix.html#method.schur

If this doesn't work and you need the classical eigen decomp then the Lapack version of the crate has: https://docs.rs/nalgebra-lapack/latest/nalgebra_lapack/struct.Eigen.html

aujxn avatar Jun 10 '23 20:06 aujxn

How can we implement this in rust though, seems a bit too complex tbh

yyadav1120 avatar Jun 22 '23 18:06 yyadav1120

I believe the current implementation of the Schur Decomposition uses an Arnodli iteration so it is already implemented in Nalgebra! This Schur decomposition is a generalization of the Spectral Decomposition.

I believe the Lapack implementation uses Householder transformations? The Lapack API and algorithm is well studied and documented if you want to RWIR. The Nalgebra-Lapack bindings are quite simple to use, though, and works great for integrating Rust code with existing numerical implementations.

I see that @sebcrozet added some tags to this issue, I'm wondering what kind of algorithms or enhancements you would like to see in this area? Maybe a simpler API for generalized eigen problems using the available tools?

aujxn avatar Jun 24 '23 00:06 aujxn

Nalgebra already has a method to calculate the eigen values of a matrix. Why not use it to find the eigenvalues and then create the diagonal matrix corresponding and return it or there is some problem?

spoloxs avatar Jun 25 '23 08:06 spoloxs

I would prefer a pure Rust verrsion for compatibility with Wasm. For now, I'm able to use nalgebra_lapack but I get slightly different results between nalgebra_lapack, ndarray-linalg and custom code to retrieve the vectors. I'm not super familiar with the different algorithms so I'm not sure if those are just precision issues or errors.

demurgos avatar Sep 02 '23 16:09 demurgos

get slightly different results between nalgebra_lapack, ndarray-linalg

This is likely a difference in implementation and probably isn't a bug, especially if the results are similar. Also note, if you have repeated eigenvalues the decomposition isn't unique, so this could be why you get different results.

There are a lot of things to consider here when choosing an API for your task. For instance, you mention that you want the entire decomposition for a non-symmetric matrix, but not all non-symmetric matrices are diagonalizable unless you allow for the complex values.

Computing eigenvectors for ill-conditioned non-symmetric matrices is a complicated task and there are many specialized algorithms and parameters built into lapack for it, see this and this from the lapack guide.

ndarray-linalg probably uses _geev (docs), I don't think there exists a pure rust implementation but the tools to build a naive version definitely exist in these libraries. It would be cool to try and put together an implementation in nalgebra based off what is documented above. If anyone wants help with that I can help when I have time. If no one has done it yet, it is likely because generalizations such as the SVD or Schur Decomposition can often be used instead and have consistent, stable, and fast algorithms. If you provide more context for your problem I might be able to suggest a way to use these existing tools instead.

If your matrices are small, do the easy thing taught in into linear by computing the eigenvalues, subtract from diagonal, and then solve for Av=0. Using the available solve and eigs functions you could stay in pure rust this way, just be aware of the footguns like:

  • needing complex entries
  • could have multiple eigenvectors for an eigenvalue
  • if you have repeated eigenvalues then the decomposition is not unique, as any orthogonal basis of the associated eigenspace are valid eigenvectors

aujxn avatar Sep 02 '23 19:09 aujxn