Hyperspectral unmixing decomposes a hyperspectral imagery (HSI) into a number of constituent materials and associated proportions. Recently, nonnegative tensor factorization (NTF)-based methods have been proposed for hyperspectral unmixing thanks to their capability in representing an HSI without any information loss. However, tensor factorization-based HSI processing approaches often suffer from low-signal-to-noise ratio condition of HSI and nonuniqueness of the solution. This problem can be effectively alleviated by introducing various spatial constraints into tensor factorization to suppress the noise and decrease the number of extreme, stationary, and saddle points. On the other hand, total variation (TV) adaptively promotes piecewise smoothness while preserving edges. In this paper, we propose a TV regularized matrix-vector NTF method. It takes advantage of tensor factorization in preserving global spectral-spatial information and the merits of TV in exploiting local spatial information, thus generating smooth abundance maps with preserved edges. Experimental results on synthetic and real-world data show that the proposed method outperforms the state-of-the-art methods.