The performance of branch-and-bound algorithms for solving non-convex optimization problems greatly depends on convex relaxation techniques. They generate convex regions which are used for improving the bounds of variable domains. In particular, convex polyhedral regions can be represented by a linear system A.x=b. Then, bounds of variable domains can be improved by minimizing and maximizing variables in the linear system. Reducing or contracting optimally variable domains in linear systems, however, is an expensive task. It requires solving up to two linear programs for each variable (one for each variable bound). Suboptimal strategies, such as preconditioning, may offer satisfactory approximations of the optimal reduction at a lower cost. In non-square linear systems, a preconditioner P can be chosen such that P.A is close to a diagonal matrix. Thus, the projection of the equivalent system P.A.x=P.b over x, by using an iterative method such as Gauss–Seidel, can significantly improve the contraction. In this paper, we show how to generate an optimal preconditioner, i.e., a preconditioner that helps the Gauss–Seidel method to optimally reduce the variable domains. Despite the cost of generating the preconditioner, it can be re-used in sub-regions of the search space without losing too much effectiveness. Experimental results show that, when used for reducing domains in non-square linear systems, the approach is significantly more effective than Gauss-based elimination techniques. Finally, the approach also shows promising results when used as a component of a solver for non-convex optimization problems.