Generalized Eigenvalue Decomposition in C#

/ cesarsouza / 16 Comments
5f7de24799f03b87862a0cdad7b0ed9e

For quite some time I have been trying to find a Java or C# implementation of the Generalized Eigenvalue Decomposition. I just wanted something simple, that didn’t depend on any external libraries such as BLAS or LAPACK. As I could not find any simple and understandable implementation of this decomposition, I decided to implement one myself.

Introduction

As Wikipedia says, a generalized Eigen value problem is the problem of finding a vector v that obeys

Amathbf{v} = lambda B mathbf{v} quad quad

where A and B are matrices. If v obeys this equation, with some λ, then we call v the generalized eigenvector of A and B, and λ is called the generalized eigenvalue of A and B which corresponds to the generalized eigenvector v. The possible values of λ must obey the following equation

det(A - lambda B)=0.,

An interesting feature of the generalized eigenvalue decomposition is that it finds the eigenvectors of the matrix B-1A even if B is singular and does not have an inverse. If B is nonsingular, the problem could be solved by reducing it to a standard eigenvalue problem of the form B-1Ax=λx. However, because B can be singular, an alternative algorithm, called the QZ method, is necessary.

Solving the problem using the QZ method is equivalent to computing the Eigen decomposition of B-1A without the need of inverting B, which could be impossible or ill-conditioned if B is singular or near-singular. It is also computationally simpler as it saves the time and memory needed for inverting and storing the inverse of B.

EISPACK

EISPACK is a software library for numerical computation of eigenvalues and eigenvectors of matrices, written in FORTRAN. It was originally written around 1972–1973, based heavily on algorithms originally implemented in ALGOL. Because those algorithms are much more human-readable than their LAPACK counterparts, I decided to port the related functions from EISPACK.

The functions which implement the Generalized Eigenvalue Decomposition in EISPACK are called QZHES, QZIT, QZVAL and QZVEC. As their name implies, they use the QZ method for finding the generalized eigenvalues of a matrix pair (A,B).

Source code

The source code is available in the download link in the upper part of this article. It can also be found in the latest versions of the Accord.NET Framework. The algorithm presented in the code is equivalent to the eig(A,B) command of MATLAB.

Because the code has been ported from FORTRAN, many sections of the algorithms have stayed as they appear in their original functions. Please be aware that it may contain lots of labels and goto’s, as most structured programming features we take for granted today were not commonplace at the time.

Using the code

Using the code is much simpler than dealing with LAPACK functions. Because the functions have been implemented directly in C#, the code is somewhat more human readable and has also been wrapped in the same nice object-oriented approach as the other matrix decompositions that can be found in MAPACK and JAMA.

Using the given sources, the following MATLAB code:

Becomes equivalent to the C# code:

in the sense that both satisfy the identity A*V = B*V*D .

References

16 Comments

  2. Hi Miguel,

    Sorry about that. The code is part of the Accord.NET Framework, and the Special class is contained in the Accord.Math namespace to host special math functions. When I posted the code I must have forgotten about those dependencies. Anyway, the code for the GeneralizedEigenvalueDecomposition and the Special classes can also be found on Google Code in the aforementioned links. You may also download the Accord.NET Framework if you wish.

    Hope it helps!

    Best regards,
    César

    Reply

  4. Hi Anonymous,

    You can try putting the Identity matrix. However, in this case, it would be better to use the standard Eigenvalue decomposition instead. The standard EVD is available in Accord.NET or in the Mapack for .NET from Lutz Roeder.

    Best regards,
    César

    Reply

  5. Hi Cesar,

    Accord.NET looks great! One question — is there a way in Accord to do an iterative eigenvalue decomposition? Specifically, I am just looking for the N-largest eigenvalues (and eigenvectors). I am only seeking this for symmetric matrices. I think there is a method called the power iteration, which can be much faster than doing the global decomposition, when N is small enough compared to the size of the matrix. I guess I could use Accord to write this myself, but I thought I’d ask you first. Thanks!

    -Dave

    Reply

  6. Hi dsv,

    Thanks for the interest in the library! Currently there are no iterative algorithms for finding Eigenvectors, though their addition would be really nice. If you are looking for the N-largest Eigenvectors, I would suggest using Lanczsos algorithm. The power method may seems simpler and works nice for finding the largest vector, but I am under the impression it is often tough to converge, specially if you have two equally largest Eigenvectors.

    One of the most promising variations of the Lanczos algorithm is the “Implicitly Restarted Lanczos Method” which is available in ARPACK. Since ARPACK is licensed under BSD, it would be possible to translate the method from ARPACK’s DSEUPD function (if you are really feeling adventurous with Fortran!). Otherwise Wikipedia’s page for the Lanczos method also have some pseudocode for the algorithm.

    Regards,
    Cesar

    Reply

  9. Dear Cesar

    Thanks for quick reply last time, I have been using the code fairly successfully but I still get slightly different results compared to the matlab method. The matrices are quite large (300 x 300) would this cause a problem.

    Thanks
    Guy

    Reply

  10. hi i am new to c, and i need the code for genralized eigen values in C++ , or C#, plz tell me how can i run ur code in my project.

    Regards

    Reply
  12. Pingback: Matrix manipulation using Accord.NET | ~/cesarsouza/blog

Leave a Reply to flora Cancel reply

Your email address will not be published. Required fields are marked *