In our work, we consider the linear least squares problem for m × n-systems of linear equations Ax = b, m ≥ n, such that the matrix A and right-hand side vector b can vary within an interval m × n-matrix A and an interval m-vector 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 ∈A and b ∈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.