The least-squares finite element method (LSFEM) has drawn much attention with desirable properties, such as always symmetric positive-definite stiffness matrix and approximation of primal and non-primal variables with equivalent accuracy. Despite a condition number comparable to that of Galerkin FEM, of O(h−2), it is sometimes found that LSFEM with low-order elements does not show a satisfactory accuracy of solution. In this report, we propose an easy-to-implement technique that optimizes the condition number of the stiffness matrix. The least-squares functional is applied to the first-order governing equations and a linear combination of the weak form domain and boundary equations, with normalization (i.e., weighting) parameters is constructed. The impact of these parameters on the condition number of the resulting global stiffness matrix is proved to be independent of mesh size or shape function order. By use of a global optimization scheme that tunes these parameters, we observe notable decrease in the condition number compared to cases where the least-squares principle is directly applied to the governing equations, and hence convergence is significantly improved. It is further shown that the technique can be applied to the elemental stiffness matrix for a computationally efficient means of determining normalization parameters. To show the general applicability of the method, we also apply this technique to solve for a two-dimensional (2D) lid-driven flow problem and to calculate the electromagnetic field distribution in a 2D rectangular resonant cavity. In both problems, the LSFEM solution accuracy is poor prior to optimization of weighting parameters, and accuracy does not increase much by mesh refinement alone. Moreover, the process can be generalized to Galerkin FEM and other methods to improve convergence properties.