EIGEN Call :: SAS/IML(R) 13.1 User's Guide

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}$ :

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

$a_{max} = \max _{i,j=1,\ldots ,n} |a_{i,j}|$
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.
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.