It is well-known that geometrical variations due to manufacturing tolerances can degrade the performance of optical devices. In recent literature, polynomial chaos expansion (PCE) methods were proposed to model this statistical behavior. Nonetheless, traditional PCE solvers require a lot of memory and their computational complexity leads to prohibitively long simulation times, making these methods non-tractable for large optical systems. The uncertainty quantification (UQ) of various types of large, two-dimensional lens systems is presented in this paper, leveraging a novel parallelized Multilevel Fast Multipole Method (MLFMM) based Stochastic Galerkin Method (SGM). It is demonstrated that this technique can handle large optical structures in reasonable time, e.g., a stochastic lens system with more than 10 million unknowns was solved in less than an hour by using 3 compute nodes. The SGM, which is an intrusive PCE method, guarantees the accuracy of the method. By conjunction with MLFMM, usage of a preconditioner and by constructing and implementing a parallelized algorithm, a high efficiency is achieved. This is demonstrated with parallel scalability graphs. The novel approach is illustrated for different types of lens system and numerical results are validated against a collocation method, which relies on reusing a traditional deterministic solver. The last example concerns a Cassegrain system with five random variables, for which a speed-up of more than 12× compared to a collocation method is achieved.