The standard non-negative matrix factorization focuses on batch learning assuming that the fixed global latent parameters completely describe the observations. Many online extensions assume rigid constraints and smooth continuity in observations. However, the more complex time series processes can have multivariate distributions switch between a finite number of states or regimes. In this paper we proposes a regime-switching model for non-negative matrix factorization and present a method of forecasting in this lower-dimensional regime-dependent space. The time dependent observations are partitioned into regimes to enhance factors' interpretability inherent in non-negative matrix factorization. We use weighted non-negative matrix factorization to handle missing values and to avoid needless contamination of observed structure. Finally, we propose a method of forecasting from the regime components via threshold autoregressive model and projecting the forecasts back to the original target space. The computation speed is improved by parallelizing weighted non-negative matrix factorization over multiple CPUs. We apply our model to hourly air quality measurements by building regimes from deterministically identified day and night observations. Air pollutants are then partitioned, factorized and forecasted, mostly outperforming the results standard non-negative matrix factorization with respect of the Frobenius norm of the error. We also discuss the shortcomings of the new model.