Forward modeling approaches in cosmology seek to reconstruct the initial conditions at the beginning of the Universe from the observed survey data.However the high dimensionality of the parameter space poses a challenge to explore the full posterior with traditional algorithms such as Hamiltonian Monte Carlo (HMC) and variational inference (VI).Here we develop a hybrid scheme called variational self-boosted sampling (VBS)that learns a variational approximation for the proposal distribution of HMC with samples generated on the fly, and in turn generates independent samples as proposals for MCMC chain to reduce their auto-correlation length. We use a normalizing flow with Fourier space convolutions as our variational distribution to scale to high dimensions of interest.We show that after a short initial warm-up and training phase, VBS generates better quality of samples than simple VI and reduces the correlation length in the sampling phase by a factor of 10-50 over using only HMC.