I don't quite understand why this solution is preferable to a solution which simply evaluates S(P_j) for each j and returns the minimum in O(n^2) time. I'm also not sure why I would prefer it to a solution that is faster.
One O(n lg n) solution requires that we convert the "Minimum sum of distances at a house for houses on lattice points in 2-space under the square metric" problem to an equivalent "Minimum sum of distances at a house for houses on lattice points in 2-space under the taxicab metric" problem, at which point we can consider the two dimensions independently. With some trickery we can find the sum of distances to other houses for all houses in the 1-dimensional version of the problem in linear time plus the amount of time required to sort the houses.
To do this, we begin at the house with the lowest coordinate and calculate its sum of distances. From the the sum of distances for one house, we can calculate the sum of distances for the next house in constant time: S(P_j) = S(P_j-1) + |P_j - P_j-1|(j) - |P_j - P_j-1|(n-j). This is not the most compact form and is an abuse of the notation considering that it refers only to the 1-dimensional problem, but I think it conveys its meaning well: Once the houses are sorted, the sum of distances for a house other than the first is the sum of distances for the previous house plus the distance between them for each house with a lower coordinate than this house, minus the distance between them for each house with a higher coordinate than this house. The phrasing of that is not quite right for houses that have the same coordinate, but that doesn't particularly matter since the distance between them will be 0.
After that we can calculate S(P_j) as the the sum of distances for P_j along one dimension and along the other. Finally, since we've calculated the sum of distances for each house, we can just print the least of these. We could also solve the problem in linear time if we were willing to use radix sort, but we probably wouldn't want to do that.