Eigenvalue extraction

The application of the subspace iteration method to eigenproblems arising from finite element models of the dynamic behavior of structures has been discussed by a number of authors—see Ramaswami (1979) for references. The basic idea is a simultaneous inverse power iteration. A small set of base vectors is created, thus defining a “subspace”: this “subspace” is then transformed, by iteration, into the space containing the lowest few eigenvectors of the overall system.

Assume that, at the i-th iteration, a set of m vectors, ${\left\{V\right\}}_{\left(i\right)}^n$, $n=1,\mathrm{\dots },m$, exists, where m is less than the number of variables in the finite element model. These are considered as the “base vectors” that define the m-dimensional subspace out of the n dimensions defined by the variables in the finite element model. We arrange these vectors as the columns in the matrix ${\left[V\right]}_{\left(i\right)}$. The number of rows of ${\left[V\right]}_{\left(i\right)}$ is, thus, the size of the complete set of equations (the number of variables in the finite element model, N), and the number of columns is the dimension of the subspace, m.

The first step in the algorithm is to define a new set of base vectors by solving

This operation, a generalized “inverse power sweep” with m vectors, involves the solution of the complete set of linear stiffness equations for several right-hand-side vectors (with, after the first iteration of the method, a previously decomposed matrix).

The stiffness and mass matrices of the structure are then projected onto the subspace by

thus defining a mass matrix, $\left[M^{*}\right]$, and a stiffness matrix, $\left[K^{*}\right]$, in the subspace. These matrices are of dimension m by m. The eigenproblem

is now solved completely in the subspace using the Householder and Q-R methods, which are discussed below.

The eigenvectors have now been defined in the reduced space and can be transformed back to the full space of the structural problem to define $\left[V\right]$ for the next iteration:

where $\left[\varphi \right]=\left[{\left\{\varphi \right\}}_1\right|{\left\{\varphi \right\}}_2\left|\mathrm{\cdots }\right|{\left\{\varphi \right\}}_m]$. This completes an iteration.

With $m=1$ (a one-dimensional subspace), the method reduces to the simple “inverse power sweep” method (see Wilkinson, 1965, or Strang, 1976).

The advantage of the subspace method is the extraction of the eigenvalues in reduced space, which will cause a rapid convergence to the eigenvectors in full space. The number of base vectors carried in the iterations and the choice of initial base vectors are, therefore, important for an economical solution. The convergence rate of a particular eigenvalue is proportional to ${\omega }_i^2/{\omega }_{\left(m+1\right)}^2$, where ${\omega }_{\left(m+1\right)}^2$ is the eigenvalue corresponding to the next vector, $m+1$ beyond the m vectors used to define the subspace. The default value of m used in Abaqus is $m=\min \left(2p,p+8\right)$, where p is the number of eigenvectors requested. If more vectors are used, the number of required iterations is reduced, but each iteration takes longer because of the greater number of right-hand sides. Increasing m can sometimes improve the performance of the algorithm significantly.

The choice of starting vectors is also important. There is no requirement that these should be close to the eigenvectors for rapid convergence. The Householder and Q-R steps rotate the vectors into the eigenvectors as long as the base vectors span the space of the eigenvectors. The approach used is to choose initial vectors that span this space as completely as possible. The algorithm used in Abaqus/Standard for this purpose is that of Bathe and Wilson (1972). They recommend using the diagonal mass terms as one vector: the other vectors are unit vectors, providing a set of single unit entries at the nodes and the degrees of freedom with the largest mass terms. A check is made to ensure that all degrees of freedom are represented in these additional vectors.

This completes the description of the basic subspace iteration algorithm as it is implemented in Abaqus/Standard. The small eigenproblem, for which all eigenpairs must be found, is solved by the Householder and Q-R algorithms. These algorithms are described at the end of this section.

The implementation of the Lanczos eigensolver as a powerful tool for extraction of the extreme eigenvalues and the corresponding eigenvectors of a sparse symmetric generalized eigenproblem has been discussed by a number of authors; see, for example, Parlett (1980), Parlett and Nour-Omid (1989), Simon (1984), and Ericsson and Ruhe (1980). A shifted block Lanczos algorithm was developed and described in detail by Grimes, Lewis, and Simon (1994).

