Summary: In our work, we consider the linear least squares problem for \(m \times n\)-systems of linear equations \(Ax = b\), \(m \geq n\), such that the matrix \(A\) and right-hand side vector \(b\) can vary within an interval \(m \times n\)-matrix \(\boldsymbol{A}\) and an interval \(m\)-vector \(\boldsymbol{b}\), respectively. We have to compute, with a prescribed accuracy, outer coordinate-wise estimates of the set of all least squares solutions to \(Ax = b\) for \(A \in \boldsymbol{A}\) and \(b \in \boldsymbol{b}\). Our article is devoted to the development of the so-called PPS-methods (based on partitioning of the parameter set) to solve the above problem. We reduce the normal equation system, associated with the linear lest squares problem, to a special extended matrix form and produce a symmetric interval system of linear equations that is equivalent to the interval least squares problem under solution. To solve such symmetric system, we propose a new construction of PPS-methods, called ILSQ-PPS, which estimates the enclosure of the solution set with practical efficiency. To demonstrate the capabilities of the ILSQ-PPS-method, we present a number of numerical tests and compare their results with those obtained by other methods.