We propose an algorithm for computing a class of least squares polynomials on polygonal regions of the complex plane. An important application of this technique to solving large sparse linear systems is considered. The advantage of using general polygonal regions instead of ellipses as was done in previous work, is that elliptic regions may fail to accurately represent the convex hull of the spectrum of the matrix A. An attractive feature of the algorithm is that it does not explicitly require any numerical integration. Numerical experiments show that the least-squares based methods for solving linear systems are competitive with the Chebyshev based methods and are more reliable.