The Lanczos procedure in Abaqus/Standard consists of a set of Lanczos “runs,” in each of which a set of iterations called steps is performed. For each Lanczos run the following spectral transformation is applied:

where $\sigma$ is the shift, $\theta$ is the eigenvalue, and $\left\{\varphi \right\}$ is the eigenvector. This transformation allows rapid convergence to the desired eigenvalues. The eigenvectors of the symmetrized problem (Equation 2) and the transformed problem (Equation 5) are identical, while the eigenvalues of the original problem and the transformed problem are related in the following manner:

A Lanczos run will be terminated when its continuation is estimated to be inefficient. In general, only tens of eigenvalues (closest to the shift value) are computed in a single Lanczos run. The possibility of computing many eigenmodes by carrying out several runs with different shift values is an important feature of the Lanczos eigensolver.

Within each run a sequence of Krylov subspaces is created, and the best possible approximation of the eigenvectors on each subspace is computed in a series of “steps.” In each Lanczos step the dimension of the subspace grows, allowing better approximation of the desired eigenvectors. This is in contrast to the subspace iteration method, in which the dimension of the subspace used to approximate the eigenvectors is fixed.

In theory the basic Lanczos process (in the assumption of “exact” computations without taking into account round-off errors) is able to determine only simple eigenvalues. The shifting strategy (and the Sturm sequence check as a part of it) detects missing modes and enforces computation of all the modes during the subsequent Lanczos runs. However, this strategy is expensive if the multiplicity of certain eigenvalues is high. Therefore, a “blocked” version of the Lanczos algorithm is implemented in Abaqus/Standard. The idea is to start with a block of orthogonal vectors and to increase the dimension of the Krylov subspaces by the block size at each Lanczos step. This approach allows automatic computation of all multiple eigenvalues if the largest multiplicity does not exceed the block size. Another important advantage of the blocked Lanczos method is that it allows efficient implementation of expensive computational kernels such as matrix-blocked vector multiplications, blocked back substitutions, and blocked vector products.

As discussed above, the Lanczos process consists of several Lanczos runs. Each Lanczos run is associated with some shift value that remains constant during the run. The initial shift value, ${\sigma }_1$, is determined (by a heuristic approach) using the geometric mean of the centers of the Gershgorin circles, termed the “problem scale.” At the beginning of each Lanczos run a factorization of the shifted matrix $\left[K\right]-{\sigma }_p\left[M\right]$ (where p is the run number) is carried out, after which a sequence of Lanczos steps is performed.

Each block Lanczos step i is implemented in Abaqus/Standard in the following manner:

1. Solve the system of linear equations with b right-hand sides (b is the Lanczos block size) using the factorized shifted matrix:

   $$\left(\left[K\right]-{\sigma }_p\left[M\right]\right){\left[\widehat{V}\right]}_{\left(i+1\right)}=\left[M\right]{\left[V\right]}_{\left(i\right)},$$

   where ${\left[V\right]}_{\left(i\right)}$ is a block of Lanczos vectors. The number of rows of ${\left[V\right]}_{\left(i\right)}$ is the number of variables in the finite element model, and the number of columns is the Lanczos block size b. The initial block of vectors ${\left[V\right]}_{\left(1\right)}$ is a set of random vectors orthonormalized using the procedure: $\left[I\right]={\left[V\right]}_{\left(1\right)}^T\left[M\right]{\left[V\right]}_{\left(1\right)},$ where $\left[I\right]$ is the identity matrix.

2. Set the auxiliary block of vectors:

   $${\left[U\right]}_{\left(i+1\right)}={\left[\widehat{V}\right]}_{\left(i+1\right)}-{\left[V\right]}_{\left(i-1\right)}{\beta }_i,$$

   where ${\beta }_i$ is a $b\times b$ upper triangular matrix. ${\beta }_1=\left[0\right]$.

3. Compute the symmetric matrix ${\alpha }_i$ of size b:

   $${\alpha }_i={\left[U\right]}_{\left(i+1\right)}^T\left[M\right]{\left[V\right]}_{\left(i\right)}.$$

4. Formulate the Lanczos residual:

   $${\left[R\right]}_{\left(i+1\right)}={\left[U\right]}_{\left(i+1\right)}-{\left[V\right]}_{\left(i\right)}{\alpha }_i.$$

