We consider the problem of computing the cannonical polyadic decomposition (CPD) for large-scale dense tensors. This work is a combination of alternating least squares and fiber sampling. Data sparsity can be leveraged to handle large tensor CPD, but this route is not feasible for dense data. Inspired by stochastic optimization's low memory consuming and per iteration complexity, we propose a stochastic optimization framework for CPD which combines the insights of block coordinate descent and stochastic proximal gradient. At the same time, we also provide an analysis for the convergence property which are unclear to many state-of-art. In addition, our algorithm can handle many frequently used regularizers and constraints, and thus is more flexible compared to many existing algorithms. Applications to hyperspectral images are considered in the experiment analysis part.