<p style='text-indent:20px;'>A primary interest in dynamic inverse problems is to identify the underlying temporal behaviour of the system from outside measurements. In this work, we consider the case, where the target can be represented by a decomposition of spatial and temporal basis functions and hence can be efficiently represented by a low-rank decomposition. We then propose a joint reconstruction and low-rank decomposition method based on the Nonnegative Matrix Factorisation to obtain the unknown from highly undersampled dynamic measurement data. The proposed framework allows for flexible incorporation of separate regularisers for spatial and temporal features. For the special case of a stationary operator, we can effectively use the decomposition to reduce the computational complexity and obtain a substantial speed-up. The proposed methods are evaluated for three simulated phantoms and we compare the obtained results to a separate low-rank reconstruction and subsequent decomposition approach based on the widely used principal component analysis.</p>