5. Normalize the residual. Compute a block of vectors, ${\left[V\right]}_{\left(i+1\right)}$, such that

   $${\left[V\right]}_{\left(i+1\right)}{\beta }_{i+1}={\left[R\right]}_{\left(i+1\right)}$$

   and

   $$\left[I\right]={\left[V\right]}_{\left(i+1\right)}^T\left[M\right]{\left[V\right]}_{\left(i+1\right)},$$

   where ${\beta }_{i+1}$ is a $b\times b$ upper triangular matrix.

6. Estimate the loss of orthogonality between ${\left[V\right]}_{\left(i+1\right)}$ and ${\left[V\right]}_{\left(j\right)},$ for $j=1,2,\mathrm{\dots },i-1$, using the partial reorthogonalization technique; and perform reorthogonalization if necessary (see Simon, 1984; Grimes, Lewis, and Simon, 1994).

7. Perform “local reorthogonalization”: recompute ${\left[V\right]}_{\left(i+1\right)}$ to provide

   $$\left[0\right]={\left[V\right]}_{\left(i+1\right)}^T\left[M\right]{\left[V\right]}_{\left(i\right)}.$$

8. Solve the following reduced eigenvalue problem for the band matrix $\left[T\right]$:

   $$\left[T\right]\left[S\right]=\left[S\right]\left[\theta \right],$$

   where

   $$\left[T\right]=\left[\begin{array}{cccccc}\hfill {\alpha }_1\hfill & \hfill {\beta }_2^T\hfill & & & & \\ \hfill {\beta }_2\hfill & \hfill {\alpha }_2\hfill & \hfill {\beta }_3^T\hfill & & & \\ & \hfill {\beta }_3\hfill & \hfill {\alpha }_3\hfill & \hfill {\beta }_4^T\hfill & & \\ & & \hfill {\beta }_4\hfill & \hfill {\alpha }_4\hfill & \hfill \mathrm{\dots }\hfill & \\ & & & \hfill \mathrm{\dots }\hfill & \hfill \mathrm{\dots }\hfill & \hfill {\beta }_i^T\hfill \\ & & & & \hfill {\beta }_i\hfill & \hfill {\alpha }_i\hfill \end{array}\right]$$

   and $\left[S\right]$ and $\left[\theta \right]$ are, respectively, the matrices containing the eigenvectors and eigenvalues of the reduced eigenproblem. This problem is solved by the Householder and Q-R algorithms, which are discussed below.

9. Determine the error bounds in eigenvalue approximation for the symmetrized eigenvalue problem (Equation 2) (see Parlett, 1980; Grimes, Lewis, and Simon, 1994). Check the termination conditions of the Lanczos run.

The Lanczos run terminates if one of the following conditions is satisfied:

• All the eigenvalues required for the current run are extracted.

• The triangular matrix ${\beta }_{i+1}$ is singular or ill-conditioned.

• The number of Lanczos steps exceeds the maximum number allowed.

• Continuation of the current run is evaluated to be inefficient. This decision is based on the estimation of the “cost per eigenvalue” over the next few steps (the Lanczos run is continued as long as the cost per eigenvalue is decreasing).

After termination of each Lanczos run, the converged eigenvectors of the symmetrized problem (Equation 2) are recovered using the blocks of vectors ${\left[V\right]}_{\left(j\right)},j=1,2,\mathrm{\dots },i$ and the eigenvectors $\left[S\right]$.

Once the Lanczos run for the shift ${\sigma }_p$ is completed, the shift value ${\sigma }_{p+1}$ is computed using the results of all the previous Lanczos runs. To describe the shifting strategy, the following concepts are introduced:

The “computational interval” is the interval between the minimum and maximum eigenvalues of interest. This interval can be finite (both ends are finite), semi-infinite (only one end is finite), or infinite (both ends are infinite). The “center point” is the point in the computational interval nearest to the desired eigenvalues.

The Sturm sequence number, $\nu ([A])$, of a real nonsingular symmetric matrix $\left[A\right]$ is the number of negative eigenvalues. This number is equal to the number of negative terms in the diagonal matrix $\left[D\right]$ of the Cholesky decomposition $\left[A\right]=\left[L\right]\left[D\right]{\left[L\right]}^T$ and is, therefore, available after the Cholesky decomposition is completed.

