  We propose inertial versions of block coordinate descent methods for solving
non-convex non-smooth composite optimization problems. Our methods possess
three main advantages compared to current state-of-the-art accelerated
first-order methods: (1) they allow using two different extrapolation points to
evaluate the gradients and to add the inertial force (we will empirically show
that it is more efficient than using a single extrapolation point), (2) they
allow to randomly picking the block of variables to update, and (3) they do not
require a restarting step. We prove the subsequential convergence of the
generated sequence under mild assumptions, prove the global convergence under
some additional assumptions, and provide convergence rates. We deploy the
proposed methods to solve non-negative matrix factorization (NMF) and show that
they compete favorably with the state-of-the-art NMF algorithms. Additional
experiments on non-negative approximate canonical polyadic decomposition, also
known as non-negative tensor factorization, are also provided.
