The Khatri-Rao product is extensively used in array processing, tensor decomposition, and multi-way data analysis. Many applications require a least-squares (LS) Khatri-Rao factorization. In broadband sensor array problems, polynomial matrices effectively model frequency-dependent behaviors, necessitating extensions of conventional linear algebra techniques. This paper generalizes LS Khatri-Rao factorization from ordinary to polynomial matrices by applying it to the discrete Fourier transform (DFT) samples of polynomial matrices. Phase coherence across bin-wise Khatri-Rao factors is ensured via a phasesmoothing algorithm. The proposed method is validated through broadband angle-of-arrival (AoA) estimation for uniform planar arrays (UPAs), where the steering matrix is a polynomial matrix, which can be represented as a Khatri-Rao product between steering matrix in azimuth and elevation directions.