  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
and an interval $m$-vector respectively. We have to compute, as sharp as
possible, an interval enclosure of the set of all least squares solutions to
$Ax = b$ when $A$ and $b$ independently vary within their interval bounds. 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 original interval least squares problem. 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.
