bug-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] Re: [Help-gsl] error in Schur decomposition of a non-symmetric


From: Patrick Alken
Subject: [Bug-gsl] Re: [Help-gsl] error in Schur decomposition of a non-symmetric matrix ?
Date: Mon, 22 Mar 2010 22:09:44 -0600
User-agent: Thunderbird 2.0.0.24 (X11/20100317)


Francesco Abbate wrote:
Hi all,

it seems that I've found an error in GSL manual in the section about
eigenvalues/vectors determination of non-symmetric real matrices. The
manual says:

14.3 Real Nonsymmetric Matrices
===============================

The solution of the real nonsymmetric eigensystem problem for a matrix
A involves computing the Schur decomposition

    A = Z T Z^T

  where Z is an orthogonal matrix of Schur vectors and T, the Schur
form, is quasi upper triangular with diagonal 1-by-1 blocks which are
real eigenvalues of A, and diagonal 2-by-2 blocks whose eigenvalues are
complex conjugate eigenvalues of A. The algorithm used is the
double-shift Francis method.
-----------------------

and after about the function gsl_eigen_nonsymmv and gsl_eigen_nonsymmv_Z :

-- 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 T, and
    Schur vectors. Then it finds eigenvectors of T 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 unit magnitude. On
    output, the upper portion of A contains the Schur form T. 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 that it
    also saves the Schur vectors into Z.
------------------------------

The problems seems to be about the Schur decomposition formula. What
seems to be true is that:

A = Z T Z^(-1)

and not that: A = Z T Z^T as asserted in the documentation. I've found
that empirical guesses.

Hi and thanks for your report. I did some quick tests and it seems you're right, Z may not be orthogonal in some cases, and so it does appear there is a bug somewhere. It does seem to be true that A = Z T Z^-1, and the eigen test code tests for: A Z = Z T which would be true in that case. I will have to look into it more to find out what is breaking orthogonality in the computation of Z.
Otherwise I'm perplexed about the statement in the gsl_eigen_nonsymmv:
"On output, the upper portion of A contains the Schur form T". I don't
understand this statement because before it was stated that: "T, the
Schur form, is quasi upper triangular with diagonal 1-by-1 blocks
which are real eigenvalues of A, and diagonal 2-by-2 blocks whose
eigenvalues are complex conjugate eigenvalues of A". So the T matrix
is not really upper triangular.

Yes, T is not strictly upper triangular, and the standard terminology seems to be "quasi upper triangle" since complex eigenvalues will have a non-zero subdiagonal element in T. So T is upper triangle with some possible non-zero subdiagonal elements corresponding to complex eigenvalues.
I've made some tests and it seems that the function gsl_eigen_nonsymmv
set the whole matrix A to the Schur form T and not only the upper part.

I don't believe the lower elements of A are zero'd out and could contain garbage values. I wouldn't rely on these values being 0, but you could call gsl_linalg_hessenberg_set_zero to force them to be 0.
It would be nice if someone could check because what I've seen is
based on empirical tests using the functions.

Best regards,
Francesco


_______________________________________________
Help-gsl mailing list
address@hidden
http://lists.gnu.org/mailman/listinfo/help-gsl





reply via email to

[Prev in Thread] Current Thread [Next in Thread]