The joint singular value decomposition of multiple rectangular matrices is formulated as a
Riemannian optimization problem on the product of two Stiefel manifolds. In this paper, the
geometry of the objective function and the Riemannian manifold for this problem are studied
to develop a Riemannian trust-region algorithm. The proposed algorithm globally and locally
quadratically converges, and our numerical experiments demonstrate that it performs much
better than the steepest descent method.