Nonnegative Tucker Decomposition (NTD) is a powerful tool to extract nonnegative parts-based and physically meaningful latent components from high-dimensional tensor data, while providing natural multiway representations. However, as the data tensor often has multiple modes and is large-scale, existing NTD algorithms suffer from very high computational complexity in terms of both storage and computation time, which has been one major obstacle for practical applications of NTD. To overcome these disadvantages, in this paper we show how low (multilinear) rank approximation (LRA) of tensors is able to significantly simplify the computation of the gradients of the cost function, upon which a family of efficient first-order NTD algorithms are developed. Besides dramatically reducing the storage complexity and running time, the new algorithms are quite flexible and robust to noise because any well-established LRA approaches can be easily applied. We also show how nonnegativity incorporating sparsity substantially improves the uniqueness property and partially alleviates the curse of dimensionality of the Tucker decompositions. Simulation results on synthetic and real-world data justify the validity and high efficiency of the proposed NTD algorithms. ; Comment: 13 pages, 8 figures