Generalized Symmetric Definite Eigenproblems

This section is concerned with the solution of the generalized eigenvalue problems $A z = \lambda B z$, $A B z = \lambda z$, and $B A z = \lambda z$, where A and B are real symmetric or complex Hermitian and B is positive definite. Each of these problems can be reduced to a standard symmetric eigenvalue problem, using a Cholesky factorization of B as either B=LL^T or B=U^TU (LL^H or U^HU in the Hermitian case). In the case $Ax = \lambda Bz$, if A and B are banded then this may also be exploited to get a faster algorithm.

With B = LL^T, we have

$$Az=\lambda Bz \quad \Rightarrow \quad (L ^{-1}AL ^{-T})(L^Tz)= \lambda(L^Tz).$$

Hence the eigenvalues of $A z = \lambda B z$ are those of $Cy=\lambda y$, where C is the symmetric matrix C = L^-1 A L^-T and y = L^T z. In the complex case C is Hermitian with C = L^-1 A L^-H and y = L^H z.

Table 2.13 summarizes how each of the three types of problem may be reduced to standard form $Cy=\lambda y$, and how the eigenvectors z of the original problem may be recovered from the eigenvectors y of the reduced problem. The table applies to real problems; for complex problems, transposed matrices must be replaced by conjugate-transposes.

Table 2.13: Reduction of generalized symmetric definite eigenproblems to standard problems

	Type of	Factorization	Reduction	Recovery of
	problem	of B		eigenvectors
1.	$A z = \lambda B z$	B = LLT	C = L-1 A L-T	z = L-T y
		B = UTU	C = U-T A U-1	z = U-1 y
2.	$A B z = \lambda z$	B = LLT	C = LT A L	z = L-T y
		B = UTU	C = U A UT	z = U-1 y
3.	$B A z = \lambda z$	B = LLT	C = LT A L	z = L y
		B = UTU	C = U A UT	z = UT y

Given A and a Cholesky factorization of B, the routines xyyGST overwrite A with the matrix C of the corresponding standard problem $Cy=\lambda y$ (see Table 2.14). This may then be solved using the routines described in subsection 2.4.4. No special routines are needed to recover the eigenvectors z of the generalized problem from the eigenvectors y of the standard problem, because these computations are simple applications of Level 2 or Level 3 BLAS.

If the problem is $A z = \lambda B z$ and the matrices A and B are banded, the matrix C as defined above is, in general, full. We can reduce the problem to a banded standard problem by modifying the definition of C thus:

$$C = X^T A X, \quad \mbox{ where } X = U^{-1} Q \quad \mbox{ or } \quad L^{-T} Q,$$

where Q is an orthogonal matrix chosen to ensure that C has bandwidth no greater than that of A. Q is determined as a product of Givens rotations. This is known as Crawford's algorithm (see Crawford [19]). If X is required, it must be formed explicitly by the reduction routine.

A further refinement is possible when A and B are banded, which halves the amount of work required to form C (see Wilkinson [104]). Instead of the standard Cholesky factorization of B as U^T U or L L^T, we use a ``split Cholesky'' factorization B = S^T S (S^H S if B is complex), where:

$$S = \left( \begin{array}{cc} U_{11} & \\ M_{21} & L_{22} \\ \end{array} \right)$$

with U₁₁ upper triangular and L₂₂ lower triangular of order approximately n/2; S has the same bandwidth as B. After B has been factorized in this way by the routine xPBSTF, the reduction of the banded generalized problem $A z = \lambda B z$ to a banded standard problem $Cy=\lambda y$ is performed by the routine xSBGST (or xHBGST for complex matrices). This routine implements a vectorizable form of the algorithm, suggested by Kaufman [77].

Table 2.14: Computational routines for the generalized symmetric definite eigenproblem

Type of matrix	Operation	Single precision		Double precision
and storage scheme		real	complex	real	complex
symmetric/Hermitian	reduction	SSYGST	CHEGST	DSYGST	ZHEGST
symmetric/Hermitian	reduction	SSPGST	CHPGST	DSPGST	ZHPGST
(packed storage)
symmetric/Hermitian	split Cholesky	SPBSTF	CPBSTF	DPBSTF	ZPBSTF
banded	factorization
	reduction	SSBGST	DSBGST	CHBGST	ZHBGST