Massive MIMO is a variant of multiuser MIMO in which the number of antennas at the base station (BS) $M$ is very large and typically much larger than the number of served users (data streams) $K$. Recent research has illustrated the system-level advantages of such a system and in particular the beneficial effect of increasing the number of antennas $M$. These benefits, however, come at the cost of dramatic increase in hardware and computational complexity. This is partly due to the fact that the BS needs to compute suitable beamforming vectors in order to coherently transmit/receive data to/from each user, where the resulting complexity grows proportionally to the number of antennas $M$ and the number of served users $K$. Recently, different algorithms based on tools from random matrix theory in the asymptotic regime of $M,K \to \infty$ with $\frac{K}{M} \to ρ\in (0,1)$ have been proposed to reduce such complexity. The underlying assumption in all these techniques, however, is that the exact statistics (covariance matrix) of the channel vectors of the users is a priori known. This is far from being realistic, especially that in the high-dim regime of $M\to \infty$, estimation of the underlying covariance matrices is well known to be a very challenging problem. In this paper, we propose a novel technique for designing beamforming vectors in a massive MIMO system. Our method is based on the randomized Kaczmarz algorithm and does not require knowledge of the statistics of the users channel vectors. We analyze the performance of our proposed algorithm theoretically and compare its performance with that of other competitive techniques via numerical simulations. Our results indicate that our proposed technique has a comparable performance while it does not require the knowledge of the statistics of the users channel vectors.