14. Eigensystems

This chapter describes functions for computing eigenvalues and eigenvectors of matrices. There are routines for real symmetric, real nonsymmetric, and complex hermitian matrices. Eigenvalues can be computed with or without eigenvectors. The hermitian matrix algorithms used are symmetric bidiagonalization followed by QR reduction. The nonsymmetric algorithm is the Francis QR double-shift.

These routines are intended for “small” systems where simple algorithms are acceptable. Anyone interested in finding eigenvalues and eigenvectors of large matrices will want to use the sophisticated routines found in LAPACK. The Fortran version of LAPACK is recommended as the standard package for large-scale linear algebra.

The functions described in this chapter are declared in the header file ‘gsl_eigen.h’.

14.1 Real Symmetric Matrices

Function: gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n)

This function allocates a workspace for computing eigenvalues of n-by-n real symmetric matrices. The size of the workspace is .

Function: void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w)

This function frees the memory associated with the workspace w.

Function: int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w)

This function computes the eigenvalues of the real symmetric matrix A. Additional workspace of the appropriate size must be provided in w. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The eigenvalues are stored in the vector eval and are unordered.

Function: gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n)

This function allocates a workspace for computing eigenvalues and eigenvectors of n-by-n real symmetric matrices. The size of the workspace is .

Function: void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w)

This function frees the memory associated with the workspace w.

Function: int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w)

This function computes the eigenvalues and eigenvectors of the real symmetric matrix A. Additional workspace of the appropriate size must be provided in w. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The eigenvalues are stored in the vector eval and are unordered. The corresponding eigenvectors are stored in the columns of the matrix evec. For example, the eigenvector in the first column corresponds to the first eigenvalue. The eigenvectors are guaranteed to be mutually orthogonal and normalised to unit magnitude.

14.2 Complex Hermitian Matrices

Function: gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n)

This function allocates a workspace for computing eigenvalues of n-by-n complex hermitian matrices. The size of the workspace is .

Function: void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w)

This function frees the memory associated with the workspace w.

Function: int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval, gsl_eigen_herm_workspace * w)

This function computes the eigenvalues of the complex hermitian matrix A. Additional workspace of the appropriate size must be provided in w. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The imaginary parts of the diagonal are assumed to be zero and are not referenced. The eigenvalues are stored in the vector eval and are unordered.

Function: gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n)

This function allocates a workspace for computing eigenvalues and eigenvectors of n-by-n complex hermitian matrices. The size of the workspace is .

Function: void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w)

This function frees the memory associated with the workspace w.

Function: int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval, gsl_matrix_complex * evec, gsl_eigen_hermv_workspace * w)

This function computes the eigenvalues and eigenvectors of the complex hermitian matrix A. Additional workspace of the appropriate size must be provided in w. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The imaginary parts of the diagonal are assumed to be zero and are not referenced. The eigenvalues are stored in the vector eval and are unordered. The corresponding complex eigenvectors are stored in the columns of the matrix evec. For example, the eigenvector in the first column corresponds to the first eigenvalue. The eigenvectors are guaranteed to be mutually orthogonal and normalised to unit magnitude.

14.3 Real Nonsymmetric Matrices

The solution of the real nonsymmetric eigensystem problem for a matrix involves computing the Schur decomposition where is an orthogonal matrix of Schur vectors and , the Schur form, is quasi upper triangular with diagonal -by- blocks which are real eigenvalues of , and diagonal -by- blocks whose eigenvalues are complex conjugate eigenvalues of . The algorithm used is the double shift Francis method.

Function: gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n)

This function allocates a workspace for computing eigenvalues of n-by-n real nonsymmetric matrices. The size of the workspace is .

Function: void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w)

This function frees the memory associated with the workspace w.

Function: void gsl_eigen_nonsymm_params (const int compute_t, const int balance, gsl_eigen_nonsymm_workspace * w)

This function sets some parameters which determine how the eigenvalue problem is solved in subsequent calls to gsl_eigen_nonsymm.

If compute_t is set to 1, the full Schur form will be computed by gsl_eigen_nonsymm. If it is set to 0, will not be computed (this is the default setting). Computing the full Schur form requires approximately 1.5-2 times the number of flops.

If balance is set to 1, a balancing transformation is applied to the matrix prior to computing eigenvalues. This transformation is designed to make the rows and columns of the matrix have comparable norms, and can result in more accurate eigenvalues for matrices whose entries vary widely in magnitude. See Balancing for more information. Note that the balancing transformation does not preserve the orthogonality of the Schur vectors, so if you wish to compute the Schur vectors with gsl_eigen_nonsymm_Z you will obtain the Schur vectors of the balanced matrix instead of the original matrix. The relationship will be where Q is the matrix of Schur vectors for the balanced matrix, and D is the balancing transformation. Then gsl_eigen_nonsymm_Z will compute a matrix Z which satisfies with . Note that Z will not be orthogonal. For this reason, balancing is not performed by default.

