An LR type algorithm is derived to compute the singular values of a product of many matrices. A bound on the backward error on each of the matrices is derived and it is shown how hyperbolic rotations can be used to apply a shift to accelerate convergence. It is shown by numerical examples how the algorithm is useful to establish stability of matrix recurrences approximating dynamical systems.