We extend our low-scaling variational Monte Carlo (VMC) algorithm to optimize the symmetry-projected Jastrow mean-field (SJMF) wave functions. These wave functions consist of a symmetry-projected product of a Jastrow and a general broken-symmetry mean-field reference. Examples include Jastrow antisymmetrized geminal power, Jastrow Pfaffians, and resonating valence bond states, among others, all of which can be treated with our algorithm. We will demonstrate using benchmark systems including the nitrogen molecule, a chain of hydrogen atoms, and 2-D Hubbard model that a significant amount of correlation can be obtained by optimizing the energy of the SJMF wave function. This can be achieved at a relatively small cost when multiple symmetries including spin, particle number, and complex conjugation are simultaneously broken and projected. We also show that reduced density matrices can be calculated using the optimized wave functions, which allows us to calculate other observables such as correlation functions and will enable us to embed the VMC algorithm in a complete active space self-consistent field calculation.