A tangent algorithm for computing the generalized singular value decomposition

We present two new algorithms for floating-point computation of the generalized singular values of a real pair (A, B) of full column rank matrices and for floating-point solution of the generalized eigenvalue problem Hx = lambda Mx with symmetric, positive definite matrices H and M. The pair (A, B) is replaced with an equivalent pair (A', B'), and the generalized singular values are computed as the singular values of the explicitly computed matrix F = A'B'(?1). The singular values of F are computed using the Jacobi method. The relative accuracy of the computed singular value approximations does not depend on column scalings of A and B; that is, the accuracy is nearly the same for all pairs (AD(1); BD2), with D-1, D-2 arbitrary diagonal, nonsingular matrices. Similarly, the pencil H - lambda M is replaced with an equivalent pencil H' - lambda M', and the eigenvalues of H - lambda M are computed as the squares of the singular values of G = L-H L-M(-1), where L-H, L-M are the Cholesky factors of H', M', respectively, and the matrix G is explicitly computed as the solution of a linear system of equations. For the computed approximation lambda + delta lambda of any exact eigenvalue lambda, the relative error \delta lambda\/lambda is of order p(n)epsilon max{min(Delta is an element of D) (kappa 2)(Delta H Delta); min(Delta is an element of D) (kappa 2) (Delta M Delta)}, where p(n) is a modestly growing polynomial of the dimension of the problem, epsilon is the round-off unit of floating-point arithmetic, D denotes the set of diagonal nonsingular matrices, and kappa(2)(.) is the spectral condition number. Furthermore, ?oating-point computation corresponds to an exact computation with H + delta H, M + delta M, where, for all i, j, \delta H-ij\/root HiiHjj and \delta M-ij\/root MiiMjj are of order of epsilon times a modest function of n.