Let ${\sigma }_1$ and ${\sigma }_2$ be two shift values such that $\left({\sigma }_1<{\sigma }_2\right)$. Then the number of eigenvalues of the symmetrized problem (Equation 2, assuming that $\left[M\right]$ is positive definite or positive semidefinite) in the interval $\left[{\sigma }_1,{\sigma }_2\right]$ is equal to $\nu \left(\left[K\right]-{\sigma }_2\left[M\right]\right)-\nu \left(\left[K\right]-{\sigma }_1\left[M\right]\right)$. If this number is equal to the number of eigenvalues actually determined by the Lanczos algorithm, the interval $\left[{\sigma }_1,{\sigma }_2\right]$ is called “a trust interval.” The trust interval containing the center point is referred to as “a primary trust interval” (denoted by the key “+” in the trust intervals list printed in the message file). The shifting strategy is aimed at constructing the primary trust interval inside the computational interval containing the required number of eigenvalues closest to the center point.

One of the most important features in the formulation of the shift update strategy is the concept of a “sentinel.” The sentinels are the endpoints of the intervals containing exclusively converged eigenvalues closest to the current shift value (in each direction). The sentinels are computed during each Lanczos run and are updated at the end of each step after the eigenvalue analysis of the reduced matrix $\left[T\right]$ is completed. The basic assumption is that there are no missing eigenvalues inside the sentinels; therefore, they are excluded from the computation intervals for the upcoming Lanczos runs. This assumption is later verified on the basis of the Sturm sequence check, and a special procedure is activated if some eigenvalues are missing. If no missing eigenvalues are detected, the sentinels transform directly into the corresponding trust intervals.

The new shift values are selected on the basis of the nonconverged eigenvalue approximations after the Lanczos run is terminated. Assuming the same convergence properties for the upcoming runs, the new shift values are selected in such a way that the number of eigenvalues expected to be found in the upcoming runs will be close to the number of eigenvalues found inside the corresponding sentinels on the previous run.

The Householder method is used to reduce a general matrix to a symmetric tridiagonal form. A tridiagonal matrix is one whose only nonzero entries are on or immediately adjacent to the diagonal. The first step is to transform Equation 4 to the form

where $\left[I\right]$ is the identity matrix.

This is done by using a Cholesky decomposition of $\left[M^{*}\right]$ and then premultiplying and postmultiplying $\left[K^{*}\right]$ by the inverse of the lower and upper triangular matrices that are, thus, obtained:

where

To preserve symmetry in $\left[A\right]$, it is necessary that the decomposition of $\left[M^{*}\right]$ result in two matrices that are the transpose of each other. A Cholesky decomposition produces the desired result but adds the requirement that the matrix $\left[M^{*}\right]$ be positive definite. $\left[M^{*}\right]$ should always be positive definite (because $\left[M\right]$ is positive definite in all of the problems considered here) if the base vectors defining the subspace are not linearly dependent across $\left[M\right]$. In practice the choices made for starting values of the base vectors usually satisfy this requirement, although cases can arise when this is not true. Then Abaqus/Standard will reduce the dimensionality of the subspace to obtain a positive definite $\left[M^{*}\right]$.

The Householder transformation starts with the matrix $\left[A\right]$ and, proceeding one column at a time, reduces all the entries outside the tridiagonal part of the matrix to zero by the transformation

The matrix $\left[U\right]$ is of the form

For the first transformation $\left[H\right]$ is of the same order as $\left[A\right]$ and $\left[U\right]$; that is, $m\times m$, the dimension of the subspace. However, each transformation reduces the order of $\left[H\right]$ by one, the leading parts of $\left[U\right]$ being 1 on the diagonal and 0 outside the diagonal, as illustrated above.

The matrix $\left[H\right]$ is obtained by the following algorithm:

1. If this is the kth iteration, write the kth column of ${\left[A\right]}_k$ below the diagonal as

   $$\left\{X\right\}={\lfloor a_{k+1,k},a_{k+2,k}\mathrm{\dots }{a}_{n,k}\rfloor }^T.$$

2. Calculate the norm $\left|X\right|=\sqrt{\lfloor X\rfloor \left\{X\right\}}.$

3. Define a vector,

   $$\lfloor V\rfloor =\lfloor a_{k+1,k}\rfloor +\left|X\right|{\lfloor a_{k+1,k}\mathrm{\dots }{a}_{n,k}\rfloor }^T.$$

