Language Reference

EIGEN Call

CALL EIGEN (evals, evecs, A<VECL=vl>);

The EIGEN subroutine computes eigenvalues and eigenvectors of an arbitrary square numeric matrix.

The A argument is the input argument to the EIGEN subroutine. The EIGEN call returns the following values:

evals

names a matrix to contain the eigenvalues of A.

evecs

names a matrix to contain the right eigenvectors of A.

vl

is an optional $n \times n$ matrix that contains the left eigenvectors of A in the same manner that evecs contains the right eigenvectors.

The EIGEN subroutine computes evals, a matrix that contains the eigenvalues of A. If A is symmetric, evals is the $n \times 1$ vector that contains the n real eigenvalues of A. If A is not symmetric (as determined by the criteria in the symmetry test described later), evals is an $n \times 2$ matrix. The first column of evals contains the real parts, $\mbox{Re}(\lambda )$, and the second column contains the imaginary parts, $\mbox{Im}(\lambda )$. Each row represents one eigenvalue, $\mbox{Re}(\lambda ) + i\mbox{Im}(\lambda )$.

If A is symmetric, the eigenvalues are arranged in descending order. Otherwise, the eigenvalues are sorted first by their real parts, then by the magnitude of their imaginary parts. Complex conjugate eigenvalues, $\mbox{Re}(\lambda )\pm i\mbox{Im}(\lambda )$, are stored in standard order; that is, the eigenvalue of the pair with a positive imaginary part is followed by the eigenvalue of the pair with the negative imaginary part.

The EIGEN subroutine also computes evecs, a matrix that contains the orthonormal column eigenvectors that correspond to evals. If A is symmetric, then the first column of evecs is the eigenvector that corresponds to the largest eigenvalue, and so forth. If A is not symmetric, then evecs is an $n \times n$ matrix that contains the right eigenvectors of A. If the eigenvalue in row i of evals is real, then column i of evecs contains the corresponding real eigenvector. If rows i and $i+1$ of evals contain complex conjugate eigenvalues $\mbox{Re}(\lambda )\pm i\mbox{Im}(\lambda )$, then columns i and $i+1$ of evecs contain the real part, $\mb{u}$, and imaginary part, $\mb{v}$, of the two corresponding eigenvectors $\mathbf{u}\pm i\mathbf{v}$.

The following paragraphs present some properties of eigenvalues and eigenvectors. Let $\mathbf{A}$ be a general $n\times n$ matrix. The eigenvalues of $\mathbf{A}$ are the roots of the characteristic polynomial, which is defined as $p(z) = \det (z\mathbf{I} - \mathbf{A})$. The spectrum, denoted by $\lambda (A)$, is the set of eigenvalues of the matrix A. If $\lambda (\mathbf{A})=\{ \lambda _1,\ldots ,\lambda _ n\} $, then $\det (\mathbf{A}) = \lambda _1 \lambda _2 \cdots \lambda _ n$.

The trace of $\mathbf{A}$ is defined by

\[ \mr{tr}(\bA ) = \sum _{i=1}^ n a_{ii} \]

and tr$(\bA )= \lambda _1 + \ldots + \lambda _ n$.

An eigenvector is a nonzero vector, $\mb{x}$, that satisfies $\mathbf{A}\mathbf{x} = \lambda \mathbf{x}$ for $\lambda \in \lambda (\mathbf{A})$. Right eigenvectors satisfy $\mathbf{A}\mathbf{x} = \lambda \mathbf{x}$, and left eigenvectors satisfy $\mathbf{x}^{H}\mathbf{A} = \lambda \mathbf{x}^{H}$, where $\mathbf{x}^{H}$ is the complex conjugate transpose of $\mathbf{x}$. Taking the conjugate transpose of both sides shows that left eigenvectors also satisfy $\mathbf{A}^{\prime }\mathbf{x} = \bar{\lambda }\mathbf{x}$.

The following are properties of the unsymmetric real eigenvalue problem, in which the real matrix $\mathbf{A}$ is square but not necessarily symmetric:

  • The eigenvalues of an unsymmetric matrix $\mathbf{A}$ can be complex. If $\mathbf{A}$ has a complex eigenvalue, $\mbox{Re}(\lambda )+i\mbox{Im}(\lambda )$, then the conjugate complex value $\mbox{Re}(\lambda )-i\mbox{Im}(\lambda )$ is also an eigenvalue of $\mathbf{A}$.

  • The right and left eigenvectors that correspond to a real eigenvalue of $\mathbf{A}$ are real. The right and left eigenvectors that correspond to conjugate complex eigenvalues of $\mathbf{A}$ are also conjugate complex.

  • The left eigenvectors of $\mathbf{A}$ are the same as the complex conjugate right eigenvectors of $\mathbf{A}^{\prime }$.

The three routines, EIGEN, EIGVAL, and EIGVEC, use the following test of symmetry for a square argument matrix $\mathbf{A}$:

  1. Select the entry of $\mathbf{A}$ with the largest magnitude:

    \[ a_{max} = \max _{i,j=1,\ldots ,n} |a_{i,j}| \]

  2. Multiply the value of $a_{max}$ by the square root of the machine precision, $\epsilon $. The value of $\epsilon $ is the largest value stored in double precision that, when added to 1 in double precision, still results in 1.

  3. The matrix $\mathbf{A}$ is considered unsymmetric if there exists at least one pair of symmetric entries that differs in more than $a_{max}\sqrt {\epsilon }$:

    \[ |a_{i,j}-a_{j,i}| > a_{max} \sqrt {\epsilon } \]

