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=λ 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'⁻¹. 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₁,BD₂), with D₁, D₂ arbitrary diagonal, nonsingular matrices. Similarly, the pencil H-λ M is replaced with an equivalent pencil H'-λ M', and the eigenvalues of H-λ M are computed as the squares of the singular values of G=LH LM⁻¹, where LH, LM 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 λ+δλ of any exact eigenvalue λ, the relative error |δλ|/λ is of order p(n)min_Δ∈ Dκ₂(Δ HΔ), min_Δ∈ Dκ₂(Δ MΔ)\, where $p(n)$ is a modestly growing polynomial of the dimension of the problem, {ε} is the round-off unit of floating-point arithmetic, D denotes the set of diagonal nonsingular matrices, and κ₂(·) is the spectral condition number. Furthermore, floating-point computation corresponds to an exact computation with H+δ H, M+δ M, where, for all i, j, |δ Hᵢⱼ|/√HᵢᵢHⱼⱼ and |δ Mᵢⱼ|/√MᵢᵢMⱼⱼ are of order of {ε} times a modest function of n.