Superfast and stable structured solvers for Toeplitz least squares via randomized sampling

We present some superfast (O((m + n) log²(m + n)) complexity) and stable structured direct solvers for m× n Toeplitz least squares problems. Based on the displacement equation, a Toeplitz matrix T is first transformed into a Cauchy-like matrix C, which can be shown to have small off-diagonal numerical ranks when the diagonal blocks are rectangular. We generalize standard hierarchically semiseparable (HSS) matrix representations to rectangular ones, and construct a rectangular HSS approximation to C in nearly linear complexity with randomized sampling and fast multiplications of C with vectors. A new URV HSS factorization and a URV HSS solution are designed for the least squares solution. We also present two structured normal equation methods. Systematic error and stability analysis for our HSS methods is given, which is also useful for studying other HSS and rank structured methods. We derive the growth factors and the backward error bounds in the HSS factorizations, and show that the stability results are generally much better than those in dense LU factorizations with partial pivoting. Such analysis has not been done before for HSS matrices. The solvers are tested on various classical Toeplitz examples ranging from well-conditioned to highly ill-conditioned ones. Comparisons with some recent fast and superfast solvers are given. Our new methods are generally much faster, and give better (or at least comparable) accuracies, especially for ill-conditioned problems.