If $\mathbf{A}$ is a symmetric matrix and $\mb{M}$ and $\mb{E}$ are the eigenvalues and eigenvectors, respectively, of $\mb{A}$, then the matrices have the following properties:

\begin{eqnarray*} \mb{A}*\mb{E} & = & \mb{E}*\mbox{diag}(\mb{M}) \\ \mb{E}^{\prime }*\mb{E} & = & \mb{I} \end{eqnarray*}

These properties imply the following:

\[ \mb{E}^{\prime } = \mbox{inv}(\mb{E}) \]
\[ \mb{A} = \mb{E}*\mbox{diag}(\mb{M})*\mb{E}^{\prime } \]

The QL method is used to compute the eigenvalues (Wilkinson and Reinsch, 1971).

In statistical applications, nonsymmetric matrices for which eigenvalues are desired are usually of the form $\mb{E}^{-1} \mb{H}$, where $\mb{E}$ and $\mb{H}$ are symmetric. The eigenvalues $\mb{L}$ and eigenvectors $\mb{V}$ of $\mb{E}^{-1}\mb{H}$ can be obtained by using the GENEIG subroutine , or by using the following statements: 

F = root(einv);
A = F*H*F`;
call eigen(L, W, A);
V = F`*W;

The computation can be checked by forming the residuals, r, as shown in the following statement: 

r = einv*H*V - V*diag(L);

The values in r should be of the order of rounding error.

The following statements compute the eigenvalues and left and right eigenvectors of a nonsymmetric matrix with four real and four complex eigenvalues: 

A = {-1  2  0       0       0       0       0  0,
     -2 -1  0       0       0       0       0  0,
      0  0  0.2379  0.5145  0.1201  0.1275  0  0,
      0  0  0.1943  0.4954  0.1230  0.1873  0  0,
      0  0  0.1827  0.4955  0.1350  0.1868  0  0,
      0  0  0.1084  0.4218  0.1045  0.3653  0  0,
      0  0  0       0       0       0       2  2,
      0  0  0       0       0       0      -2  0 };
call eigen(val, rvec, A) vecl="lvec";
print val;

The sorted eigenvalues of the A matrix are shown in Figure 24.117.

Figure 24.117: Complex Eigenvalues of a Nonsymmetric Matrix

val
1 1.7320508
1 -1.732051
1 0
0.2087788 0
0.0222025 0
0.0026187 0
-1 2
-1 -2 




            

            


            
You can verify the correctness of the left and right eigenvector computation by using the following statements: 


/* verify that the right eigenvectors are correct */
vec = rvec;
do j = 1 to ncol(vec);
  /* if eigenvalue is real */
  if val[j,2] = 0. then do;
    v = A * vec[,j] - val[j,1] * vec[,j];
    if any( abs(v) > 1e-12 ) then
      badVectors = badVectors || j;
    end;
  /* if eigenvalue is complex with positive imaginary part */
  else if val[j,2] > 0. then do;
    /* the real part */
    rp = val[j,1] * vec[,j] - val[j,2] * vec[,j+1];
    v = A * vec[,j] - rp;
    /* the imaginary part */
    ip = val[j,1] * vec[,j+1] + val[j,2] * vec[,j];
    u = A * vec[,j+1] - ip;
    if any( abs(u) > 1e-12 ) | any( abs(v) > 1e-12 ) then
      badVectors = badVectors || j || j+1;
    end;
  end;
if ncol( badVectors ) > 0 then
  print "Incorrect right eigenvectors:" badVectors;
else print "All right eigenvectors are correct";



Similar statements can be written to verify the left eigenvectors. The statements use the fact that the left eigenvectors
               of $\mathbf{A}$ are the same as the complex conjugate right eigenvectors of $\mathbf{A}^{\prime }$: 
            


/* verify that the left eigenvectors are correct */
vec = lvec;
do j = 1 to ncol(vec);
  /* if eigenvalue is real */
  if val[j,2] = 0. then do;
    v = A` * vec[,j] - val[j,1] * vec[,j];
    if any( abs(v) > 1e-12 ) then
      badVectors = badVectors || j;
    end;
  /* if eigenvalue is complex with positive imaginary part */
  else if val[j,2] > 0. then do;
    /* Note the use of complex conjugation */
    /* the real part */
    rp = val[j,1] * vec[,j] + val[j,2] * vec[,j+1];
    v = A` * vec[,j] - rp;
    /* the imaginary part */
    ip = val[j,1] * vec[,j+1] - val[j,2] * vec[,j];
    u = A` * vec[,j+1] - ip;
    if any( abs(u) > 1e-12 ) | any( abs(v) > 1e-12 ) then
      badVectors = badVectors || j || j+1;
    end;
  end;
if ncol( badVectors ) > 0 then
  print "Incorrect left eigenvectors:" badVectors;
else print "All left eigenvectors are correct";



The EIGEN call performs most of its computations in the memory allocated for returning the eigenvectors.