Function: int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval, gsl_eigen_nonsymm_workspace * w)

This function computes the eigenvalues of the real nonsymmetric matrix A and stores them in the vector eval. If is desired, it is stored in the upper portion of A on output. Otherwise, on output, the diagonal of A will contain the -by- real eigenvalues and -by- complex conjugate eigenvalue systems, and the rest of A is destroyed. In rare cases, this function will fail to find all eigenvalues. If this happens, an error code is returned and the number of converged eigenvalues is stored in w->n_evals. The converged eigenvalues are stored in the beginning of eval.

Function: int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval, gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w)

This function is identical to gsl_eigen_nonsymm except it also computes the Schur vectors and stores them into Z.

Function: gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n)

This function allocates a workspace for computing eigenvalues and eigenvectors of n-by-n real nonsymmetric matrices. The size of the workspace is .

Function: void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w)

This function frees the memory associated with the workspace w.

Function: int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval, gsl_matrix_complex * evec, gsl_eigen_nonsymmv_workspace * w)

This function computes eigenvalues and right eigenvectors of the n-by-n real nonsymmetric matrix A. It first calls gsl_eigen_nonsymm to compute the eigenvalues, Schur form , and Schur vectors. Then it finds eigenvectors of and backtransforms them using the Schur vectors. The Schur vectors are destroyed in the process, but can be saved by using gsl_eigen_nonsymmv_Z. The computed eigenvectors are normalized to have Euclidean norm 1. On output, the upper portion of A contains the Schur form . If gsl_eigen_nonsymm fails, no eigenvectors are computed, and an error code is returned.

Function: int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval, gsl_matrix_complex * evec, gsl_matrix * Z, gsl_eigen_nonsymmv_workspace * w)

This function is identical to gsl_eigen_nonsymmv except it also saves the Schur vectors into Z.

14.4 Sorting Eigenvalues and Eigenvectors

Function: int gsl_eigen_symmv_sort (gsl_vector * eval, gsl_matrix * evec, gsl_eigen_sort_t sort_type)

This function simultaneously sorts the eigenvalues stored in the vector eval and the corresponding real eigenvectors stored in the columns of the matrix evec into ascending or descending order according to the value of the parameter sort_type,

GSL_EIGEN_SORT_VAL_ASC

ascending order in numerical value

GSL_EIGEN_SORT_VAL_DESC

descending order in numerical value

GSL_EIGEN_SORT_ABS_ASC

ascending order in magnitude

GSL_EIGEN_SORT_ABS_DESC

descending order in magnitude

Function: int gsl_eigen_hermv_sort (gsl_vector * eval, gsl_matrix_complex * evec, gsl_eigen_sort_t sort_type)

This function simultaneously sorts the eigenvalues stored in the vector eval and the corresponding complex eigenvectors stored in the columns of the matrix evec into ascending or descending order according to the value of the parameter sort_type as shown above.

Function: int gsl_eigen_nonsymmv_sort (gsl_vector_complex * eval, gsl_matrix_complex * evec, gsl_eigen_sort_t sort_type)

This function simultaneously sorts the eigenvalues stored in the vector eval and the corresponding complex eigenvectors stored in the columns of the matrix evec into ascending or descending order according to the value of the parameter sort_type as shown above. Only GSL_EIGEN_SORT_ABS_ASC and GSL_EIGEN_SORT_ABS_DESC are supported due to the eigenvalues being complex.

14.5 Examples

The following program computes the eigenvalues and eigenvectors of the 4-th order Hilbert matrix, .

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int
main (void)
{
  double data[] = { 1.0  , 1/2.0, 1/3.0, 1/4.0,
                    1/2.0, 1/3.0, 1/4.0, 1/5.0,
                    1/3.0, 1/4.0, 1/5.0, 1/6.0,
                    1/4.0, 1/5.0, 1/6.0, 1/7.0 };

  gsl_matrix_view m 
    = gsl_matrix_view_array (data, 4, 4);

  gsl_vector *eval = gsl_vector_alloc (4);
  gsl_matrix *evec = gsl_matrix_alloc (4, 4);

  gsl_eigen_symmv_workspace * w = 
    gsl_eigen_symmv_alloc (4);
  
  gsl_eigen_symmv (&m.matrix, eval, evec, w);

  gsl_eigen_symmv_free (w);

  gsl_eigen_symmv_sort (eval, evec, 
                        GSL_EIGEN_SORT_ABS_ASC);
  
  {
    int i;

    for (i = 0; i < 4; i++)
      {
        double eval_i 
           = gsl_vector_get (eval, i);
        gsl_vector_view evec_i 
           = gsl_matrix_column (evec, i);

        printf ("eigenvalue = %g\n", eval_i);
        printf ("eigenvector = \n");
        gsl_vector_fprintf (stdout, 
                            &evec_i.vector, "%g");
      }
  }

  gsl_vector_free (eval);
  gsl_matrix_free (evec);

  return 0;
}

