Dimensionality reduction and matrix factorization techniques are important and useful machine-learning techniques in many fields. Nonnegative matrix factorization (NMF) is particularly useful for spectral analysis and image processing in astronomy. I present the vectorized update rules and an independent proof of their convergence for NMF with heteroscedastic measurements and missing data. I release a Python implementation of the rules and use an optical spectroscopic dataset of extragalactic sources as an example for demonstration. A future paper will present results of applying the technique to image processing of planetary disks.