A new method is presented in this article to determine loads on all rollers in a cylindrical roller bearing. By introducing a pair of nonlinear springs for contact of each roller with its inner and outer races, the equations of equilibrium of the multi-component system are established by means of the virtual work principle. A set of linear complementary equations are deduced and solved using the Lemke algorithm for gaps and contact forces between all potentially engaged components. An iterative scheme is employed to effectively deal with the nonlinearity of the Hertzian contact between non-conforming bodies. Numerical results for a 19-roller cylindrical bearing, having various combinations of roller sizes, show that the proposed method is convergent and accurate. With this method, effects of uneven roller sizes caused by manufacturing errors, on load distributions can be accurately and efficiently determined. The proposed method can be extended to deal with other types of manufacturing errors such as uneven roller angular spacing, eccentricity, ovality, friction, etc., which are of significant interest to the bearing manufacturers.