The minimization of the relativistic total energy with respect to the linear and nonlinear parameters in the wave function is a notoriously cumbersome task due to the phenomenon usually referred to as "variational collapse." This property of the Dirac matrix equation is generally attributed to the fact that the Dirac operator is not bounded from below. Several methods have been developed to avoid these unphysical solutions. Here we examine the possibility of finding the solutions to the Dirac equation by minimizing the variance of the energy. Monte Carlo integration techniques are used to evaluate the variance functional. The results are very encouraging. No spurious roots are found for any of the basis sets used in the calculations. However, the solutions are not strictly variational and hence are not upper bounds to the energy.