4. Then,

   $$\left[H\right]=\left[I\right]-2\frac{1}{\lfloor V\rfloor \left\{V\right\}}\left\{V\right\}\lfloor V\rfloor .$$

This algorithm is developed in detail in Strang's (1976) book.

The Householder algorithm produces a symmetric tridiagonal matrix, which has the same eigenvalues as the original matrix, because the transformation (Equation 6) does not alter the eigenvalues.

The next step is to calculate the eigenvalues of the tridiagonal matrix. This is done by the Q-R algorithm. In this method the matrix $\left[A\right]$, which is now tridiagonal (although this is not a requirement for the method to work), is factored into $\left[Q\right]\left[R\right]$, where $\left[Q\right]$ is an orthonormal matrix (that is, ${\left[Q\right]}^T\left[Q\right]=\left[I\right]$) and $\left[R\right]$ is an upper-triangular matrix (that is, all the terms in $\left[R\right]$ below the diagonal are zero).

The next matrix in the iteration is obtained by premultiplying and postmultiplying $\left[A\right]$ by $\left[Q\right]$:

Because $\left[Q\right]$ is orthonormal, this will not affect the eigenvalues. This process gradually reduces the off-diagonal terms of $\left[A\right]$ so that the diagonal terms approach the eigenvalues. To speed up convergence, a method of shifting is introduced, leading to the following iterative loop:

The shift cannot be used until the iteration is converging toward an eigenvalue; but once that is obvious, the shift value ${\alpha }_i$ can be set to the expected eigenvalue. This will lead to a significant improvement in the convergence rate. If the shift is done too early, with a value that is not close to an eigenvalue, it is quite possible that the process will converge to an incorrect number. The Q-R process will converge to the eigenvalues in ascending order; and as soon as an eigenvalue is obtained, the order of the matrix can be reduced by one.

The final step after the eigenvalues have been obtained is to calculate the eigenvectors by using the inverse power method to solve for an eigenvector, given any right-hand side:

where $\lambda$ is the eigenvalue just obtained and $\left\{V\right\}$ is any vector. Because the left-hand-side matrix is singular in the direction of the eigenvector $\left\{\varphi \right\}$, this vector will be obtained regardless of the right-hand-side vector $\left\{V\right\}$, as long as $\left\{V\right\}$ is not orthogonal to the eigenvector. To ensure that consecutive vectors are orthogonal, especially in the case of multiple eigenvalues, $\left\{V\right\}$ is always chosen to be orthogonal to the previously extracted eigenvectors. Since $\left[\left[K^{*}\right]-\lambda \left[M^{*}\right]\right]$ is singular, a slight numerical shift must be included to decompose it and, thus, solve for $\left\{\varphi \right\}$. In the standard notation used in this guide, the matrix of eigenvectors $\left[\varphi \right]$ is written as ${\varphi }_{\alpha }^N$, where N refers to the nodal variable in the problem and $\alpha =1,2,\mathrm{\dots },p$ is the mode number.

In the subspace projection method the original eigensystem (Equation 1) is projected onto a subspace spanned by the eigenvectors of the undamped, symmetric system (Equation 2). Thus, the symmetrized eigenproblem must be solved prior to the complex eigenvalue extraction procedure to create the subspace onto which the original system will be projected. Next, the original mass, damping, and stiffness matrices are projected onto the subspace of N eigenvectors:

Then, we arrive at the following eigenproblem:

Typically, the number of eigenvectors is relatively small; a few hundred is common. This small complex eigenvalue system is solved using the standard QZ method for generalized nonsymmetric eigenproblems. Complex eigenvalues of the projected system are the approximation of the eigenvalues of the original system (Equation 1), but the right eigenvectors of the original system need to be recovered:

where ${\left\{\varphi \right\}}^k$ is the approximation of k-th right eigenvector of the original system.

The only energy density calculated in the eigenvalue extraction procedures is the elastic strain energy density, which for the complex eigenvalue extraction is defined as follows:

where ${\mathit{\sigma }}_r,{\mathit{\sigma }}_i$ are real and imaginary components of the eigenstress tensor and ${\mathit{\epsilon }}_r,{\mathit{\epsilon }}_i$ are real and imaginary components of the eigenstrain tensor.