An Implicitly Restarted Joint Bidiagonalization Algorithm for Large GSVD Computations

Abstract. The joint bidiagonalization (JBD) process of a regular matrix pair [Formula: see text] is mathematically equivalent to two joint Lanczos bidiagonalization processes of the upper and lower parts of the Q-factor of QR factorization of the stacked matrix [Formula: see text] when their starting vectors are closely related in a specific way. The resulting JBD method for computing extreme generalized singular values and corresponding left and right generalized singular vectors of [Formula: see text] realizes the standard Rayleigh–Ritz projection of the generalized singular value decomposition (GSVD) problem of [Formula: see text] onto the two left and one right subspaces generated by the JBD process. In this paper, the implicit restarting technique is nontrivially extended to the JBD process, and an implicitly restarted JBD (IRJBD) algorithm is developed with the shifts proposed and a few key implementation details addressed in finite precision arithmetic. Compact upper bounds are established for the residual norm of an approximate GSVD component in both exact and finite precision arithmetic, which are used to design efficient and reliable stopping criteria without explicitly computing approximate left and right generalized singular vectors at each iteration before convergence. Numerical experiments illustrate that IRJBD performs well and is at least competitive with the thick-restart JBD algorithm in terms of restarts and robustness. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SIMAX and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/jiazhongxiao/JBDmethod . [Formula: see text]