From: Francesco Abbate <francesco.bbt@gmail.com>
To: gsl-discuss@sourceware.org
Subject: error in Schur decomposition of a non-symmetric matrix ?
Date: Tue, 23 Mar 2010 15:33:00 -0000 [thread overview]
Message-ID: <4ff79edb1003220246p7ef776b2o78d7f4d7dcc5622d@mail.gmail.com> (raw)
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
it is 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 empirically.
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.
I've made some tests and it seems that the function gsl_eigen_nonsymmv
set the matrix A to the Schur form T and not only the upper part.
It would be nice if someone could check because what I've seen is
based on empirical tests using the functions.
Best regards,
Francesco
reply other threads:[~2010-03-23 15:33 UTC|newest]
Thread overview: [no followups] expand[flat|nested] mbox.gz Atom feed
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=4ff79edb1003220246p7ef776b2o78d7f4d7dcc5622d@mail.gmail.com \
--to=francesco.bbt@gmail.com \
--cc=gsl-discuss@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).