We devise an optimal method to measure the temporal power spectrum of the lensing and intrinsic fluctuations of multiply-imaged strongly lensed quasar light curves , along with the associated time delays . The method is based on a Monte-Carlo Markov Chain ( MCMC ) sampling of a putative gaussian likelihood , and accurately recovers the input properties of simulated light curves , as well as the “ Time Delay Challenge ” . We apply this method to constrain the dimensionless cosmological ( non-linear ) matter power spectrum on milliparsec scales ( comparable to the size of the solar system ) , to \Delta _ { NL } ^ { 2 } < 4 \times 10 ^ { 7 } at k _ { NL } \sim 10 ^ { 3 } { pc } ^ { -1 } . Using a semi-analytic nonlinear clustering model which is calibrated to simulations , the corresponding constraint on the primordial ( linear ) scalar power spectrum is { \cal P } _ { \cal R } < 3 \times 10 ^ { -9 } at k _ { L } \sim 3 pc ^ { -1 } . This is the strongest constraint on primordial power spectrum at these scales , and is within an order of magnitude from the standard \Lambda CDM prediction . We also report measurements of temporal spectra for intrinsic variabilities of quasar light curves , which can be used to constrain the size of the emitting region in accretion disks . Future cadenced optical imaging surveys , such as LSST , should increase the number of observed strongly lensed quasars by 3 orders of magnitude and significantly improve these measurements , even though improvements in modelling quasar accretion and stellar microlensing are necessary .