This paper considers the problem of recovering the covariance matrix of a stream of high-dimensional data instances from a minimal number of stored measurements. We develop a quadratic random sampling method based on rank-one measurements of the covariance matrix, which serves as an efficient covariance sketching scheme for processing data streams. This also allows modeling of phaseless measurements that arise in high-frequency wireless communication and signal processing applications. We propose to recover the covariance matrix from the above quadratic measurements via convex relaxation with respect to the presumed parsimonious covariance structure. We show that in the absence of noise, exact and universal recovery of low-rank or Toeplitz low-rank covariance matrices can be achieved as soon as the number of stored measurements exceeds the fundamental sampling limit. The convex programs are also robust to noise and imperfect structural assumptions. Our analysis is established upon a novel notion called the mixed-norm restricted isometry property (RIP- ℓ2/ℓ1), as well as the conventional RIP-ℓ2/ℓ2 for near-isotropic and bounded measurements. Our results improve upon best-known phase retrieval performance guarantees with a significantly simpler approach. Numerical results are provided to demonstrate the practical applicability of our technique.