  We present a technique for optimizing hundreds of thousands of variational
parameters in variational quantum Monte Carlo. By introducing iterative Krylov
subspace solvers and by multiplying by the Hamiltonian and overlap matrices as
they are sampled, we remove the need to construct and store these matrices and
thus bypass the most expensive steps of the stochastic reconfiguration and
linear method optimization techniques. We demonstrate the effectiveness of this
approach by using stochastic reconfiguration to optimize a correlator product
state wavefunction with a pfaffian reference for four example systems. In two
examples on the two dimensional Hubbard model, we study 16 and 64 site
lattices, recovering energies accurate to 1% in the smaller lattice and
predicting particle-hole phase separation in the larger. In two examples
involving an ab initio Hamiltonian, we investigate the potential energy curve
of a symmetrically dissociated 4x4 hydrogen lattice as well as the
singlet-triplet gap in free base porphin. In the hydrogen system we recover 98%
or more of the correlation energy at all geometries, while for porphin we
compute the gap in a 24 orbital active space to within 0.02eV of the exact
result. The numbers of variational parameters in these examples range from
4x10^3 to 5x10^5, demonstrating an ability to go far beyond the reach of
previous formulations of stochastic reconfiguration.
