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_{1"0" SRC="img>b = U2^T U1^T b}. U₁^T b can be computed using xORMBR, and xBDSQR has an option to multiply a vector by U₂^T.

Table 2.12: Computational routines for the singular value decomposition

Type of matrix	Operation	Single precision		Double precision
and storage scheme		real	complex	real	complex
general	bidiagonal reduction	SGEBRD	CGEBRD	DGEBRD	ZGEBRD
general band	bidiagonal reduction	SGBBRD	CGBBRD	DGBBRD	ZGBBRD
orthogonal/unitary	generate matrix after	SORGBR	CUNGBR	DORGBR	ZUNGBR
	bidiagonal reduction
	multiply matrix after	SORMBR	CUNMBR	DORMBR	ZUNMBR
	bidiagonal reduction
bidiagonal	SVD using	SBDSQR	CBDSQR	DBDSQR	ZBDSQR
	QR or dqds
	SVD using	SBDSDC		DBDSDC
	divide-and-conquer