Arnoldi Procedure in GEMV Form

To begin our discussion, we shall look at the Arnoldi procedure from a slightly different perspective. The Arnoldi procedure was derived previously in §7.5 and has been presented in modified Gram-Schmidt algorithmic form in Algorithm 7.3. In this section, we describe it in matrix-vector (or GEMV) form using the classical Gram-Schmidt process together with the orthogonalization refinement [96]. The Arnoldi relation between the matrix $A$, the basis matrix $V_k$, and the residual vector $f_k$ is of the form

$$A V_k = V_k H_k + f_k e_k^{\ast},$$

where $V_k \in {\cal C}^{n \times k}$ has orthonormal columns, $V_k^{\ast} f_k = 0$, and $H_k \in {\cal C}^{k\times k}$ is upper Hessenberg with nonnegative subdiagonal elements. We shall call this a k-step Arnoldi factorization of $A$. It is easily seen from the construction that $H_k = V_k^{\ast} A V_k$ is upper Hessenberg. When $A$ is Hermitian, this implies that $H_k$ is real symmetric and tridiagonal. The columns of $V_k$ are referred to as the Arnoldi vectors.

A template for computing a $k$-step Arnoldi factorization is given in Algorithm 7.6.

We will now describe some implementation details, referring to the respective phases in Algorithm 7.6.

(1)

Choose the initial starting vector $v$ and normalize. Ideally, for eigenvalue calculations, one should attempt to construct a $v$ that is dominant in eigenvector directions of interest. In the absence of any other considerations, a random vector is a reasonable choice.

(2)

Initialize the Arnoldi procedure. In the subsequent step, $f_1$ will become the new Arnoldi basis vector. This step should include a reorthogonalization just as all of the others do.

(7)

Normalize $f_j$ to get the new basis vector $v_{j+1}$, partially update $H_j$ to get the $(j+1) \times j$ matrix $\hat{H}_{j}$, and compute the new information $w \leftarrow Av_{j+1}$. The norm calculation $\beta_j = \Vert f_j \Vert$ need not be recomputed here. It has already been computed at step (12) and can be saved for reuse here.

(10)

This is the classical Gram-Schmidt step to orthogonalize $w$ with respect to the columns of $V_{j+1}$. This formulation allows level 2 matrix-vector operations. A more detailed discussion is given below.

(12)

The sequence of statements in this if-clause assure that the new direction $f_{j+1}$ is numerically orthogonal to the previously computed directions, i.e., to the columns of $V_{j+1}$. The parameter $\eta$ must be specified ( $0 \le \eta < \infty$). The test asks the question, ``Is $w = Av_{j+1}$ nearly in the space that already has been constructed?" The ratio $\frac{\Vert f_{j+1}\Vert}{ \Vert h \Vert} =\frac{\sin(\theta)}{\cos(\theta)}$, where $\theta$ is the angle that the vector $w$ makes with $\mbox{span}(V_{j+1})$ (i.e., the angle $\theta$ between $w$ and its projection $V_{j+1} h$). A larger value of $\eta$ is a more stringent orthogonality requirement. With $\eta = 0$ the algorithm reverts to classical Gram-Schmidt (CGS). An additional detail is needed to completely specify this process. If, after updating $f_{j+1}$ and $h$, the relation $\Vert f_{j+1} \Vert < \eta \Vert h \Vert$ still holds, then numerically we have $w \in \mbox{span}(V_{j+1})$. In this case, the Arnoldi procedure should either be terminated or a randomly generated vector should be orthogonalized with respect to the columns of $V_{j+1}$ to replace $f_{j+1}$, and $\beta_{j+1}$ should be set to zero to continue the factorization. The value of $\eta$ used in the ARPACK implementation is $1/\sqrt{2}$ (see §7.6.9). Additional discussion of this orthogonalization device is given below.

(17)

$H_{j+1}$ is now a $(j+1) \times (j+1)$ upper Hessenberg matrix.

The dense matrix-vector products $V_{j+1}^{\ast} w$ and $V_{j+1} h$ may be expressed with the Level 2 BLAS operation GEMV. This provides a significant performance advantage on virtually every platform from workstation to supercomputer. Moreover, considerable effort has been made within the ScaLAPACK project [52] and also by several high-performance computer vendors to optimize these kernels for a variety of parallel machines.

The mechanism used to orthogonalize the new information $Av_k$ against the existing basis $V_k$ is the CGS. It is notoriously unstable and will fail miserably in this setting without modification. One remedy is to use the modified Gram-Schmidt process as used in Algorithm 7.3. Unfortunately, this will also fail to produce orthogonal vectors in the restarting situation we are about to discuss and it cannot be expressed with Level 2 BLAS in this setting. At step (5), we have included an iterative refinement technique proposed by Daniel, Gragg, Kaufman, and Stewart (DGKS) in [96]. This scheme provides an excellent way to construct a vector $f_{j+1}$ that is numerically orthogonal to $V_{j+1}$. The simple test specified at step (5) is used to avoid this DGKS correction if it is not needed.

This mechanism maintains orthogonality to full working precision at very reasonable cost. The special situation imposed by the restarting scheme we are about to discuss makes this modification essential for obtaining accurate eigenvalues and numerically orthogonal Schur vectors (eigenvectors in the Hermitian case). Schemes based on MGS are often sufficient for solving linear systems because they do construct a basis set that is linearly independent even though it might not be numerically orthogonal. The quality of MGS orthogonality is dependent upon the condition number (linear independence) of the original set of vectors. The implicit restarting mechanism we are about to present will be less effective and may even fail if numerical orthogonality is not maintained.

It has been well documented that failure to maintain orthogonality leads to several numerical difficulties within the Lanczos or Arnoldi procedure. In the Hermitian case, Paige [347] showed that the loss of orthogonality occurs precisely when an eigenvalue of $H_j$ is close to an eigenvalue of $A.$ In fact, the Arnoldi vectors lose orthogonality in the direction of the associated approximate eigenvector. Moreover, failure to maintain orthogonality results in spurious copies of the approximate eigenvalue produced by the Arnoldi method. In the Hermitian case, several remedies short of full reorthogonalization have been proposed. These are discussed in §4.4 (p. ).