Here is the beginning of the output from the program,

$ ./a.out 
eigenvalue = 9.67023e-05
eigenvector = 
-0.0291933
0.328712
-0.791411
0.514553
...

This can be compared with the corresponding output from GNU OCTAVE,

octave> [v,d] = eig(hilb(4));
octave> diag(d)  
ans =

   9.6702e-05
   6.7383e-03
   1.6914e-01
   1.5002e+00

octave> v 
v =

   0.029193   0.179186  -0.582076   0.792608
  -0.328712  -0.741918   0.370502   0.451923
   0.791411   0.100228   0.509579   0.322416
  -0.514553   0.638283   0.514048   0.252161

Note that the eigenvectors can differ by a change of sign, since the sign of an eigenvector is arbitrary.

The following program illustrates the use of the nonsymmetric eigensolver, by computing the eigenvalues and eigenvectors of the Vandermonde matrix with .

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int
main (void)
{
  double data[] = { -1.0, 1.0, -1.0, 1.0,
                    -8.0, 4.0, -2.0, 1.0,
                    27.0, 9.0, 3.0, 1.0,
                    64.0, 16.0, 4.0, 1.0 };

  gsl_matrix_view m 
    = gsl_matrix_view_array (data, 4, 4);

  gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
  gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);

  gsl_eigen_nonsymmv_workspace * w = 
    gsl_eigen_nonsymmv_alloc (4);
  
  gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);

  gsl_eigen_nonsymmv_free (w);

  gsl_eigen_nonsymmv_sort (eval, evec, 
                           GSL_EIGEN_SORT_ABS_DESC);
  
  {
    int i, j;

    for (i = 0; i < 4; i++)
      {
        gsl_complex eval_i 
           = gsl_vector_complex_get (eval, i);
        gsl_vector_complex_view evec_i 
           = gsl_matrix_complex_column (evec, i);

        printf ("eigenvalue = %g + %gi\n",
                GSL_REAL(eval_i), GSL_IMAG(eval_i));
        printf ("eigenvector = \n");
        for (j = 0; j < 4; ++j)
          {
            gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
            printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
          }
      }
  }

  gsl_vector_complex_free(eval);
  gsl_matrix_complex_free(evec);

  return 0;
}

Here is the beginning of the output from the program,

$ ./a.out 
eigenvalue = -6.41391 + 0i
eigenvector = 
-0.0998822 + 0i
-0.111251 + 0i
0.292501 + 0i
0.944505 + 0i
eigenvalue = 5.54555 + 3.08545i
eigenvector = 
-0.043487 + -0.0076308i
0.0642377 + -0.142127i
-0.515253 + 0.0405118i
-0.840592 + -0.00148565i
...

This can be compared with the corresponding output from GNU OCTAVE,

octave> [v,d] = eig(vander([-1 -2 3 4]));
octave> diag(d)
ans =

  -6.4139 + 0.0000i
   5.5456 + 3.0854i
   5.5456 - 3.0854i
   2.3228 + 0.0000i

octave> v
v =

 Columns 1 through 3:

  -0.09988 + 0.00000i  -0.04350 - 0.00755i  -0.04350 + 0.00755i
  -0.11125 + 0.00000i   0.06399 - 0.14224i   0.06399 + 0.14224i
   0.29250 + 0.00000i  -0.51518 + 0.04142i  -0.51518 - 0.04142i
   0.94451 + 0.00000i  -0.84059 + 0.00000i  -0.84059 - 0.00000i

 Column 4:

  -0.14493 + 0.00000i
   0.35660 + 0.00000i
   0.91937 + 0.00000i
   0.08118 + 0.00000i

Note that the eigenvectors corresponding to the eigenvalue are slightly different. This is because they differ by the multiplicative constant which has magnitude 1.

14.6 References and Further Reading

Further information on the algorithms described in this section can be found in the following book,

• G. H. Golub, C. F. Van Loan, Matrix Computations (3rd Ed, 1996), Johns Hopkins University Press, ISBN 0-8018-5414-8.

The LAPACK library is described in,

• LAPACK Users' Guide (Third Edition, 1999), Published by SIAM, ISBN 0-89871-447-8.

  http://www.netlib.org/lapack

The LAPACK source code can be found at the website above along with an online copy of the users guide.

This document was generated by Autobuild on September, 26 2007 using texi2html 1.78.