Language Reference

EIGEN Call

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

This subroutine is supported by the IML procedure and the iml action.

The EIGEN subroutine computes eigenvalues and eigenvectors of an arbitrary square numeric matrix. The EIGEN subroutine will use vendor-supplied eigenvalue routines if they are available on your system. (An example is the Intel Math Kernel Library (MKL), which is tuned to provide optimal performance for a given Intel processor.) Because eigenvectors are not unique, the results of eigenvector computations that use vendor-supplied routines are not necessarily identical to the results from earlier releases. Use the RESET EIGEN93 statement to prevent SAS/IML from using vendor-supplied routines.

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, Re left-parenthesis lamda right-parenthesis, and the second column contains the imaginary parts, Im left-parenthesis lamda right-parenthesis. Each row represents one eigenvalue, Re left-parenthesis lamda right-parenthesis plus i Im left-parenthesis lamda right-parenthesis.

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, Re left-parenthesis lamda right-parenthesis plus-or-minus i Im left-parenthesis lamda right-parenthesis, 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 plus 1 of evals contain complex conjugate eigenvalues Re left-parenthesis lamda right-parenthesis plus-or-minus i Im left-parenthesis lamda right-parenthesis, then columns i and i plus 1 of evecs contain the real part, bold u, and imaginary part, bold v, of the two corresponding eigenvectors bold u plus-or-minus i bold v.

The following paragraphs present some properties of eigenvalues and eigenvectors. Let bold upper A be a general n times n matrix. The eigenvalues of bold upper A are the roots of the characteristic polynomial, which is defined as p left-parenthesis z right-parenthesis equals det left-parenthesis z bold upper I minus bold upper A right-parenthesis. The spectrum, denoted by lamda left-parenthesis upper A right-parenthesis, is the set of eigenvalues of the matrix A. If lamda left-parenthesis bold upper A right-parenthesis equals StartSet lamda 1 comma ellipsis comma lamda Subscript n Baseline EndSet, then det left-parenthesis bold upper A right-parenthesis equals lamda 1 lamda 2 midline-horizontal-ellipsis lamda Subscript n.

The trace of bold upper A is defined by

normal t normal r left-parenthesis bold upper A right-parenthesis equals sigma-summation Underscript i equals 1 Overscript n Endscripts a Subscript i i

and trleft-parenthesis bold upper A right-parenthesis equals lamda 1 plus ellipsis plus lamda Subscript n.

An eigenvector is a nonzero vector, bold x, that satisfies bold upper A bold x equals lamda bold x for lamda element-of lamda left-parenthesis bold upper A right-parenthesis. Right eigenvectors satisfy bold upper A bold x equals lamda bold x, and left eigenvectors satisfy bold x Superscript upper H Baseline bold upper A equals lamda bold x Superscript upper H, where bold x Superscript upper H is the complex conjugate transpose of bold x. Taking the conjugate transpose of both sides shows that left eigenvectors also satisfy bold upper A prime bold x equals lamda overbar bold x.

The following are properties of the unsymmetric real eigenvalue problem, in which the real matrix bold upper A is square but not necessarily symmetric:

  • The eigenvalues of an unsymmetric matrix bold upper A can be complex. If bold upper A has a complex eigenvalue, Re left-parenthesis lamda right-parenthesis plus i Im left-parenthesis lamda right-parenthesis, then the conjugate complex value Re left-parenthesis lamda right-parenthesis minus i Im left-parenthesis lamda right-parenthesis is also an eigenvalue of bold upper A.

  • The right and left eigenvectors that correspond to a real eigenvalue of bold upper A are real. The right and left eigenvectors that correspond to conjugate complex eigenvalues of bold upper A are also conjugate complex.

  • The left eigenvectors of bold upper A are the same as the complex conjugate right eigenvectors of bold upper A prime.

The three routines, EIGEN, EIGVAL, and EIGVEC, use the following test of symmetry for a square argument matrix bold upper A:

  1. Select the entry of bold upper A with the largest magnitude:

    a Subscript max Baseline equals max Underscript i comma j equals 1 comma ellipsis comma n Endscripts StartAbsoluteValue a Subscript i comma j Baseline EndAbsoluteValue
  2. Multiply the value of a Subscript normal m normal a normal x 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 bold upper A is considered unsymmetric if there exists at least one pair of symmetric entries that differs in more than a Subscript normal m normal a normal x Baseline StartRoot epsilon EndRoot:

    StartAbsoluteValue a Subscript i comma j Baseline minus a Subscript j comma i Baseline EndAbsoluteValue greater-than a Subscript normal m normal a normal x Baseline StartRoot epsilon EndRoot

If bold upper A is a symmetric matrix and bold upper M and bold upper E are the eigenvalues and eigenvectors, respectively, of bold upper A, then the matrices have the following properties:

StartLayout 1st Row 1st Column bold upper A asterisk bold upper E 2nd Column equals 3rd Column bold upper E asterisk diag left-parenthesis bold upper M right-parenthesis 2nd Row 1st Column bold upper E prime asterisk bold upper E 2nd Column equals 3rd Column bold upper I EndLayout

These properties imply the following:

bold upper E prime equals inv left-parenthesis bold upper E right-parenthesis
bold upper A equals bold upper E asterisk diag left-parenthesis bold upper M right-parenthesis asterisk bold upper 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 bold upper E Superscript negative 1 Baseline bold upper H, where bold upper E and bold upper H are symmetric. The eigenvalues bold upper L and eigenvectors bold upper V of bold upper E Superscript negative 1 Baseline bold upper 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 137.

Figure 137: Complex Eigenvalues of a Nonsymmetric Matrix

val
10
11.7320508
1-1.732051
0.20877880
0.02220250
0.00261870
-12
-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 bold upper A are the same as the complex conjugate right eigenvectors of bold upper 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.

Last updated: July 20, 2026