diff options
Diffstat (limited to 'gsl-1.9/doc/eigen.texi')
| -rw-r--r-- | gsl-1.9/doc/eigen.texi | 462 |
1 files changed, 462 insertions, 0 deletions
diff --git a/gsl-1.9/doc/eigen.texi b/gsl-1.9/doc/eigen.texi new file mode 100644 index 0000000..e2c291a --- /dev/null +++ b/gsl-1.9/doc/eigen.texi @@ -0,0 +1,462 @@ +@cindex eigenvalues and eigenvectors +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. + +@cindex LAPACK, recommended for linear algebra +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 +@sc{lapack}. The Fortran version of @sc{lapack} is recommended as the +standard package for large-scale linear algebra. + +The functions described in this chapter are declared in the header file +@file{gsl_eigen.h}. + +@menu +* Real Symmetric Matrices:: +* Complex Hermitian Matrices:: +* Real Nonsymmetric Matrices:: +* Sorting Eigenvalues and Eigenvectors:: +* Eigenvalue and Eigenvector Examples:: +* Eigenvalue and Eigenvector References:: +@end menu + +@node Real Symmetric Matrices +@section Real Symmetric Matrices +@cindex symmetric matrix, real, eigensystem +@cindex real symmetric matrix, eigensystem + +@deftypefun {gsl_eigen_symm_workspace *} gsl_eigen_symm_alloc (const size_t @var{n}) +This function allocates a workspace for computing eigenvalues of +@var{n}-by-@var{n} real symmetric matrices. The size of the workspace +is @math{O(2n)}. +@end deftypefun + +@deftypefun void gsl_eigen_symm_free (gsl_eigen_symm_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun int gsl_eigen_symm (gsl_matrix * @var{A}, gsl_vector * @var{eval}, gsl_eigen_symm_workspace * @var{w}) +This function computes the eigenvalues of the real symmetric matrix +@var{A}. Additional workspace of the appropriate size must be provided +in @var{w}. The diagonal and lower triangular part of @var{A} are +destroyed during the computation, but the strict upper triangular part +is not referenced. The eigenvalues are stored in the vector @var{eval} +and are unordered. +@end deftypefun + +@deftypefun {gsl_eigen_symmv_workspace *} gsl_eigen_symmv_alloc (const size_t @var{n}) +This function allocates a workspace for computing eigenvalues and +eigenvectors of @var{n}-by-@var{n} real symmetric matrices. The size of +the workspace is @math{O(4n)}. +@end deftypefun + +@deftypefun void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun int gsl_eigen_symmv (gsl_matrix * @var{A}, gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_symmv_workspace * @var{w}) +This function computes the eigenvalues and eigenvectors of the real +symmetric matrix @var{A}. Additional workspace of the appropriate size +must be provided in @var{w}. The diagonal and lower triangular part of +@var{A} are destroyed during the computation, but the strict upper +triangular part is not referenced. The eigenvalues are stored in the +vector @var{eval} and are unordered. The corresponding eigenvectors are +stored in the columns of the matrix @var{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. +@end deftypefun + +@node Complex Hermitian Matrices +@section Complex Hermitian Matrices + +@cindex hermitian matrix, complex, eigensystem +@cindex complex hermitian matrix, eigensystem + +@deftypefun {gsl_eigen_herm_workspace *} gsl_eigen_herm_alloc (const size_t @var{n}) +This function allocates a workspace for computing eigenvalues of +@var{n}-by-@var{n} complex hermitian matrices. The size of the workspace +is @math{O(3n)}. +@end deftypefun + +@deftypefun void gsl_eigen_herm_free (gsl_eigen_herm_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun int gsl_eigen_herm (gsl_matrix_complex * @var{A}, gsl_vector * @var{eval}, gsl_eigen_herm_workspace * @var{w}) +This function computes the eigenvalues of the complex hermitian matrix +@var{A}. Additional workspace of the appropriate size must be provided +in @var{w}. The diagonal and lower triangular part of @var{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 +@var{eval} and are unordered. +@end deftypefun + +@deftypefun {gsl_eigen_hermv_workspace *} gsl_eigen_hermv_alloc (const size_t @var{n}) +This function allocates a workspace for computing eigenvalues and +eigenvectors of @var{n}-by-@var{n} complex hermitian matrices. The size of +the workspace is @math{O(5n)}. +@end deftypefun + +@deftypefun void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun int gsl_eigen_hermv (gsl_matrix_complex * @var{A}, gsl_vector * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_hermv_workspace * @var{w}) +This function computes the eigenvalues and eigenvectors of the complex +hermitian matrix @var{A}. Additional workspace of the appropriate size +must be provided in @var{w}. The diagonal and lower triangular part of +@var{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 @var{eval} and are unordered. The corresponding +complex eigenvectors are stored in the columns of the matrix @var{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. +@end deftypefun + +@node Real Nonsymmetric Matrices +@section Real Nonsymmetric Matrices +@cindex nonsymmetric matrix, real, eigensystem +@cindex real nonsymmetric matrix, eigensystem + +The solution of the real nonsymmetric eigensystem problem for a +matrix @math{A} involves computing the Schur decomposition +@tex +\beforedisplay +$$ +A = Z T Z^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = Z T Z^T +@end example + +@end ifinfo +where @math{Z} is an orthogonal matrix of Schur vectors and @math{T}, +the Schur form, is quasi upper triangular with diagonal +@math{1}-by-@math{1} blocks which are real eigenvalues of @math{A}, and +diagonal @math{2}-by-@math{2} blocks whose eigenvalues are complex +conjugate eigenvalues of @math{A}. The algorithm used is the double +shift Francis method. + +@deftypefun {gsl_eigen_nonsymm_workspace *} gsl_eigen_nonsymm_alloc (const size_t @var{n}) +This function allocates a workspace for computing eigenvalues of +@var{n}-by-@var{n} real nonsymmetric matrices. The size of the workspace +is @math{O(2n)}. +@end deftypefun + +@deftypefun void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun void gsl_eigen_nonsymm_params (const int @var{compute_t}, const int @var{balance}, gsl_eigen_nonsymm_workspace * @var{w}) +This function sets some parameters which determine how the eigenvalue +problem is solved in subsequent calls to @code{gsl_eigen_nonsymm}. + +If @var{compute_t} is set to 1, the full Schur form @math{T} will be +computed by @code{gsl_eigen_nonsymm}. If it is set to 0, +@math{T} will not be computed (this is the default setting). Computing +the full Schur form @math{T} requires approximately 1.5-2 times the +number of flops. + +If @var{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 @ref{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 @code{gsl_eigen_nonsymm_Z} you will obtain the Schur +vectors of the balanced matrix instead of the original matrix. The +relationship will be +@tex +\beforedisplay +$$ +T = Q^t D^{-1} A D Q +$$ +\afterdisplay +@end tex +@ifinfo + +@example +T = Q^t D^(-1) A D Q +@end example + +@end ifinfo +@noindent +where @var{Q} is the matrix of Schur vectors for the balanced matrix, and +@var{D} is the balancing transformation. Then @code{gsl_eigen_nonsymm_Z} +will compute a matrix @var{Z} which satisfies +@tex +\beforedisplay +$$ +T = Z^{-1} A Z +$$ +\afterdisplay +@end tex +@ifinfo + +@example +T = Z^(-1) A Z +@end example + +@end ifinfo +@noindent +with @math{Z = D Q}. Note that @var{Z} will not be orthogonal. For +this reason, balancing is not performed by default. +@end deftypefun + +@deftypefun int gsl_eigen_nonsymm (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_eigen_nonsymm_workspace * @var{w}) +This function computes the eigenvalues of the real nonsymmetric matrix +@var{A} and stores them in the vector @var{eval}. If @math{T} is +desired, it is stored in the upper portion of @var{A} on output. +Otherwise, on output, the diagonal of @var{A} will contain the +@math{1}-by-@math{1} real eigenvalues and @math{2}-by-@math{2} +complex conjugate eigenvalue systems, and the rest of @var{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 @code{w->n_evals}. +The converged eigenvalues are stored in the beginning of @var{eval}. +@end deftypefun + +@deftypefun int gsl_eigen_nonsymm_Z (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix * @var{Z}, gsl_eigen_nonsymm_workspace * @var{w}) +This function is identical to @code{gsl_eigen_nonsymm} except it also +computes the Schur vectors and stores them into @var{Z}. +@end deftypefun + +@deftypefun {gsl_eigen_nonsymmv_workspace *} gsl_eigen_nonsymmv_alloc (const size_t @var{n}) +This function allocates a workspace for computing eigenvalues and +eigenvectors of @var{n}-by-@var{n} real nonsymmetric matrices. The +size of the workspace is @math{O(5n)}. +@end deftypefun + +@deftypefun void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun int gsl_eigen_nonsymmv (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_nonsymmv_workspace * @var{w}) +This function computes eigenvalues and right eigenvectors of the +@var{n}-by-@var{n} real nonsymmetric matrix @var{A}. It first calls +@code{gsl_eigen_nonsymm} to compute the eigenvalues, Schur form @math{T}, and +Schur vectors. Then it finds eigenvectors of @math{T} and backtransforms +them using the Schur vectors. The Schur vectors are destroyed in the +process, but can be saved by using @code{gsl_eigen_nonsymmv_Z}. The +computed eigenvectors are normalized to have Euclidean norm 1. On +output, the upper portion of @var{A} contains the Schur form +@math{T}. If @code{gsl_eigen_nonsymm} fails, no eigenvectors are +computed, and an error code is returned. +@end deftypefun + +@deftypefun int gsl_eigen_nonsymmv_Z (gsl_matrix * @var{A}, gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_matrix * @var{Z}, gsl_eigen_nonsymmv_workspace * @var{w}) +This function is identical to @code{gsl_eigen_nonsymmv} except it also saves +the Schur vectors into @var{Z}. +@end deftypefun + +@node Sorting Eigenvalues and Eigenvectors +@section Sorting Eigenvalues and Eigenvectors +@cindex sorting eigenvalues and eigenvectors + +@deftypefun int gsl_eigen_symmv_sort (gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_sort_t @var{sort_type}) +This function simultaneously sorts the eigenvalues stored in the vector +@var{eval} and the corresponding real eigenvectors stored in the columns +of the matrix @var{evec} into ascending or descending order according to +the value of the parameter @var{sort_type}, + +@table @code +@item GSL_EIGEN_SORT_VAL_ASC +ascending order in numerical value +@item GSL_EIGEN_SORT_VAL_DESC +descending order in numerical value +@item GSL_EIGEN_SORT_ABS_ASC +ascending order in magnitude +@item GSL_EIGEN_SORT_ABS_DESC +descending order in magnitude +@end table + +@end deftypefun + +@deftypefun int gsl_eigen_hermv_sort (gsl_vector * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type}) +This function simultaneously sorts the eigenvalues stored in the vector +@var{eval} and the corresponding complex eigenvectors stored in the +columns of the matrix @var{evec} into ascending or descending order +according to the value of the parameter @var{sort_type} as shown above. +@end deftypefun + +@deftypefun int gsl_eigen_nonsymmv_sort (gsl_vector_complex * @var{eval}, gsl_matrix_complex * @var{evec}, gsl_eigen_sort_t @var{sort_type}) +This function simultaneously sorts the eigenvalues stored in the vector +@var{eval} and the corresponding complex eigenvectors stored in the +columns of the matrix @var{evec} into ascending or descending order +according to the value of the parameter @var{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. +@end deftypefun + + +@comment @deftypefun int gsl_eigen_jacobi (gsl_matrix * @var{matrix}, gsl_vector * @var{eval}, gsl_matrix * @var{evec}, unsigned int @var{max_rot}, unsigned int * @var{nrot}) +@comment This function finds the eigenvectors and eigenvalues of a real symmetric +@comment matrix by Jacobi iteration. The data in the input matrix is destroyed. +@comment @end deftypefun + +@comment @deftypefun int gsl_la_invert_jacobi (const gsl_matrix * @var{matrix}, gsl_matrix * @var{ainv}, unsigned int @var{max_rot}) +@comment Invert a matrix by Jacobi iteration. +@comment @end deftypefun + +@comment @deftypefun int gsl_eigen_sort (gsl_vector * @var{eval}, gsl_matrix * @var{evec}, gsl_eigen_sort_t @var{sort_type}) +@comment This functions sorts the eigensystem results based on eigenvalues. +@comment Sorts in order of increasing value or increasing +@comment absolute value, depending on the value of +@comment @var{sort_type}, which can be @code{GSL_EIGEN_SORT_VALUE} +@comment or @code{GSL_EIGEN_SORT_ABSVALUE}. +@comment @end deftypefun + +@node Eigenvalue and Eigenvector Examples +@section Examples + +The following program computes the eigenvalues and eigenvectors of the 4-th order Hilbert matrix, @math{H(i,j) = 1/(i + j + 1)}. + +@example +@verbatiminclude examples/eigen.c +@end example + +@noindent +Here is the beginning of the output from the program, + +@example +$ ./a.out +eigenvalue = 9.67023e-05 +eigenvector = +-0.0291933 +0.328712 +-0.791411 +0.514553 +... +@end example + +@noindent +This can be compared with the corresponding output from @sc{gnu octave}, + +@example +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 +@end example + +@noindent +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 @math{V(x;i,j) = x_i^{n - j}} with +@math{x = (-1,-2,3,4)}. + +@example +@verbatiminclude examples/eigen_nonsymm.c +@end example + +@noindent +Here is the beginning of the output from the program, + +@example +$ ./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 +... +@end example + +@noindent +This can be compared with the corresponding output from @sc{gnu octave}, + +@example +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 + +@end example +Note that the eigenvectors corresponding to the eigenvalue +@math{5.54555 + 3.08545i} are slightly different. This is because +they differ by the multiplicative constant +@math{0.9999984 + 0.0017674i} which has magnitude 1. + +@node Eigenvalue and Eigenvector References +@section References and Further Reading + +Further information on the algorithms described in this section can be +found in the following book, + +@itemize @asis +@item +G. H. Golub, C. F. Van Loan, @cite{Matrix Computations} (3rd Ed, 1996), +Johns Hopkins University Press, ISBN 0-8018-5414-8. +@end itemize + +@noindent +The @sc{lapack} library is described in, + +@itemize @asis +@item +@cite{LAPACK Users' Guide} (Third Edition, 1999), Published by SIAM, +ISBN 0-89871-447-8. + +@uref{http://www.netlib.org/lapack} +@end itemize + +@noindent +The @sc{lapack} source code can be found at the website above along with +an online copy of the users guide. |